CKrige: Spatial interpolation (2D) by constrained, covariance-matching constrained and universal, global Point or Block Kriging

Description

Function for constrained, covariance-matching constrained and universal (external drift kriging) point or block (of any shape) kriging in a global neighbourhood and for isotropic covariance models.

Usage

CKrige( formula, data, locations, object, ...)
## S3 method for class 'formula,data.frame,formula,preCKrigePolygons':
CKrige(formula, data, locations, object, method = 2, ex.out = F)
## S3 method for class 'formula,data.frame,formula,preCKrigePoints':
CKrige(formula, data, locations, object, method = 2, ex.out = F)

Arguments

formula

formula of the linear regression model in the form response ~ terms of covariates, for ordinary kriging use the formula response ~ 1

data

a data frame with the values of the covariates, the names of the covariates used in the formula object must match the column names of data.

locations

a formula object that describe the coordinates of the data locations (e.g ~ x+y)

object

either an object of the class preCKrigePolygons for block kriging or of the class preCKrigePoints for point kriging. In general the output object of the p

...

two further arguments to control the spatial interpolation method and the output

method

numeric value to choose the kriging method 1 universal (external drift), method = 2 constrained, method = 3 covariance-matching constrained kriging. By default, method = 2

ex.out

logical value, if ex.out is set TRUE CKrige returns an extended output with additional information, see details for more informations, by default ex.out = F

Value

By default, CKrige returns an object of the class SpatialPointsDataFrame or SpatialPolygonsDataFrame depending whether the input object for the object argument is of the class preCKrigePoints or preCKrigePolygons. The data frame of the returned object contains the following columns independent of the selected kriging method:
predictionnumeric vector with the kriging prediction of the chosen method
prediction.senumeric vector with the root mean square error (kriging standard error)
The data frame contains 3 additional columns with constrained kriging parameters, if the argument method = 2 of the CKrige function:
sqrt.Pnumeric vector with sqrt( Var[ target point or block ] - Var[ fitted values ] )
sqrt.Qnumeric vector with sqrt( Var[ universal kriging predictor ] - Var[ fitted values ] )
Knumeric vector with sqrt.P / sqrt.Q
The data frame contains 3 additional columns with covariance-matching constrained kriging parameters, if the argument method = 3 of the CKrige function:
P1.11numeric vector, first element of the matrix P1 = ( Cov[target point or block] - Cov[fitted values] )^(1/2)
Q1.11numeric vector, first element of the matrix Q1 = ( Cov[universal kriging predictor] - Cov[fitted values] )^(1/2)
K.11numeric vector, first element of the matrix K = O1^{-1}P1[1,1]
The CKrige function returns a list with the following components if the argument ex.out = T and the argument method is either 1 or 2:
objecteither an object of the class SpatialPolygonsDataFrame or SpatialPointsDataFrame as described above
krig.methodnumeric scalar, number of the chosen kriging method 1, 2 or 3.
parameterlist with 2 components. First component beta.coef is the vector with the Generalized Least Square coefficients of the linear regression and the second component cov.beta contains the covariance matrix of the Generalized Least Square coefficients.
sk.weightsif argument method = 1 or method = 2 sk.weights is a matrix with the simple kriging weights. The ith column contains the simple kriging weights of the ith prediction target object. If the argument method = 3 the list component sk.weights is a list. Each list component contains the matrix with the simple kriging weights of the prediction target and its defined neighbours.
inv.Sigmamatrix, inverse covariance matrix of the data
residualsnumeric vector with the Generalized Least Square residuals of the linear regression.
The list of the extended output contains the additional component CMCK.par if the argument method = 3. The CMCK.par component is a list of lists with CMCK parameters, in particular P1 list of the P1 matrices, Q1 list of the Q1 matrices and K list of the K matrices.

Details

The CKrige function depends always on a preCKrige output object that contains the parameter of the isotropic covariance model as well as the covariates of the prediction targets.

References

See main help page of the constrainedKriging package.

Examples

Run this code

# load data
data(meuse,meuse.blocks)

# approximation of block variance
# pixel area = 75m x 75m
# exponential covariance function with measurement error = 0, nugget = 0.05,
# part. sill =  0.15 and range parameter = 192.5
preCK <- preCKrige(newdata = meuse.blocks,
         model = covmodel("exponential", 0, 0.05, 0.15, scale = 192.5),
	 pwidth = 75, pheight = 75)
    
# block prediction by constrained kriging
ck <- CKrige(formula = log(zinc) ~ sqrt(dist),
    data = meuse, locations = ~ x + y, object = preCK, ex.out = TRUE)
    
# backtransformation to the original scale for the CK prediction
beta <- ck$parameter$beta.coef
c1 <-  0.2
c2  <- beta[2]^2 * meuse.blocks@data$M
c3  <- unlist( preCK@covmat )
ck$object@data$Zn <- exp(ck$object@data$prediction + 0.5*(c1 + c2 - c3))

# factor: the factor multiplied with the corresponding block mean
# prediction gives the upper bound of the 97.5% confidence interval
# of the prediction
ck$object@data$factor <- exp(ck$object@data$prediction
                         + 1.96*ck$object@data$prediction.se)/ck$object@data$Zn

# plots with spplot, generic function in the sp package

# plot of the constrained block kriging prediction
# function ck.colors(n) create a vector of n colors
ticks <- seq(0, 1850, by = 185) 
spplot(ck$object, zcol ="Zn" , col.regions = ck.colors(10), at = ticks,
    colorkey= list(labels = list(at = ticks, labels = ticks)))
    
# plot of the factor to get the upper bound of the 97.5% confidence interval
spplot(ck$object, zcol = "factor", col.regions = ck.colors(10), at = seq(1,3, by = 0.2))

Run the code above in your browser using DataLab