It is desired to estimate a real valued function F on the unit square having observed F with error at N points in the square. F is assumed to be drawn from a particular Gaussian process and measured with independent Gaussian errors. The proposed estimate is the Bayes estimate of F given the data. The roughness penalty corresponding to the prior is derived and it is shown how the Bayesian technique can be regarded as a generalisation of variance components analysis. The proposed estimate is shown to be consistent in the sense that the expected squared error averaged over the data points converges to zero as $N\rightarrow\infty$. Upper bounds on the order of magnitude of magnitude of the expected average squared error are calculated. The proposed technique is compared with existing spline techniques in a simulation study. Generalisations to higher dimensions are discussed.