We propose a multivariate smoothing method based on products of localized orthogonal polynomial series estimators for a smooth regression surface in the fixed-design regression model. The estimation of partial derivatives is included. The proposed method provides for automatic and efficient boundary modifications near the edges of the surface, assuming that the boundary of the support of the regression function satisfies some regularity conditions. By allowing for a preaveraging step, the corresponding algorithms are speeded up considerably and are easy to implement. Computation of special boundary kernels, as required by the kernel method to avoid edge effects, is not necessary. It is shown that under sufficient smoothness assumptions, the global average mean squared error has the same optimal rate of convergence as the mean squared error at an interior point; that is, the boundary correction is asymptotically effective. The method depends on two smoothing parameters, one determining the amount of preaveraging and the other determining the amount of smoothing after preaveraging. Theoretical and practical bounds for the choice of these parameters are discussed. A Monte Carlo study based on a bivariate Gaussian surface indicates that increasing the preaveraging parameterδhas a negative effect on the average mean squared error, which is not unexpected. On the other hand, larger values ofδare computationally more economical. The effects of boundary correction as compared to noncorrected estimates are investigated for the example of a quadratic surface. The numerical complexity of the proposed method is discussed. The methods are demonstrated and compared to kriging for two data sets, one on nonuniformly measured groundwater levels in Arizona and the other on cover-clay thickness data from Iran measured on a regular mesh. The two data analyses include regular and irregular designs and supports; they seem to indicate that the method works well, particularly when compared to kriging.