---
title: "Kriging from scratch -- Spatial Classes version"
author: "D G Rossiter"
date: "`r Sys.Date()`"
output:
  html_document:
    number_sections: TRUE
    theme: "lumen"
    code_folding: show
    fig_height: 4.5
    fig_width: 4.5
    fig_align: 'center'
    toc: TRUE
    toc_float: TRUE
---

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

# 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.

# 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.

```{r}
require(sp)
require(gstat)
data(meuse)
coordinates(meuse) <- c("x", "y")
class(meuse)
```

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.

```{r}
bbox(meuse)
x0 <- data.frame(x=median(bbox(meuse)["x",]),
                 y=median(bbox(meuse)["y",]))
coordinates(x0) <- c("x", "y")
print(x0)
```

# Build the OK system matrices

Build a matrix of the distances between known points:

```{r}
dim(coordinates(meuse))
dm.A <- spDists(coordinates(meuse))
dim(dm.A)
dm.A[1:5, 1:5]
```

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:

```{r}
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)))
str(A <- variogramLine(object=vmf, dist_vector=dm.A))
A[1:5, 1:5]
```

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

```{r}
dim(A)
A <- cbind(A, 1)
dim(A)
A <- rbind(A, 1)
dim(A)
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.

```{r}
dm <- spDists(rbind(coordinates(x0), coordinates(meuse)))
dim(dm)
# 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)
```

# 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.

```{r}
lambda <- drop(solve(A) %*% b)
head(lambda)
tail(lambda, 1) # LaGrange multiplier
```

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:

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

# Compare to the `krige` function

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

```{r}
(k.pt <- krige(logZn ~1, loc=meuse, newdata=x0,model=vmf))
```

As promised!

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