Automated Modelling of Digital Elevation Models for
Predictive Ecosystem Mapping in GIS.
By
Nancy Alexander
BSc, University of Calgary, 1993

THESIS SUBMITTED IN PARTIAL FULFILMENT OF
THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE
IN
NATURAL RESOURCES MANAGEMENT

© Nancy Alexander, 2001
UNIVERSITY OF NORTHERN BRITISH COLUMBIA
June,2001
All rights reserved. This work may not be
reproduced in whole or part, by photocopy
or other means, without the permission of the author.
BRITISH COlUMBIA
LI Br.~ ARY

Prl ce Gcorj'~, BC

Table of Contents
Table of Contents ...... .......... ................... ........... ...... .... ..... ........ .. ... ...................................................... ... ...... .... ....... i
Table ofFigures ................ .. ...................... ......... ...... ................ .... .... ......... ..... ........ ............................................... . ii
Table of Tables ......... .... ................. ... ................................. .......................... ....... ................................. .... ............. iii
Abstract ......... ... ............. .. ..................... ..... ... ........... ... ...... .... .. ... ... ........... ......................... .. .. .. ........ ... ......... .... ....... iv
Chapter 1 Introduction and Rationale .... .. .. .... .... .. ........... .. ...... ............. ..... ............... ... .. .. ..... ....... ...... .. .. .... .... ... .. ...... 1
1.1 Introduction.. .. .. .... ... .... .................... ................... .... .. ......... ... .. .... ........... .... .... .... ..... .... ... ... ............ ..... .. .......... 1
1.2 TRIM Elevation Data ..... .... .......... .. ..... ..... .................. .......................... ............................ .. .. .. ... .... ........ ... ..... 3
1.3 History ........................................................................ ...... ............ ... ... ....... .. ........... ......... .... ......... ...... ... ....... 4
1.4 Statement of the Problem and Objectives ................. ........................... .. ................... .. ............... ............... ...... 5
1.5 Thesis Organization ......................... ... .. ....... .... .... ........... ..................... ...................... .. .. ... ............. ... ... .. ..... ... 5
Chapter 2 Literature Review .. ... .... ........ ............................................ ... ...... ..... ... ................... ... ......... .... .... ......... ..... .?
2.1 Classification Paradigms ......................... ... ..... .. ....................... ................ .......... .. .. ..... ......................... .... ..... 7
2. 2 Parameterization ofFeatures ...... .................................. ....... ....... ..... .. ..................... .... .... ..... .. .. .. ... .. ... ............. 9
2.2.1 Slope Gradient and Aspect- The First Derivatives of Elevation ...... .... ............ ........................... .... ...... ... 9
2.2.2 Curvature -The Second Derivative of Elevation ................................ .. .......... .. ....................................... 9
2.2.3 Slope Normal ........................................................................... .. ............. ............ ............ .. .. ........ .. .. ...... 10
2.3 Operational Definitions: Case Studies Classifying Topographic Elements .................................................... 10
2.4 Data Model Selection ....................................................................................................... ..... .. .... ... ............. 21
2.4.1 Regularly-Gridded Elevation Values: The Digital Elevation Model (DEM) ............................................ 21
2.4.2 Triangular Irregular Networks (TIN) ..................................................................................................... 22
2.5 Accuracy Assessment ............................ ... ................. ...... ....... ...... ........ ..................................... .. .... .. .......... 24
2.5.1 Classification Accuracy Assessment. ...... .. .... ........ ............................. .... ..... .... ... ......... ... ......... .... ........... 24
2.5.2 DEM Accuracy Assessment ..................... .... .... .... ...... ... .......... .... .... ....... .. ...... .... ..... ... ..... ... ........ ... ........ 26
2.5.3 Data Filtering for Error Reduction .............................. ...... ... ............................. .... .. ............... .... ............ 30
2.6 Interpolation ................ ... ... ............ .... .................................... ......................... .. ........................... .... ............ 32
2.6.1 Deterministic Models .... ..... ..... ............. ........................... ... .......................................... .. ... .................... 3 3
Linear Interpolation and Inverse Distance Weighting ........ .... ..................................................................... 33
Thin Plate Splines ........... ............................ .......... .. ........ ....... .... .... ... .. ........ ................ .... .............. .. ........... 35
2.6.2 Optimal Interpolation ..................... ..... ..................... .. ..... ...... ... .. ............. ..... ... ... .... ....... ...... .. .... ....... ..... 36
Kriging........................................................................ ...... ............. .... .. .. ........................... .. .. .... ... ............. 36
2. 7 Interpolation Implementation and Accuracy Assessment .............................. .. ..................... ........ ................ .40
Chapter 3 Methods ............ ..... .................... ........... ..... ..... ... ............................... .... ................ .... ...................... ...... 43
3.0 Introduction .... ..... ............................ .. ....... ..................... ... .... ... ......... ......... ...... ................... ...... ..... .. .... .. ...... 43
3.1 Study Area ....................................... ................ .... .................... ....... .. ...... ... ... ..... .. .. ................ .... .. ... .... ... ..... 44
3.2 DEM Generation and Accuracy Assessment ......................................... .......... ...... .. .. .. .............. .. ................ .44
3.3 Elevation Data Preparation: Blunder Assessment and Control Point Removal... .................... .......... .. ............ 45
3.4 Interpolation Techniques and Accuracy Assessment. .......... ................................................... ...... ....... ..... .. ... 46
3.5 Topographic Discrimination ...................................... .... .... .. .... ..... ........... ...... ..... ......... .... ....... ... .. ................ 47
3 .5.1 Mapping Objectives: Selecting Topographic Features and Map Scale ........ .......................... .. .. .. ............ 4 7
3.6 Parameterization................................... ......... ...... .. ................. ...... .. .. ............................................... .... ........ 48
3.8 Accuracy Assessment ..... ........... .... ......... ... .. ..... .... .... ............................... .................... ....... ........... ........ ... ... 52
Chapter 4 Results and Discussion ... .... ................. .. ....... ... ... .......... .... .. .. ..... ........... ... ...... ... ....... .................. ..... ....... 54
4.1 DEM Interpolation and Error Assessment ................................ ...... .. .. .. .... ....... ... ......................... ... .... ....... ... 54
4.2 Future Research for DEM Interpolation .......................... ................................................... ........ .. ..... ........ ... 58
4.3 Feature Extraction ................................................................................................................. ........... ........... 59
4.3 .1 Curvature ..................... ... .... ... ....... .. ... .... .... .... ........ .... ...... ... .. ....... ...... .. .......... .... .. .. .... ..... .. .. .. ... ...... .. ..... 59
4.3 .2 Slope Gradient and Aspect ........ ...... ......................... .. ......................................... .. ................................ 60
4. 4 Accuracy Assessment ....................... .. ............................. ....................................... ........ ....... .... .................. 64
4.5 Discussion ................................................. .. ... .. ... ................. ...... .. .. ... .. .... .... .......................... .... ... ... .. .. ........ 67
4.5.1 Implementations ...................................... .. ...................................... .... .... .... ............ ... ............. ... .... ..... ..... 69
4.5.2 Limitations and Usefulness ....................... ... ... .. .... ..... ... ...... ............. .. .. .... .... .... .... ... ........................... ....... 70
References ................................................... .. ...................................................................................................... .7 4
Appendix 1 Figures ................................... .... .... .... .... .... ......... ........................... ... .... ................... .............. ............ 81
Appendix ll ............................................................ ....... .. ...... ...... .... .... ...................................... .... .. ... ..... .. ......... 117

ii

Table of Figures

Figure 1: Map ofBritish Columbia and Study Area Location .... ...... ......... .... ... ...... ...... .... .......... 82
Figure 2: Digital Elevation Profile ...................................... .. .. ............. ........ .. .... .... .. .. .. .. ...... ... ... 83
Figure 3: Slope Normal ............................... ..... ...................... ..... .... ... ........................ .... .. ...... .. .84
Figure 4: Curvature Classified using Eigenvector Analysis ..... .... .... ........................ .... ......... ...... 85
Figure 5: Ridges, Passes, Peaks and Valleys .................... .............. ...................................... ...... 87
Figure 6: Experimental Variogram .. ........... .... .... ............................. ..... ...... .. ................ ............. 88
Figure 7: Semivariogram ...................................... .. ................. ....... .. ... ........ ............... ... .. .......... 89
Figure 8: Randomly Generated Check Points and Residual Surfaces ......................... ..... .. ... ...... 90
Figure 9: Shaded ReliefModel .. .... .. ........................................... ..... ........ ... ... .. ....... .... ..... .. ... ..... 92
Figure 10: Airphoto ofValley ................................. ............ .... .. ........ .. .... ... .. ............ ............ ..... 93
Figure 11 : Elevation Profile and Observer Perspective .... ............................................ ....... .. .. ... 94
Figure 12: Perspective and Planimetric Views of Shaded ReliefModels .... ... .... .... .......... ... .... .... 95
Figure 13 : Elevation Profiles ............ ............ ..... .......................... ... ............. ... ......... ......... ......... 98
Figure 14: Slope Profiles .................... ................................... ..... .... .......... ...... ...................... ..... 99
Figure 15 : Flow Accumulation ............... .. ................... ........................ ... .. .... ... ........... .. .... ....... 101
Figure 16: Sampling Transect for Elevation and Profile Curvature ......................... .. ........ ....... 102
Figure 17: Curvature ......... ................................... ... ........ ..... ........................ .. ...... ... .. ... ...... ..... 103
Figure 18 : Profile Curvature ...... ... ..... .. .. ..... ... ............. ............ ................................. ................ 106
Figure 19: Ridge and Gully ................... .... ........... .. ........ ... .. .. .... ....... ... .... .................. ....... ...... . 107
Figure 20 : Slope Classes ....... .... ..................................................................................... .. ....... 108
Figure 21 : Slope Normal Vector Magnitudes ............ ..... ... ...... .... ............................................. 109
Figure 22: Relief .................... ... ... ............ .. ... .......................................................................... 110
Figure 23 : Cumulative Slope Length ...... .. .................. ............ ................. ..... ... .. ...... .. ........ ... ... 111
Figure 24 : Bioterrain Layers ........................... ... ... .. ......... .... .... ........ ................ .... ............. ....... 112
Figure 25 : Query ofBioterrain Layers ..... .. ..... ........................................................ .. ........... .... 115
Figure 26 : Soil Moisture Regime ............. .. .... ........... ....... ..................... ............... ..... .... ..... .. .... 116

iii

Table of Tables
Table 1: Topographic Map Elements ........................ .. ... ... .. .. ....................................... ... .... .. ...... .2
Table 2: Steps in Modelling Land Surface Data into a Map Classification System ........ .. ....... ..... .7
Table 3: Class Cover Changes Due to DEM Filtering ................... .. .. .. ... ...... ..................... ..... .... 20
Table 4: Overview oflnterpolation Methods .... .... .... ... ..... ... .... .. .... ............... .... .... ........ ... ..... .. ... 32
Table 5: Interpolation Root Mean Square Error .............. ..... ... ...... ......... ....... .... .. .............. ... .. .... 54
Table 6: Root Mean Square Error Stratified by Slope Class ..................... ...... ................... ......... 55
Table 7: Variation with Respect to the Mean Residual Stratified by Slope Class .... ... ... ... .. .... ..... 55
Table 8: Interpolation Results ......... .... ...... .... .... ...... ....... .... .... .... .. ... .. ..................... .... ........ ... ..... 56
Table 9: Error Matrix ...... .... ........................... ......... ....... ... .... ..... .... ... ....... ........ .. .. .... .... ...... ....... 65
Table 10: Errors ofComission and Omission .. ......... ........................ ......................... ........ ...... .. 65
Table 11 : Assessment of Topographic Elements .. ......... .... .... .... .... .... .... ......................... ...... ...... 67

iv

Abstract
This thesis is an exploratory analysis of automated mapping protocols that can be used to support
Terrestrial Ecosystem Mapping and Predictive Ecosystem Mapping in British Columbia. This
thesis employs neighbourhood analysis of elevation and its derivatives to discriminate the
bioterrain elements defined by Terrestrial Ecosystem Mapping standards. In achieving these
standards, discrimination beyond the basic topographic forms presented in current research is
explored.
The method developed strives to be
•

easily implemented by mapping projects employing standard GIS software

•

flexible so that the extracted topographic forms can be tailored to varying project objectives

•

compatible with the hierarchical procedure employed in Terrestrial Ecosystem Mapping

•

efficient and accurate in that the process is advantageous over manual mapping methods .

The effect of data quality is addressed through an assessment ofDEM data interpolation
techniques and classification accuracy. Random and systematic artifacts of the DEM that
influence the quality of the derivatives are explored. The issue of scale-dependent shape is
addressed by the constraints of objective-based mapping in which a map scale is specified and
the most basic shape elements are aggregated into contiguous classes by a roving neighbourhood
window.
The results indicate that basic topographic elements are mapable from relief as well as first and
second order elevation derivatives. These results give preliminary accuracy of 80% based on the
three classes tested. The procedure requires decisions at every step, but it is felt that this
complements the traditional mapping process in that it is hierarchical, and requires a synthesis of
extensive knowledge of vegetation and landscape across many scales.
Key Words: elevation, digital elevation model, topography, slope, aspect, curvature, Terrestrial
Ecosystem Mapping, Predictive Ecosystem Mapping, scale, random, systematic error.

Chapter 1 Introduction and Rationale
1. 1 Introduction
In the province ofBritish Columbia, the Resource Inventory Committee provides provincial
standards for Terrestrial Ecosystem Mapping with applications in vegetation, biodiversity and
ecological mapping. These ecosystem maps consist of landscape units, or polygons delineated
by comparatively homogeneous physical and biotic characteristics (Jones et al., 1999).
Terrestrial Ecosystem Mapping, as defined by the Province of British Columbia (Ecosystems
Working Group, 1998, page 1), is "the stratification of a landscape into map units, according to a
combination of ecological features, primarily climate, physiography, surficial material, bedrock
geology, soil and vegetation" . The resulting maps are valuable tools for resource management
and wilderness preservation within the 59 million hectares of crown land in the Province of
British Columbia.
Fieldwork and airphoto interpretation, the foundations ofterrestrial mapping, are costly, timeconsuming, and labour intensive to undertake across a large land base (Dikau, 1989 and Jones et
al. , 1999). In order to address the expense of mapping, the government ofBritish Columbia
assembled a Terrestrial Ecosystem Mapping Alternatives Task Force to develop standards for
Predictive Ecosystem Mapping. To reduce costs, Predicitive Ecosystem Mapping is a digital
modelling approach incorporating Terrestrial Ecosystem Mapping standards, pre-existing spatial
and attribute digital data sources and remotely acquired data sources (for example, vegetation
maps resulting from classification of remotely sensed data) into a model to predict a final
ecosystem map product.
The classification scheme to be used in this thesis is described in the guidelines of Standards for

Terrestrial Ecosystem Mapping in British Columbia (Ecosystems Working Group, 1998). In

2

the Terrestrial Ecosystem Mapping guidelines, a distinction is made between terrain and
topography. Terrain reflects the landform and parent material as well as any active and/or recent
geomorphic processes. Topography describes landscape position, shape, aspect, and slope that
are collectively referred to as bioterrain units. The latter will be the focus of this thesis.
Bioterrain unit definitions from Vegetation Resources Inventory Photo Interpretation
Procedures (Resource Inventory Branch, 1999) are given in Table 1.
Table 1: Topographic Map Elements
Bioterrain Unit
flat (plain)
gentle slope
moderate slope
moderately steep slope
steep slope
toe slope
lower slope
middle slope
upper slope
ridge
depressions

warm aspect
cool aspect

Description
level or very gently sloping <=5% (<= 30), unidirectional
(planar), local surface irregularities generally have a relief < 1m
unidirectional (planar), >5%, <=26% (>30, <= 15°), smooth
profile either straight, concave, or convex, local irregularities
generally < 1 m relief
unidirectional (planar), >26% <=49% (> 15° <=26°) smooth
'
'
'
longitudinal profile either straight, concave or convex; local
surface irregularities generally < 1 metre relief
unidirectional (planar), >49%, <=70%, (>26° ,<=3 5 °) smooth,
longitudinal profile that is either straight, concave or convex,
local surface irregularities generally < 1 m relief
unidirectional (planar), >70%, (> 35 °) smooth, longitudinal
profile that is either straight, concave or convex, local surface
irregularities generally < 1 m relief
area differentiated from the lower slope of the hill, generally
concave surface, specific aspect
between valley floor and midslope, significant break in slope
falling to moderate gradients, generally concave, specific aspect
steep to moderately steep consistent slopes and steep rock walls
with talus along the base, involve large areas, specific aspect
directly below steep bedrock summits and above rnidslope,
generally convex below crest, specific aspect
generally convex uppermost portion of a hill, convex in all
directions, no distinct aspect, also a crest
circular or irregular area of lower elevation (hollow) than the
surrounding terrain marked by an abrupt break in slope, greater
than 2 m in depth, generally at foot of meso scale hill or
generally level area
135°-285°, slope>15° (26%)
285° - 135°, slope>15° (26%)

Adapted from Resource Inventory Branch, 1999; Howes and Kenk, 1997; Ryder, 1994; and
Luttmerding, et al. , 1990.

3

According to the Standard for Terrestrial Ecosystem Mapping in British Columbia (Ecosystem
Working Group, 1998), mapping classification follows a hierarchical system. Pre-typing
involves delineating alpine tundra, parkland and other biogeoclimatic zones (Meidinger and
Pojar, 1991 ). Prior to more detailed mapping, it is advised that the pre-typing involve
considering the study area in terms of broad landscape areas of similar, repeatable units such as
alpine, mountain slopes, valley floors, plateaus and plains. To this framework, more detailed
mapping information is added to the site series level (i.e. , vegetation composition, soil and site
description) (Ecosystem Working Group, 1998). The bioterrain elements considered in this
thesis are part ofthe site description and, as illustrated in Table 1, are summarized as landscape
position, shape, aspect and slope gradient.
Since fieldwork and airphoto interpretation are costly, time-consuming, and labour intensive to
undertake across a large land base, any supporting method that automates, in keeping with
British Columbia provincial RIC mapping standards (Resource Inventory Branch, 1999), all or
part of the mapping process will be a valuable, cost-saving tool. This thesis investigates the
usefulness of digital elevation data in mapping the topographic component of the Terrestrial
Ecosystem Mapping process employing a digital elevation model of a portion of mountainous
landscape in northern British Columbia (Figure 1).

1.2 TRIM Elevation Data
The government ofBritish Columbia' s Terrain Resource Information Management (TRIM)
Program provides digital planimetric maps at 1:20,000 scale. Of interest to this thesis is the
TRIM Digital Elevation Model (DEM) data layer photogrammetrically captured under the
direction of British Columbia Surveys and Resource Mapping Branch (Surveys and Resource
Mapping Branch, 1990). ADEM is a set of points representing an elevation surface in explicit x,
y and z values (in the case of the TRIM data, the coordinates are in Universal Transverse

4

Mercator (UTM) and the elevation in metres). These data will be used to create an elevation
relief surface from which the bioterrain units will be extracted.

1.3 History
Prior to the advent of readily-available computers, limitations in processing large quantities of
data resulted in a focus on manual techniques emphasizing the analysis of elevation profiles
(Evans, 1972; Embleton, 1990). While these techniques were sophisticated, and in fact surface
modelling pre-dates the electronic age, manual computation of even the smallest number matrix
is an extremely time consuming process. As a result, statistical analysis for gridded elevation
data was not used extensively prior to the 1970s. By the mid 1970s, Mark's (1975) literature
review noted few attempts in which researchers employed the computational power and data
handling capabilities of computer software to analyze topography quantitatively as a tessellated
model. Computer processing power was greatly improved by the late 1980s. However,
researchers had not yet taken advantage of the ease with which morphometric parameters could
be modelled with a digital elevation model (DEM). O'Neill and Mark (1987) believed the gap in
the literature was a result of difficulties in constructing DEMs from analogue (paper) sources.
The early 1990s saw the advent of numerous, readily-available, government-generated DEMs.
Embleton and Liedtke (1990) noted that by the mid 1980s the literature began to show evidence
ofincreasingly challenging attempts at data reduction with respect to DEMs. Pike (1988) wrote
that the 'enabling technology' of computer-aided analysis would be useful in the automation of
shape analysis. Looking back through two decades of research, by the mid 1990s Gong (1994)
noted that rapid advancements in computer and software technology had led to a revolution in
surveying and mapping in the digital realm. However, the potential for digital elevation data
analysis as a landform classification tool that automatically subdivided landscape into its
constituent landforms had yet to be fully realized by the mid-1990s (Dymond et al. , 1995).

5

Within a year of Dymond's review of contemporary literature, there appeared two studies
(Pelligini, 1995 and Wood, 1996) presenting methodologies that classified aDEM into
topographic groupings.

1.4 Statement of the Problem and Objectives
The overall objective is to develop and test procedures for using a digital elevation model and its
derivatives to meet Terrestrial Ecosystem Mapping standards for the classification of the
topographic elements described in Table 1. Specific objectives are:
•

to develop and test a method to generalize topography from the raster DEM and its
mathematical derivatives, based on Terrestrial Ecosystem Mapping standards

•

to assess the validity of the method by comparing its results with independently derived field
and ground truth data.

It is intended that the method be considered in terms of:
•

ease of execution and further use employing standard GIS software

•

flexibility, so that the extracted topographic forms can be tailored to changing project
objectives

•

compatibility with the hierarchical procedure employed in Terrestrial Ecosystem Mapping

•

efficiency and accuracy in comparison with manual mapping methods.

1.5 Thesis Organization
This thesis is organized into six chapters:
In Chapter llntroduction and Rationale I introduce and define Terrestrial Ecosystem Mapping
and its importance to resource management in British Columbia as well as discussing the merit
of an automated procedure that supports Terrestrial Ecosystem Mapping known as Predictive
Ecosystem Mapping. TRIM (Terrain Resource Inventory Mapping), (Surveys and Resource
Mapping Branch, 1990) elevation data are introduced as a tool to support mapping objectives.

6

The classification terminology is established, a brief history of landform mapping from elevation
data is given, and the thesis objectives are stated.
In Chapter 2 Literature Review I describe automated classification paradigms and compare them
with traditional manual approaches. The discrimination of features is discussed in terms of
elevation derivatives. A literature review is presented based on six sections:

•

classification paradigms

•

feature parameterization

• topographic elements
•

data model selection

•

accuracy assessment

•

interpolation ofDEM

Chapter 3 Methods I outline the method for DEM interpolation accuracy assessment, and
neighbourhood analysis of elevation and DEM derivatives for bioterrain element classification.
Chapter 4 Results and Discussion I discuss the strengths and weakness of the methods tested and
give an accuracy assessment.

7

Chapter 2 Literature Review
2. 1 Classification Paradigms
Two views ofthe process of automated map classification are summarized in Table 2.

Table 2: Steps in Modelling Land Surface Data into a Map Classification System

Step 1
Step 2
Step 3

Townshend, (1981)
Feature selection
Definition of rules delimiting classes
Extraction of features from images
Assigning features of the image to classes

Brassel & Weibel, (1988)
Structural recognition
Process recognition
Process modelling
Process execution

Brassel and Weibel's (1988) term "structural recognition" aims to identify the desired map
objects and is controlled by the objectives of the mapping procedure, the quality of the database,
map scale and communication rules (graphic and perceptual limits).
Townshend's (1981) equivalent term for structural recognition is feature selection. In this stage
the target structures are identified. Townshend discusses the traditional (human) contribution at
this stage as being a process whereby the interpreter draws upon keys, analogues (possibly
existing only in the interpreter's imagination) and deductive reasoning to classify the geographic
data. Brassel and Weibel (1988) describe this phase as an intellectual evaluation. In both
discussions, the authors emphasize the importance of the decisions of the human interpreter at
the onset of the automation process. Dikau ( 1989) notes that relief forms are subjective in that
their definition is pragmatic and often influenced by the knowledge of the interpreter.
Deductive reasoning involves conclusions reached by the interpreter based on evidence in the
image supported by ancillary information. The central component of human interpretation is the
use of contextual information. The interpreter sees more than slope angle, for example, and can

8

integrate a vast set of topographic features at different scales. This is a process described by
various researchers as synthetic: a synthesis of multiple scales, multiple views and the human
experience (Townshend, 1981 ; Pike, 1988; Dikau, 1989). In this context, computer automation
of feature selection can be limited by the simplistic linear processing capabilities of the computer
language. Brassel and Weibel (1988) consider the process recognition stage a difficult stage to
automate. Artificial Intelligence, for example, neural network analysis, is a an analysis
proceedure resembling that of the thought processes of the human brain and is capable of pattern
recognition. Artificial Intelligence in the form of neural network analysis is not considered here
for two reasons. Firstly, it does not meet the criteria set out in this thesis that the method be
easily executed and employed in standard GIS software. Secondly, as Gong (1995) notes, neural
network analysis remains poorly understood in terms of its usefulness .
Once the target structure parameters are established, these are modelled as a sequence of
operational steps viewed by Brassel and Weibel (1988) as "process modelling" . "Process
recognition" and "process modelling" in Brassel and Weibel's description of map generalization
are described by Townshend (1981) as the "definition of rules delimiting classes" . Pike (1988) t
viewed errain features as part of a geometric continuum. The challenge, in his view, is to
categorized discretely this continuum. The compilation of rules outlined in these operational
steps relies on the derivatives of elevation and a statistical method (Pike, 1988).
The running ofthe model is termed "process execution" and the results are viewed as a map.
Brassel and Weibel's process execution is described by Townshend in two steps, "extraction of
features from images" and "assigning parts of the image to classes" .
The strengths of an automated approach are computer speed and consistency during the feature
extraction stage, which contrasts with potential human error due to fatigue, variation between

9

interpreters and the slower pace of work. A weakness of an automated approach is the inability
to reproduce the synthetic results of the highly experienced human interpreter.

2.2 Parameterization of Features
Elevation derivatives of slope, slope gradient, slope between slope reversals, slope length, aspect
and curvature are illustrated in Figure 2 and are discussed in the following sections.

2.2.1 Slope Gradient and Aspect- The First Derivatives of Elevation
Slope is the first derivative of elevation (the rate of change of elevation in any direction) and is
defined as angle between a tangent to the surface and a horizontal (geoid-parallel) surface.
Gradient is the slope in the direction of maximum slope (Evans, 1972; Evans, 1980). While
gradient is given as an instantaneous value, it is measured over an area (Hodgson, 1995). The
definition of that area in the field can be arbitrary (in raster DEMs it is by convention a distance
of three cells). Aspect is the direction of the slope gradient for any given location and is
measured in degrees from grid north.
The slope gradient measure is derived from a mathematical surface. Florinsky (1998) states that
true landscape is not smooth, therefore a comparison ofDEM derivatives to a reference (ground
truth data) is not valid. He states that aspect errors are found in flat areas while slope gradient
errors are predominantly found on steep slopes. The effect ofDEM resolution and accuracy are
a significant factor in derivative accuracy. Florinsky further proposed assessing the accuracy of
the derivative calculation technique along with assessing the accuracy of the DEM.

2.2.2 Curvature- The Second Derivative of Elevation
Curvature is the rate of change of slope (Evans, 1972; Evans 1980) and can be measured as
profile curvature (the rate of change of gradient) or planimetric curvature (the rate of change of
aspect) or a mean curvature (average ofboth). Curvature can be expressed as convex, concave
or both (such as a saddle which is assigned a value of zero).

10

2.2.3 Slope Normal
The slope normal is a vector oriented at 90° to a plane tangent to any point on the surface of the
DEM. It is a unit vector defining the orientation of the cell surface in terms of slope gradient and
aspect. Figure 3 illustrates the three angles that define the unit vector in relation to the
orthogonal axes x, y and z. Characterization of a neighbourhood can then be undertaken using
directional statistics (Hodgson and Gaile, 1999). For example, unidirectional (planar), convex
and concave are terms used to describe the majority of the physiographic units in Terrestrial
Ecosystem Mapping. In the planar class, the slope normals are parallel or sub-parallel and
dispersion is zero. Dispersion indicates the variability or spread of the unit vectors in space and
is similar to the standard deviation of the normal distribution. In the concave class the
dispersion of the slope normals is minimized, approaching the value of the resultant normal
(mean normal for an area of interest). In the convex class, the slope normals approach maximum
dispersion away from the resultant normal.

2.3 Operational Definitions: Case Studies Classifying Topographic
Elements
Aggregation of an elevation surface into classes requires the ability to describe desired shape
elements using meaningful parameters. The shape of any part of an elevation surface is
recognized as a spatial phenomenon. Its identification is inherent in how elevation varies across
the surface of the shape of interest. Roughness describes the variation amongst elevation values.
Roughness is defined by vertical (relief) and horizontal (texture and grain) criteria. The primary
elevation derivatives of gradient and aspect, along with the secondary derivative of curvature,
describe roughness in terms of the interaction of vertical and horizontal criteria. The use of
gradient and aspect as succinct descriptors of roughness finds its origins in an overview paper on
research in geornorphornetrics by Evans (1972) .

11

Having found variable terminology and methods throughout his literature review, Evans (1972)
set out to critique and clarify operational definitions and sampling procedures related to
geomorphometrics. Evans's critique noted the use of as many variables as possible in
multivariate analysis resulted in the potential contribution of irrelevant variables. Gong ( 1994)
wrote that the low accuracy reported in many research initiatives he reviewed was a result of a
lack of knowledge of suitable variables. Gong stated "we have a long way to go in studying
algorithms for integrated analysis of multi-source data before an optimal selection of algorithms
for certain tasks can be definitely made without much work on a trial-and-error basis" (p. 353).
Evans critiqued the variation in terrain parameters used by researchers as well as what he viewed
as an erroneous quest for one single unifying parameter which could discriminate landform.
Evans proposed a unifying suite ofvariables (slope, aspect, profile curvature and plan curvature)
to describe elevation in a non-redundant manner. He considered slope as arguably the most
important descriptor of surface form .
Evans also noted a distinction between general and specific geomorphometrics. General
geomorphometry measured and analyzed landform characteristics applicable to any continuous
rough surface. Specific geomorphometry dealt with the measurement of specific landforms (e.g.
cirques, drumlins, or stream channels.) in a manner that was able to distinguish a specific
landform from another. Specific geomorphometry suggested landform processes as described by
terrain definitions. Mark' s (1975) research continued with a focus on landscape shape
independent of geomorphic process (general geomorphometry).
The relation between the horizontal and vertical dimensions of relief was noted by Mark (1975)
as being described by gradient and aspect. He hypothesized that all form could be found in
surface roughness (local relief) . He noted that roughness operates at two scales, grain (longest

12

wavelength - usually on the order of ridges between valleys) and texture (based on Nyquist
sampling theory, the shortest wavelength is twice the cell size).
Evans ( 1972) was critical of the analysis of gradient and aspect together as producing
undesirable results, and recommends that they be treated separately. Mark supported Evans'
view as valid when the tessellation model is composed of irregularly distributed surface points,
but noted that with a regularly spaced model (raster DEM or regularly sized triangles), or
weighted, irregular tessellation, the use of unit vectors to characterize slope gradient and aspect
may be valid.
Valleys might not form a regularly-spaced pattern and, in reality, tend to converge up-valley.
Evans (1972) noted that, for this reason, research attempts to use spectral (Fourier) analysis of
wavelengths as a model for valleys failed . Spectral analysis assumes a stationary signal such that
the mean, variance etc. should be independent of location. This is not necessarily true of
elevation data. Fourier transforms use sine functions as the basic functions into which sample
data are decomposed . Sine functions may not necessarily represent the underlying landscape,.
nor do non-sinusoidal shapes give harmonics that represent the underlying landscape. Focal
analysis, in which a floating window performs an analysis on each data cell in turn, is useful in
Fourier analysis in that it helps to overcome the problem of non-stationary data. However, the
size of the window may not reflect the scale of the underlying landscape feature .
Chui (1992) addressed this issue by employing wavelet analysis that can be translated and dilated
to reflect the underlying position and local wavelength of elevation. Wavelet analysis produces a
power spectrum at every cell. This creates a large output dataset requiring further processing in
order to obtain meaningful results.

13

Pike (1988) employed non-redundant topographic attributes to quantifY a geometric signature
(i.e., a set of measurements capable of discriminating landform units) that described the
continuous topography of landslide types. The signature was derived from pre-mapped landslide
types for which statistical estimates of central tendency and dispersion for altitude, altitude
variance spectrum, slope between slope reversals, slope gradient, and curvature at fixed slope
lengths were calculated.
Using multivariate analysis Pike statistically discriminated two extreme landslide types: very
hard and soft. The very hard landslide type consisted of regularly spaced, straight ribs between
sharply incised flutes, sharp crests and steep slopes. "Soft" consisted of no flutes, irregularlyspaced and poorly-incised drainages, broadly-rounded, largely gentle slopes with some steeper
than 26°. He was not able to discriminate any of the intermediate types.
Based on field engineering studies, Dikau (1990) hypothesized a high correlation between
landslide location, geological formation, and geometrical slope properties of gradient, aspect and
curvature. He hypothesized these attributes to be predictors of landslide hazard. He aimed to
analyze the frequency distribution of these variables in order to provide weighting factors for
landslide-susceptibility maps. It was a system intended to predict the potential for landslides in
any given area and not to delineate a present landslide as a landform unit. Dikau's (1989)
detailed description of landform analysis involved 20 variables, 200m cell resolution and filter
windows 9.8 km x 9.8 km in size. He viewed his discriminant units as reliefforms composed of
form elements, which in turn are composed of form facets, which in turn can be grouped by userdefined criteria.
Schmid-McGibbon ( 1993) tested a pre-existing terrain map for unit homogeneity by employing a
neighbourhood filter for geomorphometric parameters of slope gradient and down-slope

14

curvature. She tested the results using relative variance and intra-class correlation. Her
assumption was that the variance within map units is smaller than variance between map units
for units of homogeneous terrain characteristics. As with Pike's work, the values of slope
gradient and curvature were grouped using boundaries from a pre-existing vector map of the
desired landform types and the summary statistics derived.
At the onset of Schmid-McGibbon's thesis, she noted that the qualitative selection of filter size
"does not consider accurately the level of detail that may be lost in the generalization process"
(p. 27). She concluded that a reliable quantitative method for filter size selection cannot be
determined and that neighbourhood sizes should be determined in conjunction with a visual
interpretation of the maps, one that is left to human judgment and should consider the objective
ofthe mapping project in the decision process.
Her conclusions were:
•

The selection of a larger neighbourhood size for filtering the slope gradient map was
appropriate to represent the dominant slope gradient classes within units of homogeneous
landform patterns, and to eliminate the slope gradient differences among individual landform
elements.

•

A by-product of the generalization- the short-wave map of the slope gradient- is of
geomorphological utility, as the map shows considerable detail ofthe specific landform
elements.

•

A by-product of the generalization - the filtered curvature map - shows sufficient detail with
which to extract individual landform elements and patterns.

15

Blaszczynski (1997) undertook a two-stage automated landform classification. The first stage
classified the DEM according to a cell being in a focal neighbourhood of general convexity,
concavity and flatness . The second stage relied upon hydrologic terrain analysis capabilities
derived from methods introduced by Jenson and Domingue (1988) in order to classify the terrain
into ridge and valley lines. His final product was a map that included crests, concave areas,
convex areas and troughs.
Blaszczynski employed a 3 x 3 floating neighbourhood window in order to evaluate elevation
using a modified average percent slope gradient equation by removing the absolute value for
calculations of rise. Positive, negative or zero values were assigned to the centre cell based on
whether the equation resulted in a predominantly concave, convex or flat neighbourhood. The
use of local relief to assign concavity or convexity for ridge and ravine discrimination as
described by Blaszczynski previously appeared in Jenson's (1985) work on hydrologic basins.
In Blaszczynski's methods, flat areas could not be determined as sloping or non-sloping. He
suggested that the application of a slope gradient algorithm to flat areas would then discriminate
horizontal from sloping areas. In areas of saddle points, convexity canceled out concavity and
the cell was assigned a value of zero.
Further results were derived by Blaszczynski using the watershed analysis of Jenson and
Domingue (1988):
•

fill single-cell depressions (considered error in the models)

•

calculate the flow directions for the filled DEM (the flow direction for each cell is the
direction water will flow out of the cell along the steepest path that is encoded in a flow
direction map)

16

•

assign integers to cells representing how many other cells are flowing into it, where zero
represents crests of ridges.

This process is more computationally intensive by many orders of magnitude compared, for
example, to a single pass of a 5x5 low pass filter over the same data set. The results delineate
ridges as a single cell-width line of zero flow accumulation. This single line becomes chaotic at
the down drainage end of the ridge.
An examination ofBlaszczynski's map reveals terrain that is gently curving, repeating, linear
ridges separated by gently curving, evenly spaced drainages. It is questionable whether his
procedure is repeatable in rough mountain terrain. He noted that the method to identify ridge
crests was no more than adequate. Extracting ridges and gullies from DEMs has been attempted
by many researchers in the area of hydrologic modelling (for example Riazanoff and Cervelle,
1988; Jensen and Domingue, 1988; Mark, 1984).
Hobson ( 1972) investigated the use of slope normals as a means to describe surface roughness
related to topography. He noted that vector magnitude ranged in value from zero (no preferred
orientation to the normals, i.e., random terrain) to a value of one (identical orientation, i.e., a flat
surface, sloping or horizontal) . He concluded that magnitude was usually high and vector
dispersion low in areas characterized by similar elevations or equal rates of elevation change. In
areas of non-systematic elevation changes, values of low vector strength and high vector
dispersion were found.
Hobson's use of vectors has an extensive basis in geostatistics. Hodgson and Gaile (1999)
identified certain problems in using vector fields in GIS : an absence of data type and/or operators
in GIS or Remote Sensing systems, and a lack of directional statistical capabilities in GIS ,
Remote Sensing and statistical packages.

17

Brabyn ( 1997) used neighbourhood analysis for macro landform classification on a 200m
resolution DEM. No justification for the neighbourhood shape or size was used in his analysis.
He grouped slope gradient as flat ( <8%) or sloping. A focal mean with a radius of 5600 m was
applied to the slope gradient classes. Relative relief was calculated within the 5600 m
neighbourhood and then classed into intervals. Brabyn then classed flat areas as upland or
lowland. He noted an undesirable result of this process to be banding between plains and
mountains. This progressive zonation was considered undesirable and was probably a result of
the neighbourhood window incorporating flat and sloped areas in progressively changing ratios.
The resulting map legend contained the following map elements:

•

flat or nearly flat plains

•

smooth plains, some local relief

• tablelands with moderate relief
•

plains with hills

•

plains with high hills

•

plains with low mountains

•

plains with high mountains

•

open high hills

•

open low hills

•

open low mountains

•

open high mountains

•

low mountains

• high mountains .
Because of coarse cell size, 5 km neighbourhood filter size and differences in legends, the
Brabyn (1997) method does meet TEM standards. In an example ofEvans' specific

18

geomorphometry, Tribe (1992) outlined a multi-step, detailed, decision-based procedure to
extract valley heads and valley bottoms from aDEM. The results of this detailed process were
partially successful and were implemented in a hydrologic model on an artificial data set.
In Dymond' s (1995) work on mapping land components, he commented that no research had
appeared that automatically subdivided the whole landscape into landform units. He generalized
aspect by constraining the generalization process using streams with ridges defined by flow
accumulation. Next, he split the aspect map into slope position by employing elevation
thresholds to derive seven slope position classes. The flow accumulation data are, as has already
been pointed out, computationally intensive to derive. The accuracy assessment of this
procedure was a visual comparison to airphotos. Dymond's procedure would not meet
Terrestrial Ecosystem Mapping standards for slope position. The break in slope used to define
toe slope and valley bottom does fall consistently an incremental distance between river and
ridge crest. Elevation thresholds for slope position would not be effective in all but the most
idealized situation (very low gradient along valleys and ridges).
Pellengrini (1995) and Wood (1995) published theses that automatically subdivided DEMs into
landform units. Pellengrini makes the distinction between continuous and discontinuous analysis
in order to classify shape from aDEM. Discontinuous analysis, in his view, refers to the
neighbourhood analysis (floating window) used to identify ridges, valley lines and cumulative
flow used in watershed analysis. Continuous analysis involves slope gradient and curvature
analysis.
Pellengrini further made a distinction between the analysis situation in which the DEM has a set
resolution and the analysis is based in part on testing the results ofvarying sized window.
Because the researcher seeks a set of unambiguous signatures for shape, the filter size is

19

important to that result. His alternative employed a fixed window size and altered the data set
until it was sufficiently smooth to match the ideal characteristics of the filter itself
He employed Fourier analysis as a tool to reduce noise by filtering in the frequency domain. The
resulting surface was then analyzed employing a Hessian 2x2 matrix created from the second
derivative of a partial quartic equation. From this matrix, equations and unknowns are derived
which correspond to maximum and minimum curvature. The determination of the roots of these
equations employs an area of linear algebra theory that uses algebraic eigenvalues. His resulting
surface is discriminated into the following legend:
•

Peak (cell higher than neighbours)

•

Pit (reverse of peak)

•

Ridge (cell higher that surface points on either side)

•

Valley (opposite of ridge)

•

Hillside (not flat, exclusive of previous elements and sorted into convex, concave, saddle,
sloping and inflection).

Figure 4 a and b illustrate an implementation ofPellengrini's work on the study area to be used in
this thesis. The map is visually difficult to interpret and illustrates that the classes are
fragmented into small units. The map could be simplified based on cartographic scale
considerations. Based on project objectives, there is a minimum polygon size usually based on
map output constraints. At a scale of 1:20000 the minimum polygon size for this map would be
32 contiguous raster cells at a cell resolution of 25m. A more detailed description ofthe class
definitions would simplify this map. For example, are ridges taken to mean any raster cell above
its neighbouring cells to the side? The definition can be expanded to include that the ridge is a
long, linear feature, in an area of specific relief Ridges in low-lying areas and of small areal
extent would be excluded. If the input DEM were smoothed further utilizing frequency domain

20

filtering, Figure 4b results. Table 3 illustrates the percentage of class cover types resulting from
two levels of filtering. It can be seen that the degree of filtering has a profound effect on the
percentage of class coverage.

Table 3: Class Cover Changes Due to DEM Filtering
Highly Filtered 3 b
Shape
Class Value Frequency Percent
0.02
Peak
1
48
Ridge
13746
5.24
2
Saddle
3621
1.38
3
3.12
Flat
4
8172
Ravine
13813
5.27
5
Pit
6
69
0.03
13.43
Convex Hillside
7
35208
91672
34.97
Saddle Hillside
8
8.76
Slope Hillside
9
22970
14.06
Concave Hillside
10
36856
Inflection Hillside
13.45
11
35269
Unknown Hillside
12
700
0.27

Moderately Filtered 3 a
Shape
Class Value Frequency Percent
Peak
1
218
0.08
Ridge
2
27262
10.4
Saddle
3
4381
1.67
Flat
4
4057
1.55
Ravine
5
27345
10.43
Pit
6
278
0.11
Convex Hillside
7
29237
11 .15
Saddle Hillside
8
108973
41.57
Slope Hillside
9
6727
2.57
Concave Hillside
10
29092
11 .1
Inflection Hillside
11
19773
7.54
Unknown Hillside
12
4801
1.83

Table 3 illustrates that the degree offiltering produces two different landscapes from the same
DEM.
Pellingini's work is an excellent tool which can be used to group curvature classes. The question
arises as to whether the technique could be modified to fit project objectives (e.g. Terrestrial
Ecosystem Mapping) and as to how to evaluate the degree of filtering required, and whether it is
necessary to undertake eigenvector calculations to achieve project goals. This last question
illustrates a valid need to consider computation intensity as one criterion by which to evaluate
methods. The gains of faster computer speed and larger storage capacity are being taken up by
research and applications covering huge land areas, add to this computations which are process
intensive and data layers modelled as 32-bit, and computer system processing speeds and storage
capacities are quickly stretched to the limit of effective processing.

21

One other researcher to successfully classify landform using fitted surfaces is Wood (1995).
Wood employed first and second derivatives of quadratic surfaces fitted over multiple scales to
discriminates pits, passes, peaks, ridges, and valley lines. Further to his PhD work in this area,
Wood ( 1998) has recently employed the semi-axis of conic sections to model surface form at the
ideal scale for that feature. Figure 5 illustrates the results ofWood's reseaerch as applied to the
inset area of interest (Figure 1). It is seen from this illustration that ridge, valley lines, passes
and peaks have been identified. This technique does not address the features required for
Terrestrial Ecosystem Mapping (Table 1).

An ideal system would identify the discrete landform desired. It is known that these forms can
be identified over many scales. While Wood addresses these issues, and creates a system that is
superior for the extraction of pits, passes, peaks, ridges, and valley lines, this thesis intends to
address a method that is objective based: that is, one that incorporates techniques which
discriminate the features required for Terrestrial Ecosystem Mapping. This requires more than
the items discriminated by Pellengrini (1995) and Wood (1996), and is a technique that aims to
be usable by the community employing Terrestrial Ecosystem Mapping as a tool.

2.4 Data Model Selection
The Digital Elevation Model and Triangulated Irregular Network Model are discussed in terms
fo their strengths and weaknessed with respect to data accuracy and facilitity of use in modelling
elevation derivatives .

2.4.1 Regularly-Gridded Elevation Values: The Digital Elevation Model
{OEM)
Samples of elevation data derived from photogrammetry and tied to benchmark elevation values
can be interpolated into a regular grid of elevation points. The topology of the model is explicit

22

in the raster file format provided that the row and column dimensions are known. As a result, z
values are stored for each point. Only minimum and maximum x and y coordinates of the matrix
are stored and coordinate values for any elevation point are derived as needed.
The implicit topology and resulting neighbourhood relations of regularly-gridded elevation cells
facilitate matrix analysis. Peuquet's ( 1984) view is that while the array structure is built into
such lower order programming languages as FORTRAN, higher order languages enable the user
to program for more complex data structures. Data redundancy occurring in flat areas can
contribute to large file size. Another disadvantage of the grid model, in Peuquet's view, is the
radial asymmetry that results in diagonal neighbours further apart than nearest neighbours. This
can be addressed by applying a weighting to the diagonal values.
Stocks and Heywood ( 1994) observed that the re-sampling of a random sample of elevation data
to the regular grid introduces elevation errors that may be unacceptable. Gridding errors are
reduced if the input elevations (control points) are more closely spaced and sufficiently dense to
adequately represent the desired features . In variable terrain, data redundancy on a regular grid
is minimized. Stocks and Heywood regard the grid structure ofDEMs, at suitable resolution, to
be an effective terrain model in mountainous regions

2.4.2 Triangular Irregular Networks {TIN)
A TIN is an irregular array of elevation points assembled as a network of edges forming planar,
space-filling, non-overlapping triangles. Significant landform features can be included in the
original irregular elevation points as breaklines which can be set to influence or constrain the
interpolation of sample elevation points to reflect such features as ridges, rivers or lakes, etc.
(Kumler 1994; Douglas, 1986; Stocks and Heywood, 1994).

23

An advantage of this elevation model is its variable resolution that allows for few, large triangles
in subdued landscapes and numerous, small triangles in high relief landscapes. In order to take
full advantage of variable resolution, only significant data points are extracted from the data set
prior to the interpolation. The inclusion of significant landforms as breaklines and breakpoints
allows for placement of these landforms as precisely as the original sampling protocol and the
storage precision ofthe model allow.
Each point is stored explicitly as x, y and z values. The topology of the points, triangle edges
and triangle adjacencies are also stored explicitly in files . As a result TINs are process intensive
because the tabular records must be searched in every analysis step as well as when viewing the

TIN . TIN derivatives of elevation, slope gradient and aspect for any given location must be
calculated and the resulting polygon coverage converted to a raster grid for analysis beyond
queries for slope gradient, aspect and elevation based on Boolean logic. The derivative values
can be stored for each triangle facet as part of the topological files . The input control points are
explicitly used as the triangle vertices and, as a result, the surface on any given triangle facet is
modelled as a flat surface which may not represent the true shape.
Kumler ( 1994) tested the accuracy of linear versus inverse-distance interpolators for a digital
elevation model and found that the linearly interpolated DEM had smaller root mean square
(RMS) error. He stated that no research existed that established the superiority of the TIN versus
DEM model. He concluded that a TIN derived from points extracted from a DEM of adequate
resolution was inferior to the DEM in modelling terrain based on an assessment of RMS error. A
more useful method would have been to create the TIN directly from control points prior to
comparison to the DEM.

24

While the variable resolution of TINs is advantageous with respect to minimizing data
redundancy, other errors can be problematic. With a variety oftriangulation algorithms to select
from, an ideal model should approximate the land surface with triangles approaching equilateral.
Where there is rapid transition from a generally flat area to an area of rough topography, the
triangles become long and narrow artifacts within the terrain model. Based on the work of
nineteenth century mathematician Herman Amanders Schwartz, the triangles became
increasingly long and thin, and the sum of the areas of these triangles diverges from the area of
the true surface.
The regular grid is an efficient method for retrieval and processing while the TIN is the more
efficient method for minimizing file size as well as being an effective model for view shed
analysis (Kumler, 1994; Stocks and Heywood, 1994). The gridded data structure is chosen in
this work for maximum accuracy in representing rough terrain and to facilitate computation.

2.5 Accuracy Assessment
2.5.1 Classification Accuracy Assessment
,

The aggregation of surface derivatives into classes of ecologically significant surface forms will
require an accuracy assessment. The classification will be compared against results determined
by airphoto interpretation.
According to Jensen (1996), Berry (1995) and Story and Congelton (1986), categorical data
require a different approach for accuracy assessment than discussed in the subsequent section on
DEM accuracy. The derived cartographic product can be tested against reference data known as
ground truth. Ideally, ground truth should rely on field site visits, but due to costs, views from
aircraft, airphotos and map products can be used as supporting data. Map products can be
considered ground truth if the accuracy is known and can be shown to be of higher quality than

25

the expected accuracy ofthe mapping project in process. While a single value such as percent
can be given for accuracy, a more thorough means of evaluating the accuracy of the map product
can be found in the error matrix.

;
•.

An error matrix is a way of comparing ground truth collected in the field or from airphotos

against predicted categories in the classified map . The two data sets are then tabulated as
normalized frequencies in a co-occurrence matrix. The rows and columns are labeled in the
same order with the class names. The column headings are assigned the ground truth class for
any given point while the row headings are assigned the class type found from the classification
method. The cell at which the same column and row number intersect, for example for ridge
position, would indicate how many of the ground truthing points known to be landform type
"ridge" were classified accurately. An examination of other cells in the same column gives
information regarding what class of landform ground truthing locations of ridge positions were
misclassified.
Overall accuracy is calculated by dividing the total number of correct ground truthing points
(the sum of the values along the matrix diagonal) by the total number of ground truthing
locations overall. In the error matrix, each map class can be evaluated for errors of omission
(producer's accuracy) and commission (user's accuracy). Errors of omission are evaluated as the
number of locations correctly classified for a landform type divided by the total number (all the
column values) of ground truthing points for that landform type. However, as with the running
example, each misclassified ground truth point for a ridge position represents an error for the
misclassified class in question. Such an error is an error of commission. Errors of commission
are calculated as the total number of ground truthing points correctly classified divided by all
ground truthing point locations (total row values) that were classified (correctly or incorrectly) as
that type. Values of the probability of chance agreement are also given for each class.

··.

'·:

...

26

In testing the significance of the results, the null hypothesis states that the predicted classification
differs from the expected classification for a given cell in the matrix due to chance sampling. In
order to obtain a measure of chance agreement, Congalton and Mead ( 1983) recommend KAPPA
analysis (Cohen, 1960). Equation 1 is the KAPPA coefficient of agreement :

Equation 1

where Po is the observed proportion of agreement (overall accuracy) and Pe is the expected
proportion of agreement due to chance (Rosenfield and Fitzpartick-Lins, 1986). A KAPPA
value incorporates row and column products and can be used to determine if the correct results of
the classification are a result of chance agreement, because it is a measure of the proportion of
agreement determined once chance agreement is removed (Foody, 1992).

2.5.2 OEM Accuracy Assessment
The quality of the cartographic results of this thesis or any DEM derivative will depend in part
on the quality of the elevation model from which they were derived . DEM error can propagate
through the derived products and influence the user's ability to successfully extract desired
landforms (Felicisimo, 1994; Florinsky, 1998; Lee et al., 1992; Bolstand and Stowe, 1994). In
order to understand the nature of any error and its possible sources (which includes the impact
of implementation of the procedures put forth in this thesis in practical applications), part of the
accuracy assessment should involve assessing the faithfulness of the DEM in modelling the land
surface.
Errors in DEMs can result from two sources: the error inherent in the raw data points and error
produced by the interpolation procedure with which the DEM is built (Brown and Bara, 1994;
Desmet, 1997). The effects of various interpolators are discussed in a subsequent section.

27

The accuracy of the source data and the interpolator will have a significant effect on the accuracy
ofthe DEM derivatives (e.g. slope gradient, aspect and curvature) (Gao, 1998). The error can
have three sources: blunders, system specific and random (USGS, 1998). Blunders, consisting of
misinterpretations, careless observations and recording failures, should be easily identified
removed from the data. System-specific error occurs in a predictable pattern with a constant
magnitude. This type of error can be caused by photogrammetric system errors such as
deformation in the camera system, deficiencies in ground control points and stereo model
orientation, and distortions introduced by the scanner. Systematic error can appear as a vertical
elevation shift in all or parts of the data, as fictitious features, for example striations, or as
improper readings caused by shadows, trees, or buildings. Random error is unpredictable in
origin, having a normal distribution with small error more frequent than large error.
Besides the following methods of error detection, visual verification is recommended (USGS,
1998) using display techniques such as hypsometric tints, stereoscopic viewing using anaglyphic
filters and shaded relief enhancement. These are particularly useful for blunders and for some
forms of systematic error.
Human error contributes to error in the raw data points. The photogrammetric stereo plotting of
the data points are gathered along lines. The operator moves south to north, for example, then
north to south to gather the adjacent line. As operators are evaluating the stereo image and keep
pace with the progress of the automated plotter, the operator tends to underestimate elevation
values as points are gathered as the line travels up slope, and overestimate elevation values as
points are gathered as the line travels down slope gradient. As Hunter and Goodchild (1995)
note, this results in striping between adjacent columns in the DEM.

28

In accuracy assessment, it is desirable to compare two independent data sets, one with higher
accuracy that describes the same elevation surface. The need for a second data set not used in
the interpolation results from the fact that most interpolators are exact in that the input elevation
values appear again in the interpolated surface. The test points should be selected in order of
preference from field control points (points sampled and surveyed on the ground),
areotriangulated test points, and spot elevations (USGS, 1998). Edge consistency of each map
sheet should also be checked.
Root Mean Square Error (RMSE) (Equation 2) is the most widely used technique when
comparing two data sets (e.g. Desmet, 1997). RMSE is a dispersion measure approximately
equivalent to the standard deviation between two data sets where the larger the RMSE, the larger
the difference between the two data sets.

Equation 2
where Z obs and Z int represent the observed and interpolated sample point and N represents the
total number of sample points.
RMSE does not provide a description ofthe variation of mean deviation across the DEM. This
assumption of zero mean deviation may not be valid as it implies zero systematic error (Desmet,
1997). The point density in the input data to the DEM is homogeneous, but the underlying
landforms may be heterogeneous and it is known that error is often dependant on the type of
landscape being modelled stereoscopically. As a result, the stated RMSE for the DEM cannot be
considered a stationary random function (Lee et al. , 1992). As a summary value, RMSE may
hide variation in error across a DEM.

29

Li ( 1994) proceeded with RMSE and since the accuracy of the check points was greater than that
of the input points. However, spot heights available in TRIM data are not randomly distributed
as they are selected from significant landforms (often ridges and peaks). Li (1988) suggests the
following indices in addition to RMSE: the mean absolute difference between interpolated and
true values (Equation 3) and the standard deviation of these differences (Equation 4).
l"lz = Li Z obs - Zint I
N

~

L(i Z obs - Z int \-l"lz)
N-1

Equation 3
2

Equation 4

These equations allow for the modelling of error across the map by providing a measure of the
variation of the error. The residuals can be mapped and a visual display of spatial distribution
and magnitude of the error can be viewed. It may appear that error in the data points could be
found clustered in one area, along the edges of the map, or associated with particular landforms.
Based on the spatial variability of a surface modelled by the shape of the semivariogram2 of that
surface, Brown and Bara (1994) used the estimates ofthe parameters that adequately modelled
the semivariogram to characterize the shape patterns inherent in the DEM surface. These authors
concluded that the fractal dimension of the surface as a function of the log-log transform of the
semivariogram could detect DEM striping from human error in the photogrammetric elevation
point acquisition phase. Systematic error in the DEM is enhanced in the derivatives (e.g. slope
gradient, aspect) . (Brown and Bara, 1994; Walsh et al. , 1987; Lanter and Veregin, 1992). Such
observations underscore the need to evaluate the quality and accuracy of the DEM prior to
modelling.

A semivariogram is estimated by the separation distance (lag), and is the average squared
difference in an elevation value between pairs of input sample points separated by the lag.

2

30

As a mathematical function, a DEM is smooth in comparison to the ground it models. Florinsky
(1998) notes that derivatives ofthe DEM are measured from the smooth mathematical function
and not a measurement of the true ground surface. It would not be valid to compare the values or
derivatives (e.g . slope gradient) ofthe mathematical function to measurements taken on the
ground. In order to assess the accuracy of the derivatives (slope gradient, aspect, horizontal and
vertical curvature), Florinsky proposes that the user consider the accuracy of the DEMand the
precision of the calculation technique used to extract the derivative. It is obvious then that
effective interpolation of sample points into the DEM is a critical step that can have profound
impact on subsequent modelling.

2.5.3 Data Filtering for Error Reduction
Much of the literature reviewed included preparation of the DEM for classification by filtering to
remove random and systematic error. The techniques included mean filtering via a moving
window (Jenson and Dorninque, 1988), filtering in the frequency domain (Pellegrini, 1995), and
pit and peak removal (Tarboton et al. , 1991).
It has been previously noted that the process of using a stereo plotter to create the elevation point
data can introduce random and systematic error. A random error point can appear as an extreme
value when a window is passed over the data at any given neighbourhood locatio. Such an
outlier would affect the result of neighbourhood analysis of shape within the neighbourhood.
DEM error is magnified in first and second derivatives of elevation. Brown (1994) notes that
systematic error not evident by error analysis of the DEM can become evident in the derivative
surfaces. However, it should be noted that random and systematic error may fall within
acceptable limits of accuracy for the DEM, as well as within acceptable limits for the scale of the
mapping product based on project objectives.

31

Brown ( 1994) states that averaging filters in the spatial domain can reduce RMSE accuracy
while attempting to reduce systematic error. The choice of error reduction can be made based on
the objectives of the DEM analysis. If the analysis requires an accurate model of the elevation in
aerial extent, as is the case in shape analysis, then the reduction of systematic error should be
examined. Systematic error has been studied using fractal analysis (Brown and Bara, 1994).
Once it is assumed that the fractal analysis has detected systematic error, the authors
recommended mean filtering of the spatial data to remove the banding error. They recommended
a filter window of 5x7 cells aligned in the direction of the banding.
Pits (single cells significantly lower than neighbouring elevations) are an impediment to accurate
modelling of neighbourhood shape and are especially problematic for hydrologic modelling. Pits
can be removed by assessing neighbouring elevations, setting a threshold for acceptable
differences and filling the pit by averaging the neighbouring elevations (Tarboton et al ., 1991).
An erroneous single control point will cause an artifact of areal extent in the DEM in that it will

influence its neighbours during the interpolation phase. An ideal treatment would be either a) a
criteria that evaluates the control points prior to interpolation b) an identification of erroneous
locations in the DEM, then elimination of control points at those locations and a re-interplolation
oftheDEM.
Frequency filtering is an improvement over the drawbacks of spatial domain mean filtering. In
frequency filtering, short spatial wavelengths are defined as noise while longer spatial
wavelengths are defined as representing real features in the data. Pelligrini ( 1995) notes that
filtering in the frequency domain can remove outliers without modifying the entire surface. The
magnitude and frequency of the error can be visualized in the frequency domain.

32

2. 6 Interpolation
Elevation data can be modelled as a mathematical function using one of a variety of spatial
interpolation techniques. Prior to interpolating the elevation data points into a digital elevation
model, an understanding ofthe effects ofthe interpolator on the data should be gained. In order
to gain such an understanding this section will compare point interpolation methods as applied to
a scalar field z = f(x, y) (where z is the elevation and x, yare the Cartesian coordinates).
With the exception of regionalized variable theory, interpolation is based on the assumption that
the function (surface) is smoothly varying. With this basic assumption in mind, interpolation can
fall into the three basic categories (Table 4, summarized from Burrough and McDonnell, 1998).
Table 4: Overview of Interpolation Methods
Method:
Basic
assumption:

Global
a global spatial structure
imposed on
interpolation while
short-range, location
variation assumed to be
random, unstructured
nmse

Local
value at unsampled
point similar to values
measured close by
using a deterministic
model

Interpolation
Technique:

• trend surface

•

•

•

•
•

analysis
classification
regressiOn
spectral analysis

linear and inverse
distance weighting
thin plate splines

Geostatistical
uses regionalized
variable theory and
optimal
interpolation
aspects of spatial
variation when
estimating
interpolation points
• kriging

In the case of high quality elevation point data characterized by minimal error, the user can
assume that the majority of local, short-range variation is not noise but a meaningful model of
the terrain. This discussion will focus on deterministic and stochastic models of local and
geostatistical interpolators.

33

2.6.1 Deterministic Models
Burrough and McDonnell (1998) summarize the process oflocal interpolation as follows :
•

defining a neighbourhood around the point to be predicted, or establishing the number of
sample points to be included in the prediction;

•

finding the data points within this neighbourhood;

•

choosing a mathematical function to represent the variation over this limited number of
points;

•

evaluating the function for the point on a regular grid.

Linear Interpolation and Inverse Distance Weighting
Linear interpolation (bilinear interpolation in a map plane) is based on the a priori assumption
that elevation values change in a linear fashion (DeMers, 1997). An estimation ofthe value of
an attribute at unsampled points is based on the measurements made at surrounding sampled
points.
In the inverse-distance method, the elevation at a point p, Zp , is interpolated using

Equation 5

where /i is the pre-determined distance within which the control point and the grid point fall
(Unwin, 1981).
It is possible for the operator to change the resolution of the grid, altern (the number of
neighbours) and alter the weighting of 1;. The cell resolution to be used is not based on any
numerical rules and is a function of operator choice, usually chosen to approximate the average
point distance. While offering no solution to the selection of cell size, Zevenbergen and Thome

34

(1987) note "the choice of an appropriate grid mesh distance is clearly an important
consideration" (p 55). While Clarke (1995) suggests statistics such as nearest neighbour, mean,
minimum and maximum point separation, can support the selection of cell size, Lam ( 1983)
notes that the nearest neighbour statistic will give the same result for specific examples of
random and non-random patterns.
Increasing the number of neighbours would have the effect of increasing the search distance
resulting in a larger denominator in equation one and a smoother map .

If li is weighted as

J: then values of b greater than one decrease the relative effect of distant
1

points, resulting in a bumpier surface. Values less than one increase the importance of distant
points and result in a smoother map . Commonly, b=2, the inverse distance squared weighting is
used (Unwin, 1981).
Distance li is zero when the point on the grid to be interpolated is coincident with a sample
point. In this situation, the unsmoothed value of the sample point is copied directly to the output
grid. Since the interpolated points are a local average (local being the search radius) and the
sample points are copied directly to the output grid, the surface reaches a maximum and
minimum at the sample points. This sensitivity to sample points can result in a surface with
unnatural bumps and hollows at the sample sites. Another result of using the local average is the
estimated values are never outside the range of the sample point values within the search radius.
(Berry, 1995; Unwin, 1981 ; Burrough and McDonnell, 1998). Sample points may or may not be
coincidental with the vertices of the grid, (i.e., li may not be zero), so the sample values may not
in all cases be copied across to the DEM (Clarke, 1995).

35

Thin Plate Splines

A spline is a polynomial surface used to represent spatial variation smoothly. Splining is based
on the piece-wise fitting of a surface to a neighbourhood of sample points as exactly as possible,
while at the same time ensuring that the transition from one neighbourhood to the next is as
smooth as possible. Spline methods of surface interpolation assume that an approximation
function should pass as close as possible to data points while being as smooth as possible.
Sample points may contain values representing either errors or natural variation such that an
exact spline could result in local artifacts of unnaturally high or low values. The thin plate spline
attempts to overcome this problem by replacing the exact spline surface with a locally smoothed
average. This is shown by:

Equation 6

where z is a sample point at location Xi and epsilon is the associated random error. The spline
function p(x) should pass not too far from the data values where the smoothing spline is the
function

J that minimizes

A<f ) + f w:[f Cx)- Y(x)]

2

Equation 7

I= I

A( J ) represents the smoothness of the function J and the second term represents the

proximity or fidelity of the data. The weights w ~ are chosen to be inversely proportional to the
error vanance;

2

2

Wi =pI Var [ & (xi)]= p/ s i

Equation 8

36

where the value of p reflects the relative importance given by the user to each characteristic of
the smoothing spline. Thin plate splines are used to interpolate digital elevation models when it
is necessary to interpolate large areas quickly and efficiently (Burrough and McDonnell, 1998).
The piece-wise interpolation approach results in quick calculation. Surfaces are very close to the
values being interpolated only as long as the measurement errors associated with the data are
small. In contrast to weighted averages, splines will retain small-scale features . Splines are
aesthetically pleasing and quickly produce a clear overview of the data. The smoothness of
splines means that mathematical derivatives can easily be calculated for direct analysis of surface
geometry and topology. There are no direct estimates of errors associated with spline
interpolation (Burrough and McDonnell, 1998)

2.6.2 Optimal Interpolation
Kriging

Kriging quantifies the spatial variation of a regionalized variable by employing the semi
variogram, which is analogous to the autocorrelation function in time series analysis. The value
(i.e., elevation) at each grid point is then expressed as an unbiased weighted average of nearby
measured elevations. The weights are chosen to minimize the RMS error over the area of the
grid. In this sense, kriging is an "optimum interpolation method". Kriging also permits
estimation of the error at each grid point associated with the variability of the data and the
distances from the grid points to the measured elevations. This does not include, however,
additional error that may be associated with the measured elevations.
In contrast, linear, inverse-distance and spline interpolation methods lack:
•

direct estimates of the accuracy of the predictions

37

•

an a priori method of determining if the best values have been chosen for the weighting
parameters

•

an appropriate method of determining the size of the search neighbourhood

•

a technique to determine the number of points needed to compute the local average

•

a technique to determine the size, orientation and shape of the neighbourhood of interpolation

•

a technique to determine the interpolation weights other than as a simple function of distance

kriging attempts to address these considerations (Burrough and McDonnell, 1998).
Geostatistical methods such as kriging are based on the observation that spatial variation of a
surface can be too irregular to be modelled by a simple, smooth mathematical function.
Regionalized variable theory assumes that the spatial variation of any variable can be expressed
as the sum of three major components. The components are trend, a spatially correlated
component and random noise.
The trend is a structural component having a constant mean. It is possible to eliminate trend also referred to as "drift"- in some kriging methods, (e.g. , universal kriging) . Universal kriging
does assume a deterministic trend (Lam, 1983). The choice of method is based on an evaluation
ofthe semivariogram. The second major component is a spatially correlated component known
as the variation of the regionalized variable. Random noise is the third major component and is
considered spatially uncorrelated or a residual error term. Kriging assumes that the measured
data points are exact. By computing a semivariogram, however, it is possible to estimate the
random noise in the data from the "nugget effect" (at zero lag) .
In the case of the assumption of a trend, these three components are described at a surface point
location, in this case known as X where the value of a random variable z at x is given by:

Z(x) = m (x) +

i (x) + & .

38

Equation 9

where m(x) is a deterministic function describing the structural component ofZ at x (the trend
and not the elevation value z), E is the stochastic, locally varying but spatially dependent
residuals from m(x) (i.e., the regionalized variable), and E is a residual, spatially independent
2

Gaussian noise term having zero mean and variance C5 . In kriging, h, usually denotes lag.
Logical increments for lag would be multiples of the grid cell size. Therefore, h ranges as 0, la,
2a, ..... na, where a is the grid cell dimension.

In ordinary kriging the drift is not taken into account, although it is not necessarily zero. Where
no drift is assumed, m(x) equals the mean value in the sampling area and the average or expected
difference between any two places x and x+h separated by a distance vector h will be zero.

E[Z(x) - Z(x+h)] = 0

Equation 10

where Z(x), Z(x+h) are the values of regionalized variable Z at locations x, x+h . It is assumed
that the variance of differences depends on h, so that

Yz *E[{Z(x)- Z(x+h)} 2 ] = Y (h)

where

Equation 11

r (h) is the sernivariance at lag h.

The assumption underlying kriging is that the data are stationary (i.e. , that the statistical
properties of the data are constant spatially). While hardly any data are strictly stationary, it is
sufficient that the data is nearly stationary. Semivariance can be estimated from the sample data:

39

1

r (h)= -L:
{z(x)- z(x i + h)}
2n
n

2

Equation 12

i= l

Equation 12 is the equation for semivariance, where n is the number of points of variable z
separated by distance h. A plot of semivariance against lag is known as the experimental
semivariogram. Figure 6 shows an example of an experimental variogram. The curve fitted
through the derived data points can be described by the sill, range and nugget. At large values
of the lag, h, the curve levels off at the sill. Values of the lag equal or greater than the range
exhibit no spatial dependence.
The curve reaches the sill at some value of y (h) known as the range. The range is the maximum
lag for which the data can be considered spatially correlated. Within the range, the smaller the
lag, the more similar the values of the regionalized variable (i.e., elevation) are likely to be. The
range can be used to define the neighbourhood of interpolation since points separated by greater
distances (lags) are uncorrelated. Measured elevations at distances father than the range can,
however, be used in the kriging equations for given grid points, if necessary. They are
"downweighted" because oftheir great distance.
The model variogram intersects the y-axis at a positive value of h representing the "nugget
effect" and is an estimate of E . For a regionalized variable characterized by noise, the nugget
effect is the square of the standard error of a large number of repeat measurements at the same
geographic location. That is, it is the square of the noise level in the data. It must be inferred by
extrapolating the variogram to zero lag, and cannot be computed directly.
Abbass et al., (1990) compared surface fitting algorithms for geophysical data and concluded
that no one method of interpolation is universal. The authors felt that a choice of interpolator

40

should be made with respect to the distribution of points and the expected quality of the resulting
surface.

2. 7 Interpolation Implementation and Accuracy Assessment
Varekamp et al. (1996) recommended that a spatial interpolation be carried out in the following
steps:
•

exploratory data analyses

•

parameterization

•

validation

•

interpolation and display.

During exploratory data analyses outliers are identified. Varekamp et al. recommend using the
coefficient of variation which is a dimensionless quantity measuring the amount of variation with
respect to the mean. Global outliers can then be removed from the data by examining the
histogram and determining the value of any outliers that can then be removed from the data set.
Local outliers would require filtering of the file of control points prior to interpolation. Some
assumption of maximum difference between neighbouring points would need to be made and
any points exceeding this value can be removed from the file.
The parameterization step involves determining the parameters required for implementing the
interpolation procedure. As already discussed this involves selecting an appropriate cell
resolution, neighbourhood size, weighting value or the spline smoothing parameter, for
deterministic models, or determining the initial lag and computing an experimental variogram in
order to select the best model variogram for kriging.
The validation step involves estimating the accuracy of the interpolation procedure. A second
set of sample points not used in the interpolation is compared to the surface values. The

41

interpolated surface is then checked against the sample data points kept aside for accuracy
assessment using, for example, root mean square. The best interpolator is then selected.
In the interpolation stage, the best interpolation method is implemented with optimal parameter
values .. The results are then displayed for visual inspection.
Desmet (1997) reviewed the literature on error in DEMs and notes that the effects of density and
distribution of sample points, data sampling techniques and terrain characteristics are covered to
some extent, but information about interpolation error is sparse.
Data management and analysis include (Dikau, 1989):
•

storage of digitized spatial data;

•

evaluation of these data for the applied purpose;

•

simplification of the mapping process by modelling;

•

modelling of geomorphic processes based on quantified landforms.

Assessing the accuracy of aDEM and its derivatives involves an understanding of the effect of
the interpolator on the elevation data. Desmet (1992) concluded that a spline interpolator
produced the most accurate interpolation. He notes that his results contradict those found in the
literature and that his results are a product of the smoothness of the actual land surface for which
his digital elevation model was built.
Cressie's (1993) review ofthe literature showed that kriging did as well or better than other
techniques. The smoothing effects of splines may not be desirable in the case of high quality,
high-resolution elevation data as detail may be lost. The smoothing effect may be aesthetically
pleasing but may not reflect the underlying landscape. Splines also require the user to select
breakpoint positions in the absence of unbiased guidelines. The arrangment of the input elevation

42

data can profoundly influence the interpolate. It is desirable that the points be as randomly
arranged as possible and include points places at important features . An inverse distance
weighting interpolator produces a local average within the search radius resulting in interpolated
points always lying within the range of the surrounding data. The output is influenced by the
values of the weights and the search radius. The weights and the search radius are user-defined,
lacking an unbiased method of selection. Kriging addresses most of the weakness of the
previous interpolators in that it provides an a priori method of selecting the optimum weights; the
size of the search neighbourhoods, and the nature of the error. In order to use kriging, the spatial
variation of the data must be quantified by computing a variogram. Local and regional trends
within elevation data can exist, and can be taken into account using "universal kriging" (UK).
UK can cause problems, however, if the data are widely spaced or if extrapolating beyond the
data domain, because the interpolated (or extrapolated) points do not necessarily lie within the
elevation range of the data.
The effect of the selected interpolator on the data can be managed using the implementation
procedure of Varecamp et al. (1996), in which the input data quality, error removal, subsequent
selection of interpolator, (including parameters), validation and final interpolation and display
are carried out. This step-wise approach to managing the selection and use of an interpolator
allows the user to consider the influence of the interpolator on the data with respect to a control
set of elevation data points.

43

Chapter 3 Methods
3.0 Introduction
Two steps are involved in computer-automated mapping oftopographic features :
•

creation of aDEM, and

•

discrimination of topographic features from the DEM.

Inaccuracies in the mapping product can be interpreted as the result of either DEM error, or poor
discrimination techniques. Both random and systematic errors should be accounted for. As
previously discussed, error in geographic information system (GIS) data layers will propagate
through the analysis steps into subsequent data layers and be magnified in some cases.
These two general steps are executed as follows :
Creation of a DEM:
•

define required accuracy level

•

assess the input points to identify outliers representing "blunders"

•

delete a random set of points from the control points to be used for accuracy assessment

•

grid the elevation data

•

assess the gridded data through a) visual comparison with airphotos and shaded relief of
DEM and b) statistical comparison with the assessment control points.

Discrimination of topographic features from a DEM:
•

establish the mapping objectives through a) a statement and definition of the topographic
features to be used, and b) the map scale

•

establish discrimination techniques based on neighbourhood analysis and relief, slope
gradient, aspect, slope normal, average curvature, plan curvature and profile curvature

44

•

test the neighbourhood analysis with respect to neighbourhood size, DEM filtering, and a
hierarchy of discrimination

•

grade the results using accuracy assessment.

The following software was used : ESRI, Redlands California, ARC Version 7.2.1 and GRID
Version 7.2.1; Dr. J. Wood's LANDSERF; Dr. R. Eastman's IDRISI32, Clark University, MA,
USA and PCI Enterprises, EASIPACEIIMAGEWORKS, Version 6.3, Richmond Hill, Canada.

3. 1 Study Area
2

The study area encompasses 300 km in the Omineca Mountains of north-central British
Columbia (Figure 1). Nearly seventeen kilometres across each side, the area is centred on NTS
1:20,000 map sheets 093N086, 093N087, 093N076 and 093N077 and can be defined by the
following Universal Transverse Mercator coordinates: 379600 E. , 396700 E. , 6179200 N. ,
6197000 N ., (Zone 10 U, Datum NAD83).

3.2 DEM Generation and Accuracy Assessment
The DEM will be interpolated from a set of elevation points photogramrnetrically captured under
the direction ofBritish Columbia Surveys and Resource Mapping Branch (Surveys and Resource
Mapping Branch, 1990) representing the topographic surface in explicit x, y and z values
(Universal Transverse Mercator- UTM coordinates and elevation in metres) stored in ASCII
text files . The elevation points were captured under the British Columbia Terrain Resource
Inventory Mapping !initiative (TRIM) . TRIM mapping data is available from the Government of
British Columbia Ministry of Environment Lands and Parks.
The data is spatially partitioned into files according to 1 :20 000 scale map sheet boundaries. The
data accuracy meets the standards of the North American Treaty Organization (NATO) Standard
Agreement (ST ANAG) for maps. The error for the set of elevation points is given as the mean
of the linear standard error in Z where 90 percent of all elevation control points are accurate

45

within 5 m of their true elevation, and linear standard error is greater than or equal to mean
standard error. The error for data interpolated from those points is within 10m of their true
elevation 90% ofthe time (Surveys and Resource Mapping Branch, 1990).
The elevation data set does not contain the original control points as separate entities identifiable
by a unique code. This creates the situation in which I am not able to extract a subset of points
of higher accuracy. Since both elevation point types are grouped together in the data set, the
assumption of vertical accuracy is 10 m. Permanently marked control points (accurate within a
few centimetres) are available in the TRIM data, however, none of the markers are present in the
study area. A check for their location revealed that in order to include 40 of these markers, the
study area would have to be 90 x 90 kilometers in size producing a dataset 3600 x3600 pixels, at
25 m per pixel. This proved prohibitive in terms of computer processing.
A set of 500 easting and northern UTM coordinate numbers were randomly generated. These
coordinates are used to select and remove data points from the input data set. These points are
considered the check points. These points are of equal accuracy, matching the error for the
permanent control points. While a check point set of higher accuracy is desirable, it is also
desirable that the check points be randomly generated. The permanent control points are known
not to be randomly selected as they represent heights of land, prominent ridges, bridges and road
intersections.

3.3 Elevation Data Preparation: Blunder Assessment and Control

Point Removal

For the interpolator to be chosen it must be within 10 m vertical accuracy with respect to the
check points. A second level of evaluation will be made utilizing a shaded relief model of the
DEM. A visual comparison will be made between the airphotos and the shaded relief to ensure
the absence of artificial forms (artifacts) not evaluated by the check point method.

46

The elevation data are stored as a table, partitioned on an individual map sheet basis (BCGS
1:20,000). The files will be extracted and translated from government file format into ASCII
format. The files will be appended and sorted based on the elevation value. Blunders, as
defined by Survey and Mapping Branch, are errors in excess of 18.3 m (99.73% probability).
While these were removed by Surveys and Resource Mapping Branch, blunders have been
noted to appear in DEMs. The expected minimum and maximum elevation values will be
established for the data based on an examination of the histrogram and any obvious blunders or
zero values will be removed. It is assumed that any outliers to the histogram are blunders. These
values will be confirmed by reviewed the 1:50,000 NTS maps for the expected maximum and
minimum values for the area. The text files will be modified to match the required input format
for ARC/INFO GRID.

3.4 Interpolation Techniques and Accuracy Assessment
Four methods of interpolation will be employed: kriging, spline, linear inverse distance
weighting, linear and quintic interpolation of a Triangulated Irregular Network. The rationale
behind the fourth choice is as follows . The TIN allows for the incorporation of breaklines, which
are linear features used to enforce the TIN to retain important features,( i.e., a river, or a ridge) .
Areal features such as lakes can also be used. The TIN is forced to place the triangle edge along
the river, ridgeline, or lake shorelines. It is then possible to extract the grid DEM from this
surface using a linear or a fifth order polynomial (quintic) interpolator.
The following ARC and GRID commands were used to interpolate the elevation surfaces: ARC
KRIGING and CREATETIN commands and GRID IDW and SPLINE commands. Default
parameters, 25 m cell resolution, and a search radius of 15 points were used where applicable.
The ARC TINLATTICE command was used to convert the resulting TINs into raster format.

47

Once the DEMs have been interpolated, the locational coordinates of the check points will be
used to extract the elevation value for check locations within each of the four DEMs. The error
for each DEM will be assessed using the following equations.

Equation 13

fu =

Li

Z obs- Z int

N

I

Equation 14
Equation 15

All three equations refer to the check point data set and the corresponding interpolated elevation
points. In these three equations, z obs is an elevation point in the check point data set (observed),

. is an interpolated elevation point, N is the number of data points, & is the mean absolute
Z mt
difference between observed and interpolated points in the domain, and p is the standard
deviation.

3.5 Topographic Discrimination
3.5.1 Mapping Objectives: Selecting Topographic Features and Map Scale
As referenced by the Terrestrial Ecosystem Mapping manual, the bioterrain units found in
Vegetation Resources Inventory Photo Interpretation Procedures (Resource Inventory Branch,
1999) are defined in Table 1 and are the features to be discriminated from the DEM.
Valley line (gully) will be added to this list as it is the complement to ridge and is expected to be
useful in deriving slope units.
While a fully automated approach is desirable, based on the literature review it is expected that a
method by which the features will be discriminated in an heuristic and hierarchical manner will

48

be required in order to compliment the mapping standards as outlined by Resource Inventory
Branch, (1999). In keeping with the standards, the initial mapping phase involves the derivation
ofbroad landscape units . These include alpine, mountain slopes, and valley floors . In the next
step, the topography is described in terms of landscape position and shape, slope, and aspect.
The features will be mapped to 1:20,000 and, according to Terrestrial Ecosystem Mapping
standards, this requires a minimum polygon size of2 hectares (ha). Therefore, with a cell size of
25 m, the minimum polygon size is 32 cells.
The goal of the thesis is to produce a method that discriminates shape from aDEM that has the
flexibility to meet stated project objectives. In this context, Terrestrial Ecosystem Mapping
dictates landforms to be extracted. The scale has been stated as 1:20,000, which controls
minimum class size.

3. 6 Parameterization
Can this method work at all map scales and at all landform sizes? The dilatation wavelets of
Chui ( 1992) and conic sections of Wood ( 1998) address the issue that a specific landform can
exist in varying sizes across the landscape. Pellengrini (1995) filtered the DEM until it matched
an ideal fixed-sized filter with which to extract a fitted surface.
The method selected here employs a neighbourhood window that will vary in size. The rationale
is that the map units to be discriminated are assumed to exhibit a strong degree of homogeneity.
A filter window will evaluate in a local neighbourhood the nature of the central cell with respect
to its immediate neighbours. This creates a basic form element that meets the desired criteria
(e.g. convex), which is taken to be an attribute of a specific legend class. Spatially, these
elements will aggregate to one degree or another into contiguous areas that meet the minimum
class size based on the map scale. If the aggregation does not meet the minimum class size, then

49

the unit is dissolved into the class that forms the greatest enclosing neighbour. Such an approach
meets Terrestrial Ecosystem Mapping standards of ignoring local irregularities and basic
mapping techniques of generalization to match a given scale.
The fundamental approach incorporates Tomlin's (1990) neighbourhood analysis and employs a
floating window to evaluate each DEM cell based on its neighbourhood values. The roving
neighbourhood window will require a standard notation system. It will be based on the grid cell
notation of the data layer known as G. Each cell is assigned a unique location Gi,j where i=
1,2,3 ..... ,m; j=1,2,3, ...,n and i, j represent the column and row address with respect to the central
cell of the window m x m cells in size.
The parameters to be tested in order to discriminate the surface forms from Table 1 are relief,
slope, aspect, slope normal, average curvature, plan curvature and profile curvature. The slope
gradient (8), aspect(¢), and curvature (K) will be derived from the DEM using the following
equations based on the partial quartic equation cited in Zevenbergen and Thorne (1987), and
derived by Evans (1972). They employ a finite difference approach in order to find the
derviatives of the raster surface by fitting a local quadratic surface to each point. The full
quadratic surface is given in Equation 16.

Equation 16
where Z is the elevation and A through I are surface coefficients and X, Y are coordinates with
respect to Z 5 (the central point). Given the neighbourhood window:
Zt

z2 z3

z4

Zs

z6

Z1

Zs

Z9

50

these nine elevations are used to solve Equation 16. The coefficents for the surface function in
Equation 16 are solved by:

•

'

.. '

;

·..

H = (Z2 - Z&) I 2L

I= Zs
where Lis the distance of the cell size (25 m) .
Slope Gradient= G cos()+ H sin() , at the origin cos() = -G I (G2 + H 2 ) 0 ·5

Equation 17
By convention, the negative sign in Equation 17 is ignored.
Aspect ( ¢) = arctan (-H I -G)

Equation 18

Profile curavture is solved by employing the partial quadratic equation:

51

PrC = - 2(DG 2 + EH 2 +FGH)/ (G 2 + H 2 ) Equation 19
Plan curvature is solved employing the partial quadratic equation by:
PIC= 2(DH 2 +EG 2 -FGH)I(G 2 +H 2 )

Equation 20

The values for curvature are given as radians per 1/100 metre where positive values are concave
and negative values are convex. While Equations 19 and 20 do not give true curvature,
Zevenbergen and Thome ( 1987) consider them to be directional derivatives that measure the
same topographic properties as plan and profile curvature.
Slope and aspect will be used to calculate the directional cosines that will describe each surface
normal. The first step in calculating a surface normal is to determine the rectangular coordinates
based on spherical coordinates p, 8 , ¢ (unit vector length 1, slope (in radians), and aspect (in
radians), respectively). Figure 3 illustrates the angular descfiption of the slope normal with
respect to orthagonal axes (James, 1994).
y= p sin ¢ cos

e

X=

p sin ¢ sin 8

z = p cos 8 Equation 21

In which x, y, z are the rectangular coordinates on an orthogonal set of axes used to define a
slope normal with its origin at 0,0. x, yare axes coincident with Cartesian coordinates and are
oriented along east and north respectively. The z-axis is coincident with the vertical component
of elevation (Figure 3 ).
For a sample of slope normals within a neighbourhood window, a resultant normal, Rh , of the
sample can be determined by:

Equation 22

52

where Rh is the resultant normal for the sample, and Xi, Yi and Zi are the z, y and z components of
the normal vectors from Equation 21 .

An estimate of the hemispherical variance, Sh , of the sample of slope normals around Rh can be
determined by:

Equation 23

where N is the number of samples within the neighbourhood window (Hodgson and Gaile, 1999)

An estimate of the hemispherical standard deviation, s h , of the sample of slope normals around
Rh can be determined by:

Equation 24

(Hodgson and Gaile, 1999)
A resultant normal ~ will be used to characterize the degree of dispersion of a sample of slope
normals with Equations 23 and 24 .
The results of the calculations of slope gradient, aspect, relief and slope normal will be placed
into GIS layers and tested in various combinations in order to determine their suitability for
discriminating the bioterrain units identified in Table 1. Thresholds based on the mean, variance
and standard deviation for topographic legend elements will be interactively tested .

3.8 Accuracy Assessment
The resulting classification will be tested against reference ground truth derived from airphoto
interpretation of the area. The two data sets will be tabulated as normalized frequencies in a co-

53

occurrence matrix. Each map class will be evaluated for errors of omission and commission. It
is theorized that the classification differs from the observed values for a given cell in the matrix
due to chance sampling (Williams, 1984). Therefore, a KAPPA coefficient will be used to
determine if the correct results of the classification are a result of chance agreement by removing
the proportion of accuracy due to chance agreement.

Kappa= (Po- Pe)/(1-pe)

Equation 25

Where Po is the observed proportion of agreement (overall accuracy) and Pe is the expected
proportion of agreement due to chance.

L p (i) * p /i)
n

Pe =

r

,

Equation 26

i= l

where n is the number of classes , r represents the row marginals and c represents the column
marginals (Foody, 1992). The value for KAPPA ranges from 0 to 1, and is usually expressed as
a percentage.

;·

54

Chapter 4 Results and Discussion
4. 1 DEM Interpolation and Error Assessment
Interpolators run on the raw elevation points were universal kriging, linear and quintic TIN,
spline and inverse weighted distance. Initial attempts to derive a semivariogram from the data
produced the curve shown in Figure 7. The parabolic curve is indicative of non-stationary data
(contains local trends) . Universal kriging, which does not assume local stationarity, was used as
the kriging interpolator both with linear and quintic drift.
The RMSE results of these interpolations found by comparing the 500 randomly removed
control points with the DEM surface are shown in Table 5. The results for Equation 13, RMSE,
for each of the tested interpolators are given in the second column. The results for Equation 14,
mean absolute difference for each of the interpolators tested is given in the third column.
Table 5: Interpolation Root Mean Square Error
Interpolator for DEM:
Universal Kriging- linear
Universal Kriging- quintic
Linear Interpolated from TIN
Quintic Interpolated from TIN
Spline
Inverse Weighted Distance

RMSE{m)
5.6
5.8
15 .3
12.5
5.4
7.6

Mean AbsoluteDifTerence (m)
3.9
4.2
14.3
11.3
3.9
5.4

An illustration of the randomly generated check points is shown in Figure 8a. The error
residuals calculated at these points are shown in Figure 8b. The discrete, random error points are
shown in the figure as an interpolated surface in order to enhance visualization. Only the largest
of the error residuals are seen. Figure 8b illustrates the magnitude of the residuals across the
DEM. Since a TIN is the single most commonly used interpolator for TRIM data in British
Columbia, and kriging is highly recommended in the literature as the most accurate interpolator,
these two are chosen for illustration. RMSE does not display the spatial distribution of random

55

error. The three dimensional representation of residual error shown in Figure 8b illustrates that
the highest residuals for TIN are along the edge of the map extent. This is due to the inability of
the convex hull to constrain the interpolation. Secondly, when the residuals are stratified by
slope gradient class, regardless of the high spikes along the map edge, high residuals occur in
moderate and steep slope areas. Table 6 lists the RMSE by slope class, while Table 7 lists the
variation with respect to the mean residual by slope class.
Table 6: Root Mean Square Error Stratified by Slope Class.
Slope
Classes
flat
gentle
moderate
moderate
steep
steep

0

/o

iwd

tin I

tinq

spline

krigl

krigq

0-5
5-26
26-49
49-70

3.2
6.1
7.9
10.7

7.5
20.8
9.0
9.5

8.0
15 .7
9.1
8.9

2.5
4.0
5.7
7.3

2.5
4.4
6.1
7.4

2.7
4.2
6.2
7.2

70+

15 .0

12.0 13 .1

12.0

11.8

13 .8

Table 7: Variation with Respect to the Mean Residual Stratified by Slope Class.
Slope
Classes
flat
gentle
moderate
moderate
steep
steep

%

iwd

tin I

tinq

spline

krigl

krigq

0-5
5-26
26-49
49-70

4.0
4.0
5.2
7.3

7.3
20.0
6.7
6.7

7.4
14.9
6.8
6.3

3.3
2.8
3.7
4.8

3.2
2.8
3.8
4.9

3.2
2.8
3.8
4.8

70+

11.7

7.7

9.7

9.0

8.9

11.0

The high residuals along the TIN map are not apparent for the kriged surface shown in Figure 8b.
The high error points for TIN along the edge of the map are removed by clipping the DEM to a
slightly smaller extent that excluded the edge errors, and the RMSE is recalculated and listed in
Table 8.

56

Table 8: Interpolation Results
Interpolator for DEM:
Universal Kriging- linear
Universal Kriging- quintic
Linear Interpolated from TIN
Quintic Interpolated from TIN
Spline
Inverse Weighted Distance

RMSE (m)
5.6
5.8
7.1
7.0
5.4
7.5

Mean AbsoluteDifference (m)
3.9
4.2
5.4
5.2
3.9
5.4

All interpolators are within the stated limit ofRMSE <10m. Kriging and splining display the
two lowest measures of random error.
RMSE is not sensitive to systematic error in the form of artificial landforms since the basic
assumption of RMSE is that the error it is testing is random. Among the bioterrain units to be
discriminate are ridges. A second set of points is selected to focus on this bioterrain unit. The
points removed from the data set were selected by Geodata B.C. for illustrating significant
landforms (i.e., ridges and peaks), for cartographic output purposes. The reason for this step is
that one of the stated advantages of a TIN model is the introduction of significant landform
features via the breakline.
Is the TIN better than kriging for representing ridges? These non-random points were also coincident with the TRIM breakline data. The results of the non-random ridge and hill points show
that RMSE for kriging was 6 m. For a raster interpolated from a TIN that incorporated
breaklines, the RMSE was 4 m, and when breaklines were not used, the RMSE for the raster
from TIN was 6 m. In this case, neither interpolator was significantly superior. Assuming the
breakline data are accurate, its use in the construction of a TIN is superior in retaining an
accurate representation of the breakline features .

57

An inspection was made of the DEM using shaded relief in order to identify visual signs of
systematic error (Figure 9). Visual comparisons were made between the shaded relief and the
stereo pair of airphotos for the area. As Figure 9 illustrates, pronounced striping in the form of
parallel ridges is seen in the lower left quadrant of the DEM. These are examples of the
systematic error not revealed by the RMSE values. It is not possible to remove the stripping
directly from the raster surface, as filtering will alter adjacent valid data as well. For this reason,
any attemp to remove the stripping should be executing on the raw elevation point files prior to
DEM interpolation.
Figure 10 illustrates a subsection ofthe airphoto centred on a mountain valley. This valley is
representative of the mountainous terrain from forested valley to alpine ridges and peaks. Figure
11 illustrates the same valley as contour lines. Subsequent discussion focuses on a planimetric
and perspective view of the shaded relief of this valley. The legend key in Figure 11 shows the
field of view and location of observer for the perspective views.
The perspective view (on the left of the figure) and shaded relief (on the right of the figure)
samples are given in Figure 12 a-d. The inverse weighted distance produced strong systematic
error in the form of regular bumps (Figure 12 c). The rasters extracted from TINs in Figure 12a
and b show bright mounds across the slopes that are not apparent in the airphoto. The kriged
surfaces and the splined surfaces illustrate the most realistic surfaces of the data set based on a
visual assessment of shaded relief and a comparison with stereo airphotos.
The transect 1 in Figure 11 was used to sample elevation and slope gradient data from the
various surfaces. Figure 13 compared the results in graph format. Mid-slope, as seen in Graph
B, the elevation values for all surfaces are similar. At a slope crest, as shown in Figure 13,
Graph A, elevation differences are noticed. Graph B illustrates a mid-slope section of the

58

elevation profile in which the surfaces are graphed as being close in value. Graph C illustrates a
portion of the sampling transect which passes over an area of the TIN containing breaklines.
TIN surfaces emphasized features constrained by breaklines.
Since the literature states that elevation derivatives emphasize error in the DEM, transect 1 from
Figure 11 was next used to sample the slope gradients derived from each surface. The results are
shown in Figure 14 a-b . Ideally, the slope transect would be compared to slope values acquired
in a more accurate manner. This has been shown to be a difficult process to undertake since
some authors find it unacceptable to compare slope gradient derived from a mathematical
function to slope gradient derived from the ground surface. It would be expected that such data
would have some acceptable variation with respect to the slope gradient derived from the DEM,
and it is not known how to reasonably compare the two in a precise fashion. The line transect for
slope gradient in Figure 14 a and 14 b illustrates that elevation interpolation can significantly
effect the derivatives, in this case slope gradient.
The kriged DEM with linear drift was selected as the DEM. It has a low RMSE and the lowest
systematic error based on visual inspection.

4.2 Future Research for DEM Interpolation
Future research can focus on quantifying random and systematic errors and their effective
removal and/or reduction from the DEM. A more thorough understanding ofthe effect of a
DEM interpolator on derivatives and, in turn, the impact on project goals would be a possible
research undertaking. In the area of neighbourhood shape analysis, it may prove significant to
understand how elevation derivatives are affected by changes in surface form as produced by the
interpolator and by attempts at reducing random and systematic error.

59

These considerations would need to be taken in context with the intended use of the DEM.
Viewshed (line of site) analysis, hydrologic analysis and shape analysis are three examples of
DEM implementations that may prove to have individual requirements with respect to accuracy
assessment of elevation and its derivatives.

4.3 Feature Extraction
The parameters of curvature, slope, and aspect are iteratively classified employing varying
thresholds, which were assessed on-screen in order to compare their merits to discriminate shape.
The comparison was made with airphotos and to a simultaneous cursor position on a shaded
relief model, which was considered the visual analogue of shape.
Is it possible that pre-existing methods could be used to extract ridges and valley lines? This
question is based on the availability of methods used in hydrologic watershed analysis. As noted
in the literature review, part of the process of creating watersheds involves extracting raster cells
that are cells of zero flow accumulation (cells which are of higher elevation than their
neighbours) . Figure 15 illustrates cells of zero flow accumulation in white. The study area is
shown as a bounding white rectangle. Based on visual inspection, this output has no advantage
over relief and curvature as a tool for ridge and valley line extraction. All three discriminate
ridge and valley lines of various sizes and at all elevations. It may be desirable to be able to
discriminate the top ridge as separate from a side slope ridge. The zero accumulation surface is
part of a suite of output for watershed analysis that is computationally intensive and could prove
prohibitive for processing large areas.

4.3.1 Curvature
Transect 2, shown in Figure 11 , is used to sample profile curvature at 25m intervals. Would it be
possible to recognize the ridge and slope positions (especially toe slope) from the values of
curvature? Figure 16 illustrates the results . The top graph illustrates the elevation along the

60

transect while the bottom graph illustrates profile curvature along the transect. The only
identifiable signature in the bottom graph is for ridge and valley (in this deeply-incised mountain
terrain) . Wood (1995) and Pellengrini (1995) were also able to extract ridges and valley lines
employing their techniques, as can watershed analysis protocals (e.g. Jenson and Domingue
(1988).
Three types of curvature were reclassified according to convexity and concavity and the results
for the valley of focus are shown in Figure 17 a-c. The light gray indicates areas that are convex,
and darker gray areas that are concave and the black areas in which the neighbouring cells are
equally convex and concave or flat. Average and plan curvature both capture ridges and valley
lines. Figure 18 illustrates how a curvature map enhances systematic error in aDEM. The area
of focus is shown in the figure as a black rectangle. The north-south ridges evident in the shaded
relief of the DEM are seen in the curvature map illustrated in the figure.
A further consideration is made based on computational efficiency. Curvature derivation creates
raster layers stored in 32-bit format. To compare efficiency, an 8-bit format was assessed for its
ability to discriminate ridges and valley lines in comparison to curvature. Appendix II, Filter 1,
contains the algorithm for the filter used in this step. After applying the filter, a value greater
than zero represents a cell in a neighbourhood that is concave and a value that is negative
represents an area that is convex. The results, shown in Figure 19, were visually comparable to
the convexity layers, and can be stored as an 8-bit layer in comparison to a 32-bit curvature
layer.

4.3.2 Slope Gradient and Aspect
In a landscape such as the one used for this thesis, the highly varying relief and associated
structure result in a slope class map containing much detail. Minor perturbations can be removed
employing either a modal filter or dissolving criteria based on polygon GIS practices (available

61

in ARC/INFO and other commercial GIS) to match the desired scale. The criteria by which the
contiguous cell values (analogous to polygons) are dissolves are discussed. A threshold is
selected for the minimum polygon size for the resulting map. When a polygon is located below
the threshold size, it is assigned the value of the largest adjacent neighbouring class. The results
are shown in Figure 20. It is reasonable in this landscape that slope gradient values over
extended areas would increase up the mountain sides. Lower slopes are found on the mountains
ridges but are not seen in the scale ofFigure 20.
Ridge crests are coincident with areas of maximum convexity in terms of profile curvature. The
ridge descends down-valley and retains profile curvature convexity along its sides. Along-ridge
shape may, at the same time, be concave. Along the ridge, landscape areas as narrow as single
cells may cover areas of zero or low slope gradient. However, the slope gradient derivative is
calculated from a 3x3 cell neighbourhood and such small landscape areas are not discernable in
the DEM. The slope gradient value assigned to any given cell is the gradient determined from
eight neighbouring cells. Zero or low slope gradient cells do not, as a result, consistently
coincide with the ridge line in this case. Further investigation of the DEM over a greater area
also indicated that low relief, broad areas several hundred metres in extent, are found in high
elevation, alpine areas. This also hinders the effectiveness of using low slope values at high
elevations to discriminate ridges.
This pattern may not hold in other landscape areas. Other areas may exhibit low relief that in
turn is generally smooth, or may exhibit rapid changes in slope gradient. Such variation supports
the necessity that topographic analysis be considered in context with knowledge of the landscape
in question and a clear vision of the intended map use (in terms ofboth scale and legend items).
Due to the underlying geologic structure that exerts a significant control on the local topography,
the slope gradients in the mountains of this area may be assessed as behaving differently from

62

another area. The clean grouping of slope classes may not repeat in other areas. This
underscores the need for knowledge of the area to be mapped and, in terms of automating
mapping support, the need for a suite of tools that is flexible and adjustable, based on an overall
assessment of the area in question.
These slope groupings can be interpreted in context with slope position. Upper mountain
positions include high mountain shoulders, steep rock slopes and exposed ridges. Mid-slope
positions include steep to moderately steep slopes and steep rock walls with talus along the base.
Mid-slope positions cover large areas. Lower slope positions include a significant break in slope
gradient that falls to moderate gradient. Lower slopes are areas of deposition. These
descriptions, along with their respective ecological implications, illustrate the complexity of a
straightforward discrimination of topographic form . The exceptions can make the discrimination
more challenging. For example, even though the lower slopes are moderate in gradient and are
areas of accumulation, they also could contain canyons, cliffs, and deeply-incised river channels.
These examples of local extreme relief, in an area otherwise characterized as having moderate
slope gradient, illustrate the complexity in classifYing form.
Table 1 defines slope gradient groupings as being unidirectional (planar) with a smooth
longitudinal profile that is straight, concave or convex, and local surface irregularities generally
less than one metre relief This description describes an area of generally similar aspect. To
evaluate the possibility of groupings being made based on the stated criteria a tool was required
to quantifY slope gradient and aspect together.
The three data layers given in equation 21 were derived that describe the slope normal. The
resultant vector R - a measure of the magnitude of the sum of normal vectors within a given
neighbourhood - was derived (illustrated in Figure 21 ). Visual inspection of the magnitude of R

63

indicates that it is analogous to roughness. High values indicate areas of systematic elevation
change (flat to relatively smooth), and low values indicate elevation change that is chaotic or
abrupt. For example, ridges and valley lines have the same range of values for R. One possible
use ofR may be quantification of planar slope units which form part of the bioterrain unit
descriptions listed in Table 1. The units are not expected to be truly planar, but to be relatively
smooth and to fall within relatively narrow aspect ranges.
The first method stratified the study area by overall relief: mountains, moderate relief and low
relief This was achieved by testing a relief analysis neighbourhood of increasing size. Relief is
the difference between the maximum elevation value in the analysis window and the minimum
elevation value in the analysis window. The relief output from each filter was interactively
queried to determine the relief value at the break in slope between mountainous and nonmountainous areas. This value was then used as a threshold to classify the relief layer into a
binary map that is mountainous and non-mountainous. This map is then compared to the shaded
relief model and the airphoto to assess the results. It was determined that a smaller filter of size
is required to discriminate medium relief and low relief (in non-mountainous areas) . Successful
relief filtering requires an iterative and hierarchical approach. Medium/low relief discrimination
requires a different window size than mountainous/non-mountainous. The relief threshold value
will be unique for each landscape. It was necessary to derive mountainous and non-mountainous
relief first. Once mountainous areas were identified, they were masked out of futher filtering.
Only non-mountainous areas were filtered to determine moderate from low relief The results
are shown in Figure 22. This layer can be used to discriminate ridges in mountains versus ridges
in non-mountainous areas.
The second method identified slope length. Cumulative slope length, shown in Figure 23 , was
calculated using RUSLE code written by Hickey et al. (1994), and modified by Van Remortel et

64

al. (2001). This process is computationally intensive: a 500 x 600 raster was processed in four
hours on a workstation with two Gb of memory. Each convergent catchment along the length of
the whole slope contributes to the cumulative slope length. When the values are queried along
the river at the bottom of the slope a value of3000 can be next to a value of 12000, reflecting the
contributions of convergent and/or divergent catchments up-flow from each of the neighbouring
cells. Cumulative slope length has a direct application in soil and moisture flow modelling but
does not appear to be useful in modelling slope position. Figure 23 illustrates the convergence
and divergence that occur along the gradient length. The red linework represents the watersheds
derived using Jensen and Dominigue' s (1988) watershed analysis.

An approach for future research that might prove more effective would be to employ Pike ' s
slope between slope reversals (Pike, 1988). This approach would eliminate the effect of
divergence and convergence contributing to slope length. As a result, the length of the slope
from valley bottom to the ridge is determined as one unit that could be scaled to a specified
length. Once all slope runs from ridge to valley bottom are determined and scaled to an equal
length, then the resulting raster can be partitioned into slope position units.

4.4 Accuracy Assessment
Table 9 contains the error matrix for a preliminary assessment of the accuracy of the ridge and
gully classification. These observed data points were collected through a visual assessment of
the stereo airphotos. A series of points were listed in a digital file and assigned a class value:
gully, slope and ridge. These three classes were selected to illustrate the accuracy assessment
process. The remaining classes listed in Table 1 were not incorporated due to the absence of
independantly acquired ground truth. Correctly classified forms fall on the diagonal. See section
2.5.1 for details on accuracy assessment.

65

Table 9: Error Matrix
Observed
Gully Slope Ridge Row Totals
Predicted
Gully

13

0

0

13

Slope

7

10

3

20

Ridge

0

2

31

33

Column Totals

20

12

34

66

The overall accuracy is the number correct divided by the total number of ground truth points.
Overall accuracy is 80%.
The error of omission (producer's accuracy) and the error of commission (user's accuracy) are
listed in Table 10.
Table 10: Errors of Comission and Omission
Gully

Omission
Commission
65%
100%

Slope

83%

50%

Ridge

94%

94%

Chance agreement is revealed by an analysis of the marginal values in the matrix, namely those
values that are rnisclassified. When the proportion of accuracy due to change agreement is
removed, the overall KAPPA statistic is 70% (the overall accuracy of the classification with
chance agreement removed) . The KAPPA statistic for each class are: gully 100%, slope 39%
and ridge 88%. The low value for the slope class is due to the ambiguity of the break between a
ridge and the slope, and for the gully to the slope transition. This is due to the ground truth
including a broader definition of slope and a narrower definition of ridge and gully (e.g. the ridge
in the automated process extends further down slope than the ground truth).
Based on the maps of ridge and valley line (gully) classification, the areas that are classified as
valley lines also incorporate some gullies up the sides of the steep slopes. This over-

66

classification causes the ground control points that are slopes to include ground points labeled as
the gully class. As a result, the error of commission for the slope class is 50% and, once chance
agreement is removed, the KAPPA statistic for this class is 40%.
Preliminary results indicate that it is difficult to differentiate between valley side gullies and
valley lines. As noted previously, it may be possible for valley lines to be incorporated as a tool
to use in concert with the ridge lines in order to differentiate slopes Also, the ridge class
incorporates the main ridges and subsidiary ridges that descend the valley sides. It is necessary
under certain conditions to differentiate the two types of ridges.
Differentiation of ridges could be carried out employing criteria of identifying local elevation
maxima providing the filter window is of complementary (i.e., appropriate) size. Vegetation
maps could also be employed to differentiate alpine ridges from forested ridges. Further
investigation of more precise differentiation of ridges and valley lines and gullies could be
carried out by incorporating watershed delineation. The drawback of this technique is that the
watershed size is extracted based on a user-determined minimum size. The valley of interest is
subdivided into five watersheds when a threshold size is set to capture two narrow watersheds to
the east (Figure 23). For the purpose of this classification, the valley of focus was considered
one watershed. The watershed boundaries are shown in red. There is no technique known in the
literature that allows for watershed extraction based on variable size.
The slope normals (Figure 21) did not contain information about neighbourhood shape.
However, the resultant normal surface does contain information about roughness. Extensive
research has been done on roughness and texture analysis (Haralick et. al., 1973; Haralick, 1979
and Conners, et. al. , 1980; Wang L. and D.C. He, 1990). It is not known ifthe slope normals
contain information that improves on texture analysis.

67

4.5 Discussion
The accuracy ofDEMs and the derivative algorithms are cited in the literature as being of great
importance in terms of mapping accuracy. Kriging proved in this case to be the most effective
interpolator in terms of random and systematic error. Choosing interpolator parameters, and the
selection of appropriate data spacing and cell size for the interpolators, are all contributing
factors to the effective interpolation of an elevation surface. An understanding of the methods by
which the derivatives are calculated with respect to neighbourhood size and their relation to the
natural landscape surface are significant in terms of effective use of these data in terrestrial
ecosystem mapping.
The methods tested in this study and summarized in Table 11 achieve discrimination of features
of complexity comparable to the complexity of the topographies addressed in other research
initiatives cited in this thesis. These methods are decision-based, so further investigations should
be undertaken in conjunction with specific Terrestrial Ecosystem Mapping objectives.

Table 11: Assessment of Topographic Elements
Bioterrain
Unit
flat (plain)

I Description

METHOD

Modal filter to remove local irregularities; slope classes in
conjunction with roughness (R) to determine planar areas.

gentle
slope

unidirectional (planar), >5%, <=26% (>3°, <=15°), smooth profile either
straight, concave, or convex, local irregularities generally < 1 m relief

METHOD

Modal filter to remove local irregularities; slope classes in
conjunction with roughness (R) to determine planar areas.

level or very gently sloping <=5% (<= 3°), unidirectional (planar), local
surface irregularities generally have a relief < 1m

68

Table 11 Continued
unidirectional (planar), >26 %,<=49%,(>15°,<=26°) smooth longitudinal
moderate
profile either straight, concave or convex; local surface irregularities
slope
generally< 1 metre relief

METHOD

Modal filter to remove local irregularities; slope classes in
conjunction with roughness (R) to determine planar areas.

moderately
steep slope

unidirectional (planar), >49 %,<=70%, (>26° ,<=35 °) smooth, longitudinal
profile that IS either straight, concave or convex, local surface
irregularities generally< 1 m relief

METHOD

Modal filter to remove local irregularities; slope classes in
conjunction with roughness (R) to determine planar areas.

steep slope

unidirectional (planar), >70%, (> 35 °) smooth, longitudinal profile that is
either straight, concave or convex, local surface irregularities generally <
1m relief

METHOD*

Modal filter to remove local irregularities; slope classes in
conjunction with roughness (R) to determine planar areas.

toe slope

area differentiated from the lower slope of the hill , generally concave
surface, specific aspect

METHOD*

Slope between slope reversal I slope length and distance from valley
line; assess for rapid change in profile curvature.
between valley floor and rnidslope, significant break in slope falling to
moderate gradients, generally concave, specific aspect

lower slope

METHOD*
middle slope

Slope between slope reversal I slope length (scaled) divided into 3
slope position classes.
steep to moderately steep consistent slopes and steep rock walls with talus
along the base, involve large areas, specific aspect

METHOD*

Slope between slope reversal I slope length (scaled) divided into 3
slope position classes.

upper slope

directly below steep bedrock summits and above rnidslope, generally
convex below crest, specific aspect

METHOD*

Slope between slope reversal I slope length (scaled) divided into 3
slope position classes.
* These methods are suggested, but were not implemented, and require further research.

69

Table 11 Continued
generally convex uppermost portion of a hill, convex in all directions, no
ridge
distinct aspect, also a ridge
METHOD
depressions

Manipulate thresholds for profile curvature.
circular or irregular area of lower elevation (hollow) than the surrounding
terrain marked by an abrupt break in slope, greater than 2 m in depth,
generally at foot of meso scale hill or generally level area

METHOD
warm aspect

Watershed prepared DEM difference with original DEM, or
watershed pit analysis techniques.
135°-285°, slope>15° (26%)

METHOD
cool aspect

Aspect and slope class groupings.
285° - 135°, slope>15° (26%)

METHOD

Aspect and slope class groupings.

Definitions adapted from Resource Inventory Branch, 1999; Howes and Kenk, 1997; Ryder,
1994; and Luttmerding, et al., 1990. Methods summarized from thesis.

The methods listed in Table 11 illustrate that it is possible to generalize topographic elements
from aDEM, based on Terrestrial Ecosystem Mapping standards. For slope position, the
modification of cumulative slope length to Pike's slope between slope reversal is proposed, but
was not executed. The way in which these methods are integrated into a Predictive Ecosystem
Model will depend on the design of that model. This in tum is dependent upon the project
objectives, map scale, data quality and the resulting legend design. The methods listed in Table
11 constitute a flexible suite of tools that can be used to meet highly diverse mapping objectives.

In Figure 24 the graphic output of the bioterrain raster layers described in Table 11 are shown.
The relief layer is not listed in Table 11 but is useful as the TEM standards recommend a pretyping phase during which the area is stratified into broad landscape classes.
The inability of cumulative slope (Figure 23) to be useful in discriminating slope position and
modification of cumulative slope to Pike's slope between slope reversals as a solution to define
slope position was discussed in section 4.3.2.

70

4.5. 1 Implementations
Figure 25 illustrates another example of the usefulness of these bioterrain data layers in a
vegetation mapping model. Predictive Ecosystem Mapping was developed in part with the
objective of increased use of remote sensing to acquire data layers. For example, once a
vegetation layer is classified from remotely sensed data, the site attributes of each vegetation
object can be determined by a direct query of the bioterrain data layers determined using the
methods discussed in this thesis. The attributes of the polygon resulting from the query are
stored in an attribute database table that is linked to the final digital map. In the event that the
polygon being queried covers multiple classes in any given bioterrain layer, a statistical summary
is made of the resulting classes and the attribute is stored as a compound of those classes.
Figure 26 illustrates a soil moisture regime model (Resource Inventory Branch, 1999) for the
coastal rain forest of British Columbia. The bioterrain units as determined by the methods
presented here can be incorporated into this model as five raster layers: slope position, slope
classes, warm and cool aspect, ridge and valley line and depressions . As indicated in the figure,
these layers are integrated into the binary decisions of the model where circled. As required in a
Predictive Ecosystem Model, other data layers are needed to represent soil, bedrock and moisture
conditions.

4.5.2 Limitations and Usefulness
The usefulness of these bioterrain layers is dependant upon the decisions during the production
phase. This requires a thorough understanding of the landscape being interpreted in order to
make effective decisions regarding the classes and thresholds needed to derive the bioterrain
layers. This digital approach requires sufficient ground truth for each bioterrain class in order for
an effective evaluation of the accuracy of the mapping products. Trained and knowledgeable
landscape interpreters and the acquisition of quality ground truth for accuracy assessment are

71

required in order for effective implementation of these digital bioterrain layers. The absence of
expert knowledge and accurate data can limit the effectiveness of the automated approach . Since
the automated approach will always produce a data layer, it is important for the process to be
assessed using quality ground truth .
In comparing the digital approach to manual mapping, it can be stated that both approaches
require high quality ground truth. In order for the ground truth to be statistically valid, a
stratified random approach to selecting the ground truth sites is required. This can be time
consuming and expensive. Once the data layers are established, then the mapping process is fast,
consistent, updateable and flexible as opposed to human interpreters who introduce variation in
mapping standards, human error and the time-consuming process of analog mapping. The digital
approach is scale flexible in that the maps can be generalized as required. The objects can be
treated as contiguous raster cells or converted to polygon line work allowing for further
flexibility of the map product.
Neighbourhood analysis at a fixed window size can be critiqued in that it may not reflect the
underlying scale of the bioterrain unit. What is demonstrated in this thesis is that in conjunction
with carefully designed project objectives (the identification of the objects -bioterrain units- to
be extracted), neighbourhood analysis is an effective analysis tool. While the window analysis
assigns a single cell to a class, the bioterrain unit is considered to be defined by the aggregation
or contiguity of the cells defined as a specific unit type.
The minimum resolvable unit of the DEM, as controlled by the cell size, is known. In addition,
the mapping criteria (the definition of the bioterrain unit, and its expected minimum size, median
size and maximum size) are controlled by the map scale and the Terrestrial Ecosystem Mapping
specifications. Prior knowledge on the part of the analyst required by Terrestrial Ecosystem

72

Mapping standards is expected for these methods to be effectively used. The methods are
limited when the natural variation of the landscape within each bioterrain unit is not known. In
this case, it is not known whether the DEM, in terms of scale and accuracy, is an effective model
of that landscape variation, nor if the methods would be effective in capturing that variation .
As stated, these methods are effective, in part, because they are easy to determine from a gridded
DEM. This fulfills the criteria of Terrestrial Ecosystem Mapping standards in part, in that the
methods are required to be accessible and non-proprietary. These methods do not fulfill all
possible legend criteria in mapping with respect to all potential landscape features. There is a
sense that such measures have been developed because they are easy to measure from a gridded
DEM, not because they are geomorphologically significant, and they certainly do not pick up
important features . As required, such features as drainage pattern or river profiles can be
extracted based on other criteria and data processing approaches not cited in this thesis .
These methods are flexible and allow for the integration of expert knowledge on the part of the
analyst. These methods can be used in context with the hierarchical approach used in terrestrial
ecosystem mapping. The overall discrimination of the landscape into broad units is supported by
selecting threshold values from a relief layer. This approach allowed for the discrimination of
the landscape into mountainous, non-mountainous, medium relief and low relief areas.
Depending on the project criteria, the analyst may then choose to further discriminate these
broad classes based on the methods listed in Table 11.
It should be noted that the standards that govern Terrestrial Ecosystem Mapping are an initiative

aimed at producing accurate, timely, and cost-effective map products. These automated methods
incorporate expert knowledge at various steps. Once assessments are made for threshold values,
the procedures can process large data sets in a consistent manner as quickly as the computer can

73

manage. This is an advantage over traditional methods in that interpretations between people
vary slightly, and human error can propagate throughout the mapping process.
These processes are dynamic in that it is expected that they can be adjusted for various
landscapes. The processes are programmable and digital as are the resulting layers that lend
themselves to digital revision. This is an advantage over traditional methods which result in
analogue layers which are extremely time-consuming to revise. The traditional mapping
approach is also rigid in that the map product is based on set polygon boundaries. The approach
developed in this thesis result in data layers that are continuously varying from cell to cell. As a
result, the final map product of polygon-analogous raster cell groupings can be adjusted to revise
cell groupings to give up-dated polygon results reflecting changing project objectives or the
incorporation of additional information.
As per the requirements of Terrestrial Ecosystem Mapping and the stated objectives of this
thesis, the methods produced should be easily employed in standard Geographic Information
System (automated, georeferenced, spatial analysis software). The criteria developed by the
committee are open to further definition that will continue to incorporate the developments made
available through research in digital modelling.

74

References
Abbass, T.E., et al. 1990. A Comparison of Surface Fitting Algorithms for Geophysical Data.
Terra Research. Vol2.
Berry, J.K. 1995 . Spatial Resoning for Effective GIS . GIS World Books. Ft. Collins, Colorado.
Blaszczynski, J. S. 1997. Landform Classification with Geographic Information Systems.
Photogrammetric Engineering and Remote Sensing. Vol. 63 . No.2.
Bolstad, P.V. and T.Stowe. 1994. An Evaluation ofDEM Accuracy: Elevation, Slope and
Aspect. Photogrammetric Engineering and Remote Sensing. Vol. 60. No. 11.
Brabyn, L. 1997. Classification of Macro Landforms using GIS. ITC Journal. Vol. 1.
Brassel, K.E. and R. Weibel. 1988. A Review and Conceptual Framework of Automated Map
Generalization. International Journal of Geographical Information Systems. Vol. 2. No . 3.
Brown, D.G. 1994. Anisotropy in Derivative Surfaces as an Indication of Systematic Errors in
DEMS . International Symposium on the Spatial Accuracy ofNatural Resource Data Bases. 1620 May, 1994, Williamsburg, Virginia. American Society for Photogrammetry and Remote
Sensing.
Brown, D.G. and T.J.Bara. 1994. Recognition and Reduction of Systematic Error in Elevation
and Derivative Surfaces from 7.5-Minute DEMs. Photogrammetric Engineering and Remote
Sensing. Vol. 60. No. 2.
Burrough, P. A and R. A McDonnell. 1998. Principles of Geographical Information Systems.
Spatial Information Systems and Geostatistics. Oxford University Press. New York.
Chui, C.K. 1992. An Introduction to Wavelets. Academic Press. Boston.
Clarke, K. C. 1995 . Terrain Analysis. In: Analytical and Computer Cartography. 2nd edition.
Prentice Hall, New Jersey.
Cohen, J. 1960. A Coefficient of Agreement for Nominal Scales. Education and Psychological
Measurement. Vol 20 .

75

Conners R.W. and C. A Harlow. 1980. A Theoretical Comparison ofTexture Algorithms .
IEEE Tr. on Pattern Analysis and Machine Intelligence, Vol PAMI-2, No. 3, May 1980.
Congalton, R.G. and R.A. Mead. 1983 . A Quantitative Method to Test for Consistency and
Correctness in Photointerpretation. Photogrammetric Engineering and Remote Sensing. Vol. 49.
No.1.
Cressie, N .A.C. 1993 . Statistics for Spatial Data. John Wiley & Sons, Inc. New York.
DeMers, M. N . 1997. Fundamentals of Geographic Information Systems. John Wiley &
Inc. New York.

Sons,

DeSmet, P.J.J. 1997. Effects oflnterpolation Errors on the Analysis ofDEMs. Earth Surface
Processes and Landforms. Vol 22.
Dikau, R. 1989. The Application ofDigital ReliefModels in Landform Analysis, In: J. Raper,
(ed.). 1989. Three Dimensional Applications in Geographic Information Systems. Taylor &
Francis, London.
Dikau, R. 1990. Derivatives from Detailed Geoscientific Maps using Computer Methods.
Zietschrift fur Geomorphologie. Suppl.-Bd. 80.
Dymond, J. , et.al. 1995. Automated Mapping ofLand Components from Digital Elevation
Data. Earth Surface Processes and Landforms. Vol. 20.
Ecosystems Working Group . 1998. Standards for Terrestrial Ecosystem Mapping in British
Columbia. Resources Inventory Committee. Province ofBritish Columbia.
Embleton, C. and H. Liedtke. 1990. Geomorphological Mapping, Remote Sensing and Terrain
Models. Zeitschrift fur Geomorphologie. Suppl.-Bd. 80.
Evans, I. S. 1972. General Geomorphology, Derivatives of Altitude, and Descriptive Statistics.
In: R. J. Chorley, (ed.). Spatial Analysis in Geomorphology. British Geomorphological
Research Group. Harper & Row. New York.
Evans, I.S . 1980. An Integrated System ofTerrain Analysis and Slope Mapping. Zietschrift fur
Geomorphologie. Suppl.-Bd. 36.

76

Felicisimo, A.M. 1994. Parametric Statistical Method for Error Detection in Digital Elevation
Models. ISPRS Journal ofPhotogrammetry and Remote Sensing. Vol. 49. No.4.
Florinsky, I.V. 1998. Accuracy ofLocal Topographic Variables Derived from Digital Elevation
Models. International Journal of Geographical Information Science. Vol. 12. No. 1.
Foody, G.M . 1992. On the Compensation for Chance Agreement in Image Classification
Accuracy Assessment. Photogrammetric Engineering and Remote Sensing. Vol. 58 . No.10
Gao, J. 1998. Impact of Sampling Intervals on the Reliability of Topographic Variables Mapped
from Grid DEMs at a Micro-Scale. International Journal of Geographical Information Science.
Vol. 12. No. 8.
Gong, P. 1994. Integrated Analysis of Spatial Data from Multiple Sources: An Overview.
Canadian Journal ofRemote Sensing. Vol. 20. No. 4.
Haralick, R.M, K. Shanmugan and I. Dinstein. 1973 . Textural Features for Image Classification.
IEEE Tr. on Systems, Man, and Cybernetics. Vol SMC-3 , No. 6, November 1973 .
Haralick, R.M. 1979. Statistical and Structural Approaches to Texture. Proceedings ofthe
IEEE, Vol. 67, No. 5, May 1979.
Hickey, R, A. Smith, and P. Jankowski, 1994. Slope length calculations from aDEM within
ARC/INFO GRID: Computers, Environment and Urban Systems. Vol. 18. No 5.
Hobson, R.D. 1972. Surface Roughness in Topography: Quantitative Approach. In: R. J.
Chorley, (ed.). Spatial Analysis in Geomorphology. British Geomorphological Research Group.
Harper & Row. New York.
Hodgson, M.E. 1995 . What Cell Size Does the Computed Slope/Aspect Angle Represent?
Photogrammetric Engineering and Remote Sensing. Vol. 61. No. 5.
Hodgson, M.E. and G.L. Gaile. 1999. A Cartographic Modeling Approach for Surface
Orientation-Related Applications. Photogrammetric Engineering and Remote Sensing. Vol. 65 .
No. 1.

77

Howes, D .E . and E . Kenk (eds.). 1997. Terrain Classification System for British Columbia
(Version 2). Publication 258 . Surveys and Resource Mapping Branch. British Columbia
Ministry of Environment.
Hunter, G.J. and M F. Goodchild. 1995. Dealing with Error in Spatial Databases: A Simple
Case Study. Photogrammetric Engineering and Remote Sensing. Vol. 61. No . 5.
James, S. 1994. Multivariable Calculus. Brooks/Cole Publishing Company. Pacific Grove.
Jensen, J.R. (1996) Thematic Information Extraction: Image Classification. In: J.R. Jensen.
Introduction Digital Image Processing: A Remote Sensing Perspective. 2 nd . edition. Prentice
Hall. Upper Saddle River. New Jersey.
Jenson, S.K. and J.O . Dominigue. 1988. Extracting Topographic Structure from Digital
Elevation Data for Geographic Information System Analysis. Photogrammetric Engineering and
Remote Sensing. Vol. 54. No . 11.
Jenson, S.K. and C.M. Trautwein. Year Unknown. Methods and Applications in Surface
Depression Analysis. U.S . Geological Survey Contract no 14-08-0001-20129.
Jones, K. , D . Medidinger, D . Clark and F. Schultz. 1999. Towards the Establishment of
Predictive Ecosystem Mapping Standards: A White Paper. Terrestrial Ecosystem Mapping
Alternatives Task Force; Resource Inventory Committee.
Kumler, M .P . 1994. An Intensive Comparison of Triangulated Irregular Networks (TINs) and
Digital Elevation Models (DEMs). Cartographica. Vol. 31.
Lam, N.S . 1983 . Spatial Interpolation Methods: A Review. The American Cartographer. Vol.
10. No . 2.
Lanter, D.P. and H. Veregein. 1992. A Research Paradigm for Propagating Error in Layer-Based
GIS . Photogrammetric Engineering & Remote Sensing. Vol. 58 . No . 6.
Lee, J. , P . K. Snyder and P . F. Fisher. 1992. Modeling the Effect of Data Errors on Feature
Extraction from Digital Elevation Models. Photogrammetric Engineering & Remote Sensing.
Vol. 58 . No . 10.
Li, Z. 1988. On the Measure ofDigital Terrain Model Accuracy. Photogrammetric Record, Vol.
12 (72).

78

Li, Z. 1994. A Comparative Study ofthe Accuracy ofDigital Terrain Models (DTMs) Based on
Various Data Models. ISPRS Journal ofPhotogrammetry and Remote Sensing. Vol49. No. 1.
Luttmerding, H.A. et al. , (eds). 1990. Describing Ecosystems in the Field. MOE manual 11. 2nd
Edition. Ministry ofEnviornment, Lands and Parks.
Mark, D.M. 1975 . Geomorphometric Parameters: A Review and Evaluation. Geografiska
Annaler. Vol. 57. No. 3-4.
Mark, D.M. 1984.Automated Detection ofDrainage Networks from Digital Elevation Models.
Cartographica. Vol. 21.
Meidinger D. and J. Pojar, 1991 . Ecosystems of British Columbia. Research Branch, British
Columbia Ministry of Forests. Special Report Series, No. 6.
O'Neill, M.P. and D.M. Mark. 1987. On the Frequency Distribution ofLand Slope. Earth
Surface Processes and Landforms. Vol 12.
Pellegrini, G.J. 1995 . Terrain Shape Classification ofDigital Elevation Models Using
Eigenvectors and Fourier Transforms. Ph.D. State University ofNewYork College of
Environmental Science and Forestry. Syracuse, New York.
Peququet, D.J. 1984. A Conceptual Framework and Comparison of Spatial Data Models.
Cartographica. Vol. 21. No. 4.
Pike, R.J. 1988. The Geometric Signature: Quantifying Landslide- Terrain Types from Digital
Elevation Models. Mathematical Geology. (20) (5).
Resource Inventory Branch. 1999. Vegetation Resources Inventory: Photo Interpreation
Procedures. Ministry afForests. Terrestrial Ecosystem Task Force. Vegetation Resources
Inventory Committee. Version 2.2.
Riazanoff, S. et al. 1988. Ridge and Valley Line Extraction from Digital Terrain Models.
International Journal ofRemote Sensing. Vol. 9. No. 6.
Rosenfield, G.H. and K. Fitzpatrick-Lins, 1986. A Coefficient of Agreement as a Measure of
Thematic Classification Accuracy. Photogrammetric Engineering and Remote Sensing. Vol. 52.

79

Ryder, J.M. 1994. Guidelines and Standards to Terrain Geology Mapping in British Columbia.
Terrain Geology Task Group, Earth Science Task Force, Resource Inventory Committee,
Government ofBritish Columbia.
Schmid-McGibbon, G. 1993 . Landform Mapping Analysis and Classification using Digital
Terrain Models. Ph.D . Thesis. Department of Geography. University of Alberta.
Stocks, A.M. and D. I. Heywood. 1994. Terrain Modelling for Mountains. In Mountain
Environments and Geographic Information Systems. Martin F. Price and D . Ian Heywood, eds.
Taylor and Francis. London.
Story, M. and R.G. Congalton. 1986. Accuracy Assessment: A User's Perspective.
Photogrammetric Engineering and Remote Sensing. Vol. 52. No. 3.
Surveys and Resource Mapping Branch. 1990. British Columbia Specifications and Guidelines
for Geomatics. Content Series Volume 3: Digital Baseline Mapping at 1:20,000. Release 1.0.
Ministry of Crown Lands. Province of British Columbia.
Tarboton, D.G. et al. 1991 . On the Extraction of Channel Networks from Digital Elevation Data.
Hydrologic Processess. Vol. 5.
Tomlin, C.D. 1990. Geographic Information Systems and Cartographic Modeling. Prentice Hall.
Englewood Cliffs, N .J.
Townshend, J.R.G. 1981. An Introduction to the Study of Terrain. In: J.R.G. Townshend (ed.).
Terrain Analysis and Remote Sensing. George Allen and Unwin Ltd. London.
Tribe, A 1992. Problems in Automated Recognition of Valley Features from Digital Elevation
Models and a New Method Toward their Resolution. Earth Surface Processes and Landforms.
Vol. 17.
Unwin, D. 1981. Introductory Spatial Analysis. Methuen. New York. 1981.
USGS, 1998. National Mapping Program Technical instructions part 2: Specifications, Standards
for Digital Elevation Models. U.S . Department of the Interior. U.S . Geological Service.

80

Van Remortel, R., M. Hamilton, and R. Hickey, 2001 , Estimating the LS factor for RUSLE
Trough Iterative Slope Length Processing of Digital Elevation Data. Cartography. June (in
press).
Varekamp, C. , A.K. Skidmore, P.A.B . Burrough. 1996. Using public Domain Geostatistical and
GIS Software for Spatial Interpolation. Photogrammetric Engineering & Remote Sensing. Vol.
62, No. 7.
Walsh, S.J. et al. 1987. Recognition and Assessment ofError in Geographic Information
Systems. Photogrammetric Engineering and Remote Sensing. Vol. 53 . No . 10.
Wang L. and D.C. He. 1990. A new statistical approach to texture analysis. Photogrammetric
Engineering and Remote Sensing. Vol. 56(1).
Williams, R.B .G. 1984. Introduction to Statistics for Geographers and Earth Scientists.
University of Sussex, Brighton. MacMillan.
Wood, J. 1996. The Geomorphological Characterisation ofDigitial Elevation Models. Ph.D.
Thesis. Department of Geography. University of Leicester.
Wood, J. 1998. Modelling the Continuity of Surface Form Using Digital Elevation Models. 8th
International Symposium on Spatial Data Handling. Vancouver, Canada. July 11-15th, 1998.
Zevenbergen, L.W. and C. Thorne. 1987. Quantitative Analysis ofLand Surface Topography.
Earth Surface Processes and Landforms. Vol. 12.

81

Appendix 1 Figures

82

Figure 1: Map of British Columbia and Study Area Location

Study Area

5km
The black line within the study area
represents an inset area of interest
discussed in detail in the thesis.

Study

Area
0

':>

~.' .,
• Prince
\..-,
George
~

'\_

'""'"""
500 km

83

Figure 2: Digital Elevation Profile

Resolution (distance between elevation points)
Slope between slope reversals \

:_,.

,., •

Range (relief)

I
Dots and circles represent elevation values.

•\

I

•
I
•'•-' •

Adapted from Pike, 1988, figure 2.

Figure 2 illustrates an elevation profile with graphical representation
of relief and the geometry of first and second derivatives of elevation.
The circles represent grid cell spacing on an elevation surface.

J

•

84

Figure 3: Slope Normal

z

X

The thick black line in Figure 3 represents the slope normal in relation to the
orthagonal axes (x, y represent planar coordinates and axis z represents elevation).
See section 2.2.3 for further information.

85

Figure 4a: Curvature Classified using Eigenvector Analysis

Legend

•
••
•
••
D

D

peak
ridge
saddle
flat
ravine
pit
convex hillside
saddle hillside
slope hillside
concave hillside
inflection hillside
unknown hillside

2 km

Figure 4a, top, illustrates the inset area of interest
(see Figure 1) classified using Idrisi Toposhape,
employing Pellengrini's Eigenvector analysis.
The DEM is unfiltered.
A comparison can be made to the
airphoto (bottom).

2 km

86

Figure 4b: Curvature Classified using Eigenvector Analysis

Legend

•
••
•
••
D

D

peak
ridge
saddle
flat
ravine
pit
convex hillside
saddle hillside
slope hillside
concave hillside
inflection hillside
unknown hillside

2 km

Figure 4b illustrates the inset area of interest (see
Figure 1) classified using Idrisi Toposhape,
employing Pellengrini's Eigenvector analysis (top).
The DEM has been altered following Pellengrini's
recommended frequency filtering. See Table 3 and
accompanying text for discussion of the effect of
DEM frequency filtering on the Eigenvector
analysis.
The airphoto (bottom) is for comparison.

2 km

87

Figure 5: Ridges, Passes, Peaks and Valleys

Figure 5 illustrates Wood's
discrimination of passes peaks,
ridges and valleys employing
quadratic surfaces

2 km

••
•

Ridge
Valley
Peak
Pass

The top example of the study area
illustrates Wood's results in raster
format draped on a shaded relief.
The middle example illustrates the
same results in vector format
draped on a shaded relief.
The bottom airphoto is for
comparison.

2 km

88

Figure 6: Experimental Variogram

l

~

~

~

....

•

•ge a

•

•

~

~

nugget
Lag(h)

From Burrough and McDonnell, 1998, figure 6.2, p.135

Figure 6 is an example of an experimental variogram
representing a simple transitional variogram with range,
nugget and silL

89

Figure 7: Semivariogram

80000
70000
60000
(J)

u
c
0
0

50000

·;:::

.E>

40000

(J)

(/)

30000
20000

Gaussian Curve

10000

0

5000

10000

15000

20000

Distance h (metres) between sample points.

Figure 7 illustrates the parabolic concave upward shape of the estimated
semivariogram at the origin. This is indicative oflocally changing linear drift
in the data points. As a result, universal kriging with a first order polynomial
to approximate the drift is employed to interpolate the DEM. The curve is
interpreted to half the domain size, in this case, 8500 metres.

90

Figure 8a: Randomly Generated Check Points

+

*

++ +

+

+
~

\ t+
+

+

+

+

+

+

+

++

++

+ +
+

+ +

+

+

+

+
+

.t

+ +
+ ++ +

+

+

+

+

+ ++ + +
+ ++

+

+
#

+

+

++

++ +

Skm

Figure 8a illustrates the randomly generated points selected
for removal from the input control point data. These points
were subsequently used to assess the accuracy of the DEMs.

91

Figure 8b : Residual Surface

Kriging

Figure 8b illustrates a three-dimensinal perspective view of the residual error
interpolated as a surface. At this scale, only the largest residuals are seen and represents
a visualization of the spatial arrangement of the error show by the values in Table 5.
In this case, it can be seen that significant TIN error is located at the edge of the surface.

92

Figure 9: Shaded Relief

5km

Figure 9 illustrates systematic error evident
in a shaded relief model ofthe DEM surface.
Examples of stripping (evident as thin north
south ridges) are seen within the oval.

93

Figure 10: Airphoto ofValley

BC87097 No. 056

2km

The airphoto illustrated in this figure represents the
inset area shown in Figure 1. The inset area
is discussed in the thesis with respect to
Figures 11 , 12 and 13.

94

Figure 11: Elevation Profile and Observer Perspective

Legend

Sampling Transect

0
Location of Observer

Direction of View

0

2km

Figure 11 illustrates the orientation of the perspective view
and the sampling transects employed in Figures 12, 13
and 14 respectively.

95

Figure 12a

A perspective and plan view of shaded relief
representing a raster from TIN (Linear) for the
purpose of visual quality assessment..

Figure 12b
A perspective and plan view of shaded relief
representing a raster from TIN (Quintic) for the
purpose of visual quality assessment.

96

Figure 12c
A perspective and plan view of shaded relief
representing a raster (Inverse Weighted Distance)
for the purpose of visual quality assessment.

Figure 12d
A perspective and plan view of shaded relief
representing a raster (Kriging - Liner Drift)
for the purpose of visual quality assessment.

2km

97

Figure 12e

A perspective and plan view of shaded relief
representing a raster (kriging- quinitic drift)
for the purpose of visual quality assessment.

Figure 12f
A perspective and plan view of shaded relief
representing a raster (Spline) for the purpose
of visual quality assessment.

98

Figure 13: Elevation Profiles

1650 .
Graph C

1550

Graph B

,...,

.."".

1450

Graph A

....I" 1350

•0
·=..
1250
ii:i"
If

1150
1050

0

500

1000

1500

2000

2500

3000

Distance (metres)

-

Weighted Distance

I~

Krige -Linear
Krige-Quintic

-

Spline
Tin - Linear

-

Tin - Quintic

Figure 13 illustrates elevation values sampled from Transect 1 (Figure 11 ).
Details of the graph discussed in the text refer to the boxes labelled
Graph A, Graph B and Graph C.

3500

99

Figure 14a: Slope Profiles
Slope Profile

uo ~

,

Graph C

90 +--------------------------------------f===========T-----i

80

=
t

G;
~

i

r ~ ' ~~~~

o ~~~
~~

6o

I

~

~

,

t

Ill

I

J~

~~~ ~ ~~ ~A ~

~ 3o

~

~

J

'Q ,

~

~

fY

~~

~

v~

~

~~~~

~~
~
~

~

~

~

~

~

~

~

~

uoo

500

0

-

~

~~ r~

,;

~

0

~L ~

~ ~~~ . ~~
~

\W

,A~~~~~~

~
,,~ r~
. A~~
~~ ~~~~~~~,~~

~,~~~~~~~~
~ ~~~~~ ~ ~
~,~~~~~~,~~r

o

~

~~

~ ~ Ill j
~~A~A, ~~~~~~~~,~~~~ ~ ~~~
~~r

~

~~ ~ ~~

Graph B

~

~

r

,

,

1500

2000

2500

3000

~

3500

Distance (metres)
-

Weighted Distance

-

90 +-------------------------------------------- -

Krige -Linear
Krige-Quintic
Spline
Tin - Linear

-

Tin - Quintic

~

,

Slope- Percent Graph A

uo ~
80

x

~

._.

l

Q

~~~~

o ~~

~

~
~

~~~~~~~

~~
~

G~~~j
~

~~

~

~~

~

~~~

~

u

~

0
200

~

,

,

,

,

300

400

500

600

700

~

Distance (metres)

Figure 14a, top, illustrates slope values sampled from Transect 1 (Figure 11).
The bottom of Figure 14a illustrates the detail of the transect shown in Graph A.

100

Figure 14b: Slope Profiles

Slope- Percent Graph B
100
90

--

80
70

X

60

a,o

50

Q.

0

tii

40
30
20
10
0
1400

1500

1600

1700

1800

1900

2000

Distance (metres)

-

L--------------------------------1 -

2100

Weighted Distance
Krige -Linear
Krige-Quintic

~

Slope- Percent Graph C

90 ~
85

=
=
-l
a,o
~

0

~

~~

.~~
~

70

~~~~~~~

65 ~~~

~~~~~~

~~

~~~~

~~~~

r

~

~~

~

~~~

~

~~

~r ~~~~

~

~r

~

~~~

~~~

~

~~~~

50

~

45

~~

2250

~ c~

~~~~

~~~~~

60 ~~~~~~~L ~~~~

40

~~

~~~

75

55

Tin- Linear

~~

80

..t

-

~

~

~

~

~

,

,

,

,

,

,

,

2350

2450

2550

2650

2750

2850

2950

~

Distance (metres)

Figure 14b, top, illustrates a subsection of the slope values sampled from Transect 1
(Figure 11), shown in Graph B, Figure 14a. The bottom ofFigure 14b illustrates a
subsection of the slope values sampled from Transect 1, shown in Graph C, Figure 14a.

101

Figure 15: Flow Accumulation

Skm

Figure 15 illustrates, in white, raster cells of zero flow accumulation.
For reference, the inset area of interest in shown as a white box.

102

Figure 16: Sampling Transect for Elevation and Profile Curvature

I 8 0 0

~

I 7 0 0

s

I 6 0 0

<1.)
.....

~

<1.)

"-'

§

""§
5)
~

I 50 0

I 4 0 0
I J 0 0

I 2 0 0
I I 0 0

0

I 0 0 0

:

2 o:o 0
:

Distance (metres)
.

~

0 . 6
0 . 4

0 . 2

et _o .

0
2

~

- 0 . 4

:s

- 0 . 6

2:
u

~
t;:::

- 0 . 8

p..

-I

8

- I . 2

0

I 0 0 0

2 0 0 0

3 0 0 0

Distance (metres)

Transect 2, Figure 11 was used to sample elevation (top) and profile curvature (bottom)
in order to quantify a threshold value to discriminate ridge and valley lines.
The arrows represent the threshold values. Each of the broken, vertical lines contect
points of equal distance along Transect 2.

103

Figure 17a: Profile Curvature

2km

Figure 17a represents profile curvature (the rate of change of gradient) to which a threshold
has been applied in order to stratify the profile curvature data into a binary map with convex
and concave values. Dark gray values correspond to raster cells pf convex value. Light gray
values correspond to raster cells of concave values. Black cells are flat.
A visual comparison is made between Figure 17a, band c in order to determine which form
of curvature best represents ridges and valley lines. Contours are overlain in order
to facilitate the process.

104

Figure 17b: Plan Curvature

2km

Figure 17b represents plan curvature (the rate of change of aspect) to which a threshold
has been applied in order to stratify the plan curvature data into a binary map with convex
and concave values. Dark gray values correspond to raster cells pf convex value. Light gray
values correspond to raster cells of concave values. Black cells are flat.
A visual comparison is made between Figure 17a, band c in order to determine which form
of curvature best illustrates ridges and valley lines. Contours are overlain in order to
facilitate the process.

105

Figure 17c: Average Curvature

2km

Figure 17c represents average curvature to which a threshold has been applied in order
to stratify the average curvature data into a binary map with convex and concave values.
Dark gray values correspond to raster cells pf convex value. Light gray values correspond
to raster cells of concave values. Black cells are flat.
A visual comparison between Figure 17 a, b and c is made in order to determine which
form of curvature best represents ridges and valley lines. The contours are overlain in order
to facilitate the comparison.

106

Figure 18: Profile Curvature

5km

Figure 18 illustrates profile curvature. White represents convex cells. Black represents
concave cells. Gray represents intermediate values. For reference, the inset area of
interest is shown as a black box.

107

Figure 19: Ridge and Gully

5km

White - Ridge
Dark- Gully

Figure 19 represents a map of ridge and gullies (valley lines) derived from the DEM
employing the average difference filter (filter 1, Appendix II).

108

Figure 20: Slope Classes

Original Slope Classification

5km

Filtered Slope Classes:
modal filter, 5x5 neighbourhood

Slope Classes:

Dark to light:
Flat
Gentle
Medium
Medium-steep
Steep

109

Figure 21: Slope Normal - Vector Strength

5km

2km
Figure 21 represents the vector strength of the slope normal for the study area (top)
and the inset area of interest (bottom). White represents areas of low magnitude
(flat or relatively smooth- systematic elevation change). Black represents areas of
high magnitude (elevation change is chaotic or abrupt).

II 0

Figure 22: Relief

5km

Gray - mountains
White - medium relief
Black - low relief

Ill

Figure 23: Cumulative Slope Length

5km
Figure 23 represents cumulative slope length. Increasing values
of accumulation of downslope moisture flow are illustrated as
a progression from white to black. Hydrologic basins are delineated
by the black lines.

112

Figure 24a: Bioterrain Layers

Aspect
dark gray = warm aspect
light gray = cool aspect

See Table 11 for criteria defining
warm and cool aspect and the methods
by which they are derived.

• j:

.
•

1

'

•

•

,1

-,_

~ ~ ·.

.·.;

··.r

·' ···

.I

· ,

black = depression

'

.... :·

,,•

:

... ··~

•'

. ·. ·...'"· -....
*: ~ .

\

_. , .

..

.. . ... . -= • .
,·

. ·..·

See Table 11 for criteria defining
depression and methods by which
depressions were derived.

.
.

... ..

.

'

:

.~ .

•

Depressions

,·

.. ,

I

...,•
. .,,

.' .

.,

•;

..·... ,• ..

_

·: .• ... ·. ..

..

·.... .

o

J ll

•

I

. ~ ~ ... '··. .-... .
.

• I

.,

.

~

'

,•

.

,·

'

'

\

113

Figure 24b: Bioterrain Layers

Slope
Dark to light shades:
Flat
Gentle
Medium
Medium-Steep
Steep

See Table 11 for criteria defining
slope classes and the methods
by which they are derived.

Ridges and Gullies
White - Ridge
Dark Gray - Gullies

See Table 11 for criteria defining
ridges and methods by which
they were derived.

114

Figure 24c: Bioterrain Layers

Relief
Black - Low Relief
White - Medium Relief
Gray - Mountainous Relief

These relief classes are implemented in the
broad landscape definitions recommended
by the TEM standards for the landscape
pre-typing phase.

Slope Position
Future Research

See Table 11 for criteria defining slope
position and the suggested methods
by which they may be derived.

ll5

Figure 25 : Query ofBioterrain Layers
Query: Pine, Ridge, Upper Slope, Warm , Steep Slope

Vegetation Layer

Depressions

Ridge and Gully

Slope Position

Aspect Classes

Slope Classes

Raster Layers

Attribute Type

Figure 25 illustrates how the bioterrain raster layers are queried under a vegetation
polygon to determine site attributes. For the example of the pine polygon, the dominant
bioterrain types are: ridge upper slope, warm aspect and moderate step slope. These
attributes are stored in a database table in a record item coded to the pine polygon.

11 6

Figure 26 : Soil Moisture Regime
Generalized SMR

Site Position Meso

v

L

o ~~ .

~ L__ N

D (1- (2))
D ((1)·2)
CMorS

y

.----- 11

y

N

coarse soil

-I

y
N

N

y

water table present
productive

J- 1

water table present

Y

cGEE2D

water table present

I productive

water table present

3

D (1-2)b

~

(L- 4J'

- - - CM
- - - CM

(2-4)•
(4)

s

(5-6)
(7)

y

CM

(3-5)

N

s

(5-6)

productive

s

(5-6)

low productivity

wr

(7)

productive

s

(5)

~

i
+
y

1\1

t ~ --1 low productMty
~~

t ~ --1 low productMty

s

DorCM
y
1\1

. wr

-S
CM

s

CM

y
1\1

water table at or above soil surface.
1---t-- - - - ---1 non-forested

• Generally moister if aspect.is·N orNE

(2-3)

WT

coarse soil

y

~

-

Water table near
surface low prod.

y

fine

DorCM

. --.

y

J

D (1 -2)

_

y

N

(( 4-)5)

y
DorCM (2-3)
3
N - - - CM (3-4)

N

y

(0)

D (1-2)

r-::s""'oi"'"lve-:::ry:::-,----,- Y

ridge crests
upper slopes

_Wf

-S

s

AQ

(5·6)
(2-3)

~

(4(-5))
[5-6)
4)

m

(5(-6))

(8)

Source: Figure 4.1, Resource Inventory Branch, 1999

' Generally drier ifaspect isS or SW

Figure 26 illustrates the usefulness of the bioterrain
units as discriminated by the methods developed in this
thesis and listed in Table 11 . This Soil Moisture Regime
model can incorporate the automated procedures (circled)
tested in thi s thesis.

117

Appendix II

118

Filter 1

Employing the following neighbourhood notation related to centre cell 0,0

-1 ,1 0,1
-1 ,0 0,0
-1 ,-1 0,-1

1,1
1,0
1,-1

The following filter determines the average difference between the centre analysis cell and each
of its eight neighbours. lfthis average difference is 0, then the neighbourhood ofthe cell is
either flat, or equally convex/concave. If the average difference is positive, then the
neighbourhood is generally concave. If the average difference is negative, then the
neighbourhood is generally convex.
This filter is written as ARC/INFO AML, where RESULT is the average difference layer, and
DEM is the elevation surface.
&s gri d = DEM
RESULT= (( %gr i d % (0 , 0) - %g r id % ( -1,-1 )) + ( %g ri d % (0 , 0) %g ri d % ( -1, 0)) + ( %gr i d % (0 , 0) - %g ri d % ( -1, 1 )) + ( %g ri d % (0 , 0) %grid % (0 , 1 )) + ( %g ri d % (0 , 0) - %g r id % ( 1, 1)) + ( %gr i d % (0 , 0) %g rid % ( 1, 0)) + ( %gr id % (0 , 0) - %g ri d % ( 1, -1 )) + ( %gri d % (0 , 0) %grid % (0 , - 1))) I 8

&re t urn