This tutorial shows how to compute an Ordinary Kriging prediction directly from the OK equations.

1 Ordinary Kriging equations

The prediction by Ordinary Kriging at an unknown point \(x_0\) is computed as a \(\hat{x}_o = \sum_i \lambda_i x_i\) of known points \(\hat{x}_i\), \(i = 1 \ldots n\).

These weights \(\lambda_i\) are dervied by solving the Ordinary Kriging equations: \(\mathbf{A} \lambda = b\), solved as \(\lambda = \mathbf{A}^{-1} b\). The kriging prediction variances are computed as \(b^T \lambda\), where\(A, \lambda, b\) are defined from the semivariances \(\gamma\) and the LaGrange multiplier \(\psi\) as:

\[\begin{eqnarray*} \mathbf{A} & = & \left[ \begin{array}{*5{c}} \gamma(\mathbf{x}_1,\mathbf{x}_1) & \gamma(\mathbf{x}_1,\mathbf{x}_2) & \cdots & \gamma(\mathbf{x}_1,\mathbf{x}_N) & 1 \\ \gamma(\mathbf{x}_2,\mathbf{x}_1) & \gamma(\mathbf{x}_2,\mathbf{x}_2) & \cdots & \gamma(\mathbf{x}_2,\mathbf{x}_N) & 1 \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ \gamma(\mathbf{x}_N,\mathbf{x}_1) & \gamma(\mathbf{x}_N,\mathbf{x}_2) & \cdots & \gamma(\mathbf{x}_N,\mathbf{x}_N) & 1 \\ 1 & 1 & \cdots & 1 & 0 \\ \end{array} \right] ; \mathbf{\lambda} & = & \left[ \begin{array}{c} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_N \\ \psi \end{array} \right] ; \mathbf{b} & = & \left[ \begin{array}{c} \gamma(\mathbf{x}_1,\mathbf{x}_0) \\ \gamma(\mathbf{x}_2,\mathbf{x}_0) \\ \vdots \\ \gamma(\mathbf{x}_N,\mathbf{x}_0) \\ 1 \end{array} \right] \end{eqnarray*}\]

We will build these matrices and solve explicitly for one prediction point.

2 Example data

Load the Meuse topsoil dataset from the sp package and convert to a SpatialPointsDataFrame, with both coördinates and attributes. Rhis class is defined in the sp “Classes and Methods for Spatial Data” package.

require(sp)
## Loading required package: sp
require(gstat)
## Loading required package: gstat
data(meuse)
coordinates(meuse) <- c("x", "y")
class(meuse)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

Place a point to predict somewhere in the middle of the bounding box, and convert it to a SpatialPoints object. Note it has no attributes, we only need its coordinates.

bbox(meuse)
##      min    max
## x 178605 181390
## y 329714 333611
x0 <- data.frame(x=median(bbox(meuse)["x",]),
                 y=median(bbox(meuse)["y",]))
coordinates(x0) <- c("x", "y")
print(x0)
## SpatialPoints:
##             x        y
## [1,] 179997.5 331662.5
## Coordinate Reference System (CRS) arguments: NA

3 Build the OK system matrices

Build a matrix of the distances between known points:

dim(coordinates(meuse))
## [1] 155   2
dm.A <- spDists(coordinates(meuse))
dim(dm.A)
## [1] 155 155
dm.A[1:5, 1:5]
##           [,1]      [,2]     [,3]     [,4]     [,5]
## [1,]   0.00000  70.83784 118.8486 259.2393 366.3141
## [2,]  70.83784   0.00000 141.5662 282.8516 362.6403
## [3,] 118.84864 141.56624   0.0000 143.1712 251.0239
## [4,] 259.23927 282.85155 143.1712   0.0000 154.2628
## [5,] 366.31407 362.64032 251.0239 154.2628   0.0000

Notice the diagonals are 0, this is the distance of a point to itself.

Convert these to semivariances, using a fitted variogram model, as the upper-left of the \(\mathbf{A}\) matrix:

meuse$logZn <- log10(meuse$zinc)
v <- variogram(logZn ~ 1, loc=meuse, cutoff=1300, width=90)
(vmf <- fit.variogram(v, vgm(psill=0.12, model="Sph", range=900, nugget=0.01)))
##   model      psill    range
## 1   Nug 0.01004124   0.0000
## 2   Sph 0.11525701 967.2639
str(A <- variogramLine(object=vmf, dist_vector=dm.A))
##  num [1:155, 1:155] 0 0.0227 0.0312 0.0553 0.0724 ...
A[1:5, 1:5]
##            [,1]       [,2]       [,3]       [,4]       [,5]
## [1,] 0.00000000 0.02267992 0.03117694 0.05526735 0.07238485
## [2,] 0.02267992 0.00000000 0.03516364 0.05915613 0.07182145
## [3,] 0.03117694 0.03516364 0.00000000 0.03544430 0.05390113
## [4,] 0.05526735 0.05915613 0.03544430 0.00000000 0.03737988
## [5,] 0.07238485 0.07182145 0.05390113 0.03737988 0.00000000

Extend the \(\mathbf{A}\) matrix with the column and row of 1’s (constraints), bottom-right entry is 0:

dim(A)
## [1] 155 155
A <- cbind(A, 1)
dim(A)
## [1] 155 156
A <- rbind(A, 1)
dim(A)
## [1] 156 156
n <- dim(A)[1]
A[n,n] <- 0

Build the \(b\) vector from the distances from the target point to the known points. To do this, we compute distances of an expanded matrix of known points, with the target point as the first entry, and then take out the vector that is just the distances between the target and known points.

dm <- spDists(rbind(coordinates(x0), coordinates(meuse)))
dim(dm)
## [1] 156 156
# 2nd..nth rows of the first column are these distances
dm.b <- dm[2:n, 1]
# last entry is 1
b <- c(variogramLine(vmf, dist_vector=dm.b)$gamma, 1)
length(b)
## [1] 156

4 Solve the OK system

Solve for weights and LaGrange multiplier. Note solve with a single matrix argument inverts the matrix, and %*% is the matrix multiplication operator. The drop function removes matrix dimensions with only one entry, here it will be the column, and so converts the \(\lambda\) matrix into a vector.

lambda <- drop(solve(A) %*% b)
head(lambda)
## [1]  0.0001981746 -0.0003767448 -0.0003017196 -0.0001535747 -0.0007040921
## [6] -0.0008275446
tail(lambda, 1) # LaGrange multiplier
## [1] 0.0001190512

Now that we have the \(\lambda\) vector, we use it to solve for the prediction and its variance. Here drop will convert the results to scalars:

# prediction
(ok.pred <- drop(lambda[1:n-1] %*% meuse$logZn))
## [1] 2.270603
# variance
(ok.pred.var <- drop(t(b) %*% lambda))
## [1] 0.0321583

5 Compare to the krige function

Check that we get the same results with the krige function:

(k.pt <- krige(logZn ~1, loc=meuse, newdata=x0,model=vmf))
## [using ordinary kriging]
##            coordinates var1.pred  var1.var
## 1 (179997.5, 331662.5)  2.270603 0.0321583

As promised!

Questions? Comments? mailto:d.g.rossiter@cornell.edu?subject=Comments_on_KrigingFromScratch