The quality of maps obtained by interpolation of observations of a target environmental variable at a restricted number of locations, is partly determined by the spatial pattern of the sample locations. A method is presented for optimization of the sample pattern when the environmental variable is interpolated with the help of exhaustively known covariates, which are assumed to be linearly related to the target variable. In this method the spatially averaged universal kriging variance (MUKV), which incorporates trend estimation error as well as spatial interpolation error, is minimized. The optimal pattern is obtained using simulated annealing. The method requires that the covariance function or variogram of the regression-residuals is known. The method is tested in a case study on the Mean Highest Water table in a coversand area in The Netherlands. The patterns of 25, 50 and 100 sample locations are optimized and compared with the patterns optimized with the ordinary kriging (OK) model (assuming no trend) and with the multiple linear regression (MLR) model (assuming no spatial autocorrelation of residuals). The results show that the UK-patterns are a good compromise between spreading in geographic space and feature space. The MUKV for the UK-patterns is 19% (n = 25), 7% (n = 50) and 3% (n = 100) smaller than for the OK-patterns. Compared with the MLR-patterns the reduction is 2%, 4% and 4%, respectively.