preCKrige {constrainedKriging}

Spatial Variance-Covariance Matrices for Sets of Points and Polygons

Description

The function preCKrige computes (approximated) spatial variance-covariance matrices for user-defined sets of points or polygons (blocks) of any shape for two-dimensional isotropic random fields. The areas of a set of polygons (polygon neighbourhood configuration) are approximated by pixels and the block-block covariances are approximated by averaging covariances between the pixels used to approximate the polygons.

The object returned by preCKrige is needed by CKrige for computing spatial point or block predictions by constrained, covariance-matching constrained or universal (external drift) Kriging.

Usage

preCKrige(newdata, neighbours, model, ...)
## S4 method for signature 'SpatialPoints,ANY,covmodel'
preCKrige(newdata, neighbours, model)
## S4 method for signature 'SpatialPointsDataFrame,ANY,covmodel'
preCKrige(newdata, neighbours, model)
## S4 method for signature 'SpatialPolygons,ANY,covmodel'
preCKrige(newdata, neighbours, model,
  pwidth = 0, pheight = 0, napp = 1, ncores = 1L,
  fork = !identical( .Platform[["OS.type"]], "windows"))
## S4 method for signature 'SpatialPolygonsDataFrame,ANY,covmodel'
preCKrige(newdata, neighbours,
  model, pwidth = 0, pheight = 0, napp = 1, ncores = 1L,
  fork = !identical( .Platform[["OS.type"]], "windows"))

Arguments

newdata

either an object of the class “SpatialPointsDataFrame” or “SpatialPoints” that contains the coordinates of the prediction points and optionally additional information (covariates) stored in the data slot of the SpatialPointsDataFrame, or an object of the class “SpatialPolygonsDataFrame” or “SpatialPolygons” with the coordinates of the polygons (blocks) for which predictions are computed and optionally additional information (covariates) stored in the data slot of the SpatialPolygonsDataFrame.

neighbours

a list of length n with integer vectors as components. n is equal to the number of points if newdata is an object of class “SpatialPointsDataFrame” or “SpatialPoints” or equal to number of polygons (blocks) if newdata is an object of class “SpatialPolygonsDataFrame” or “SpatialPolygons”.

The ith list component defines the neighbours of the ith point or ith polygon (block) in newdata, which form jointly with the ith point or polygon the so-called point or polygon neighbourhood configuration. If newdata is an object of class “SpatialPolygonsDataFrame” or “SpatialPolygons” the ith list component contains the indices of the neighbouring polygons for the ith polygon. If newdata is an object of class “SpatialPoints” or “SpatialPointsDataFrame” the ith list component contains the row indices of the neighbouring points in the point coordinate matrix. The ith list component is set to integer(0) if the ith polygon or ith point have no (defined) neighbours. By default, the points or polygons have no neighbours.

See the second example below where the function poly2nb of the package spdep is used to build a list of neighbours for target polygons of the data set meuse.blocks.

model

an object of class “covmodel”. The object contains the parameters of the isotropic covariance function, generated by the function covmodel.

...

further arguments if newdata is of class “SpatialPolygonsDataFrame” or
“SpatialPolygons”.

pwidth

a positive numeric scalar, defines the width of the pixels used to approximate the polygon (block) areas.

pheight

a positive numeric scalar, defines the height of the pixels used to approximate the polygon (block) areas.

napp

a positive integer scalar. napp > 1 reduces the block-block variance-covariance approximation error. By default, napp = 1, see Details.

ncores

a positive integer scalar with the number of CPUs to use for parallel computations.

fork

a logical scalar to control whether parallel computations are done by forking using mclapply (non-windows OSes) or by socket clusters using parLapply (windows OS).

Details

If the object newdata is of class “SpatialPolygonsDataFrame” or “SpatialPolygons” then
preCKrige searches the polygon neighbourhood configuration (defined by neighbours) with the largest bounding box and generates a pixel grid that completely covers the largest bounding box. Subsequently, the covariance matrix of this set of pixels is calculated by the spatialCovariance package and the polygon (block) areas of each polygon neighbourhood configuration are approximated by intersecting the polygons with the shifted pixel grid, which yields a pixel representation of the polygon neighbourhood configuration. Finally, the block-block covariances of the polygons are approximated by averaging the covariances of the pixel representation of the polygon neighbourhood configuration.

By default, napp = 1, which means that the approximation of the block-block variance-covariance matrix for each polygon neighbourhood configuration is computed just once. If napp > 1 the approximation of the block-block variance-covariance matrix for one polygon neighbourhood configuration is based on the mean of napp repetitions of the approximation to reduce the approximation error. Each of the napp block-block variance-covariance approximations are based on a new, randomly shifted pixel gird which results each time in a new pixel representation of the polygon neighbourhood configuration. Large values of the argument napp increases the computation time.

There is a plot method plot.preCKrigePolygons for preCKrige output objects of class
“preCKrigePolygons” to visually control the polygon (block) area approximation by the pixels.

Value

preCKrige returns an S4 object, either of class “preCKrigePolygons” if newdata is of class
“SpatialPolygons” or “SpatialPolygonsDataFrame” or an S4 object of class “preCKrigePoints” if newdata is of class “SpatialPoints” or “SpatialPointsDataFrame”.

Notation:

n	number of polygons or points in newdata, i = 1, ..., n
m_i	size of point or polygon neighbourhood configuration
	m_i = 1 + number of (defined) neighbours of the ith point or ith polygon
r_{\mathrm{pix}}	number of pixel grid rows
c_{\mathrm{pix}}	number of pixel grid columns
n_{\mathrm{pix}}	number of pixels in pixel grid n_{\mathrm{pix}} = r_{\mathrm{pix}} \cdot c_{\mathrm{pix}}

An object of class “preCKrigePoints” contains the following slots:

covmat

a list of length n, the ith list component contains the point-point covariance matrix of the ith prediction point and its neighbours, i.e. of the ith point neighbourhood configuration.

posindex

a list of length n, the ith list component contains a vector with the row indices of the m_i - 1 neighbours in the ith point neighbourhood configuration.

model

an object of class “covmodel” with the parameters of the used covariance function.

data

a data frame, which is the data slot of the SpatialPointsDataFrame object. This data frame is used to build the design matrix of the prediction points by the CKrige function. data is empty with dim(data) = (0, 0) if newdata is an object of class “SpatialPoints”.

coords

a matrix with dim(coords) = (n, 2) with the coordinates of the prediction points.

An object of class “preCKrigePolygons” contains the following slots:

covmat

a list of length n, the ith list component contains the approximated block-block covariance matrix of the ith polygon and its neighbours, i.e. of the ith polygon neighbourhood configuration.

se.covmat

a list of length n, the ith list component contains a matrix with the standard errors of the approximated block-block covariances of the ith polygon neighbourhood configuration. Values are equal to NaN for argument napp = 1, see Details.

pixconfig

a list of lists of length n, the ith list component contains a list with the information about the pixels used for the covariance approximation of the ith polygon neighbourhood configuration. The components of pixconfig are described below.

pixcovmat

a matrix, dim(matrix) = (n_{\mathrm{pix}}, n_{\mathrm{pix}} ), with the covariance matrix of the pixels.

model

an object of class “covmodel” with the parameters of the used covariance function.

data

a data frame which is the data slot of the SpatialPolygonsDataFrame object. This data frame is used to build the design matrix of the prediction polygons by the CKrige function. data is empty with dim(data) = (0, 0) if newdata is an object of class “SpatialPolygons”.

polygons

a SpatialPolygons object. A list of length n with the polygons of the newdata object.

The ith component of pixconfig is a list with the following 10 components:

pixcenter

a matrix with dim(pixcenter) = (n_{\mathrm{pix}}, 2) with the coordinates of the pixels centroids for the ith polygon neighbourhood configuration.

rowwidth

preCKrige input argument pheight.

colwidth

preCKrige input argument pwidth.

nrows

a numeric scalar with number of rows r_{\mathrm{pix}} of the pixel grid.

ncols

a numeric scalar with number of columns c_{\mathrm{pix}} of the pixel grid.

no.pix.in.poly

a numeric vector of length m_i, each number indicates by how many pixels a polygon of the ith polygon configuration is approximated.

sa.polygons

a logical vector of length m_i, TRUE means that the ith polygon is treated as a point because its area is smaller than the area of a pixel, and FALSE means that the polygon is approximated by pixels, see Note for more details.

polygon.centroids

a matrix with dim(polygon.centroids) = (m_i, 2) with the coordinates of the polygon centroids of the ith polygon neighbourhood configuration.

posindex

an integer vector of length m_i with indices of the ith polygon and its neighbours as defined by the argument neighbours.

pix.in.poly

is a binary matrix with dim(pix.in.poly) = (n_{\mathrm{pix}}, m_i). pix.in.poly[k, j] = 1 indicates that the centroid of the kth pixel lies in the jth polygon, and pix.in.poly[k, j] = 0 indicates that the kth pixel centroid does not lie in the jth polygon.

Note

A polygon (block) is treated as point if the polygon area is smaller than the (defined) pixel area or if all pixel centroids of the generated pixel grid lie outside the polygon (block) area. If a pixel centroid lies inside a polygon that has a smaller area than a pixel, the pixel is allocated to the polygon (block) by which it shares the largest area.

The point-point correlations are calculated via the internal function CorrelationFct (this function implements a subset of the covariance models available previously in the function CovarianceFct of the archived package RandomFields, version 2.0.71) and the point-block covariances are calculated by the C function PointRectCov of the package.

Author(s)

Christoph Hofer, [e-mail]

References

Hofer, C. and Papritz, A. (2011). constrainedKriging: an R-package for customary, constrained and covariance-matching constrained point or block Kriging. Computers & Geosciences. 37, 1562–1569, doi:10.1016/j.cageo.2011.02.009

See Also

CKrige

Examples

### first example
### load data
data(meuse, package = "sp")
data(meuse.blocks)

### plot blocks
plot(meuse.blocks)

### compute the approximated block variance of each block in meuse.blocks
### without any neighbouring blocks (default, required for in universal
### and constrained Kriging) for an exponential covariance function without
### a measurement error, a nugget  = 0.15 (micro scale white noise process),
### a partial sill variance = 0.15 and a scale parameter = 192.5
### approximation of block variance by pixel of size 75m x 75m
preCK_1 <- preCKrige(newdata = meuse.blocks, model = covmodel(modelname =
    "exponential", mev = 0, nugget = 0.05, variance = 0.15,
    scale = 192.5), pwidth = 75, pheight = 75)

### plot block approximation for block 59
plot(preCK_1, 59)

### second example
### define neighbours by using the poly2nb function
### of the spdep package
if(!requireNamespace("spdep", quietly = TRUE)){
  stop("install package spdep to run example")
}
neighbours <- spdep::poly2nb(meuse.blocks)
class(neighbours)
### neighbours should be an object of class "list"
class(neighbours) <- "list"
### compute the approximated block variance-covariance
### matrices of each block in meuse.blocks without the
### defined block neighbours
preCK_2 <- preCKrige(newdata = meuse.blocks, neighbours = neighbours,
  model = covmodel("exponential", nugget = 0.05, variance = 0.15,
    scale = 192.5), pwidth = 75, pheight = 75)

### plot block approximation of block 59 and its
### block neighbours
plot(preCK_2, 59)