krig.image {fields}R Documentation

Spatial process estimate for large irregular 2-d dats sets.


Computes the spatial predictions for large numbers of irregularly spaced observations using the standard Kriging equations. The main approximation is that the locations are discretized to a regular grid, but the field need not be observed at all grid boxes.

In Bayesian terms this function computes the posterior mean for the field given the observations under the usual Gaussian assumptions for the fields and observations. The solution is found by the iterative solution of a large linear system using the conjugate gradient algorithm (CGA). Part of the calculations rely on discretizing the spatial locations to a regular grid to make use of the FFT for fast multiplication of a covariance matrix with a vector.


krig.image(x, Y, cov.function, m=NULL, n=NULL, lambda=0, start=NULL, 
tol=1e-05, kmax=25,  cov.obj=NULL, grid=NULL, 
weights=rep(1, length(Y)), verbose=FALSE, conv.verbose=FALSE, expand=1, ...)


x A 2 column matrix of observed locations
Y Values of observed field. Missing values are omitted from computation.
cov.function An S function that multiplies the covariance matrix by a vector. Two that are part of FIELDS are exp.image.cov ( Exponential and Gaussian) and W.image.cov ( W transform covariance model)
lambda The value of the smoothing parameter. Should be nonnegative. See the notes below for more information about this parameter
m Number of grid points in the x axis. Default is to use the length of grid$x.
n Number of grid points in the y axis. Default is to use the length of grid$y.
cov.obj A covariance object that contains information to be used by the covariance function specified above. If this is not specified this object will be created within krig.image.
grid A list with components x and y that specify the grid points in the x and y directions. The default is to use the number of point specified by m and n and use the ranges from the observed locations.
start Starting values for omega2 in the iterative algorithm. Default is zero.
tol Convergence tolerance for CGA.
kmax Maximum number of iterations for CGA
weights This vector is proportional to the reciprocal variance of the measurement errors. The default is a vector of ones.
verbose If true all kinds of stuff is printed out! Default is false of course.
conv.verbose If true the convergence criterion is printed out at each iteration of the CGA. The values are scaled as the criterion divided by the tolerance. So the algorithm terminates when the values are less than one.
expand The amount the grid should be expanded beyond the range of the observed data. For example expand 1.1 will give a range that is 10 % larger on each end.
... Any extra arguments are considered as information for the covariance function and are used to create the covariance object.


From a functional point of view krig.image and supporting functions are similar to the class Krig. The main difference is that only 2-dimensional problems are considered and the solution is calculated for a fixed value of lambda. (The Krig function can estimate lambda.) For large data sets a practical way to estimate lambda is by out of sample cross-validation and the FIELDS manual gives a detailed example of this for the precip data set. Also see the manual for an explanation of the computational strategy (Conjugate Gradient) here.

Efficiency for large datasets comes with restrictions on the range of covariance functions and some other features. Currently FIELDS just has two covarince models: exponential/Gaussian and wavelet based. However, it is not difficult to modify these to other models. The default discretization is to a 64X64 grid however even 256X256 is manageable and quite likely to separate irregular locations in most cases. The user should also keep in mind that the estimate is the result of an iterative algorithm and so issues such as good starting values and whether the algorithm converged are present.

The spatial model includes a linear spatial drift and MLE estimates of the nugget variance and sill are found based on the values of lambda. If the weights are all equal to one and the covariance function is actually a correlation function, in the notation of this function, the "sill" is sigma2 + rho and the "nugget" is sigma2. Moreover sigma2 and rho are constrained so sigma2/rho =lambda. This is why lambda is the crucial parameter in this model.

Although the field is only estimated to the resolution of the grid, prediction off of the grid is supported by bilinear interpolation using the FIELDS function interp.surface.


An list object of class krig.image. An explanation of some components:

call The calling sequence
cov.function A copy of the covariance S function
na.ind logical indicating missing values in Y
xraw Passed spatial locations having removed missing values
y Observations having omitted missing values
N Length of y
weights passed weights having omitted missing cases.
grid list with components x an y indicating grid for discretization
cov.obj List object to use with cov.function
m Number of grid point in x axis
n Number of grid point in y axis
index A two column matrix indicating the indices of the closest grid point to each observed location.
x Observed locations discretized to nearest grid point
yM Observed values but with a weighted average replacing multiple values associated with the same grid point.
xM Discretized locations associated with yM
weightsM Weight vector associated with YM.
uniquerows Logical indicating which rows of x are unique.
shat.rep Pooled standard deviation among observations that fall within the same gird points
indexM A two column matrix indicating the indices of the closest grid point to each observed location, yM.
qr.T QR decomposition of the matrix of constant and linear terms at xM
multAx The S function that is used for matrix multiplication in the CGA.
omega2 Parameter vector that describes the spatial process part of the conditional mean.
converge CGA convergence information
beta Constant, and the two linear parameters for the fixed linear part of the model
delta Covariance matrix times delta give the spatial predictions.
rhohat, rho Conditional on lambda the MLE for the parameter multiplying the covariance function.
sigma2, shat.MLE Conditional on lambda the MLE for the parameter dividing the weight function.
surface A list giving the predicted surface at the grid points.
fitted.values Predicted values at true locations


Large spatial prediction problems and nonstationary fields (1998) Nychka, D., Wikle, C. and Royle, J.A.

FIELDS manual

See Also

plot.krig.image, predict.krig.image, exp.image.cov, sim.krig.image


# fit a monthly precipitation  field over the Rocky Mountains 
# grid is 64X64 
out<- krig.image( x= RMprecip$x, Y = RMprecip$y, m=64,n=64,cov.function= 
lambda=.5, theta=1, kmax=100)  

# range parameter for exponential here is .5 degree in lon and lat.  
#diagnostic plots.  
plot( out)  

# look at the surface 
image.plot( out$surface)  #or just      surface( out)

#simulate 4 realizations from the conditional distribution 
look<- sim.krig.image( out, nreps=4) 
# take a look:  plot( look)
# check out another values of lambda reusing some of the objects from the 
# first fit 

out2<- krig.image( RMprecip$x, RMprecip$y, cov.function= exp.image.cov, 
start= out$omega2,cov.obj=out$cov.obj)
# some of the obsare lumped together into a singel grid box
# find residuals among grid box means and predictions
res<- predict( out2, out2$xM) - out2$yM
#compare with sizes of out2$residuals (raw y data) 

#starting values from first fit  in out$omega2 
# covariance and grid information are 
# bundled in the cov.obj 

## fitting a thin plate spline. The default here is a linear null space 
## and  second derivative type penalty term.  
##  you will just have to try different values of lambda vary them on  
##  log scale to 
out<- krig.image( RMprecip$x, RMprecip$y, cov.function=rad.image.cov, 
lambda=1, m=64, n=64, p=2, kmax=300) 

# take a look: image.plot( out$surface) 

# check out  different values reuse some of the things to make it quicker 
# note addition of kmax argument to increase teh number of iterations

out2<- krig.image( RMprecip$x, RMprecip$y,cov.function=rad.image.cov,  
lambda=.5, start= out$omega2, cov.obj=out$cov.obj, kmax=400) 

# here is something rougher  
out3<- krig.image( RMprecip$x, RMprecip$y,cov.function=rad.image.cov, 
lambda=1e-2, start= out2$omega2, cov.obj=out$cov.obj,kmax=400,

# here is something close to an interpolation  
out4<- krig.image( RMprecip$x, RMprecip$y,cov.function=rad.image.cov, 
lambda=1e-7, start= out3$omega2, cov.obj=out$cov.obj,kmax=500, tol=1e-3)
#compare the the four surfaces: 
# but note the differences in scales ( fix zlim to make them the same) 
# take a look
# set.panel( 2,2) 
# image.plot( out$surface) 
# points( out$x, pch=".") 

# image.plot( out2$surface) 
# image.plot( out3$surface) 
# image.plot( out4$surface) 
# some diagnostic plots) 
set.panel( 4,4) 
plot( out, graphics.reset=FALSE) 
plot( out2, graphics.reset=FALSE) 
plot( out3, graphics.reset=FALSE) 
plot( out4, graphics.reset=FALSE) 

[Package fields version 3.3.1 Index]