The nearest neighbor interpolation (NNI) is a very simple 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\) that is 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 poylgon, 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)

It is fairly straight forward to apply this approach on a real data set. 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)

## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/johan/AppData/Local/R/cache/R/renv/cache/v5/R-4.2/x86_64-w64-mingw32/rgdal/1.6-7/10b777236c9e7855bc9dea8e347e30b7/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/johan/AppData/Local/R/cache/R/renv/cache/v5/R-4.2/x86_64-w64-mingw32/rgdal/1.6-7/10b777236c9e7855bc9dea8e347e30b7/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

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