This paper proposes a data-driven approach for constructing firmly nonexpansive operators. We demonstrate its applicability in Plug-and-Play (PnP) methods, where classical algorithms such as Forward-Backward splitting, Chambolle-Pock primal-dual iteration, Douglas-Rachford iteration or alternating directions method of multipliers (ADMM), are modified by replacing one proximal map by a learned firmly nonexpansive operator. We provide sound mathematical background to the problem of learning such an operator via expected and empirical risk minimization. We prove that, as the number of training points increases, the empirical risk minimization problem converges (in the sense of Gamma-convergence) to the expected risk minimization problem. Further, we derive a solution strategy that ensures firmly nonexpansive and piecewise affine operators within the convex envelope of the training set. We show that this operator converges to the best empirical solution as the number of points in the envelope increases in an appropriate way. Finally, the experimental section details practical implementations of the method and presents an application in image denoising, where we consider a novel, interpretable PnP Chambolle-Pock primal-dual iteration.