The nearest neighbour interpolation (NNI) is a straightforward interpolation approach. Recall that interpolation seeks for a given sample of observations \(\{z_1, z_2,..., z_n\}\) at locations \(\{x_1, x_2,..., x_n\}\) to estimate the value \(z\) at some new point \(\mathbf x\).

The NNI approach uses the value \(z_i\) closest to \(\mathbf x\). Thus, the algorithm seeks to find \(i\) such that \(\vert \mathbf x_i - \mathbf x \vert\) is minimized, then the estimate of \(z\) is \(z_i\).

It should be noted that the set of closest points to \(\mathbf x_i\) forms a Thiessen polygon, also known as Voronoi polygon. The technique used to compute these polygons is referred to as Delauny triangulation. Such a particular polygon corresponds to a region consisting of all points closer to \(\mathbf x_i\) than to any other \(\mathbf x\).

In R, the dismo package, which relies on the deldir package, provides the voronoi() function to create Voronoi polygons for a set of points. The input points data set is either a SpatialPoints class object or a two-column matrix with x and y coordinates.

In order to showcase the function, we generate 25 random data points, build Thiessen polygons and plot them.

library(dismo)
set.seed(3) # set seed for reproducibility

n <- 25 # number of points
x <- runif(n, -1, 1)
y <- runif(n, -1, 1)
m <- cbind(x, y)

# plotting
plot(voronoi(m), main = "Thiessen poylgons for \nrandom data points")
points(m)

Applying this approach to a real data set is fairly straightforward. Let us compute the Voronoi polygons for the weather station data set dwd_east_sp, in particular for the mean annual temperature and the mean annual rainfall in East Germany.

We use dwd_east_sp, a spatial object of class SpatialPointsDataFrame, as input.

load(url("https://userpage.fu-berlin.de/soga/data/r-data/East_Germany.RData"))
voronoi <- voronoi(dwd_east_sp)
voronoi

## class       : SpatialPolygonsDataFrame 
## features    : 71 
## extent      : 140239.2, 528937.2, 5541059, 6107404  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs 
## variables   : 7
## names       :      X, Station_id, Altitude,                 Name, Federal_State, Rainfall, Temperature 
## min values  :   3089,        164,        2,          Angermuende,        Berlin,      317,         3.5 
## max values  : 111897,       5825,     1213, Zinnwald-Georgenfeld,    Thueringen,     1879,        10.5

Once we computed the Voronoi polygons we may plot them by using the spplot() from the sp package.

library(gridExtra)
p1 <- spplot(voronoi, "Temperature",
  main = "Mean annual temperature [°C]"
) +
  latticeExtra::layer(sp.polygons(east_germany_states_sp,
    col = "white",
    lwd = 1.5
  )) +
  latticeExtra::layer(sp.points(dwd_east_sp,
    pch = 1,
    col = "black"
  ))

plot(p1)

Exercise: Plot the Voronoi polygons of the mean annual rainfall data for East Germany. When you look at the plots, which aspect of the Voronoi polygons seems to be problematic?

## Your code here...
p2 <- NULL

Show code

p2 <- spplot(voronoi, "Rainfall",
  main = "Mean annual rainfall [mm]"
) +
  latticeExtra::layer(sp.polygons(east_germany_states_sp,
    col = "white",
    lwd = 1.5
  )) +
  latticeExtra::layer(sp.points(dwd_east_sp,
    pch = 1,
    col = "black"
  ))

p2

Looking at the plots, it is obvious that a problematic aspect of this approach is that the interpolated surfaces are discontinuous at their borders. These discontinuities are artefacts of the locations of the samples (Brundson and Comber 2015).

---

Citation

The E-Learning project SOGA-R was developed at the Department of Earth Sciences by Kai Hartmann, Joachim Krois and Annette Rudolph. You can reach us via mail by soga[at]zedat.fu-berlin.de.

You may use this project freely under the Creative Commons Attribution-ShareAlike 4.0 International License.

Please cite as follow: Hartmann, K., Krois, J., Rudolph, A. (2023): Statistics and Geodata Analysis using R (SOGA-R). Department of Earth Sciences, Freie Universitaet Berlin.