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
```

```
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.*