Chapter 12 Spatial interpolation of point data

Last updated: 2023-10-17 21:44:19

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.

Figure 12.1: Spatial interpolation (Input elevation point data, Interpolated elevation surface) (http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/understanding-interpolation-analysis.htm)

Figure 12.2: Spatial interpolation (Point locations of ozone monitoring stations, Interpolated prediction surface) (http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/understanding-interpolation-analysis.htm)

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

Figure 12.3: \(CO_2\) emissions from power plants: a localized phenomenon (https://edzer.github.io/UseR2016/)

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

Figure 12.4: \(PM_{10}\) measurements: a continuous phenomenon (https://edzer.github.io/UseR2016/)

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.

Figure 12.5: Distances between predicted point and all measured points (http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-kriging-works.htm)

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

Figure 12.6: Spatial interpolation of annual rainfall using IDW with \(p=0.25\), \(p=2\) and \(p=16\)

Figure 12.7: Nearest Neighbor interpolation (left) and Voronoi polygons (right)

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.

Figure 12.8: Variogram models: spherical (left) and exponential (right) (http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-kriging-works.htm)

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

Figure 12.9: Spatial interpolation of annual rainfall using IDW, OK and UK

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):

rainfall$elev_1km = st_extract(dem, rainfall)$elev_1km

and subset those stations that coincide with the raster:

rainfall = rainfall[!is.na(rainfall$elev_1km), ]

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)

Figure 12.10: Rainfall data points and elevation raster

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.

Figure 12.11: gstat predict methods (Applied Spatial Data Analysis with R, 2013)

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:

library(gstat)
g = gstat(formula = annual ~ 1, data = rainfall)

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³⁸.

For example, in the following expression we create a formula object named f:

f = annual ~ 1
f
## annual ~ 1

Checking the class shows that f is indeed a formula object:

class(f)
## [1] "formula"

We can also convert character values to formula using the as.formula function. For example:

f = as.formula('annual ~ 1')
class(f)
## [1] "formula"

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:

z = predict(g, dem)
## [inverse distance weighted interpolation]

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':

z = z['var1.pred',,]
names(z) = '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)

Figure 12.12: Predicted annual rainfall using Inverse Distance Weighted (IDW) interpolation

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:

v_emp_ok = variogram(annual ~ 1, rainfall)

We can examine the variogram using plot (Figure 12.13):

plot(v_emp_ok)

Figure 12.13: Empirical variogram

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:

library(automap)
v_mod_ok = autofitVariogram(annual ~ 1, rainfall)

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):

plot(v_mod_ok)

Figure 12.14: Variogram model

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:

z = z['var1.pred',,]
names(z) = 'annual'

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)

Figure 12.15: Predicted annual rainfall using Ordinary Kriging

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:

z = predict(g, dem)
## [using ordinary kriging]
z = z['var1.pred',,]
names(z) = 'elevation'

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)

Figure 12.16: Ordinary Kriging prediction of elevation

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:

i = 'may'
as.formula(paste0(i, ' ~ 1'))
## may ~ 1

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:

m = c('sep', 'oct', 'nov', 'dec', 'jan', 'feb', 'mar', 'apr', 'may')
result = list()

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):

result$along = 3
result = do.call(c, result)

The interpolated montly rainfall amounts are shown in Figure 12.17:

plot(result, breaks = 'equal', col = hcl.colors(11, 'Spectral'), key.pos = 4)

Figure 12.17: Monthly rainfall predictions using Ordinary Kriging

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')

Figure 12.18: OK and UK variogram models

Next we create a gstat object, where the formula contains the covariate and the corresponding variogram model:

g = gstat(formula = annual ~ elev_1km, model = v_mod_uk$var_model, data = rainfall)

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:

z = predict(g, dem)
## [using universal kriging]

and then subset and rename:

z = z['var1.pred',,]
names(z) = 'annual'

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)

Figure 12.19: Predicted annual rainfall using Universal Kriging

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:

cv = gstat.cv(g)

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):

library(sp)
bubble(as(cv[, 'residual'], 'Spatial'))

Figure 12.20: Cross-validation residuals

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:

sqrt(sum((cv$var1.pred - cv$observed)^2) / nrow(cv))
## [1] 36.17293

---

38. The ~ 1 part means there are no independent variables other than the intercept.↩︎