Introduction

Kriging (Matheron 1963) is a popular technique for interpolating and estimating elevation values from digital terrain data. General references on the subject include David (1977), Isaaks and Srivastava (1989), Goovaerts (1997), and Journal and Huijbregts (1997). Bailey (1994, p. 32) asserts that "there is an argument for kriging to be adopted as a basic method of surface interpolation in geographic information systems (GIS) as opposed to the standard deterministic tessellation techniques that currently prevail and which can produce artificially smooth surfaces." This argument was supported by Laslett (1994) whose study gives an example of a data set for which splines are "too smooth," and kriging results in more precise estimations. Regarding the "spline vs. kriging" debate, it has been shown (Kimeldorf and Wahba 1971; Wahba 1990) that kriging is mathematically equivalent to thin-plate splines, and Almansa et al. (2002) established its place in a more general class of functions, namely, the absolutely minimizing Lipschitz extension.

While kriging is not without its critics (Philip and Watson 1986) there is no question that its use is widespread. The properties of any mathematical surface being used as a terrain model define the properties imbued to the model. The onus is on the modeler to choose the model wisely so that its properties match the desired traits of the terrain. Continuity properties are of paramount importance. Discontinuous surfaces have "holes" or "tears" in them. Continuous surfaces might not be smooth, meaning that the surface might have one or more "creases" in it.

The computational geometry literature is replete with discussions about the continuity of patches within themselves and with their neighbors (de Boor 1978; Lancaster and Salkauskas 1986; Dierckx 1993; Farin 1993), but there is relatively little attention given to surface discontinuity properties in the GIS literature. Although Meyer (1999) gives a review of the continuity properties of digital terrain models categorizing them according to their topological and continuity characteristics, Lam (1983, p. 134) set the stage for the current research when she reported concerning kriging that, "Choice of neighborhood will also affect the continuity properties of the estimates ... If the change of data points from one neighborhood to the next is too abrupt there may be discontinuities even though the actual phenomenon is continuous." The purpose of this paper is to carefully document this undesirable characteristic of kriging, namely, that kriging as it is typically used for digital terrain modeling produces a piece-wise discontinuous surface.

This paper is organized as follows. A brief explanation of the notation and the theoretical setting of digital terrain modeling in the context of this paper is presented. This is followed by a proof that kriging can produce discontinuous surfaces along with a discussion that generalizes the results. This is followed by a section presenting an analysis of the discontinuities found in a USGS 7.5' digital elevation model, which provides some idea of the magnitude and distribution of the worst discontinuities encountered in this real-world data set. The final section of the paper provides a discussion, summary, and conclusions.

Theoretical Setting

Define a surface to be a real-valued function of two variables, z = f(x,y). Let f be defined over an open, bounded, simply connected, two-dimensional region R that is a subset of three-dimensional Euclidean real space. A digital terrain model (DTM) is often constructed from a set of samples, s, taken from the area of interest, such that these samples are intended to capture the essence of the terrain's shape. Digital terrain modeling includes the problem of interpolating s to derive heights at places in R for which no sample was taken. Some surfaces are defined upon all the sample points in s. Examples include Lagrange polynomials, Fourier transformations, and kriging. Such functions are said to have global support, meaning that every point in s contributes to the formulation of f.

Global support is generally not desirable for digital terrain modeling for several reasons. It imposes heavy computational burdens for large data sets. Also, it has the counter-intuitive property that, for certain techniques such as the Lagrange polynomials, making a small change in any particular sample can produce large changes over the entire surface. This runs contrary to Tobler's Law of Geography, "Everything is related to everything else, but near things are more related than distant things" (Tobler 1970, p. 236). Also, polynomial global support surface models which interpolate all the points in s must be of an order equal to the cardinality of s, or greater. This can produce unwanted behaviors such as extreme surface departures and unrealistic undulations. From a statistical point of view, global support leads to over-fitting and yields poor validation.

Although kriging is defined with global support (Siska et a1.1997), in practice it is not typically used that way for terrain modeling. Instead, s is subdivided into neighborhoods, which are subsets of s with relatively few elements that roughly (or strictly) partition s. Then, kriging is used to interpolate over the neighborhoods in the following way. Suppose I want to interpolate a surface value at the point p = (x,y) and p is not an element of s. Let n denote a neighborhood surrounding p. Then, the height estimate at p is a weighted sum of the heights of n. The weights are related...