Many geophysical and environmental problems depend on estimating a spatial
process based on a large number of irregularly spaced observations. In
addition, such processes often have nonstationary structure.  Standard
spatial estimates are based on the solution to a large linear system that
is of the same size as the number of observations and a straight forward
solution is prohibitive for large, spatial datasets.  
<p>
In this work we
illustrate how the iterative solution of linear systems makes it possible
to find good approximate solutions to large spatial problems. One key to
this approach is the use of spatial process models that allow for
efficient multiplication of the covariance matrix by an arbitrary vector.
This is possible for stationary covariances by making connections with
spatial convolution kernels.  Because many spatial processes are
heterogeneous in space an important extension of these algorithms is to
covariances that are nonstationary.  Recent developments in
multiresolution bases such as wavelets allow one flexibility in capturing
nonstationary structure and also permit rapid evaluation of the basis
functions.  The net result is a flexible models for spatial structure that
can be computed efficiently in an interactive data analysis environment.  
<p>
As an illustration of this approach we consider monthly precipitation
recordsat approximately 800 sites in the Rocky Mountain region of the US.
Here the goal is to produce spatial predictions on a regular grid and
sample the posterior distribution to produce valid realizations of the
precipitation spatial process.