Natural neighbor interpolation creates weights for each of the input points based on their assumed “area of influence”. These areas are determined by the generation of Voronoi polygons around each input point. In principle every grid intersection created would be in one of these polygons and could be assigned the value of the point around which the polygon has been created. This would result in a step-like surface of patches. This is the kind of result that is obtained from Nearest Neighbor interpolation. A far more effective approach involves a development of this idea, as described below. The end result is a smooth surface with discontinuities at the input points (Figure 6‑31D).

The first step in the process is to create a Delaunay triangulation of the j=1, 2…62 input data points (Figure 6‑36A) as a preliminary stage in the creation of Voronoi polygons.

The second stage is to generate a set of Voronoi polygons for the study region (Figure 6‑36B). Each of the points, j, in the source dataset has its own Voronoi polygon, which has an area Aj. In order to determine the estimated value at a sample point P the point P is temporarily added to the set (so there are now 63 points) and the Voronoi polygons are re-computed. Adding point P results in a new Voronoi polygon and redefinition of those immediately surrounding it (Figure 6‑36C). This new polygon has an area Ap.

Effectively this new point has “borrowed” some of the area of influence from each of the nearby points. This can be seen in Figure 6‑36D, where the new region has been overlaid on the original set. There are k=5 original polygons that the new polygon has borrowed area from. Let us call these borrowed areas Aip, i=1,…k then the total area of P’s Voronoi polygon is:

and thus the proportion borrowed from each of the original points is:

These proportions are the weights used to compute the estimated value at P, based on the standard linear weighting equation:

If the point P coincides with one of the existing points its area of overlap with that point would be 100%, hence its weight would be 1, as required. If the Voronoi polygon for P does not overlap a region the weight associated with that region is 0. One of the main advantages of this method of interpolation is that it requires no decision-making regarding the number of points to use, the radius or direction of search, or any other parameters. It is a surprisingly effective and straightforward technique that is widely used in many disciplines other than GIS.

The standard implementation of this method provides interpolated values for the convex hull of the input point set. However, one or two software packages extend this to enable limited extrapolation to the boundary of a rectangular region (the MBR). This is achieved by using a slightly enlarged MBR (e.g. increased by 10%). Within this enlarged region temporary points are added to the original set at each corner (and/or along the region boundary). z-values for these temporary points are estimated by simple IDW or a similar technique, and then the Voronoi regions computed for all points. The procedure continues as before, but on completion the temporary points are discarded and the MBR re-drawn based on the original point set.

Figure 6‑36 Natural Neighbor interpolation ― computation of weights

A. Delaunay triangulation of source point set |
B. Voronoi polygons of source point set |

C. Revised Voronoi polygons with additional point, P |
D. Borrowed area from original polygons |