Exercise 29b--Interpolating a surface

Grids are created in many ways:

Directly, through remote sensing devices
By scanning printed images (this is a form of "remote sensing")
By interpolating discretely measured values
Through computations based on other grids
As arrays of model results

This section discusses interpolation from point data.  (Interpolation from contour data or other non-point data is performed differently.)  Point data usually arise as measurements or observations yielding numerical results at specific locations.  Examples abound: temperature measurements at weather stations, chemical concentrations measured in soil samples, river levels measured at gauging stations.

Interpolation is the process of estimating a value at a location (x, y) from assumed or measured values at other locations (x1, y1), ..., (xN, yN).  Interpolation can be used to create a grid simply by estimating values at the centers of each grid cell.  So, to a great extent, understanding how an interpolated grid is created depends on understanding how interpolation is performed at a single point.

 

One way to interpolate is to average all the measured values.  This makes good sense when the measurements are really independent of location.  Normally, however, a variable defined continuously over part of the earth's surface does not change wildly from point to point.  We can reasonably expect the values at unmeasured points to be related to the values at nearby measured points.  Local interpolators exploit this idea by using only the values at points near (x, y) for estimates.  These interpolators vary according to how the nearby points (the "neighborhood") are selected.

Neighborhood search methods

Two common approaches are (1) select all points within a fixed distance of the estimated point or (2) select the nearest K points, where K is predetermined.  In the next figure, the fixed-distance technique (using a distance equal to the radius of the circles) would interpolate the value at point A using four neighbors but would use only one neighbor to interpolate the value at point B.

If, instead, a fixed number of neighbors were used, then the search around point B would have to expand until four neighbors are found:

Each method has its advantages and drawbacks:

The fixed distance approach may fail to find any neighbors.  It will be unable to interpolate into large spatial gaps in the data support.  (The support is the set of points where measurements are available.)  In some circumstances this is good, because it helps reveal areas of undesirably low sampling; in other circumstances, this is bad because the holes in the grid may be undesirable.

The K-nearest-neighbors approach will always find neighbors, but they may be so far away from some points (as in B above) as to be essentially useless.  The interpolator may be giving misleading information.

Clearly, which method to use will depend on the data, their method of collection, the desired grid characteristics, and the nature of the decisions or analyses that will be performed with the resulting grid.

Interpolation methods

Once a set of neighbors is found, the interpolator must combine their values to produce the estimate.  The simplest combinations, such as the average or median, are often reasonable, but they have several undesirable properties:

Any statistical combination of the neighboring data perforce depends only on the neighbors, which implies the estimated values will change only when the collection of neighbors changes.

The values of the statistical combinations do not change except when the neighborhood changes.

This leads to discontinuous, step-like surfaces:

This four-nearest-neighbor average picked up a point with a value of 410 and dropped the point with value 270 when moving from A to B.  Therefore the estimated value at B is (410 - 270)/4 = 35 greater than the estimated value at A.  The change occurred exactly at the location, roughly midway between A and B, where point 410 became one of the four closest points.  The color change in the interpolated grid is apparent there.

IDW

For these reasons, analysts usually prefer to weight the averages according to the relationships among the points.  The simplest form of weighting is Inverse Distance Weighting, or IDW.  This technique creates weights according to the distances between the interpolated location (x, y) and each of its neighbors.  To be reasonable, the weight functions should be largest at zero distance and decrease as the distance increases.

Power functions provide an easy-to-compute set of weight functions.  The values of y = x-p, for some fixed value p, decrease as x increases provided p > 0.  (The minus sign in the equation was chosen to make p positive.)  Thus, IDW works in four simple steps:

  1. Find a neighborhood {(x1, y1), ..., (xN, yN)} of points around (x, y), the location where the value is to be estimated.
  2. Compute the relative weight wK of each point (xK, yK) in the neighborhood as

wk = distance((x, y), (xK, yK))-p = [(x-xK)2 + (y-yK)2](-p/2)

  1. Divide each relative weight by the sum of all the weights: this causes the weights to sum to unity.
  2. Form the weighted average of the known or assumed values at each point (xK, yK).  This is simply

z = Sum of wk * zk

where zk is the value at point (xK, yK).

It is helpful to acquire some intuition about how IDW works.  Some general observations follow immediately from the fact that IDW estimates are averages:

Adding a constant to all known values causes all interpolated values to increase by the same constant.

Multiplying all known values by a constant causes all interpolated values to multiply by the same constant.

(In other words, IDW is a linear estimator.)   Because the IDW weights are relative weights (step 3),

Changing the distance units of measurement does not change the interpolated values.

(In other words, IDW is a scale-invariant estimator.)

The simplest case consists of two data points A and B.  By virtue of the properties above, we may rescale the known values to be equal to 0.0 at A and 1.0 at B.  Consider the IDW estimates along the line joining A and B: they must rise from a value of 0.0 at A (0 percent of the way from A to B) to a value of 1.0 at B (100 percent of the way from A to B).  The following graph profiles these estimated values:

Point A is located at X = 0.0 on this graph and B is located at X = 1.0.  Points directly between A and B correspond to distances between 0.0 and 1.0.  Points "behind" A along the same line correspond to negative distances.  The different curves correspond to different values of the power p as indicated in the legend.

Many of the evident conclusions can be proven mathematically for the general situation of many data points:

For powers p > 1.0, IDW surfaces are smooth (have first derivatives) at the data locations.

For powers p <= 1.0, IDW surfaces have cusps (discontinuous or infinite first derivatives) at the data locations.

As p >= 1.0 increases, the IDW surface becomes flatter near the data points and jumps more rapidly in value between the data points.  To put it another way, the IDW estimator acts as if data have stronger spatial correlation as the value of p increases.  (At the lower limit of p = 0.0, the IDW estimator reduces to the usual unweighted average.)

Beyond the spatial extent of the data points, IDW values asymptotically approach the (usual, unweighted) average.

For these reasons, analysts usually choose powers greater than 1.0, with 2.0 being a common default.

Exercise

  1. Download the zip file containing the data and ArcView project for this exercise.  Unzip the project and data into a folder of your choice.  Open the ArcView project.

  2. You are going to create a gridded surface from the point theme of fertilizer requirements, [Fertilizer.shp].  The first thing is to tell Spatial Analyst about the grid you intend to produce: its spatial extent and cell size.  The Analysis|Properties menu item launches the dialog for this.  While you are experimenting it's a good idea to keep grid sizes small.  Remember, an interpolation method has to create an estimate at every grid cell.  Therefore, if you halve the cell size, the number of cells quadruples.  This law of diminishing returns suggests you run software trials with a coarse grid, and then use a fine grid only when you are certain of the results.

To set the spatial extent of the grid, first adjust the View display until it just contains the area you would like to grid.  Then specify "Same as Display" in the Analysis Extent box.  If you like, you may then manually adjust the "Left", "Top", "Right", and "Bottom" values so they are nice round numbers.

A good starting point for the cell size in this data set is 0.01 map units.  After typing in the value, press enter (while the cursor is still in the Cell Size text box, as illustrated) to see how many rows and columns Spatial Analyst will create.  Their product--about 21,000 here--is the number of cells and determines the computation time.  A modest PC (200 MHz) will have little trouble quickly producing grids of 100,000 to 1,000,000 cells with most Spatial Analyst operations.

  1. With the Analysis Properties set, activate the theme containing the data.  Press the Clear Selected Features button to make sure all the data will be used: interpolation, like most ArcView operations, will use only the selected features if any are selected.  Then activate the Interpolate Surface dialog by selecting Surface|Interpolate grid from the View's menu.

  2. The Interpolate Surface dialog is filled out top to bottom.  Read each item carefully.  It is easy to overlook an item and then waste your time calculating something you will not want.  The field containing the data is called [Fertilizer]; its alias is [Fertilizer Requirement].  (Notice that this dialog does not use the alias--you have to know the actual field name.)  Ignore the "Barriers" item for now.

This figure shows the default options for interpolation: IDW-2 using the 12 nearest neighbors.  You will have to specify the Z-value field.

  1. Press OK.  After the computation is complete, the grid will appear at the top of the View's Table of Contents.  Move it beneath the data theme so it does not obscure everything.  Edit the legend to suit your needs, then make the theme visible.

  1. Experiment with the interpolation options.  After each computation, take a few seconds to give the new grid an appropriate name.  Otherwise, you will lose track of how each grid was computed.  Consider typing some information into the grid theme's comments box.

  1. Note that the grids you are creating are marked "Temporary."  They are placed in your $HOME directory, often C:/Temp.  ArcView will delete them when it thinks you no longer want them--such as when you close the project.  Therefore, when you want to keep a grid, use the Theme|Save data set menu item to specify a permanent location.