The representation of three-dimensional data by contouring or surface fitting is a general task with applications well beyond the familiar depiction of elevations on a topographic map.  Contouring involves two basic tasks:    1) finding values of the variable being contoured at target points--usually on  a rectangular grid--between control points where the values are known, 2) displaying the results in a contour plot, levelplot, perspective or wireframe plot.  

Interpolation

The general problems and solutions that are involved in three-dimensional interpolation can be illustrated by some two-dimensional examples.  Linear interpolation involves the fitting what can be thought of a as straight line-segments between the control points, while spline interpolation involves the fitting of a smooth curve that passes through the individual control points.

There are some particular configurations of the control points that produce results that are less than desirable, and the choice of an appropriate interpolation method is often an overlooked task in contouring.  A couple of examples illustrate situations when basic interpolation procedures do not work well.

• two-dimensional interpolation

Three-dimensional interpolation

Three-dimensional interpolation (the interpolation of a third variable (e.g. "Z") as a function of a two-dimensional set of control points, (X's and Y's) can be considered as a basic elaboration of simple interpolation.  There are several potential approaches

• distance-weighted averaging

Distance-weighted averaging involves the calculation of weighted averages of Z at the target points, where the weights are some function of distance between a target point and nearby control points.  The number of control points to use is an issue; more points will yield a smoother surface, but one that may not "respect" the values of Z at the control points, unless specific "fixups" are done.

• fitting a plane

Another approach involves fitting a plane to three points, and evaluating the fitted function at the control points contained within the triangle formed by the points.  A plane can be fit in R using the lm (linear model) function, the one most often used in regression analysis.  One problem that arises in fitting planes (not illustrated by this example, however) is the sometimes abrupt changes in interpolated values occur as one passes from one triangular region to another.

• fitting a polynomial

Higher-order (i.e. polynomial) functions can also be fit to a small number of control points and then evaluated at the target points,  Higher-order models require more control points, and as is the case with lower-order models, at times unusual interpolated values may occur as the control points change across the region under consideration.

In applying one of these approaches in practice, it is necessary to identify the closest control-point neighbors to the target points, or to identify which target points lie in, for example, the triangular regions formed by groups of three control points.  This task can be solved by defining the "Delaunay triangulation" of a set of points (also known as Thiessen polygons).  These definition of the triangulation, and the resulting contours can be illustrated with some simple examples.  Delaunay triangulations of points can be done using the deldir function in the package deldir.

• Delaunay triangulation contour1.gif      
• contouring using a TIN contour2.gif
   
• 3-d interpolation examples

Drawing the contour lines

One approach for drawing individual contour lines (or isolines) involves inspecting the grid of interpolated values, and then determining where along the line segments connecting individual grid points a particular contour would cross.  This is ordinarily done using linear interpolation between the grid-box corners.  Once those grid-cell-edge crossing points have been identified for each cell, then line segments joining the crossing points on the edges of each cell can be drawn, which together form a contour line (made up of individual short segments).  The general smoothness of contour lines drawn in this way depends on the  resolution of the grid:  higher-resolution grids yield smoother contours.  Alternatively, the contour lines developed using a coarser grid can be smoothed to yield smoother contours.  This line-drawing approach can be illustrated using a couple of simplified cases:

• identifying the contour-crossing points contour3.gif
• drawing the contours contour4.gif

An example

R uses the function interp to interpolate data from an irregular target point array onto a rectangular grid.  The specific approach used is called Akima's method, that involves fitting quintic (fifth-order) polynomials to sets of points, with the polynomials constrained to do nice things (e.g. to be relatively smooth).  Note that in this "traditional graphics" approach in R and S-Plus, the "object" ultimately displayed consists of a matrix (z) and two vectors x and y of the X and Y coordinates.

The Trellis Graphics system in S-Plus allows more detailed control of graphic design, ultimately allowing the particular features of individual panels of the "trellises" of plots to specified.  The specification of data in the Trellis contouring and surface-depiction functions is different from that in the traditional graphics--in the Trellis plots, the data to be plotted consist of columns in a data frame, and the particular variables to be plotted are specified in an model-description equation (e.g. elev ~ longitude * latitude).

• contouring example

Other approaches

• Loess surfaces

A locally determined surface can be constructed using loess, first fitting a model that illustrates the response of the variable being contoured as a function of two (or more) location variables, and then using the predict() function to visualize the resulting surface:

• Trend surfaces

Trend surfaces are developed by fitting low-order polynomial functions of the location variables, and then evaluating the resulting model over a grid of points.  The smoothness of the resulting surface is controlled by the order of the polynomial.  This analysis can be conveniently done using the spatial library of Venebles and Ripley.

• Kriging

Kriging is related to trend surface analysis except (as Venebles and Ripley show), an explicit "covariance structure" of the residuals of the model is included in the model specification in the surf.gls() function, where the form of the covariance structure is determined by inspecting the spatial correlogram and variogram of the variable being contoured.

• examples of other approaches