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.