Compiled with assistance from Nigel M. Waters, University of
Calgary
NOTES
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
B. CLASSIFICATION OF INTERPOLATION PROCEDURES
- there are several different ways to classify spatial
interpolation procedures:
1. Point Interpolation/Areal Interpolation
- 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
2. Global/Local Interpolators
- 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
3. Exact/Approximate Interpolators
- 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
4. Stochastic/Deterministic Interpolators
- 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)
5. Gradual/Abrupt Interpolators
- 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
C. POINT BASED INTERPOLATION - EXACT METHODS
- 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
1. Proximal
- 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
2. B-splines
- 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"
3. Kriging
- 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
Variograms
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
Deriving the variogram
- 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
Computing the estimates
- 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:
- 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
4. Manual Interpolation or "eyeballing"
- 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
D. POINT BASED INTERPOLATION - APPROXIMATE METHODS
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
- in general, any cross-section of a surface of degree n
can have at most n-1 alternating maxima and minima
- 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
2. Fourier Series
- 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)
3. Moving average/distance weighted average
REFERENCES
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 370 Home Page
Back to Geography 470 Home Page
Back
to GIS & Cartography Course Information Home Page
Please send comments regarding content to: Brian
Klinkenberg
Please send comments regarding web-site problems to: The
Techmaster
Last Updated: August 30, 1997.