|
Some methods for deriving variables from digital
elevation models for the purpose of analysis, partitioning of terrain and
providing decision support for what-if scenarios
By
Peter Caccetta
CSIRO Mathematical and Information Science (CMIS)
email: Peter.Caccetta@csiro.au
24/02/1999
Background
Landform is a significant factor in agricultural production and land
management, influencing the relative growth of crops within a paddock as
well as larger-scale processes such as salinity, waterlogging and erosion.
The potential for using digital elevation models (dems) in the analysis of
landscape processes appears large.
Previous research by members of CMIS, W.A., have concentrated on
techniques for deriving variables from dems formed from the interpolation
of contour data. The contour data are readily available for much of the
south-west of W.A. in elevation sampling intervals of either 5m, 10m or
20m. The accuracy of the dems imposes limitations on the sort of variables
which can be reliably derived. Those found to be useful include estimates
of (flow) slope, upslope area and percentage upslope cleared area. These
variables have been used in the analysis of salinity risk and for
partitioning the landscape into landform units such as hilltops,
ridges and upper slopes, upper valleys, lower valleys and broad valleys.
For a review of these and other approaches, see Caccetta (1997).
The advent of relatively high resolution dems (elevation accurate to 2
metres sampled on a 10m easting/northing grid) allows a greater variety of
variables to be derived. For instance, variables such as height above
(the nearest feature) the nearest stream or salt-affected area may be
determined with a much greater accuracy than from a dem constructed from,
say, 10m contour data.
This report describes a number of algorithms, based on principles of
water flow, which are relevant to dems. The variables presented are by no
means exclusive.
The first algorithm considered preprocesses the height data to obtain hydrologically
sound dems. The algorithm is relatively simple to apply,
computationally efficient and does not require that streams have been
captured by some other means in a pre-specified format (for example in
descending order – as in ANUDEM). Not requiring streams as input is
considered an advantage as the dems often have greater spatial accuracy
than many of the streams captured at 1:50000 or 1:100000, which
constitutes much of the region of interest.
Given a hydrologically sound dem, a number of variables are described
which can be derived directly from the dem. These include water
accumulation, relief energy, average flow path length, average slope
gradient, flow slope, height above and distance from. Hypsometric
curves, stream order, stream density (km/km2)
plan and profile curvatures can also be derived, but are not
discussed here.
Next are considered variables which may be derived from combinations of
the above variables. These variables are referred to here as compound
variables, and include wetness index, stream power and universal
soil loss equation.
Finally, notes on how to use the programs to generate the variables is
given in section 5.
1 Preprocessing of DEMs to Remove Unwanted
Minimum
Typically dems have many local minima (depressions / pits), many of
which are the result of errors in the dem production process as opposed to
a few which are real (for example lakes and soaks). The erroneous minima
cause difficulty for algorithms which simulate water flowing across the
landscape, and methods for removing them need to be employed. The first
step in this process requires an expert to provide information to locate
all true depressions. An algorithm which eliminates all remaining false
depressions can then be applied. A (raster) dem with all erroneous
depressions filled and all pixels having a defined flow direction(s) is
sometimes referred to as being hydrologically sound. The process of
removing spurious minima is commonly referred to as pit filling.
As the name implies, pit filling is the process of changing the
relative elevation of the dem for an area which does not flow such that
the modified dem does flow. Generally the elevations of the pixels within
a pit are increased until they match the elevation of the lowest boundary
pixel which does flow. In this case, a flat area is produced.
Flat areas by definition do not flow, so the next step in the procedure
requires that flow directions be assigned. Again different algorithms
exist here, but perhaps the best known is that proposed by Jensen and
Domingue (1988). For this algorithm, pixels in a flat region are
recursively assigned to flow to a neighbouring pixel which does have flow
defined. A pixel at the edge of a flat region which has a flow direction
defined is used as the seed for this recursive process. In this way, the
flat region of the dem which results from pit filling is still flat;
however, flow models may be applied using the flow direction map (commonly
referred to as a flow net). A limitation of this technique, and its
variants (for example Peckham (1995) – RIVERTOOLS), is that pixels are
defined to flow in one direction only. Algorithms which route flow in one
direction are commonly referred to as D8.
One means of minimising the number of pits encountered on a flow path
is to force, during the dem interpolation process, the surface to conform
to stream line data where they exist. Here the dem is constructed such
that it has a relative minimum at the location of streams (that is streams
will sit in the bottom of valleys) and that stream elevations are ordered
such that water will run from the stream head to the stream base. This is
the approach taken in ANUDEM (see Hutchinson, (1989)). This approach
requires that each stream vector needs to be captured in the direction of
flow; stream data in W.A. do not typically have this property.
Soille and Gratin (1993) proposed an alternative approach to producing
a hydrologically sound dem, based on principles taken from the
field of mathematical morphology. Their approach combines the pit filling
and flow direction specification into a single step. The result of
applying the approach is a pit-free relative dem having no (strictly) flat
areas; areas which would normally be flat as a result from pit filling
have very small elevations added such that flow is defined.
An attractive feature of the algorithm proposed by Soille and Gratin is
that it requires no user-defined parameters and runs in O[n] (that is, it
is fast), where n is the number of pixels in the dem. After applying the
algorithm, the user is guaranteed a model which flows – there are no
cases in which it fails to define flow directions. Unlike algorithms such
as that proposed by Jensen and Domingue, the approach may also be used
with multiple flow algorithms, which have many advantages over single flow
models. For discussions relating to single flow algorithms versus multiple
flow algorithms, see for example Quinn et al. (1991) and Caccetta (1997)
pp 131-136.
Not considered in this document are approaches commonly referred to as process
modeling. The characteristic feature of these models is that they
require many parameters and boundary conditions to be specified, and in
many instances are based on the solution of differential equations. In the
context of broad scale mapping, this requirement limits the use of such
models to fairly small catchments (or very coarse grids), and typically
accurate information regarding the boundary conditions and parameters are
not available. Many algorithms based upon process modeling exist. Examples
of this form of modeling include MODFLOW (McDonald and Harbaugh, 1988) and
TOPOG.
Figure 1 (bottom) Water accumulation algorithm applied
to a dem with pits. (top) Water accumulation algorithm applied to a dem
with pits removed. Notice that the flow paths (blue) are continuous for
the pit-filled version as compared to being discontinuous for the
unfilled version.
2 Variables derived directly from the dem
Having removed spurious minima from the dem, a number of variables
based on overland flow are readily derived. Here multiple flow is assumed
and is formulated as given by Quinn et. al. 1991. This is the formulation
adopted for TOPMODEL (Beven and Kirkby, 1979); a rainfall-runoff model
which bases its predictions on catchment topography, developed by the
hydrology and fluid dynamics group in the department of environmental
sciences at Lancaster university.
The multiple flow algorithm distributes flow according to the relative
slope of downhill pixels. That is, pixels may contribute flow to one or
more (sub) catchments. Variables derived using multiple flow algorithms as
their basis are thus weighted by the proportion of flow that each
constituent pixel contributes to the (sub)catchment in question.
All figures are coloured red, indicating relatively low values ,
through orange, yellow green, blue and then purple, which signifies
relatively high values.
2.1 Upslope area (commonly called water
accumulation)
An upslope area image has assigned to each pixel a value
corresponding to the catchment area upslope of the pixel. An example of
this variable is given in figure 2. From the figure it is observed that
hilltops and ridges have relatively low values as compared with valleys
which have relatively high values.
Figure 2 Expected upslope area, or
commonly referred to as water accumulation.
2.2 Average slope of the expected upslope
catchment area (or average expected upslope slope!)
An example of this variable is depicted in figure 3. From the figure it
is observed that break aways and hilltops reminiscent of rocky
outcrops appear as purple as compared with other regions.
Figure 3 Average slope of the expected
upslope catchment area of a given pixel.
2.3 Average expected relative altitude of the upslope catchment area
(or relief potential)
Each pixel is assigned a value equal to the average upslope catchment
elevation minus the pixel elevation.
Under gravity, the potential energy of a mass is proportional to the
magnitude of its relative elevation.
If converted to kinetic energy, this equates to velocity increasing
with relative elevation. If the mass is an idealized confined water
column, this equates to column base pressure increasing with the height of
the column. It would be of interest to see if this variable is related to
water erosion (velocity) or saline seeps (pressure).
An example of this variable is given in figure 4.
Figure 4 Average expected relative
altitude of the upslope catchment area of a given pixel.
2.4 Expected flow path length
Each pixel in the image is assigned the upslope (overland) distance
that the flow has traveled. An example of this variable is given in figure
5.
Figure 5 This figure depicts the expected
flow path length. The variable may be used to identify relative position
within a catchment.
2.5 Slope in the direction of flow
The slope of the surface may be calculated in a number of different
ways, and depending on which way it is done, pixels will have different
values. The slope in the direction of flow is calculated as the weighted
average of the slope to all neighbouring pixels with elevation less than
the current pixel. In this way, slope for a pixel wide drainage line, say,
is the slope that a flow path observes as opposed to slope produced by
common difference kernals.
Figure 6 Average slope in the direction of
flow.
2.6 The number of neighbouring downhill pixels
(or outflow, convergence index)
Figure 7 The number of neighbouring
downhill pixels. This variable appears to be of limited value for the dem
considered.
2.7 Expected height above (the nearest feature
pixel)
Each pixel in the image is assigned its relative elevation as compared
to the nearest feature pixel. Here `nearest’ is defined in terms of
overland flow. An example of the variable is given in figure 8.
Flow paths were firstly extracted by stratifying the water accumulation
image, these flow paths defined the feature pixels. A height above
calculation was then performed giving an estimate of the height above the
nearest flow path pixel. Existing salt scalds could also have been used as
feature pixels, giving the height above the nearest salt scald. This
variable may be of use in salinity prediction.
Figure 8 Height above the nearest feature
pixel.
2.8 Expected distance from (the nearest
feature pixel)
Each pixel is assigned the value of its distance to the `nearest’
feature pixel. Here `nearest’ is defined in terms of overland flow. An
example is given in figure 9.
Flow paths were firstly extracted by stratifying the water accumulation
image, these flow paths defined the feature pixels. A distance from
calculation was then performed giving an estimate of the distance from the
nearest flow path pixel. Existing salt scalds could also have been used as
feature pixels, giving the distance from the nearest salt scald. This
variable may be of use in salinity prediction.
Figure 9 Distance from the feature pixel.
Here flow paths, extracted as surrogates of stream lines, were used as
features.
3 Compound Topographic Attributes
This section covers some attributes which are commonly quoted in the
literature and may be formed by combining various layers from the previous
section. Here we merely introduce them by their commonly used names
without any consideration of their validity for the concepts implied by
their names. We note that different authors use the same concepts but
define the indices differently. The use of Insitu measurements may be
useful in determining the appropriate definition.
3.1 Wetness index
The wetness index is the log of the ratio of catchment area and
flowslope. Flat areas having a large upslope area will have a relatively
high wetness index as compared with steep areas having a relatively small
catchment area. An example of this variable is given in figure 10.
Figure 10 An example of wetness index.
3.2 Stream power
Stream power is defined as the product of the upslope area and
flowslope. Steep areas having a large catchment area will have a
relatively high value for stream power as compared with flat areas having
relatively small catchments. An example is given in figure 11.
Figure 11 An example of the stream power
variable.
4 Some extensions to the above models
The water accumulation model simulates water flowing across the terrain
surface. Obviously, many factors including land cover, soil type and
contour banks will affect this flow. Schultz (1994) discusses the use of
remote sensing as input to hydrogeological monitoring and suggests that
remote sensing input can be used both as a basis for model parameter
estimation and model input.
A simple application of this is the modification of the water
accumulation algorithm to form estimates of upslope cleared area and
percentage upslope cleared area (Caccetta, 1997 pp 234-235). This is
achieved by initializing all pixels which have remnant vegetation on them
to zero (area), and all others to 1 (unit of area) at the commencement of
the algorithm; then the map that results is upslope cleared area.
Dividing this by upslope area gives the percentage upslope cleared area.
In the context of salinity mapping, catchments which have no clearing
(zero percent upslope cleared area) would tend to have little risk of
salinity, as compared to a catchment which is completely cleared (100%
upslope cleared area). It would be of interest to examine the sensitivity
of this variable for mid-range values.
We have implemented the water accumulation algorithm with the following
two options:
- The user may specify a 1-byte land cover relative equivalent area
map. This map provides an area weighting for each pixel. For
instance, in the preceding example, remnant vegetation is given a
value of zero and cleared land a value of one. There may also be other
things worth doing, for instance a pixel in a 600mm rainfall zone,
say, may be given a value of two compared to a value of one for a
300mm zone. Here we merely speculate, to introduce the concept.
- The user may specify a 1-byte land cover barrier layer. The
barrier layer specifies the percentage of flow that is lost to further
flow. For example, a dam may be specified as a 100% barrier, whereas a
contour bank may act as a barrier some what less that 100%. Again we
merely speculate to introduce the concept.
In regards to salinity prediction, if any of the variables listed above
have predictive power, then the fast algorithms which perform these
calculations may form the basis of a scenario modeling system, for
example, for examining the likely effect of planting previously cleared
paddocks to plantation, say.
5 How to generate the above results
The variables discussed in this report were generated using a number of
C coded routines. All routines assume 2 byte unsigned integer dems. All
files are assumed to have most significant byte first (MSBFirst)
format (this is the SUN default). The routines are:
- pitfill2
– his program performs the pit filling and creates a
model which is required by the routines which derive variables based on
overland flow. The program generates two files: (1) a relative elevation
file (REV) and (2) a REV support file (SUP). To note here is that the
REF file is not 1 to 1 with the original DEM: it is purely a surface
which has the required local flow characteristics.
- mmaccum –
this program performs the standard upslope area
computation. Three options are provided: (1) standard upslope area (2)
upslope area with each pixel initialized to a value specified by an
input image (3) the same as two plus the specification of an
infiltration map.
- mmaccum_plus –
this program calculates upslope slope,
upslope relative altitude and flow path length.
- habove4 –
this program generates a height above and distance
from the nearest feature pixel. Feature pixels are specified as
non-zero values in a 1 byte input image.
5.1 pitfill2
Input file: dem, 2 byte unsigned integer.
Output files: a 2 byte unsigned integer REV file; a
4 byte unsigned integer SUP file.
5.2 mmaccum
Option 1
Input files: the REV and SUP files produced by pitfill2.
Output file: 4 byte unsigned integer upslope area
file.
Option 2
Input files: the REV and SUP files produced by pitfill2, a 1 byte
unsigned char initialisation layer specifying the relative importance of
each pixel (each pixel is worth `x’ units of area).
Output file: 4 byte unsigned integer upslope area file.
Option 3
Input files: the REV and SUP files produced by pitfill2, a 1 byte
unsigned char initialisation layer specifying the relative importance of
each pixel (each pixel is worth `x’ units of area), a 1 byte unsigned
char infiltration layer (each pixel represents the percentage of flow
that leaves to pixel).
Output file: 4 byte unsigned integer upslope area file.
5.3 mmaccum_plus
Input files: the REV and SUP files produced by pitfill2.
Output file: 2 byte unsigned integer file.
5.4 habove4
Input files: the REV and SUP files produced by pitfill2, a 1 byte
unsigned char feature file. Features are assumed to be those pixels with
non-zero values.
Output file: (options 1 and 2) 2 byte unsigned integer file.
(option 3) 4 byte unsigned integer file.
References
Beven K.J. and Kirkby M.J. (1979). A physically based variable
contributing area model of basin hydrology. Hydrol. Sci. Bull., 24, 43-69.
Caccetta, P.A. (1997), Remote Sensing, Geographical Information Systems
(GIS) and Bayesian Knowledge-Based Methods for Monitoring Land Condition, PhD
thesis, School of Computing Science, Curtin University of Technology.
Hutchinson, M.F. (1989), A new method for gridding elevation and stream
line data with automatic removal of pits, Journal of Hydrology, 106:211-232.
Jensen S. K. and Domingue, J. O. (1988), Extracting topographic
structure from digital elevation data for geographic information system
analysis, Photogrammetric Engineering and Remote Sensing
54(11):1593-1600.
McDonald, M.C., and Harbaugh, A.W., 1988, A modular three-dimensional
finite-difference ground-water flow model: U.S. Geological Survey
Techniques of Water-Resources Investigations, book 6, chap. A1, 586 p.
O' Loughlin E.M (1981). Saturation regions in catchments and their
relation to soil and topographic properties. J. Hydrol., 53, 229-246.
Peckham, S. D. (1995), Appendix A: River Network Extraction from Large
DEMs, Excerpt from the doctoral thesis of Scott Dale Peckham,
University of Colorado, Boulder
Qiuinn, P., Beven K., Chevallier, P., and Planchon, O. (1991), The
prediction of hillslope flow paths for distributed hydrological modelling
using digital terrain models, Hydrological Processes, 5:59-79.
Quinn, P., Beven, K., Chevallier, P., and Planchon, O. (1991), The
prediction of hillslope flow paths for distributed hydrological modelling
using digital terrain models. Hydrological Processes, 5:59-79.
Schultz, G. A. (1994), Meso-scale modelling of runoff and water
balances using remote sensing and other GIS data, Hydrological Sciences
– Journal des Sciences Hydrologiques, 39(2).
Soille, P. and Gratin, C. (1993), Short Communication: An Efficient
Algorithm for Drainage Network Extraction on DEMs, Centre de
Morphologie Mathematique, Ecole Nationale Superieure des Mines de Paris,
35, rue Saint-Honore F-773005 Fontainebleau Cedex, France.
 To
top
last updated
September 19, 2005 11:50 AM |
