This paper deals with a shape preserving method of interpolating positive data at points of the plane in R 2 . We formulate a problem in order to define a positive interpolation variational spline in the Sobolev space H m (Ω), where Ω is a non-empty bounded set of R 2 , by minimizing the semi-norm of order m, and we discrete such problem in a finite element space. An algorithm allows us to compute the resulting function. Some convergence theorems are established. The error is of the order O(1/N) m when N tends to +~, being N the number of the interpolating points. Some numerical and graphical examples are given in order to test the validity of this method.