# Chapter 12 Spatial interpolation of point data

*Last updated: 2024-06-24 11:09:53 *

## Aims

Our aims in this chapter are:

- Calculate an empirical variogram
- Fit a variogram model
- Interpolate using three methods:
- Inverse Distance Weighted (IDW) interpolation
- Ordinary Kriging (OK)
- Universal Kriging (UK)

- Evaluate interpolation accuracy using Leave-One-Out Cross Validation

We will use the following R packages:

`sf`

`stars`

`gstat`

`automap`

`sp`

## 12.1 What is spatial interpolation?

### 12.1.1 Interpolation models

**Spatial interpolation** is the prediction of a given phenomenon in unmeasured locations (Figures 12.1–12.2). For that, we need a spatial interpolation **model**—a set of procedures to calculate **predicted** values of the variable of interest, given **calibration** data.

Calibrarion data usually include:

**Field measurements**—available for a limited number of locations, for example: rainfall data from meteorological stations**Covariates**—available for each and every location within the area of interest, for example: elevation from a DEM

Spatial interpolation models can be divided into two general categories:

**Deterministic models**—Models using arbitrary parameter values, for example: IDW**Statistical models**—Models using parameters chosen objectively based on the data, for example: Kriging

Keep in mind that data *structure* does not imply *meaning*. It is technically possible to interpolate any numeric variable measured in a set of points, however it does not always make sense to do so. For example, it does not make sense to spatially interpolate point data when they refer to a localized phenomenon, such as amount of emissions per power plant (Figure 12.3).

Converesly, it does not make sense to sum up point measurements of a continuous phenomenon (Figure 12.4).

### 12.1.2 The weighted average principle

Many of the commonly used interpolation methods, including the ones we learn about in this Chapter (Nearest Neighbor, IDW, Kriging), are based on the same principle, where a predicted value is a **weighted average** of neighboring points. Weight are usually inveresely related to distance, i.e., as distance increases the weight (importance) of the point decreases. The predicted value for a particular point is calculated as a weighted average of all measured values in measurement points (Equation (12.1)):

\[\begin{equation} \hat{Z}(s_{0})=\frac{\sum_{i=1}^{n}w(s_{i})Z(s_{i})}{\sum_{i=1}^{n}w(s_{i})} \tag{12.1} \end{equation}\]

where:

- \(\hat{Z}(s_{0})\) is the predicted value at location \(s_{0}\)
- \(w(s_{i})\) is the weight of measured point \(i\)
- \(Z(s_{i})\) is the value of measured point \(i\)

The weight \(w(s_{i})\) of each measured point is a function of distance (Figure 12.5) from the predicted point.

In IDW, the weight is the inverse of distance to the power of \(p\) (Equation (12.2)):

\[\begin{equation} w(s_{i})=\frac{1}{d(s_{0}, s_{i})^p} \tag{12.2} \end{equation}\]

where:

- \(w(s_{i})\) is the weight of measured point \(i\)
- \(d(s_{0}, s_{i})\) is the distance between predicted point \(s_{0}\) and measured point \(s_{i}\)

The default value for \(p\) is usually \(p=2\) (Equation (12.3)):

\[\begin{equation} w(s_{i})=\frac{1}{d(s_{0}, s_{i})^2} \tag{12.3} \end{equation}\]

The \(p\) parameter basically determines how steeply does weight increase with proximity. As a result, \(p\) determines whether weights are more or less equally distributed among neighbors (low \(p\)) or whether one point (the nearest) has overwhelmingly high weight and thus the predicted value will be strongly influenced by that point (high \(p\)). In other words, when \(p\) approaches zero, the predicted result will approach a uniform surface which is just an average of all measured points. When \(p\) approaches infinity, the predicted result will approach nearest **neighbor interpolation**, which is the simplest *spatial* interpolation method there is: every predicted point gets the value of the nearest measured point (Figures 12.6–12.7).

In Kriging, the weight is a particular function of distance known as the **variogram model** (Figure 12.8). The variogram model is fitted to characterize the autocorrelation structure in the measured data, based on the **empirical variogram**.

There are two frequently used kriging methods: Ordinary Kriging (OK) and Universal Kriging (UK). Adding up the Inverse Distance Weighted (IDW) interpolation, we now mentioned three interpolation methods. We are going to cover those three methods (Figure 12.9), mostly from the practical point of view, in the next three sections (Sections 12.2–12.4).

For the examples, we will load the `rainfall.csv`

file (Section 4.4.3), calculate the `annual`

column (Section 4.5) and convert it to a point layer (Section 7.4):

```
library(sf)
rainfall = read.csv('data/rainfall.csv')
rainfall = st_as_sf(rainfall, coords = c('x_utm', 'y_utm'), crs = 32636)
m = c('sep', 'oct', 'nov', 'dec', 'jan', 'feb', 'mar', 'apr', 'may')
rainfall$annual = apply(st_drop_geometry(rainfall[, m]), 1, sum)
```

We will also use a \(1\times1\) \(km^2\) DEM raster of the area of interest (Sections 9.1, 9.3 and 10.1):

```
library(stars)
dem1 = read_stars('data/srtm_43_06.tif')
dem2 = read_stars('data/srtm_44_06.tif')
dem = st_mosaic(dem1, dem2)
borders = st_read('data/israel_borders.shp', quiet = TRUE)
grid = st_as_sfc(st_bbox(borders))
grid = st_as_stars(grid, dx = 1000, dy = 1000)
dem = st_warp(src = dem, grid, method = 'average', use_gdal = TRUE, no_data_value = -9999)
dem = dem[borders]
names(dem) = 'elev_1km'
```

Finally, we subset the `dem`

to include only the area to the north of 31 degrees latitude, where meteorological station density is relatively high:

```
y = st_as_sf(dem, as_points = TRUE)
y$lat = st_coordinates(st_transform(y, 4326))[,2]
y = st_rasterize(y[, 'lat'], dem)
y[y < 31] = NA
y = st_as_sf(y, merge = TRUE)
dem = dem[y]
```

Next, we *extract* the elevation values (Section 10.7.2):

and subset those stations that coincide with the raster:

Figure 12.10 shows the `dem`

elevation raster and the `rainfall`

point layer:

```
plot(dem, breaks = 'equal', col = terrain.colors(11), reset = FALSE)
plot(st_geometry(rainfall), add = TRUE)
```

## 12.2 Inverse Distance Weighted interpolation

### 12.2.1 The `gstat`

object

To interpolate, we first need to create an object of class `gstat`

, using a function of the same name: `gstat`

. A `gstat`

object contains all necessary information to conduct spatial interpolation, namely:

- The
**model**definition - The calibration
**data**

Based on its arguments, the `gstat`

function “understands” what type of interpolation model we want to use:

- No variogram model →
**IDW** - Variogram model, no covariates →
**Ordinary Kriging** - Variogram model, with covariates →
**Universal Kriging**

The complete decision tree of `gstat`

, including several additional methods which we are not going to use, is shown in Figure 12.11.

We are going to use three parameters of the `gstat`

function:

`formula`

—The prediction**“formula”**specifying the dependent and the independent variables (covariates)`data`

—The calibration**data**`model`

—The**variogram**model

Keep in mind that we need to specify parameter names, because these three parameters are not the first three in the `gstat`

function definition.

For example, to interpolate using the **IDW** method we create the following `gstat`

object, specifying just the `formula`

(Section 12.2.2 below) and `data`

:

### 12.2.2 Working with `formula`

objects

In R, `formula`

objects are used to specify *relation* between objects, in particular—the role of different data columns in statistical models. A `formula`

object is created using the `~`

operator, which separates names of **dependent** variables (to the left of the `~`

symbol) and **independent** variables (to the right of the `~`

symbol). Writing `1`

to the right of the `~`

symbol, as in `~ 1`

, means that there are no independent variables^{38}.

For example, in the following expression we create a `formula`

object named `f`

:

Checking the class shows that `f`

is indeed a `formula`

object:

We can also convert `character`

values to `formula`

using the `as.formula`

function. For example:

The `as.formula`

function is particularly useful when we want to construct different formulas as part of a `for`

loop (Section 12.3.3).

### 12.2.3 Making predictions

Now that our model is defined, we can use the `predict`

function to actually interpolate, i.e., to calculate predicted values. The `predict`

function accepts:

- A
**raster**—`stars`

object, such as`dem`

- A
**model**—`gstat`

object, such as`g`

The raster serves for two purposes:

- Specifying the
*locations*where we want to make predictions (in all methods) - Specifying
*covariate*values (in Universal Kriging only)

For example, the following expression interpolates `annual`

values according to the model defined in `g`

and the raster template defined in `dem`

:

The resulting `stars`

object has two attributes:

`var1.pred`

—the predicted values`var1.var`

—the variance (for Kriging only)

For example:

```
z
## stars object with 2 dimensions and 2 attributes
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## var1.pred 42.8514 357.0842 495.3477 466.7833 551.8437 946.2976 20255
## var1.var NA NA NA NaN NA NA 40194
## dimension(s):
## from to offset delta refsys x/y
## x 1 154 615965 1000 WGS 84 / UTM zone 36N [x]
## y 1 261 3691819 -1000 WGS 84 / UTM zone 36N [y]
```

We can subset just the first attribute and rename it to `'annual'`

:

The interpolated annual rainfall raster, using IDW, is shown in Figure 12.12:

```
b = seq(0, 1200, 100)
plot(z, breaks = b, col = hcl.colors(length(b)-1, 'Spectral'), reset = FALSE)
plot(st_geometry(rainfall), pch = 3, add = TRUE)
contour(z, breaks = b, add = TRUE)
```

## 12.3 Ordinary Kriging

### 12.3.1 Annual rainfall example

Kriging methods require a **variogram model**. The variogram model is an objective way to quantify the autocorrelation pattern in the data, and assign weights accordingly when making predictions (Section 12.1.2).

As a first step, we can calculate and examine the **empirical variogram** using the `variogram`

function. The function requires two arguments:

`formula`

—The dependent variable and the covariates (same as in`gstat`

, see Section 12.2.1)`data`

—A point layer with the dependent variable and covariates as attributes

For example, the following expression calculates the empirical variogram of `annual`

, with no covariates:

We can examine the variogram using `plot`

(Figure 12.13):

There are several ways to fit a **variogram model** to an empirical variogram. We will use the simplest one—*automatic* fitting using function `autofitVariogram`

from package `automap`

:

The function chooses the best fitting type of model, and also fine tunes its parameters. (Use `show.vgms()`

to display variogram model types.)

The fitted model can be plotted with `plot`

(Figure 12.14):

The resulting object is actually a `list`

with several components, including the empirical variogram and the fitted variogram model. The `$var_model`

component of the resulting object contains the actual model:

```
v_mod_ok$var_model
## model psill range kappa
## 1 Nug 695.8141 0.00 0
## 2 Ste 21611.5603 30844.48 5
```

The variogram model can then be passed to the `gstat`

function, and we can carry on with the Ordinary Kriging interpolation:

```
g = gstat(formula = annual ~ 1, model = v_mod_ok$var_model, data = rainfall)
z = predict(g, dem)
## [using ordinary kriging]
```

Again, we will subset the predicted values attribute and rename it:

The Ordinary Kriging predictions are shown in Figure 12.15:

```
b = seq(0, 1200, 100)
plot(z, breaks = b, col = hcl.colors(length(b)-1, 'Spectral'), reset = FALSE)
plot(st_geometry(rainfall), pch = 3, add = TRUE)
contour(z, breaks = b, add = TRUE)
```

### 12.3.2 Elevation example

Another example: suppose that we did not have a DEM for Israel, but only the elevation measurements at the meteorological stations. How can we produce an elevation raster using Ordinary Kriging?

First, we prepare the `gstat`

object:

```
v = autofitVariogram(altitude ~ 1, rainfall)
g = gstat(formula = altitude ~ 1, model = v$var_model, data = rainfall)
```

Then, we interpolate:

The predicted elevation raster is shown in Figure 12.16:

```
b = seq(-500, 1200, 100)
plot(z, breaks = b, col = hcl.colors(length(b)-1, 'Spectral'), reset = FALSE)
plot(st_geometry(rainfall), pch = 3, add = TRUE)
contour(z, breaks = b, add = TRUE)
```

### 12.3.3 Monthly rainfall example

In the next example we use kriging inside a `for`

loop, to make a series of predictions for different variables. Specifically, we will use Ordinary Kriging to predict *monthly rainfall*, i.e., `sep`

through `may`

columns in the `rainfall`

layer.

In each `for`

loop “round”, the formula is going to be re-defined according to the current month `i`

. For example:

First, we set up a vector with the column names of the variables we wish to interpolate, and a `list`

where we “collect” the results:

Second, we specify the `for`

loop, as follows:

```
for(i in m) {
f = as.formula(paste0(i, ' ~ 1'))
v = autofitVariogram(f, rainfall)
g = gstat(formula = f, model = v$var_model, data = rainfall)
z = predict(g, dem)
z = z['var1.pred',,]
result[[i]] = z
}
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
```

Finally, we combine the `list`

of results per month into a single multi-band raster, using `do.call`

and `c`

(Section 11.3.2):

The interpolated montly rainfall amounts are shown in Figure 12.17:

## 12.4 Universal Kriging

**Universal Kriging** interpolation uses a model with one or more independent variables, i.e., covariates. The covariates need to be known for both:

- The
*point layer*, as an attribute such as`"elev_1km"`

in`rainfall`

- The
*predicted locations*, as raster values such as`dem`

values

The `formula`

now specifies the name(s) of the covariate(s) to the right of the `~`

symbol, separated by `+`

if there are more than one. Also, we are using a subset of `rainfall`

where `elev_1km`

values were present:

```
v_emp_uk = variogram(annual ~ elev_1km, rainfall)
v_mod_uk = autofitVariogram(annual ~ elev_1km, rainfall)
```

Comparing the Ordinary Kriging and Universal Kriging variogram models (Figure 12.18):

```
plot(v_emp_ok, model = v_mod_ok$var_model, ylim = c(0, 25000), main = 'OK')
plot(v_emp_uk, model = v_mod_uk$var_model, ylim = c(0, 25000), main = 'UK')
```

Next we create a `gstat`

object, where the `formula`

contains the covariate and the corresponding variogram model:

Remember that all of the variables that appear in the `formula`

need to be present in the `data`

. In this case we have two variables: a dependent variable (`annual`

) and an independent variable (`elev_1km`

).

Now we can make predictions:

and then subset and rename:

Universal Kriging predictions are shown in Figure 12.19:

```
b = seq(0, 1200, 100)
plot(z, breaks = b, col = hcl.colors(length(b)-1, 'Spectral'), reset = FALSE)
plot(st_geometry(rainfall), pch = 3, add = TRUE)
contour(z, breaks = b, add = TRUE)
```

## 12.5 Cross-validation

In Sections 12.2–12.4, we have calculated annual rainfall surfaces using three different methods: IDW, Ordinary Kriging and Universal Kriging. Although it is useful to examine and compare the results graphically (Figures 12.12, 12.15, and 12.19), we also need an objective way to evaluate interolation *accuracy*. That way, we can objectively choose the most accurate method among the various interpolation methods there are.

Plainly speaking, to evaluate prediction accuracy we need to compare the predicted values with measured data in the same location. Since measured data are often sparse and expensive to produce, it makes little sense to collect more data merely for the sake of accuracy assessment. Instead, the available data are usually *split* to two parts, called training and test data. The training data are used to fit the model, while the test data are used to calculate prediction accuracy. The procedure is called **cross validation**. A specific, commonly used, type of cross validation is **Leave-One-Out Cross Validation** where all observations consecutively take the role of test data while the remaning observations take the role of training data. The separation of training and test data is important because evaluating a model based on the same data used to fit it gives the model an “unfair” advantage and therefore overestimates accuracy.

In Leave-One-Out Cross Validation we:

**Take out**one point out of the calibration data- Make a
**prediction**for that point **Repeat**for all points

In the end, what we get is a table with an *observed* value and a *predicted* value for all points.

We can run Leave-One-Out Cross Validation using the `gstat.cv`

function, which accepts a `gstat`

object:

The `gstat.cv`

function returns an object of class `SpatialPointsDataFrame`

(package `sp`

), which we can convert to an `sf`

object with `st_as_sf`

:

```
cv = st_as_sf(cv)
cv
## Simple feature collection with 162 features and 6 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 629301.4 ymin: 3435503 xmax: 761589.2 ymax: 3681163
## Projected CRS: WGS 84 / UTM zone 36N
## First 10 features:
## var1.pred var1.var observed residual zscore fold
## 1 632.4547 959.1984 582.8 -49.6546683 -1.60326700 1
## 2 615.1599 922.4803 608.5 -6.6598763 -0.21927424 2
## 3 598.2919 914.5156 614.7 16.4080581 0.54257730 3
## 4 609.6198 781.6000 562.7 -46.9198381 -1.67827930 4
## 5 703.0883 923.2564 682.8 -20.2882997 -0.66770478 5
## 6 607.7640 747.3236 705.5 97.7360387 3.57520027 6
## 7 644.6945 831.8391 610.4 -34.2945000 -1.18906286 7
## 8 582.0312 983.3018 583.3 1.2687933 0.04046201 8
## 9 654.9796 712.8466 709.8 54.8204170 2.05326166 9
## 10 654.2158 713.4856 654.8 0.5842106 0.02187141 10
## geometry
## 1 POINT (697119.1 3656748)
## 2 POINT (696509.3 3652434)
## 3 POINT (696541.7 3641332)
## 4 POINT (697875.3 3630156)
## 5 POINT (689553.7 3626282)
## 6 POINT (694694.5 3624388)
## 7 POINT (686489.5 3619716)
## 8 POINT (683148.4 3616846)
## 9 POINT (696489.3 3610221)
## 10 POINT (693025.1 3608449)
```

The result of `gstat.cv`

has the following attributes:

`var1.pred`

—Predicted value`var1.var`

—Variance (only for Kriging)`observed`

—Observed value`residual`

—Observed-Predicted`zscore`

—Z-score (only for Kriging)`fold`

—Cross-validation ID

A **bubble plot** is convenient to examine the residuals, since it shows positive and negative values in different color. A bubble plot can be created using the `bubble`

function from package `sp`

, which requires transforming from `sf`

to an `sp`

-compatible class `SpatialPointsDataFrame`

(Figure 12.20):

Using the “predicted” and “observed” columns we can calculate prediction accuracy indices, such as the **Root Mean Square Error (RMSE)** (Equation (12.4)):

\[\begin{equation} RMSE=\sqrt{\frac{\sum_{i=1}^{n} (pred_i-obs_i)^2}{n}} \tag{12.4} \end{equation}\]

where \(pred_i\) and \(obs_i\) are *predicted* and *observed* values for point \(i\), respectively.

For example:

The

`~ 1`

part means there are no independent variables other than the*intercept*.↩︎