Estimation of the covariance structure of spatial processes is a fundamental prerequisite for problems of spatial interpolation and the design of monitoring networks. We introduce a nonparametric approach to global estimation of the spatial covariance structure of a random functionZ(x, t) observed repeatedly at timesti(i= 1, …,T) at a finite number of sampling stationsxi(i= 1, 2, …,N) in the plane. Our analyses assume temporal stationarity but do not assume spatial stationarity (or isotropy). We analyze thespatial dispersionsvar(Z(xi, t) −Z(xj, t)) as a natural metric for the spatial covariance structure and model these as a general smooth function of the geographic coordinates of station pairs (xi, xj). The model is constructed in two steps. First, using nonmetric multidimensional scaling (MDS) we compute a two-dimensional representation of the sampling stations for which a monotone function of interpoint distancesδijapproximates the spatial dispersions. MDS transforms the problem into one for which the covariance structure, expressed in terms of spatial dispersions, is stationary and isotropic. Second, we compute thin-plate splines to provide smooth mappings of the geographic representation of the sampling stations into their MDS representation. The composition of this mappingfand a monotone functiongderived from MDS yields a nonparametric estimator of var(Z(xa, t) −Z(xb, t)) for any two geographic locationsxaandxb(monitored or not) of the formg(|f(xa) −f(xb)|). By restricting the monotone functiongto a class of conditionally nonpositive definite variogram functions, we ensure that the resulting nonparametric model corresponds to a nonnegative definite covariance model. We usebiorthogonal grids, introduced by Bookstein in the field of morphometrics, to depict the thin-plate spline mappings that embody the nature of the anisotropy and nonstationarity in the sample covariance matrix. An analysis of mesoscale variability in solar radiation monitored in southwestern British Columbia demonstrates this methodology.