Proc. of Monitoring, Prediction and Mitigation of Water-Related Disasters (MPMD2005), pp. 85-90, 2005

DEM Based Multi-Directional Flow Path Mapping Using the Raft Method
1

2Roshan Shrestha , Yasuto Tachikawa and Kaoru Takara3
1 Graduate School of Urban and Environmental Engineering. Kyoto University, Japan
2 Associate Professor, DPRI, Kyoto University, Uji, Kyoto 611-0011, Japan
3 Professor, DPRI, Kyoto University, Uji, Kyoto 611-0011, Japan

ABSTRACT: Digital elevation models (DEM) are extensively used in hydrological analysis to obtain the
direction of flow on a topographical surface. Knowing the direction of flow helps in determining channel
networks, and in obtaining the distributed specific catchments, which are important hydrological attributes in a
DEM based analysis. There are certain demerits in the known methods of DEM-based flow-path mapping,
despite the ease of handling the DEM in present generation computational capacity. A major demerit is the
problem of accounting multidirectional flow path mapping efficiently. In this study, a new approach of assigning
the flow direction is proposed. This approach considers the flow as a vector component and argues that thedirection of flow can be different from the direction of steepest slope on a planar surface depending on surface
roughness and flow vectors. The direction of flow, if assigned to the direction of the maximum flow tendency,
gives a better flow path mapping. In this case, the direction is not necessarily restrict either on grid or acrossdiagonal of the grid-centers, which enables to assign the flow direction toward multiple downstream cells. Thismethod is named as the Ranked Flow Tendency (RAFT) method. This method was tested using the DEM of
Kamishiiba Reservoir Site (210 km2) as the case. The results show the RAFT method is better than theconventional D8 method and modified D8 method. The RAFT method identifies the network of channels andlakes and can incorporate the flow dispersion properly. These abilities show higher possibility of yielding abetter result while used this method in sediment and pollutant flow simulation and water quality modeling.
1

INTRODUCTION

Many studies in hydrological science use Digital Elevation Model (DEM) for describing the catchment
topography and obtaining topographical attributes. A network of flow channels, which the DEM can produce
easily, is the most widely used product of DEM in hydrological studies. In addition to that, the DEM provide thebasic data for calculating the upslope area, specific catchment, and slope driven parameters in flow routing andsediment and contaminant movements (Beven et al. 1984;Beven and Kirkby 1979;Costa-Cabral and Burges
1994;Moore et al. 1991;Tarboton et al. 1991;Wood et al. 1990). It is therefore understandable the significant
importance of the method, which calculated the flow direction from the DEM.
The earliest and simplest method for specifying the directions of flow in a catchment using the DEM is theD8 method (O'Callaghan and Mark 1984). This method assigns the direction of flow same to the direction ofsteepest downward slope. For each cell the direction is assigned toward one of its eight neighbors. This method
has been widely used in hydrology for flow direction mapping and to evaluate hydrological attributes
Several disadvantages and limitation are reported while using D8 method (Costa-Cabral and Burges
1994;Moore et al. 1991;Tarboton 1997). Alternatives to D8 method have been investigated and tried by several
researchers. Following methods can be listed as the sequential advancement in this direction. Multiple flow
direction method (Quinn et al. 1991) (usually termed as MS method) recognizes Multiple Slopes (MS) to
allocate fraction of flow proportional to slope downstream. Another similar method presented by (Freeman
1991) uses the slope to an exponent. The MS methods are criticized for drawbacks of too much dispersed flow.
Associating a probability which refers to an aspect angle as same as that of expected flow direction, Rho8
method (Fairfield and Leymarie 1991) is suggested to obtain the flow direction, which is criticized on its
disability to reproduce the result and random wiggles. Improved Rho8 method (Lea 1992), called Lea’s method,
has routed the flow as a rolling ball released on a plane from the center of each grid cell. DEMON (Costa-Cabral
and Burges 1994) has advanced the Lea’s method. Both Lea’s method and DEMON are questioned on planefitting technique that may mislead the determined flow direction.
A triangular facet based D∞ method (Tarboton 1997) calculates steepest descent for the triangular facetssame as that of triangular irregular network (TIN) surface slope (Tachikawa et al. 1994) by constructing facet
from DEM (Douglas 1986). The distribution weights are evaluated on the basis of surface angle of the steepest

contents pagedescent to shift the steepest flow line to nearby edge of triangular facet. This method is asserted as compromise
between D8 and Lea’s method.
(Shiiba et al. 1999) presented a method to trace flow direction from two sides by tracing the steepest
descending lines and steepest ascending lines. Overlaying the layers of steepest descending lines and opposite ofascending lines, the flow-diverging grid-cells are traced. The flow is distributed to multiple downstream cell
receivers where the descending line and opposite of ascending line do not overlap each other. This method canbe said as further development of the multidirectional approach (Freeman 1991) and pure stream lines(Douglas1986).
MS methods (Freeman 1991;Quinn et al. 1991) have included all the down stream cells. (Costa-Cabral andBurges 1994) commented these multiple flow direction methods for their discontinuous nature while suggesting
DEMON. (Tarboton 1997) illustrated an erroneous result from the Lea’s method with an example, which may
appear with DEMON too and argued that MS methods are free from grid bias. In D∞ method, the receiver twocells are always together. Proposed method of (Shiiba et al. 1999) gives freedom to the receiver cells such thatthey should not necessarily be connected as that of the D∞ method.
Most of the alternate methods apply strategy of recognizing multiple downstream receiver cells, if they exist,
unlike recognizing a single downstream receiver cell in D8 method. The flow direction is then utilized to obtainsome distribution factors that represent the fractions of flow from the source to receiver cells either adjacent or
diagonal. It is therefore clear that, flow representation by using the distribution factors is the generally accepted
method for using the grid DEM surfaces to determine flow direction. However, the approaches to evaluate the
distribution factors are different and have prime concern.
Demerits, those the hydrologists want to remove, are still there even after developing many methods. Searchof an efficient method that attempts to address the encountered weaknesses preserving the strengths has inspired
to think of a new method to define the flow direction.
A new method, named as ranked flow tendency (RAFT) method, to define the flow direction is presented inthis paper. This is a multi-directional method, which has attempted to consider the flow mechanism referred toland surface slope to obtain the flow directions. The proposed method is tested at Kamishiiba reservoir with 210
sq. km. catchments in Kyushu, Japan. The results are compared against conventional D8 method (O'Callaghanand Mark 1984) and modified D8 method.
2

RANKED FLOW TENDENCY METHOD

2.1. Basic Concept
While the slope drives to cause flow, the non-linear relation between flow and slope represents the tendency offlow along the direction of slope, but the magnitude of the tendency is non-linearly dependent. The flow
tendency may be quantified to represents the possibility to occur flow on each infinite strip of slope direction,
which is termed here as “Virtual Flow Tendency” (VFT). The VFT has its own magnitude and direction;
therefore, it is a vector.
An area with variable sloped surfaces, divided into multiple infinite strips of uniform sloped surfaces, hasmany VFTs for each of these strips. A resultant vector of the multiple strip VFT vectors should represent theresultant VFT vector for that area, which is given by∫φ f d φ
i

FR =

(1)

∫i

Here, fi is ith strip VFT, which is a function of surface characteristic and slope. f is spread angle of the surface
and FR is a resultant VFT. The orientation of FR vector on horizontal plane represents the gross flow direction onthat area.
Chezy's relation Q ∝ Sm may be referred to define the nonlinear dependency of slope S with discharge Q. An
exponent m = 0.5 to the slope thus represents the VFT magnitude, which may be generally applicable. However,a different m value, for example 0.55, 0.6, etc., may sometime suit to some specific topographies (Freeman
1991).
Manning has introduced the surface roughness coefficient recognizing the effect of surface characteristic onflow velocity. The surface characteristic affects the VFT too. However, in the present experiment, this part isomitted assuming uniform surface characteristic.

862.2 Flow direction within triangular facet
Using a DEM, infinite numbers of strips cannot be involved in calculation. Inside a 9 grid cells’ window, one
can easily form 8 planar triangular facets by connecting the center point of all 8-neighbor grid cells to the centerof central grid cell (see Fig.1). Adopting the triangular facets as primary regions of flow surface, their edge strip
VFTs give the resultant VFT vector FR as a cross product of the two VFTs for the triangular surface. The
orientation of resultant slope SR may differ from that of FR (see Fig.2 and 3) Thus in a 9 grid cells' window, 8
directions of FR vectors are obtained. The FR vector which when drawn outward from the center may be at anangle that lies within or outside of the facet angle of 45o range at center point. If the VFT angle q is within thefacet angle, it represents the resultant flow direction on that facet. If the VFT angle q is outside the facet, thedirection associated with that facet is taken along the steepest edge.

Figure 1. Ranked Flow Tendency search

Figure 2.

Resultant slope vector

Figure 3.

Resultant VFT vector

2.3 Magnitude of flow tendency
To quantify the magnitude of flow tendency along the direction q of FR on that plane, the following equation isused.

F = s1 cos 2θ + s 2 sin 2θ(2)Here, F is the flow tendency (see Fig.1), s1 and s2 are adjacent and diagonal slopes respectively (see Fig.2),
which are assigned zero if negative. The maximum F is then found by comparing the F of all 8 triangular facets.
The direction corresponding to the maximum F determines the direction of flow from the central cell of thewindow.

87

2.4 Multiple flow potentials
By choosing the maximum F as the flow direction acknowledges only two downstream cells, at maximum, as thereceiver cells. To include additional downstream cells, if the receiver cells exist, should yield a better result. The
significances of including more receiver cells in a diverged flow mechanism have already been displayed byQuinn et al. (1991), Freeman (1991), Tarboton (1997) and Shiiba et al. (1999).
Shiiba et al. (1999) reported that their proposed method produced insignificant change in result. From an
investigation using arbitrary DEM data, it was revealed that the method designated multiple flow directions often
at local peak, which mostly remains dry on flow routing. Ascending and descending lines mostly overlap at
off-peak wet regions and premises of major channel, which might have caused insignificant change in results.
The proposed VFT approach attempts to fulfill the need of multidirectional divergence at off-peak region. It
involves multiple higher Fs to enable tracing almost of the off-peak cells. The multiple Fs are chosen on thebasis of their rank from the maximum (see Fig.1) that designates this method as Ranked Flow Tendency (RAFT)
method. The RAFT method provides facility to limit the maximum number of receiver cells. For example, themaximum four lower cells, those connected with two triangular facets, are traced by choosing the highest andsecond highest Fs to incorporate much larger wet zones.
2.5 Fraction of flow
The flow direction information in multidirectional mode is generally referred to calculate the fraction of flow
to distribute among the receiver cells along either grid axes or diagonal. For a triangular facet, the fractions offlow represent the flow along edges 45o apart. When a resultant VFT overlaps triangular facet's one edge, theother edge should have no flow and vice versa. To evaluate fraction reflecting these behavior, the following
techniques are applied.
• For the ease, the angle is always measured from axes, not from diagonals.
• cos 2θ represents weight (w1) to axis line and sin 2θ represents weight (w2) to diagonal line.
• Multiply the weights by respective magnitude of F to accommodate multiple flow tendencies.
• The individual weights are divided by sum of weights to get fraction to corresponding lines.
2.6 Pit removal
DEM data contains pits and flats either representing the naturally existed pits or due to some error within
itself. The VFT remains negative (that represents no down slope) in such cells. The pits need preprocessing tofill up such that a leveled surface is formed by referring the lowest neighbor cell. The adjacent flat surface gridcells are linked together to behave as a flat region. An iterative checking is performed whether the flat regions
are again behaving as a pit or not. If, entire flat region is behaving as a pit, then the region is broadened, until it
detects a cell that drains to a neighbor cell, which ultimately drains to a lower elevation. The water level insidepits and flats are maintained uniform every time allowing them to flow any direction.
3 APPLICATION AND RESULTThe proposed method is applied to Kamishiiba reservoir site using a DEM of 50-m grid resolution. This data isreferred from Digital Map of Geological Survey Institute, Japan. After necessary preprocessing of pit removals,
the flow directions and corresponding F values are evaluated. Then the connection relations among the grid cells
are defined by choosing F1 and F2. The distribution weights are determined on the basis of flow direction andchosen F1 and F2 to allocate fractions of flow to the connected cells.
In order to visualize the flow network, the flow is accumulated from upstream to downstream side. This gives
a flow accumulation map and a density map of flow accumulation. The density map plots darker segment forhigher accumulation of water and vice versa resulting a flow path map. This method works for any unidirectional
method, such as the D8 method. However, this method cannot simply be used in any multidirectional methods,
including the RAFT as it produces cyclic accumulation through the diverged flow path. Instead, by employing ahydrological model, it is possible to obtain a flow density map without cyclic conflict. The method is describedin detail by (Shrestha et al. 2003). Corresponding density maps of flow movement are plotted in Figs. 4 and 5 forcomparison.

88

4 WHAT ARE THE DIFFERENCES?
Figures 4 and 5 presented here for comparison, clearly show the difference between the D8, modified D8 and the
new RAFT methods. The RAFT method (Fig. 5) has performed better to represent the topographic
characteristics and flow map that are clearly visible. Many disconnected river segments that appeared in D8
method (Fig. 4) are improved by the RAFT method (Fig. 5). RAFT method has displayed improved modeling of
the flow network along the multiple sets of adjacent off-peak cells unlike that in D8 method.

Figure 4. Density map of flow accumulation using D8 method

Figure 5.

Density map of flow using RAFT method

89

The parallel lines, which are generated on reservoir site by D8 method, are disappeared in RAFT method, butinstead, a uniform water level encircled the reservoir area properly. The D8 method with modified pit removal
algorithm has not produced parallel lines, and it still fails to give well result as that of the RAFT method. Thus,
the proposed method displays promising ability to handle the lake/reservoir. Parallel flow line representation oflake often is matter of headache to hydrologists as they produce erroneous interpretations in hydrologicalsimulation.
Only the flow direction issues, of course, not necessarily are enough to define the correct hydrologicalprocess because the hydrological model’s algorithm and parameters do matter a lot. Nevertheless, to analyze
DEM based hydrological process without correct representation of flow direction may be analogous to a blind
man’s walking without his white stick.
5

CONCLUSION

There is a remarkably great attraction of using DEM in hydrological analysis. To represent the correct flow
condition in DEM based analysis is, therefore, an important work for a hydrologist. A new RAFT-method is
proposed and tested in Kamishiiba reservoir site, Kyushu, Japan.
The RAFT method has displayed better representation of channel networks. The RAFT method is able torepresent the diverge distribution of flow and converge distribution as well. It has displayed promising ability ofdisplaying and modeling the lake/reservoir site. The RAFT method is found able to trace most of wet regions in
the tested area.
Modified D8 method gives significant indication of rapid flow down unlike the RAFT method.
Multidirectional consideration of RAFT method should represent the real water movement more correct than themodified D8 method.
REFERENCES:

Beven, K. J. and Kirkby, M. J. 1979. A physically based variable contributing area model of basin hydrology. Hydrological
Science Bulletin 24: 43-69.
Beven, K. J., Kirkby, M. J., Schofield, N. and Tagg, A. F. 1984. Testing a physically based flood forcasting model(TOPMODEL) for three UK catchments. Journal of Hydrology 69: 119-143.
Costa-Cabral, M. and Burges, S. J. 1994. Digital elevation model networks (DEMON): A model of flow over hillslopes forcomputation of contributing and dispersal areas. Water Resources Research 30: 1681-1692.
Douglas, D. H. 1986. Experiments to loacate ridges and channels to create a new type of digital elevation model.
Cartographica 23: 29-61.
Fairfield, J. and Leymarie, P. 1991. Drainage networks from grid digital elevation models. Water Resources Research 27:
709-717.
Freeman, T. G. 1991. Calculating catchment area with divergent flow based on a regular grid. Computational Geoscience 17:413-422.
Lea, N. L. 1992. An aspect driven kinematic routing algorithm. In: A. J. Parsons and A. D. Abrahams, Editors. Overland
Flow: Hydraulic and Erosion Mechanics. Chapman & Hall, New York.
Moore, I. D., Grayson, R. B. and Ladson, A. R. 1991. Digital terrain modeling; A review of hydrological, geomorphological,
and biological applications. Hydrological Processes 5: 3-30.
O'Callaghan, J. F. and Mark, D. M. 1984. The extraction of drainage networks from digital elevation data. Comput. Vision
Graphics Image Processing 28: 328-344.
Quinn, P., Beven, K., Chevallier, P. and Planchon, O. 1991. The Prediction of Hillslope Flow Paths for Distributed
Hydrological Modeling Using Digital Terrain Models. Hydrological Processes 5: 59-79.
Shiiba, M., Ichikawa, Y., Sakakibara, T. and Tachikawa, Y. 1999. A new numerical representation form of basin topography.
Jounal of Hydraulic, Coastal and Environmental Engineering 621: 1-9.
Shrestha, R. K., Tachikawa, Y. and Takara, K. 2003. A new method to define flow direction using grid digital elevation
models. Annual Journal of Hydraulic Engineering, JSCE 47: 241-246.
Tachikawa, Y., Shiiba, M. and Takasao, T. 1994. Development of A Basin Geomorphic Information-System Using A
Tin-dem Data Structure. Water Resources Bulletin 30: 9-17.
Tarboton, D. G. 1997. A new method for the determination of flow directions and upslope areas in grid digital elevation
models. Water Resources Research 33: 309-319.
Tarboton, D. G., Bras, R. L. and Rodrigueziturbe, I. 1991. On the Extraction of Channel Networks from Digital Elevation
Data. Hydrological Processes 5: 81-100.
Wood, E. F., Sivapalan, M. and Beven, K. 1990. Similarity and Scale in Catchment Storm Response. Reviews of Geophysics
28: 1-18.

90

contents page