1. health and fitness
2. disorders
3. mental disorder
4. depression

Journal of Landscape Studies 2 (2009), 57 – 68
Received: 18 December 2009; Accepted: 6 January 2010; Published online: 9 January 2010
Journal of
Landscape
Studies
How to extract river networks and catchment boundaries from DEM: a
review of digital terrain analysis techniques
Vojtěch Barták*
Department of Applied Geoinformatics and Spatial Planning, Faculty of Environmental Sciences
University of Life Sciences Prague, Kamýcká 129, 165 21, Prague, Czech Republic
Abstract
Raster digital elevation models (DEMs) are a source of data with rapidly increasing quality and availability. Although they
are used in a wide range of research and application fields, including hydrological modelling, river basin management,
geomorphological modelling, soil science, plant ecology, and landscape characterization, only a few of the methods for
processing DEMs described in the literature are accessible in commonly-used GIS software. A review of hydrologicallyoriented digital terrain analysis (DTA) techniques is provided. Apart from computation of basic terrain attributes, the main
tasks are extraction of the automated river network and catchment boundaries. The basic approach is overland flow
simulation using the predetermined flow directions from each cell to one or several of its neighbours. Such an approach
requires treatment of the problematic parts of DEM, where the flow directions cannot be determined due to the absence of
any lower neighbour (i.e. closed depressions and flat areas). Except in the case of two recently developed computer
programs (DEMETERR and TAS), the available GIS software does not offer a satisfactory range of usable methods.
Further development of programs like ArcGIS or GRASS is therefore necessary to fill this gap.
Key words: Digital terrain analysis; Raster DEMs; Overland flow algorithms; Catchment delineation; River network
extraction.
1. Introduction
Raster digital elevation models (DEMs) are a
source of data with rapidly increasing quality and
availability. They are used in a wide range of
research and application fields, particularly
hydrological modelling, river basin management,
geomorphological modelling, soil science, plant
ecology, and landscape characterization (Moore et
al., 1991). Unfortunately, only a few of the
methods for processing DEMs described in the
literature are accessible in commonly-used GIS
software. Moreover, the accessible methods are
sometimes the oldest and worst ones. Recently, two
programs with satisfactory algorithm choice have
been developed: DEMETERR (Barták, 2008) and
* Corresponding autor; E-mail: [email protected]
Available online at: www.centrumprokrajinu.cz/jls/
TAS (Lindsay, 2005), both focused on
hydrological applications.
The aim of this paper is to provide a review of
hydrologically-oriented digital terrain analysis
(DTA) techniques. Most of them are motivated by
models of catchment rainfall-runoff response and
river basin management needs, so the main tasks
are to identify the river network and the positions
of the catchment boundaries. The oldest methods
were based on local topography evaluation (i.e.
determining ridge and valley lines; see Tribe,
1992), but since the early 1980s the main approach
has been to model overland flow (Mark, 1984;
O’Callaghan et Mark, 1984; for a review, see
57
V. Barták: Journal of Landscape Studies 2 (2009), 57 – 68
Bertolo, 2000).
In the present context, the term “overland flow”
does not refer to any real outflow process, but to
the theoretical water flow over the completely
uncovered and impermeable terrain, generated by
the initial unit rainfall regularly spaced over the
entire DEM. The overland flow simulation consists
of two basic steps. The first is the assignment of
flow directions along which the water flows from a
particular cell to one or several of its eight
neighbours. The second step is flow accumulation
computation, i.e. determining the upslope area
draining through the particular cell. The results
from the latter step are the size and the position of
the catchment (often referred to as the contributing
area or upslope drainage area) corresponding to
each cell.
Since the assignment of proper flow directions is
a necessary step for the successful solution of these
tasks, those cells that have no lower neighbour are
problematic parts of the DEM. On such cells the
flow directions remain undefined, so the resulting
drainage network is discontinuous. Such
problematic DEM parts are of two types, closed
depressions and flat areas, and both must be treated
before simulating the overland flow.
2. Primary and secondary terrain attributes
Moore et al. (1991) distinguish between primary
and secondary terrain attributes derivable from a
DEM. The primary attributes are computed directly
from elevations, while the secondary attributes are
defined as a compound of the primary attributes.
Further, the primary attributes may be computed
using either the values of the eight neighbours (i.e.
the surrounding cells) for each DEM cell, or larger
cell neighbourhoods. Examples of primary
attributes are the terrain slope, the aspect, or the
horizontal and vertical curvature; the contributing
area is an example of a secondary attribute.
The slope, aspect and curvature can be calculated
for the particular cell from the elevations of its
eight neighbours, using the first and second
discrete derivatives (for details, see e.g. Gallant et
Wilson, 1996). Another way is to compute the
slopes between the cell and its neighbours (as the
ratio between the vertical and horizontal distance),
and then use the steepest slope (i.e. greatest
downhill slope), or their weighted mean, where the
weights are proportional to the slope magnitude
(Quinn et al., 1991). These methods directly
correspond to certain methods for computing
58
contributing areas, which are described in the next
section.
A well-known example of secondary terrain
attributes is the topographic index (sometimes
called the wetness index), which indicates the
tendency of the site to soil saturation (Quinn et al.,
1995). It is defined as log(a/S), where a is the
contributing area of the cell and S is the slope. For
a comprehensive list of primary and secondary
terrain attributes, see Moore et al., 1991.
3. Assigning flow directions
3.1. SFD8
The simplest method for assigning flow directions
is to determine the flow direction from a given cell
to that of its neighbours, to which the slope is
steepest (Fig. 1). Such an approach is called SFD8
(Single Flow Direction chosen from 8
possibilities), and the first references are in Mark
(1984) and O’Callaghan et Mark (1984).
Figure 1. The SFD8 algorithm. The large numbers denote the
elevations, the small numbers denote the slope from the central
cell, and the arrow represents the determined flow direction.
The slope from the cell to its neighbour is
computed as the difference in elevation divided by
the horizontal distance (cell size for cardinal
neighbours and cell size multiplied by the square
root of 2 for diagonal neighbours), and is
considered as positive downhill.
The basic feature of this approach is that the
flow is always represented as convergent. No
splitting of the flow to several directions is
allowed. This may be regarded as an advantage as
well as a disadvantage, depending on the aim of the
V. Barták: Journal of Landscape Studies 2 (2009), 57 – 68
analysis. An additional and clearly undesirable
feature is that the flow is biased towards one of the
eight possible directions, so the result is strongly
dependent on the grid orientation of the DEM.
3.2. MFD8
The MFD8 (Multiple Flow Directions chosen from
8 possibilities) algorithm was developed as a
modification of SFD8 to represent flow divergence.
The water flow from the cell is divided among all
the neighbours to which there exists a positive (i.e.
downhill) slope (Fig. 2).
Figure 2. The MFD8 algorithm. The thickness of the arrows
symbolizes the proportion of flow falling on the particular
neighbour.
et al. (1991). The generalization to p ≠ 1 (or 1.1)
was given by Holmgren (1994) and Quinn et al.
(1995). Note that the higher the values of
parameter p, the greater the similarity of the MFD8
algorithm to SFD8 (theoretically if p →∞ then both
algorithms merge). In other words, the exponent p
controls the degree of flow convergence.
3.3. SFD∞
This approach (Single Flow Direction chosen from
infinite possibilities) was proposed by Tarboton
(1997). The algorithm chooses the direction of the
steepest slope from the (theoretically) infinite
range from 0˚ to 360˚. It proceeds by constructing
eight triangles with vertices in the centres of the
cell and its two (mutually adjacent) neighbours
(Fig. 3), and subsequently identifies the steepest
slope direction on the planes determined by the
triangles. The greatest of these slopes is then
chosen and, if it does not go directly to the centre
of one of the neighbours, the flow is apportioned to
the corresponding two neighbours, so that the
amount of flow falling on the particular neighbour
is inversely proportional to the deviation of the
direction to the neighbour from the determined
steepest slope direction (angles α1 and α2 in Fig. 3).
The proportion of the flow assigned to the i’th
neighbour is computed using the ratio
S ip
,
∑ S jp
0 < j ≤8
(1)
where Si is the slope from the cell to its i’th
neighbour, the summation is over the neighbours to
which a positive slope exists, and p is a parameter
that can be considered as constant or spatially
variable (for the latter, see Quinn et al., 1995 and
Kim et Lee, 2004). The original formulation of (1)
with p = 1.1 was proposed by Freeman (1991), and
with p = 1 and slopes multiplied by the so-called
“contour length” (coefficient greater for diagonal
neighbours than for cardinal neighbours) by Quinn
Figure 3. The SFD∞ algorithm. The dots are the centres of the
cells, the arrow is the flow direction. Note that this example
does not use the topography from Figures 1 and 2.
The fact that the flow is mostly divided into two
directions has led some authors to consider this
algorithm as “multiple flow” with the
representation of divergent flow. However, the
original aim (Tarboton, 1997) was to represent
59
V. Barták: Journal of Landscape Studies 2 (2009), 57 – 68
flow convergence, but without the undesirable bias
towards cardinal and diagonal directions. Flow
divergence obviously occurs, but it is always
constrained to the two adjacent neighbours.
3.4. MFD∞
The MFD∞ method (Multiple Flow Direction
chosen from infinite possibilities), developed by
Seibert et McGlynn (2007), is a straightforward
generalization of SFD∞, in the sense that MFD8 is
a generalization of SFD8. As in SFD∞, the
algorithm seeks the steepest slope directions on the
eight triangle-generated planes. Then, however, it
does not choose the steepest slope but all the
positive and locally steepest slopes (Fig. 4). The
flow apportioning between the chosen directions is
realized by (1), and the corresponding flow part for
each such direction is subsequently divided into
two mutually adjacent neighbours, in the same
manner as in SFD∞.
appear much more realistic than those generated by
SFD8, the random nature of the algorithm,
providing different results with the same data,
together with flow lines inconsistent with the
topography, are usually considered as undesirable.
Another method is the so-called DEMON
(Digital Elevation MOdel Networks) algorithm,
developed by Costa-Cabral et Burges (1994). The
DEM values are placed in the corners of the cells
rather than in their centres, and the flow directions
are determined on the plane fitted to the four
elevations of each cell. The flow is then
represented by flow tubes of variable width, so a
low degree of flow dispersion on a divergent
topography is allowed. This algorithm represents a
quite different approach than those described
above, but its use is restricted due to the presence
of complicated topographies requiring special
treatments (Tarboton, 1997).
The idea of a spatially variable degree of flow
divergence led Pilesjö et al. (1998) to develop the
Form Based algorithm. This consists of three basic
steps. First, each cell is indicated as undisturbed
(simple), complicated (consisting of several
valleys) or flat, according to the configuration of its
lower and higher neighbours. Then, a different
flow mechanism is carried out on the each such cell
type. On complicated cells, the flow is apportioned
between the valleys, in proportion to their width.
On undisturbed cells, the terrain shape is evaluated
as convex or concave, and then the MFD8 method
or the modified SFD∞ method is chosen.
4. Flow accumulation and watershed delineation
Figure 4. The MFD∞ algorithm. The dark lines and dots
delineate the triangles on which the locally steepest slopes are
determined. The arrows represent the determined flow
directions.
3.5. Other methods
Apart from the algorithms described above, some
other flow direction assignment methods have been
proposed in the literature. In order to avoid the bias
of SFD8 toward one of the eight possible
directions, Fairfield et Laymarie (1991) suggested
a modification that introduces a random component
into the slope calculation procedure. Although their
method (named Rho8) leads to flow lines that
60
Once the flow direction has been defined for each
cell, we can simulate the flow accumulation and
compute the contributing areas. One possible
method (O’Callaghan et Mark, 1984) is as follows.
At the beginning, the flow accumulation value is
set to 1 for each cell. Then the DEM is scanned
row-by-row and, for each cell, the number of its
inflows (i.e. the neighbours whose flow direction is
directed towards it) is counted. Subsequently, the
DEM is repeatedly scanned (row-by-row) and,
during each iteration, the cells with no inflow are
determined and their flow accumulation values are
distributed among their neighbours (by one of the
methods described above). Simultaneously, the
numbers of inflows of such neighbours are
decreased by one. This implementation has O(n2)
V. Barták: Journal of Landscape Studies 2 (2009), 57 – 68
asymptotic complexity. The complexity O(n) can
be obtained when the neighbours whose number of
inflows has dropped to zero are stored in the stack
and the computation then proceeds at the first cell
on the stack (while the stack is nonempty; see
Barták, 2008). Thus an iterative row-by-row scan is
avoided. Another approach, with asymptotic
complexity O(nlogn), uses a priority queue and
searches the DEM “from highest to lowest”.
For watershed delineation, the position of the
outlet must be known. Such an outlet may be
represented by a single cell or by a group of cells.
The catchment consists of all the cells draining
through the flow directions to the outlet. Such cells
can be found in exactly the same way as the flow
accumulation described in the previous section,
with two differences. Firstly, only the directions
are traced and no flow partitioning has to be carried
out. Secondly, the output is a list of cells belonging
to the catchment rather than their count.
depression has a “bottom” (consisting of cells with
no lower neighbour), and a “catchment” (consisting
of cells from which the flow paths terminate at the
bottom). The lowest cell on the boundary of a
catchment (i.e. the point through which the water
will overflow when the depression is filling) is
often called the “pour point”, and the depression
itself can be defined as consisting of all the cells
from its catchment whose elevation is lower than
the elevation of the pour point (Fig. 5).
5. Removing closed depressions
Closed depressions (often called sinks or pits) are
cells or groups of cells completely surrounded by
cells of higher elevation. This means that no flow
direction can be assigned to such cells, since no
positive slope is determined. Most landscapes are
modelled by fluvial and erosion processes, and the
presence of closed depressions in the DEM is
usually due to data errors, such as errors from
interpolation or measurement (Florinsky, 2002).
Sometimes, however, they can reflect real terrain
features (e.g. volcanic craters). Modellers often
treat them as spurious features that must be
removed from the DEM, otherwise they cause
discontinuities in the drainage network determined
by the flow directions, and thus prevent the
successful solution of other DTA tasks (e.g.
catchment delineation or river network extraction).
An interesting attempt to develop a method for
distinguishing between actual and spurious
depressions is presented in Lindsay et Creed, 2006.
Before describing methods for removing
depressions, it is necessary to discuss their
identification in the DEM, since they must usually
be identified before they are removed, and since
various complex depressions, such as nested or
cascade-organized
depressions,
make
this
procedure
surprisingly
complicated.
Some
terminology is needed for the discussion. Each
Figure 5. The vertical section of the DEM with a depression.
The oldest methods for finding depressions are
based on the following two-phase procedure
(Jenson et Domingue, 1988; Martz et De Jong,
1988). In the first phase, the bottoms are found by
row-by-row scanning; in the second phase, the
catchment (and the pour point on its boundary) is
found for each bottom. This type of procedure has
asymptotic complexity O(n2), where n is the
number of cells in the DEM, and since they treat
each depression separately, their real complexity
depends on the number and the complexity of the
depressions.
Soille et Gratin (1994) and later (without any
reference to their paper!) also Wang et Liu (2006),
suggested a more effective procedure. It uses a
“from lower to higher” and “from the edges
inwards” search, based on the priority queue data
structure (also known as a binary heap), rather than
a row-by-row scan. Only one search is needed and
all the depressions and pour points are easily
identified, with asymptotic complexity O(nlogn)
(complexity of the priority queue).
61
V. Barták: Journal of Landscape Studies 2 (2009), 57 – 68
5.1. Filling
The simplest and probably best known method for
removing depressions from the DEM, first
proposed by Jenson et Domingue (19898) and
Martz et DeJong (1988), is to “fill” them by
increasing the elevations of the depression cells up
to the elevation of the pour point (Fig. 6). The
theoretical background is the assumption that the
depressions are caused by an underestimation of
the elevations in the depressions. The depressions
can be effortlessly filled during the priority-queuebased depression search described above (Soille et
Gratin, 1994; Wang et Liu, 2006).
Figure 7. An example of depression carving. The medium-thick
line is the original topography before carving, the thickest line
is the new topography after carving.
Although this method is sometimes adequate,
and removes the depressions by decreasing the
elevations of only few cells, in the case of small
deep depressions it can sometimes cause unrealistic
long deep “channels” consisting of the decreased
cells. In such cases, the filling approach would
probably be better.
5.3. Combinations
Figure 6. An example of depression filling. The thick line is the
original topography before filling, the dotted line is the new
topography after filling.
An obvious disadvantage of this method is that
the filled depressions create large flat areas, where
the problem of assigning the direction of flow
arises again. Such regions must be further treated
by some method for removing flats (see the next
section).
5.2. Carving
An alternative assumption about the origin of
depressions can be that they are due to an
overestimation of the elevations in the
neighbourhood of the pour point, where the
overestimated cells have formed a “dam”. The
proper solution is therefore to decrease rather than
increase the elevations of the affected cells (Fig. 7).
This approach (sometimes called “breaching”) was
proposed independently (and with some
differences) by Rieger (1998), Jones (2002) and
Soille et al. (2003). The latest implementation is
also based on priority-queue-search, and thus
shares its complexity.
62
A logical extension of the approaches described
here is to combine them, with the particular method
being chosen according to its impact (or cost). The
cost may be measured as the number of changed
cells (the area cost), as the total change in
elevations (the volume cost), or as a combination
of these two. The last of these options is adopted in
the IRA (Impact Reduction Approach) method
proposed by Lindsay et Creed (2005), in which the
cost is computed for each of the two methods
(filling as well as carving) applied to the whole
DEM, and then the method with the lower cost is
chosen.
A straightforward improvement to such a method
would be to decide between filling and carving for
each depression separately. Unfortunately, this
cannot be done easily due to the presence of
complex depressions (Barták, 2008). However,
there are two methods that solve each depression
separately, and, moreover, enable us to treat some
part of the depression by carving and the rest by
filling (Fig. 8). The first method was proposed by
Martz et Garbrecht (1999), who called it “outlet
breaching” (sometimes referred to as “constrained
breaching”). The procedure tries to decrease the
elevations of some cells in the neighbourhood of
the pour point, before filling is performed. Soille
V. Barták: Journal of Landscape Studies 2 (2009), 57 – 68
(2004) suggested a more sophisticated method
based on mathematical morphology in which the
proportion of carving and filling is set in such a
manner that the resulting cost is the lowest
possible.
Although they generate more realistic terrains,
both methods have the awkward feature that they
are usually iterative, since in raising the elevations
some new (mostly small and shallow) depressions
may occur. This leads to a relatively high cost,
containing elevation changes on the cells which
were originally not in any depression. On the other
hand, the advantage of such methods is that the flat
areas and the depressions can be treated
simultaneously by the same procedure.
6. Treatment of flat areas
Figure 8. Combination of filling (arrows up) and carving
(arrows down).
5.4. Physically-based methods
These methods try to create a terrain which is more
realistic than the flat areas generated by filling or
carving procedures. The first and very simple
approach is to fill the depressions in a manner
similar to the filling technique, but the elevations
of the cells in the depression are increased not at
the pour point elevation but slightly higher. The
method is called SDI (Spatially Distributed
Elevation Increment; Tianqi et al., 2003), which
reflects the fact that the increment added to the
pour point elevation is spatially distributed in order
to generate a gently inclined plane. An important
assumption is that the DEM is relatively simple
and is roughly inclined to one of its corners, so the
resulting slope orientation of the plane originating
from the particular depression corresponds to this
DEM inclination.
The second, truly physically-based method is
PEM4PIT (Physical Erosion Model for PIT filling;
Grimaldi et al., 2007). It simulates filling-up of the
depressions by erosion, the intensity of which is
consistent with the surrounding topography. The
elevations of the cells in the depressions are
changed so that the new topography satisfies the
mass continuity equation for steady state
topography, reflecting the equilibrium between
tectonic, fluvial and diffusion processes in the
landscape.
Flat areas (or “flats”) are cells or groups of cells
completely surrounded by cells with the same (or
higher) elevation (Fig. 9). They are often present in
DEMs due to data rounding on a relatively flat
topography, or as a result of a particular
depression-removing procedure. Apart from
methods that treat flats together with depressions
(both of the physically-based methods described
above, and the carving implementation proposed
by Jones, 2002), the literature offers two other
usable approaches, both constructing the improved
topography of the flat according to the surrounding
topography.
Figure 9. An example of a flat. The numbers denote the
elevations. The flat is bordered by the thick line. The light grey
cells have the flat elevation, the dark grey cells are higher, and
the white cell is lower than the flat elevation.
6.1. Towards the lower gradient
The first method, originally proposed by Jenson et
Domingue (1988) and improved by Soille et Gratin
(1994), is to construct a topography inclined
towards the lower neighbours of the flat (Fig. 10).
63
V. Barták: Journal of Landscape Studies 2 (2009), 57 – 68
In the first step, the cells not adjacent to any lower
cell are increased by a sufficiently small increment.
This step is then repeated, until there is no cell with
a lower neighbour. The increment must be
determined in such a way that the elevations in the
flat will not exceed the elevation range of the
surrounding topography.
mathematical morphology, which implicitly solves
such cases and therefore does not require any
additional step.
Figure 11. From the upper gradient computed on the flat from
Figure 9. The meanings of the symbols are as in Figure 10.
Figure 10. Towards the lower gradient computed on the flat
from Figure 9. In the picture on the left, the numbers denote the
number of elevation increments added to the particular cell. The
picture on the right shows the resulting flow directions
(determined by the SFD8 algorithm).
6.2. Combined gradient
The obvious result of the previous method is a
simple inclined plane, in which the subsequently
determined flow paths form non-realistic parallel
lines (Fig. 10). To represent some degree of flow
convergence, Garbrecht et Martz (1997) suggested
a method combining the gradient “towards lower”
topography with the gradient “from higher”
topography. The first is computed in exactly the
same way as in the preceding method. The
computation of the second is analogous: in the first
iteration the cells which are adjacent to some
higher cell and are not adjacent to any lower cell
are increased. Then the algorithm iteratively
increases the cells not adjacent to any lower cell
(for the result, Fig. 11). The final step is the
summation of both gradients (i.e. of the number of
increments added to each cell through any of the
preceding steps, Fig. 12).
An additional step is usually required, due to
certain special cases where new (but always small
and simple) depressions occur (Fig. 12, left). Such
a step consists of computing the “towards lower”
gradient on this new flat, with a half elevation
increment. Soille et al. (2003) proposed an
improved version of this algorithm, based on
64
7. River network extraction
Sometimes it may be desirable to estimate the
position of the river network from the DEM, even
when the “real” position (i.e. the topographic map)
is known, because the real position often does not
agree with the DEM topography. Spatially
distributed hydrological models that compute the
outflow from a catchment using flow directions
determined on the DEM cells should distinguish
between “river cells” and “hill cells”, and represent
the outflow from these two cell types in a different
way, e.g. with flow divergence (multiple flow
algorithms) on the hills and flow convergence
(single flow algorithms) in the streams. Even in
river basin management, when subcatchment
Figure 12. The final step of the combined gradient algorithm.
The numbers in the picture on the left are the sums of the
numbers from Figures 10 and 11. The picture on the right
shows the resulting flow directions. The flow direction of the
central cell was determined through an additional step, since
after summation of the total increments a new flat arose on this
cell (see Garbrecht et Martz, 1997). (i.e. the number of
increments added to the central cell will be 3.5).
V. Barták: Journal of Landscape Studies 2 (2009), 57 – 68
delineation is required, there is a need for
information about the position of the outlets (and
thus the position of the river network). In all these
tasks, the position of the river cells must be
compatible with the remaining DEM topography so
that the water will not be forced to flow uphill.
Regardless of its real nature, the river network
estimated from the DEM typically consists of lines
that are one cell in width. Thus the methods for
river network extraction must decide for each cell
whether it is a river cell or not. Such a decision
may be made using the real position of the river
network (the so-called “blue lines” in the maps), or
may be fully automatic, using the DEM as an
exclusive source of information. The first
approach, often referred to as “enforcement” or
“burning” (see e.g. Soille et al., 2003; Turcotte et
al., 2001), proceeds by overlaying the blue lines
over the DEM and changing the elevations of the
cells intersecting with blue lines. In this paper, the
latter, fully automatic, approach is discussed.
7.1. Local curvature methods
The oldest methods for fully automated river
network extraction are based on satisfying some
terrain shape criteria that are typical for river
valleys. Each cell is analyzed within the moving
window (usually 3x3 cells in size), and the shape
of the corresponding terrain is examined. One
possible way is to make several vertical cuts in
different directions and then to analyze these cuts
as to whether they are V-shaped, which is assumed
to determine a river valley. Another way is to
measure the degree of terrain convergence on each
cell by counting the number of its neighbours that
are higher than the inspected cell. The best known
example is the Peucker et Douglas (1975)
algorithm, which scans the DEM by 2x2 sized
windows, and within each group of four cells the
highest cell is marked. Then the river cells are
determined as the cells that remain unmarked.
An extensive review of local curvature methods
along with a detailed discussion is provided by
Tribe (1992), who noted that these apparently
simple methods often require an additional step,
which makes them much more complicated. The
typical
resulting
networks
are
namely
discontinuous, full of inappropriately determined
river cells, and are too wide, so that there us a need
for a step that connects them (e.g. using previously
determined flow directions) and thins them. An
interesting algorithm solving all these obstacles
was proposed by Meisels et al. (1995).
7.2. Flow accumulation methods
The basic concept of this approach is that river
cells can be determined as cells in which the flow
accumulation (i.e. the contributing area) is high
enough for the transition of hillslope processes to
fluvial processes (Mark, 1984). Therefore the cell
may be considered as a river cell if its contributing
area exceeds some predefined threshold value. If
the SFD8 algorithm is used, the obvious advantage
of such an approach will be that the resulting river
network is fully connected. The particular methods
differ in the way the threshold value is determined.
The simplest method is to consider the area
threshold as spatially constant (Mark, 1984; Jenson
et Domingue, 1988; Martz et Garbrecht, 1982;
Quinn et al., 2005; and many others). The
corresponding theoretical background was given by
Tarboton et al. (1992). Such approaches have been
criticized for insufficient density of the resulting
networks in the higher parts of a landscape, and for
inappropriate positions of river initiations (see e.g.
Montgomery et Foufoula-Georgiou, 1993). One
possible solution is to consider the area threshold
as dependent on the slope values. Generally, the
higher the slope, the lower the sufficient area
threshold.
In
practice,
several
authors
(Montgomery et Foufoula-Georgiou, 1993;
Montgomery et Dietrich, 1994; Ijjasz-Vasquez et
Bras, 1995; and others) argue that the appropriate
variable whose critical value has to be determined
is aS2, where a is contributing area and S is slope.
This threshold often leads to discontinuous
networks, so subsequent network connecting
according to flow directions is usually required.
A substantial disadvantage of this approach is
that the proper estimation of a local slope is a very
difficult task, due to its strong scale dependency.
Therefore Roth et al. (1996) suggested replacing
the slope by the so-called “relative subcatchment
elevation”, which is defined for each cell as the
mean elevation over its contributing area (i.e. its
catchment) minus the elevation of the cell. They
argue that this variable exhibits a behaviour similar
to that of the local slope multiplied by the
contributing area; however, the estimation is less
affected by uncertainty due to its summation over a
65
V. Barták: Journal of Landscape Studies 2 (2009), 57 – 68
larger area.
An alternative method for extracting networks
with a higher density in the higher parts of a
landscape was developed by Tarboton et Ames
(2001). Their area threshold value is determined by
using the “weighted contributing areas” rather than
the normal contributing areas. The weights
reflecting local curvature enter into the
computation of the contributing areas, so that the
initial amount of water is greater in the cells
located at the convergent topography. The simplest
option is to set this amount to two for the cells
determined as river cells by the Peucker&Douglas
algorithm, and to one for the rest, and then to
compute the contributing areas in exactly the same
way as described in section 2.2.
7.3. Determining threshold values
Irrespective of how the type of threshold was
chosen, the proper determination of its value is an
additional and non-trivial question. The simplest
method is to visually compare the resulting
network with the “blue lines” (e.g. Mark, 1984;
Jenson et Domingue, 1988; Martz et Garbrecht,
1992). When the blue lines are converted to an
appropriate raster, we can improve this comparison
by minimizing the spatial differences between the
extracted and mapped river cells (for such an
optimization using a genetic algorithm, see Kim et
Lee, 2004). Montgomery et Foufoula-Georgiou
(1993) suggested using information about the real
position of the river sources for calculating the
critical value of aS2.
If no such additional information is available, the
determination should use some criterion based only
on elevations. One way is to analyze the slopecontributing area plot and determine the break
points corresponding to river initiations (Tarboton
et al., 1991; Ijjasz-Vasquez et Bras, 1995);
however, there is no agreement concerning the
proper break point. A completely different
approach was developed by Tarboton et al. (1991).
It is based on fulfilling some empirical laws (often
called Horton’s laws) describing the self-similar
and scale-invariant morphology of real river
networks. One such law is called the “constant
drop law”, and can be formulated as follows: the
average elevation drop along river reaches with the
same Strahler order is approximately constant
across the Strahler orders. According to Tarboton
66
et al. (1991), the lowest threshold value leading to
a network that satisfies the constant drop law can
be considered as the proper value. The authors
called this procedure “constant drop analysis”, and
an evaluation made by Barták (2008) proved its
usability.
8. Final remarks
A brief outline of available DTA techniques has
been presented. The algorithms for particular tasks
often differ substantially, and the individual choice
depends on the aim of the analysis to be carried
out. Mostly, it is not a matter of saying which
method is generally the best or the proper one. For
example, the SFD8 algorithm, which is unable to
represent divergent flow, is more appropriate for
watershed delineation than multiple algorithms,
since flow dispersion leads to overlapping of
adjacent catchments.
Therefore, the availability of a suitable range of
methods is desirable. Unfortunately, this is not the
typical case, since just one method or a limited
number of methods for a particular task can be
accessed in widely-used commercial or noncommercial software (e.g. ArcGIS, Olivera et al.,
2002; GRASS, GRASS Development Team, 2008;
TAPES, Gallant et Wilson, 1996; etc.). Two
exceptions are DEMETERR (Digital Elevation
Models and their Elementary TERRain analysis;
Barták, 2008) and TAS (Terrain Analysis System;
Lindsay, 2005), which are however in an early
stage of development. Further development and
upgrading of extensive GIS programs, such as
ArcGIS and GRASS, is therefore more than urgent.
Acknowledgements
This work has been supported by an IGA CULS
grant. I thank Barbora Matoušková for
substantially improving the manuscript. I thank
Petra Šímová and Vladimír Puš for initializing this
work.
References
Barták, V. 2008. Algoritmy pro zpracování digitálních
modelů
terénu
s aplikacemi
v hydrologickém
modelování. Diploma thesis [online]. Czech
University of Life Sciences, Prague. URL:
http://fzp.czu.cz/~bartakv/thesis [cit. 2009-11-23].
V. Barták: Journal of Landscape Studies 2 (2009), 57 – 68
Bertolo, F. 2000. Catchment delineation and
characterisation – a review. Technical report EUR
19563 EN [online]. Space Applications Institute –
Joint Research Institute – Ispra, Italy. URL:
http://agrienv.jrc.it/publications/pdfs/CatchRev.pdf
[cit. 2008-06-11].
Costa-Cabral, M.C., Burges, S.J. 1994. Digital elevation
model networks (DEMON): A model of flow over
hillslopes for computation of contributing and
dispersal areas. Water Resources Research 30(6):
1681-1692.
Fairfield, J., Leymarie, P. 1991. Drainage network from
grid digital elevation models. Water Resources
Research 27(5): 709-717.
Florinsky, I.V. 2002. Errors of signal processing in
digital terrain modelling. International Journal of
Geographical Information Science 16(5): 475-501.
Freeman, T.G. 1991. Calculating catchment area with
divergent
flow
based
on
regular
grid.
Computers&Geosciences 17(3): 413-422.
Gallant, J.C., Wilson, J.P. 1996. TAPES-G: a grid-based
terrain analysis program for the environmental
sciences. Computers&Geosciences 22(7): 713-722.
Garbrecht, J., Martz, L.W. 1997. The assignment of
drainage direction over flat surfaces in raster digital
elevation models. Journal of Hydrology 193(1-4):
204-213.
Garbrecht, J., Martz, L.W. 1999. An overview of
TOPAZ: an automated digital landscape analysis tool
for top. evaluation, drainage identification, watershed
segmentation, and subcatchment parameterization
[online]. URL: http://homepage.usask.ca/~lwm885/
topaz/overview.html#4.9 [cit. 2008-08-10].
GRASS Development Team 2008. Raster data
processing in GRASS GIS [web page, last modified:
2008-04-26].
URL:
http://grass.itc.it/gdp/html_
grass63/rasterintro.html [cit. 2008-08-10].
Grimaldi, S., Nardi, F., Di Benedetto, F., Istanbulluoglu,
E., Bras, R.L. 2007. A physically-based method for
removing pits in digital elevation models. Advances in
Water Resources 30(10): 2151-2158.
Holmgren, P. 1994. Multiple flow direction algorithms
for runoff modelling in grid based elevation models:
an empirical evaluation. Hydrological Processes 8(4):
327-334.
Ijjasz-Vasquez, E.J., Bras, R.L. 1995. Scaling regimes of
local slope versus contributing area in digital elevation
models. Geomorphology 12(7): 299-311.
Jenson, S.K., 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.
Jones, R. 2002. Algorithms for using a DEM for
mapping catchment areas of stream sediment samples.
Computers&Geosciences 24(4): 315-323.
Kim, S., Lee, H. 2004. A digital elevation analysis: a
spatially distributed flow apportioning algorithm.
Hydrological Processes 18(10): 1777-1794.
Lindsay, J.B. 2005. The Terrain Analysis System: a tool
for hydro-geomorphic applications. Hydrological
Processes 19(5): 1123-1130.
Lindsay, J.B., Creed, I.F. 2005. Removal of artefact
depressions from digital elevation models: towards a
minimum impact approach. Hydrological Processes
19(16): 3113-3126.
Lindsay, J.B., Creed, I.F. 2006. Distinguishing actual
and artefact depressions in digital elevation data.
Computers&Geosciences 32(8): 1192-1204.
Mark, D.M. 1984. Automated detection of drainage
networks from digital elevation models. Cartographica
21(2-3): 168-178.
Martz, L.W., De Jong, E. 1988. CATCH: a FORTRAN
program for measuring catchment area from digital
elevation models. Computers&Geosciences 14(5):
627-640.
Martz, L.W., Garbrecht, J. 1992. Numerical definition of
drainage network and subcatchment areas from digital
elevation models. Computers&Geosciences 18(6):
747-761.
Martz, L.W., Garbrecht, J. 1999. An outlet breaching
algorithm for treatment of closed depressions in a
raster DEM. Computers&Geosciences 25: 835-844.
Meisels, A., Raizman, S., Karnieli, A. 1995.
Skeletonizing a DEM into a drainage network.
Computers&Geosciences 21(1): 187-196.
Montgomery, D.R., Dietrich, W.E. 1994. Landscape
dissection and drainage area-slope thresholds. In:
Kirkby M.J. (ed.), Process models and theoretical
geomorphology. Wiley, New York, USA.
Montgomery, D.R., Foufoula-Georgiou, E. 1993.
Channel network source representation using digital
elevation models. Water Resources Research 29(12):
3925-3934.
Moore, I.D., Grayson, R.B., Ladson, A.R. 1991. Digital
terrain modelling: a review of hydrological,
geomorphological, and biological applications.
Hydrological Processes 5(1): 7-34.
O’Callaghan, J.F., Mark, D.M., 1984. The extraction of
drainage networks from digital elevation data.
Computer Vision, Graphics, and Image Processing
28(3): 323-344.
Olivera, F., Furnans, J., Maidment, D.,Djokic, D., Ye, Z.
2002. Drainage systems. In: Maidment, D. (ed), Arc
Hydro: GIS for water resources. ESRI, New York,
USA, pp 55-86.
Peucker, T., Douglas, D.H. 1975. Detection of surfacespecific points by local parallel processing of discrete
terrain-elevation data. Computer Vision, Graphics, and
Image Processing 4: 375-387.
Pilesjö, P., Zhou, Q., Harrie, L. 1998. Estimating flow
distribution over Digital Elevation Models using a
67
V. Barták: Journal of Landscape Studies 2 (2009), 57 – 68
form-based algorithm. Geographic Information
Sciences 4(1-2): 44-51.
Quinn, P., Beven, K., Chevallier, P., Planchon, O. 1991.
The prediction of hillslope flow paths for distributed
hydrological modelling using digital terrain models.
Hydrological Processes 5(1): 59-79.
Quinn, P.F., Beven, K.J., Lamb, R. 1995. The ln(a/tanβ)
index: how to calculate and how to use it within the
TOPMODEL framework. In: Beven K.J. (ed.),
Distributed hydrological modelling – applications of
the TOPMODEL concept. Wiley, New York, USA.
Rieger, W. 1998. A phenomenon-based approach to
upslope contributing area and depressions in DEMs.
Hydrological Processes 12(6): 857-872.
Roth, G., La Barbera, P., Greco M. 1996. On the
description of the basin effective drainage structure.
Journal of Hydrology 187: 119-135.
Seibert, J., McGlynn, B.L. 2007. New triangular multiple
flow direction algorithm for computing upslope areas
from gridded digital elevation models [online]. Water
Resources Research 43, W04501, doi:10.1029/
2006WR005128.
Soille, P. 2004. Optimal removal of spurious pits in
digital elevation models [online]. Water Resources
Research 40, W12509, doi:10.1029/2004WR003060.
Soille, P., Gratin, Ch. 1994. An efficient algorithm for
drainage network extraction on DEMs. Journal of
Visual Communication and Image Representation
5(2): 181-189.
Soille, P., Vogt, J., Colombo, R. 2003. Carving and
adaptive drainage enforcement of grid digital elevation
models [online]. Water Resources Research 39, 1366,
doi:10.1029/2002WR001879.
Tarboton, D.G. 1997. A new method for determination
of flow directions and upslope areas in grid digital
elevation models. Water Resources Research 33(2):
309-317.
Tarboton, D.G., Bras, R.L., Rodriguez-Iturbe, I. 1991.
On the extraction of channel networks from digital
elevation data. Hydrological Processes 5(1): 81-100.
Tarboton, D.G., Bras, R.L., Rodriguez-Iturbe, I. 1992. A
physical basis for drainage density. Geomorphology
5(59-76): 59-76.
Tarboton, D.G., Ames, D.P. 2001. Advances in the
mapping of flow networks from digital elevation data.
Paper presented at World Water and Environmental
Resources Congress, Orlando, Florida [online]. URL:
http://www.neng.usu.edu/cee/faculty/dtarb/tarpubs.ht
m [cit. 2008-04-05].
Tianqi, A., Takeuchi, K., Ishidaira, H., Yoshitani, J.,
Fukami, K. 2003. Development and application of a
new algorithm for automated pit removal for grid
DEMs. Hydrological Sciences Journal 48(6): 985-997.
Tribe, A. 1992. Automated recognition of valley lines
and drainage networks from grid digital elevation
models: a review and a new method. Journal of
Hydrology 139(1-4): 263-293.
68
Turcotte, R., Fortin, J.-P., Rousseau, A.N., Massicotte,
S., Villeneuve, J.-P. 2001. Determination of the
drainage structure of a watershed using a digital
elevation model and a digital river and lake network.
Journal of Hydrology 240: 225-242.
Wang, L., Liu, H. 2006. An efficient method for
identifying and filling surface depressions in digital
elevation models for hydrologic analysis and
modelling. International Journal of Geographical
Information Science 20(2): 193-213.