UNIT 40 - SPATIAL INTERPOLATION I

• 1. Point Interpolation/Areal Interpolation
• 2. Global/Local Interpolators
• 3. Exact/Approximate Interpolators
• 4. Stochastic/Deterministic Interpolators
• 5. Gradual/Abrupt Interpolators

• 1. Proximal
• 2. B-splines
• 3. Kriging
• Variograms
• Deriving the variogram
• Computing the estimates
• 4. Manual Interpolation or "eyeballing"

• 1. Trend Surface Analysis
• 2. Fourier Series
• 3. Moving average/distance weighted average

UNIT 40 - SPATIAL INTERPOLATION I

Compiled with assistance from Nigel M. Waters, University of Calgary

A. INTRODUCTION

• spatial interpolation is the procedure of estimating the value of properties at unsampled sites within the area covered by existing observations
  • in almost all cases the property must be interval or ratio scaled

• can be thought of as the reverse of the process used to select the few points from a DEM which accurately represent the surface

• rationale behind spatial interpolation is the observation that points close together in space are more likely to have similar values than points far apart (Tobler's Law of Geography)

• spatial interpolation is a very important feature of many GISs

• spatial interpolation may be used in GISs:
  • to provide contours for displaying data graphically
  • to calculate some property of the surface at a given point
  • to change the unit of comparison when using different data structures in different layers
  • frequently is used as an aid in the spatial decision making process both in physical and human geography and in related disciplines such as mineral prospecting and hydrocarbon exploration

• many of the techniques of spatial interpolation are two- dimensional developments of the one dimensional methods originally developed for time series analysis

• this unit introduces spatial interpolation and examines point based interpolation, while the next looks at areal procedures and some applications

• there are several different ways to classify spatial interpolation procedures:

• point based
  • given a number of points whose locations and values are known, determine the values of other points at predetermined locations

  diagram
  • point interpolation is used for data which can be collected at point locations e.g. weather station readings, spot heights, oil well readings, porosity measurements
  • interpolated grid points are often used as the data input to computer contouring algorithms
    • once the grid of points has been determined, isolines (e.g. contours) can be threaded between them using a linear interpolation on the straight line between each pair of grid points
  • point to point interpolation is the most frequently performed type of spatial interpolation done in GIS

• lines to points
  • e.g. contours to elevation grids

• areal interpolation
  • given a set of data mapped on one set of source zones determine the values of the data for a different set of target zones
    • e.g. given population counts for census tracts, estimate populations for electoral districts

  diagram

• global interpolators determine a single function which is mapped across the whole region
  • a change in one input value affects the entire map

• local interpolators apply an algorithm repeatedly to a small portion of the total set of points
  • a change in an input value only affects the result within the window

• global algorithms tend to produce smoother surfaces with less abrupt changes
  • are used when there is an hypothesis about the form of the surface, e.g. a trend

• some local interpolators may be extended to include a large proportion of the data points in set, thus making them in a sense global

• the distinction between global and local interpolators is thus a continuum and not a dichotomy
  • this has led to some confusion and controversy in the literature

• exact interpolators honor the data points upon which the interpolation is based
  • the surface passes through all points whose values are known
  • honoring data points is seen as an important feature in many applications e.g. the oil industry

  • proximal interpolators, B-splines and Kriging methods all honor the given data points
    • Kriging, as discussed below, may incorporate a nugget effect and if this is the case the concept of an exact interpolator ceases to be appropriate

• approximate interpolators are used when there is some uncertainty about the given surface values
  • this utilizes the belief that in many data sets there are global trends, which vary slowly, overlain by local fluctuations, which vary rapidly and produce uncertainty (error) in the recorded values
  • the effect of smoothing will therefore be to reduce the effects of error on the resulting surface

• stochastic methods incorporate the concept of randomness
  • the interpolated surface is conceptualized as one of many that might have been observed, all of which could have produced the known data points

• stochastic interpolators include trend surface analysis, Fourier analysis and Kriging
  • procedures such as trend surface analysis allow the statistical significance of the surface and uncertainty of the predicted values to be calculated

• deterministic methods do not use probability theory (e.g. proximal)

• a typical example of a gradual interpolater is the distance weighted moving average
  • usually produces an interpolated surface with gradual changes
  • however, if the number of points used in the moving average is reduced to a small number, or even one, there would be abrupt changes in the surface

• it may be necessary to include barriers in the interpolation process
  • semipermeable, e.g. weather fronts
    • will produce quickly changing but continuous values
  • impermeable barriers, e.g. geologic faults
    • will produce abrupt changes

• Lam (1983) and Burrough (1986) describe a variety of quantitative interpolation methods suitable for computer contouring algorithms
  • in this and the next sections, these are divided into exact and approximate methods
  • this section deals with exact methods

• all values are assumed to be equal to the nearest known point
  • is a local interpolator
  • computing load is relatively light

• output data structure is Thiessen polygons with abrupt changes at boundaries
  • has ecological applications such as territories and influence zones

• best for nominal data although originally used by Thiessen for computing areal estimates from rainfall data

• is absolutely robust, always produces a result, but has no "intelligence" about the system being analyzed

• available in very few mapping packages, SYMAP is a notable exception

• uses a piecewise polynomial to provide a series of patches resulting in a surface that has continuous first and second derivatives
  • ensures continuity in:
    • elevation (zero-order continuity) - surface has no cliffs
    • slope (first-order continuity) - slopes do not change abruptly, there are no kinks in contours
    • curvature (second order continuity) - minimum curvature is achieved
  • produces a continuous surface with minimum curvature

• output data structure is points on a raster
  • note that maxima and minima do not necessarily occur at the data points

• is a local interpolator
  • can be exact or used to smooth surfaces
  • computing load is moderate

• best for very smooth surfaces
  • poor for surfaces which show marked fluctuations, this can cause wild oscillations in the spline

• are popular in general surface interpolation packages but are not common in GISs

• can be approximated by smoothing contours drawn through a TIN model

• see Burrough (1986), Davis (1986) and mathematical aspects in Lam (1983) and Hearn and Baker (1986)
  • also described in "numerical approximation theory"

• developed by Georges Matheron, as the "theory of regionalized variables", and D.G. Krige as an optimal method of interpolation for use in the mining industry

• the basis of this technique is the rate at which the variance between points changes over space
  • this is expressed in the variogram which shows how the average difference between values at points changes with distance between points

diagram
• De (vertical axis) is E(zi - zj)2, i.e. "expectation" of the difference
  • i.e. the average difference in elevation of any two points distance d apart
• d (horizontal axis) is distance between i and j

• most variograms show behavior like the diagram
  • the upper limit (asymptote) of De is called the sill
  • the distance at which this limit is reached is called the range
  • the intersection with the y axis is called the nugget
  • a non-zero nugget indicates that repeated measurements at the same point yield different values

• in developing the variogram it is necessary to make some assumptions about the nature of the observed variation on the surface:
  • simple Kriging assumes that the surface has a constant mean, no underlying trend and that all variation is statistical

  • universal Kriging assumes that there is a deterministic trend in the surface that underlies the statistical variation

• in either case, once trends have been accounted for (or assumed not to exist), all other variation is assumed to be a function of distance

• the input data for Kriging is usually an irregularly spaced sample of points

• to compute a variogram we need to determine how variance increases with distance

• begin by dividing the range of distance into a set of discrete intervals, e.g. 10 intervals between distance 0 and the maximum distance in the study area

• for every pair of points, compute distance and the squared difference in z values

• assign each pair to one of the distance ranges, and accumulate total variance in each range

• after every pair has been used (or a sample of pairs in a large dataset) compute the average variance in each distance range

• plot this value at the midpoint distance of each range

• once the variogram has been developed, it is used to estimate distance weights for interpolation
  • interpolated values are the sum of the weighted values of some number of known points where weights depend on the distance between the interpolated and known points

• weights are selected so that the estimates are:
  • unbiased (if used repeatedly, Kriging would give the correct result on average)
  • minimum variance (variation between repeated estimates is minimum

• problems with this method:
  • when the number of data points is large this technique is computationally very intensive
  • the estimation of the variogram is not simple, no one technique is best
  • since there are several crucial assumptions that must be made about the statistical nature of the

    variation, results from this technique can never be absolute

• simple Kriging routines are available in the Surface II package (Kansas Geological Survey) and Surfer (Golden Software), and in the GEOEAS package for the PC developed by the US Environmental Protection Agency

• traditionally not a highly regarded method among geographers and cartographers

• however, Dutton-Marion (1988) has shown that among geologists this is a very important procedure and that most geologists actually distrust the more sophisticated, mathematical algorithms
  • they feel that they can use their expert knowledge, modelling capabilities and experience and generate a more realistic interpolation by integrating this knowledge into the construction of the geological surface

• attempts are now being made to use knowledge engineering techniques to extract this knowledge from experts and build it into an expert system for interpolation
  • see Unit 74 for more on this topic

• characteristics of this method include:
  • procedures are local as different methods may be used by the expert on different parts of the map
  • tend to honor data points
  • abrupt changes such as faults are more easily modelled using these methods

• the surfaces are subjective and vary from expert to expert

• output data structure is usually in the form of a contour

1. Trend Surface Analysis

• surface is approximated by a polynomial

• output data structure is a polynomial function which can be used to estimate values of grid points on a raster or the value at any location

• the elevation z at any point (x,y) on the surface is given by an equation in powers of x and y
  • e.g. a linear equation (degree 1) describes a tilted plane surface:

    z = a + bx + cy

  • e.g. a quadratic equation (degree 2) describes a simple hill or valley:

    z = a + bx + cy + dx2 + exy + fy2

• in general, any cross-section of a surface of degree n can have at most n-1 alternating maxima and minima
  • e.g. a cubic surface can have one maximum and one minimum in any cross-section
  • equation for the cubic surface:

    z = a + bx + cy + dx2 + exy + fy2 + gx3 + hx2y + ixy2 + jy3

• a trend surface is a global interpolator
  • assumes the general trend of the surface is independent of random errors found at each sampled point

• computing load is relatively light

• problems
  • statistical assumptions of the model are rarely met in practice
  • edge effects may be severe
  • a polynomial model produces a rounded surface
    • this is rarely the case in many human and physical applications

• available in a great many mapping packages
  • see Davis (1973) and Sampson (1978) for non- orthogonal polynomials; Mather (1976) for orthogonal polynomials

• approximates the surface by overlaying a series of sine and cosine waves

• a global interpolator
  • computing load is moderate

• output data structure is the Fourier series which can be used to estimate grid values for a raster or at any point

• best for data sets which exhibit marked periodicity, such as ocean waves

• rarely incorporated in computing packages
  • simple program and discussion in Davis (1973)

• estimates are averages of the values at n known points:

  z = S wizi/ S wi

  where w is some function of distance, such as:

  w = 1/dk w = e-kd

• an almost infinite variety of algorithms may be used, variations include:
  • the nature of the distance function
  • varying the number of points used
  • the direction from which they are selected

• is the most widely used method

• objections to this method arise from the fact that the range of interpolated values is limited by the range of the data
  • no interpolated value will be outside the observed range of z values

  diagram

• other problems include:
  • how many points should be included in the averaging?
  • what to do about irregularly spaced points?
  • how to deal with edge effects?

Burrough, P.A., 1986. Principles of Geographical Information Systems for land Resources Assessment, Clarendon, Oxford. See Chapter 8.

Davis, J.C., 1986. Statistics and Data Analysis in Geology, 2nd edition, Wiley, New York. (Also see the first, 1973, edition for program listings.)

Dutton-Marion, K.E., 1988. Principles of Interpolation Procedures in the Display and Analysis of Spatial Data: A Comparative Analysis of Conceptual and Computer Contouring, unpublished Ph.D. Thesis, Department of Geography, University of Calgary, Calgary, Alberta.

Hearn, D., and Baker, M.P., 1986. Computer Graphics, Prentice-Hall Inc, Englewood Cliffs, N.J.

Jones, T.A., Hamilton, D.E. and Johnson, C.R., 1986. Contouring Geologic Surfaces with the Computer, Van Nostrand Reinhold, New York

Lam, N., 1983. "Spatial Interpolation Methods: A Review," The American Cartographer 10(2):129-149.

Mather, P.M., 1976. Computational Methods of Multivariate Analysis in Physical Geography, Wigley, New York.

Sampson, R.J., 1978. Surface II, revised edition, Kansas Geological Survey, Lawrence, Kansas.

Waters, N.M., 1988. "Expert Systems and Systems of Experts," Chapter 12 in W.J. Coffey, ed., Geographical Systems and Systems of Geography: Essays in Honour of William Warntz, Department of Geography, University of Western Ontario, London, Ontario.

An important class of interpolation methods is missing here - so called radial basis functions, such as multiquadrics, thin plate spline, thin plate spline with tension, regularized spline with tension and a large number of other flavours of this approach (also sometimes refered to as variational approach). These methods are available in almost every GIS, from ArcINFO, GRASS, SURFER to specialized visualization packages. The description can be found at Mitas, L., Mitasova, H., 1999, Spatial Interpolation. In: P.Longley, M.F. Goodchild, D.J. Maguire, D.W.Rhind (Eds.), Geographical Information Systems: Principles, Techniques, Management and Applications, GeoInformation International, Wiley, 481-492.

DISCUSSION AND EXAM QUESTIONS

1. Are there other techniques for surface generation? How many of the above procedures are commonly used? How would they be ranked in terms of popularity? Give examples from the literature of where they have been used.

2. How does hand contouring rate as an alternative? What did you think of it and have you changed your mind? What are the key features and processes involved in hand contouring?

3. Explain the advantages and disadvantages of manual interpolation as used in hand contouring over computer based interpolation as used in a computer contouring package.

4. Describe the different ways in which spatial interpolation algorithms can be classified.

Back to Geography 470 Home Page

Back to GIS & Cartography Course Information Home Page