Spatial process estimate of a curve or surface, "kriging" with a known covariance function

Description

This is a simple version of the Krig function that is optimized for large data sets and a clear exposition of the computations.

Usage

mKrig(x, y, weights = rep(1, nrow(x)), 
  lambda = 0, cov.function = "stationary.cov", 
    m = 2, chol.args=NULL,cov.args=NULL, ...)

## S3 method for class 'mKrig':
predict( object, xnew=NULL, derivative=0, ...)
## S3 method for class 'mKrig':
print( x, ... )

Arguments

{ If zero the surface will be evaluated. If not zero the matrix of partial derivatives will be computed.}

Details

This function is an abridged version of Krig that focuses on the computations in Krig.engine.fixed done for a fixed lambda parameter for unique spatial locations and for data without missing values. These restriction simply the code for reading. Note that also little checking is done and the spatial locations are not transformed before the estimation. Sparse matrix methods are handled through overloading the usual linear algebra functions with sparse versions.

The sparse methods rely on two types of temporary workspace. The max.points or mean.neighbor arguments used in fields.rdist.near to find nearest neighboring points and tmpmax connected with the sparse Cholesky decompostion. The default values for each of these is a small multiple of the number of rows of the full matrix. If either of these is too small reset them to something larger. This is typically the total number of nonzero elements in the data covariance matrix.

The temp space for finding near neighbors is max.points = nrow(x1)* mean.neighbor. It is better to set mean.neighbor rather than max.points because this may then appropriate for different size matrices that key off of the same distribution of points. This will be the case in creating the mKrig object and then using this result to predict to a grid.

To explicitly reset mean.neighbor as an additional argument to sparse covariance functions add the argument.

e.g. mean.neighbor =100

To reset tmpmax include a component in the chol.args list.

e.g. chol.args=list( memory=list( tmpmax=5e5) )

will reset temp space to 500,000.

See the example below for more details.

predict.mKrig will evaluate the derivatives of the estimated function if derivatives are supported in the covariance function. For example the wendland.cov function supports derivatives.

print.mKrig is a simple summary function for the object.

Value

Author(s)

Doug Nychka

See Also

Krig

Examples

#
# Midwest ozone data  'day 16' stripped of missings 
data( ozone2)
y<- ozone2$y[16,]
good<- !is.na( y)
y<-y[good]
x<- ozone2$lon.lat[good,]

# nearly interpolate using defaults (Exponential)
mKrig( x,y, theta = 2.0, lambda=.01)-> out
#
# NOTE this should be identical to 
# Krig( x,y, theta=2.0, lambda=.01) 

# interpolate using tapered version the taper scale is set to 1.5
# Default covariance is the Wendland.
# Tapering will done at a scale of 1.5 relative to the scaling 
# done through the theta  passed to the covariance function.

mKrig( x,y,cov.function="stationary.taper.cov",
       theta = 2.0, lambda=.01, Taper.args=list(theta = 1.5, k=2)
           ) -> out2

predict.surface( out2)-> out.p
surface( out.p)

# here is a bigger problem 
# using a compactly supported covariance directly
# 
# Also note increase in the temp space sie for the
# Cholesky decomposition. Default size (1e5) is too small
#

set.seed( 334)
N<- 1000
x<- matrix( 2*(runif(2*N)-.5),ncol=2)
y<- sin( 1.8*pi*x[,1])*sin( 2.5*pi*x[,2]) + rnorm( 1000)*.01
  
mKrig( x,y, cov.function="wendland.cov",k=2, theta=.1, 
            lambda=1e2)-> look2

# The following will fail for theta=.2 because tmpmax and max.points are too
# small. (Here theta controls the support of the covariance and so 
# indirectly the  number of nonzero elements in the sparse matrix

# mKrig( x,y, cov.function="wendland.cov",k=2, theta=.3, lambda=1e2)-> look2

# as a guess on the size of tmpmax this was set to mean.neighbor* nrow(x)

 mKrig( x,y, 
            cov.function="wendland.cov",k=2, theta=.3, 
            lambda=1e2, mean.neighbor=150, 
            chol.args=list( memory=list( tmpmax=150*1000)) 
             )-> look2


predict.surface( look2)-> out.p
surface( out.p)

#
#    Using mKrig for evaluating  a solution on a big grid.
#    (Thanks to Jan Klennin for motivating this example.)

x<- RMprecip$x
y<- RMprecip$y

Tps( x,y)-> obj

# make up an 80X80 grid that has ranges of observations
# use same coordinate names as the x matrix

glist<- fields.x.to.grid(x, nx=80, ny=80) # this is a cute way to get a defautl grid that covers x

# convert grid list to actual x and y values ( try plot( Bigx, pch="."))
    make.surface.grid(glist)-> Bigx 

# include actual x locations along with grid. 
    Bigx<- rbind( x, Bigx)

# evaluate the surface on this set of points (exactly)

    predict(obj, x= Bigx)-> Bigy

# theta sets range for the compact covariance function 
# this will involve  less than 20 nearest neighbors tha have
# nonzero covariance

    theta<- c( 2.5*(glist$lon[2]-glist$lon[1]), 
                 2.5*(glist$lat[2]-glist$lat[1]))

# this is an interplotation of the values using a compact 
# but thin plate spline like covariance. 
    mKrig( Bigx,Bigy, cov.function="wendland.cov",k=4, theta=theta, 
                 lambda=0)->out2 
# the big evaluation this takes about 45 seconds on a Mac G4 latop
    predict.surface( out2, nx=400, ny=400)-> look

# the nice surface
## Not run: 
    
    surface( look)
    US( add=TRUE, col="white")
## End(Not run)