- Email: [email protected]
Application of the PROMETHEE technique to determine depression outlet location and flow direction in DEM
Journal of Hydrology 287 (2004) 49–61 www.elsevier.com/locate/jhydrol
Application of the PROMETHEE technique to determine depression outlet location and flow direction in DEM Tien-Yin Choua, Wen-Tzu Linb, Chao-Yuan Linc, Wen-Chieh Choud,*, Pi-Hui Huange a
Department of Land Management, Feng-Chia University, Taichung 407, Taiwan, ROC Graduate Institute of Environmental Planning and Design, Ming-Dao University, Changhua 523, Taiwan, ROC c Department of Soil and Water Conservation, National Chung-Hsing University, Taichung 402, Taiwan, ROC d Department of Civil Engineering, Chung-Hua University 707, Sec 2 Wu-Fu Rd, Hsinchu 300, Taiwan, ROC e Geographic Information Systems Research Center, Feng-Chia University, Taichung City 407, Taiwan
b
Received 27 August 2002; revised 22 September 2003; accepted 26 September 2003
Abstract With the fast growing progress of computer technologies, spatial information on watersheds such as flow direction, watershed boundaries and the drainage network can be automatically calculated or extracted from a digital elevation model (DEM). The stubborn problem that depressions exist in DEMs has been frequently encountered while extracting the spatial information of terrain. Several filling-up methods have been proposed for solving depressions. However, their suitability for large-scale flat areas is inadequate. This study proposes a depression watershed method coupled with the Preference Ranking Organization METHod for Enrichment Evaluations (PROMETHEEs) theory to determine the optimal outlet and calculate the flow direction in depressions. Three processing procedures are used to derive the depressionless flow direction: (1) calculating the incipient flow direction; (2) establishing the depression watershed by tracing the upstream drainage area and determining the depression outlet using PROMETHEE theory; (3) calculating the depressionless flow direction. The developed method was used to delineate the Shihmen Reservoir watershed located in Northern Taiwan. The results show that the depression watershed method can effectively solve the shortcomings such as depression outlet differentiating and looped flow direction between depressions. The suitability of the proposed approach was verified. q 2003 Elsevier B.V. All rights reserved. Keywords: Digital elevation model; Preference Ranking Organization METHod for Enrichment Evaluation; Hydrological parameter; Digital terrain analysis; TOPAZ; Depression-filling
1. Introduction Digital terrain model (DTM) represents the spatial distribution of terrain characteristics (Doyle, 1978). When there is only one property, elevation, the term * Corresponding author. Fax: þ 886-3-537-2188. E-mail address: [email protected] (W.-C. Chou). 0022-1694/$ - see front matter q 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2003.09.026
‘digital elevation model’ (DEM) is used (Collins and Moon, 1981). It can be used to derive a wealth of information about the morphology of a land surface (US Geological Survey, 1987). Because depressions sometimes exist in DEMs, a stubborn problem exists involving the automatic extraction of spatial information from the terrain. Some depressions are data errors introduced in
50
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
the surface generation process, while others represent real topographic features such as quarries or natural potholes (Jenson and Domingue, 1988). Consequently, a few depression-filling methods were presented to remove depressions from a DEM. O’Callaghan and Mark (1984) attempted to remove depressions by smoothing the DEM data. The smoothing approach removes shallow depressions, but deeper depressions remain (Jenson and Domingue, 1988). Marks et al. (1984) and Jenson and Trautwein (1987) presented an algorithm for filling depressions by increasing the values of cells in each depression to the value of the cell with the lowest value on the depression boundary. There are three situations in which flow direction assignment becomes problematic: (1) at cells where there are two or more possible flow directions; (2) in flat areas, where all of the cells have the same elevation; (3) at pits or closed depressions in the DEM (Tribe, 1992). Jenson and Domingue (1988) pointed out that depression-filling algorithms will change the natural terrain, enlarge the depression, loop the depression and/or produce an outflow point in the depression while processing a flat area. It is obvious from the calculated results from various depression-processing methods that watershed delineation, drainage network extraction and hydrologic event simulation accuracy will be affected. Recently, TOPAZ method proposed by Garbercht and Martz (1997) that moves generally away from higher and towards lower terrain in natural landscape drainage, and then repeatedly increases the depression elevation for drainage towards lower terrain. This approach can actually improve the shortcomings like depression outlet differentiating better than other depression-filling algorithms. However, the qualitative analysis for this method cannot be easily used to analyze the optimal outlet and flow direction for a large-scale flat area. Martz and Garbrecht (1998) proposed other algorithms that allow breaching of depression outlets and consider distribution of both higher and lower elevations in assigning flow directions on flat areas. Applying these algorithms can produce more realistic results in application. In this study, the watershed boundary delineated using the TOPAZ method was also compared to other methods. The objective of this study is to propose a new algorithm for determining the optimal outlet and flow
direction over large-scale flat area. The first step is to derive the incipient flow direction, non-depression flow direction, using an elevation-differentiating method coupled with a surface-inclining method. The second step is to calculate the optimal outlet and flow direction in depressions using a depression watershed method with the Preference Ranking Organization METHod for Enrichment Evaluations (PROMETHEEs) theory proposed by Brans et al. (1984). The third step is to generate the depressionless flow direction. To verify the suitability of proposed method, it was applied to delineate the watershed boundary of the Shihmen Reservoir and compared with other depression-processing methods.
2. Depression watershed method Several software systems have been developed using different depression-filling algorithms to calculate the flow direction, delineate the watershed and extract the drainage network. Methods such as the HYDROLOGY module in ArcView, the DWCON and TERRAIN ANALYSIS modules in the EASI/PACE (Chang, 2002) and TOPAZ landscape analysis tools have been used. As mentioned above, these depression-filling algorithms are not suitable for large-scale flat areas. The large-scale flat area are combining for several depressions nearby, the optimal outlet and flow direction are not easily determined in these areas. Exiting methods may result false simulation, Fig. 11 in the case study shows the unsuitable reason for some popular methods. For this reason, this study will propose an improved method to determine the optimal outlet and calculate the flow direction in depressions. Three procedures are used as follows. 2.1. Calculating the incipient flow direction A 3 by 3 moving window is used over a DEM to locate the flow direction for each cell. The elevation difference ðDhÞ from the center cell of the window to one of its eight neighbors is calculated as: ( Dh ¼
z5 2 zi ;
pffiffi ðz5 2 zi Þ= 2;
i ¼ 2; 4; 6; 8 i ¼ 1; 3; 7; 9
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
If Dh is less than zero, a zero flow direction is assigned to indicate that the flow direction is undefined. If Dh is greater than or equal to zero and occurs at one neighbor, the flow direction is assigned to that neighbor. If Dh is greater than zero and occurs at more than one neighbor, the flow direction is assigned logically according to the surface-incline method, which can be derived as follows. (1) Calculating the angle of inclination. If the coordinate of z5 is defined as (0, 0), the neighbor coordinates of z1 ; z2 ; z3 ; z4 ; z6 ; z7 ; z8 ; z9 are (2 1, 1), (0, 1), (1, 1), (2 1, 0), (1, 0), (2 1, 2 1), (0, 2 1) and (1, 2 1), respectively. The plane equation for a 3 £ 3 moving window can be expressed as z^ ¼ a þ bx þ cy: The angle of inclination ðuÞ for the plane is calculated using the least square method as u ¼ tan21 ðc=bÞ (Horn, 1981). (2) Determining flow direction of the center cell. The flow direction of the center cell is determined by the parameters b, c and u, shown in Fig. 1. If b is greater than zero, the elevation at the right neighbors is higher than that at the left neighbors. If c is greater than zero, the elevation at the top neighbors is higher than that at the bottom neighbors. In Fig. 1, the elevation of the colored cells is lower than the other cells, so the flow direction is assigned to the lowest elevation cell. If the colored cells have the same elevation, then the flow direction is determined by u.
Fig. 1. Illustration of determining the z5 flow direction.
51
2.2. Determining the optimal depression outlet using the depression watershed method After calculating the incipient flow direction, pits and/or flat areas coded as zero are necessary for advanced processing. Because drainage moves generally away from higher terrain towards lower terrain, the depression surfaces that include pits and/or flat areas will form a closed watershed, shown in Fig. 2. Four main factors (h1 ; h2 ; l1 and l2 ) can be used to determine the optimal depression watershed outlet. The h1 factor is the difference in elevation from the edge cells to the depression cell. The h2 factor is the difference in elevation from the edge cell to its neighbor cells with the lowest elevation, but not in the depression watershed. The l1 factor is the distance from the edge cells to the depression cell. The l2 factor is the distance from the edge cells to their neighbor cells that have the lowest elevation, but not in the depression watershed. The depression watershed outlet function (O) can be expressed as: O ¼ f ðh1 ; h2 ; l1 ; l2 Þ The priority for determining the optimal outlet from the four factors is h1 ¼ h2 . l1 . l2 : If the h1 and/or h2 of a cell on the edge of the depression watershed is the lowest and its l1 and/or l2 is the shortest, that cell will be the optimal outlet for the depression watershed. Multiple criteria decision-making methods are suitable for analyzing this problem because each cell around the edge of a depression watershed is likely to
Fig. 2. Illustration of depression watershed.
52
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
become the optimal candidate cell for the depression watershed outlet and be evaluated with these four factors. The PROMETHEE method presented by Brans et al. (1984) is an interactive multiple criteria decision-making technique designed to handle qualitative and discrete alternatives. In this study, the PROMETHEE method was chosen to quantify and determine the priority of the candidate cells for the depression watershed outlet. The steps applied to derive the optimal outlet using the PROMETHEE method are described as follows. (1) Establishing an alternatives and criterion matrix If the edge around a depression watershed has n candidate cells, there will be n alternatives ða1 ; a2 ; …; an Þ: Each candidate cell has four evaluative factors (h1 ; h2 ; l1 and l2 ). The evaluative function can be written as f1 ð·Þ; f2 ð·Þ; f3 ð·Þ and f4 ð·Þ: Both the alternatives and criterion can be expressed as a n by 4 ~ matrix T:
Brans et al. (1984) presented the shape of the six possible types of generalized criteria to assist the decision-maker with this selection, as shown in Table 1. He pointed out that the GAUSSIAN criterion was selected most by users for practical applications especially for continuing data. The major criteria in this study are elevation and distance, because of the criteria containing continuity, the GAUSSIAN criterion was chosen for evaluation. In the six type function, HðdÞ ¼ 1 2 expð2d2 =2s2 Þ; s is defined as the threshold value between the indifferent and strict preference areas. In real terrain, if the elevation difference in h1 or h2 is more than 1 m or the distance in l1 or l2 is more than one cell size, the effect on determining the depression outlet will be critical. Therefore, s is defined as 1 and the GAUSSIAN criterion function can be written as HðdÞ ¼ 1 2 expð2d2 =2Þ: The evaluative difference in the ith candidate cell compared with the other candidate ~ i: cells can be expressed as a n using the 4 matrix H 1 2 expð2 d1 ðai ; a1 Þ 1 2 expð2 d2 ðai ; a1 Þ Þ 1 2 expð2 d3 ðai ; a1 Þ Þ 1 2 expð2 d4 ðai ; a1 Þ Þ 2 2 2 2 d ða ; a Þ d ða ; a Þ d ða ; a Þ d ða ; a Þ 1 2 expð2 1 i 2 Þ 1 2 expð2 2 i 2 Þ 1 2 expð2 3 i 2 Þ 1 2 expð2 4 i 2 Þ ~i ¼ H 2 2 2 2 ··· ··· ··· ··· 1 2 expð2 d1 ðai ; an Þ Þ 1 2 expð2 d2 ðai ; an Þ Þ 1 2 expð2 d3 ðai ; an Þ Þ 1 2 expð2 d4 ðai ; an Þ Þ 2 2 2 2 f1 ða1 Þ f2 ða1 Þ f3 ða1 Þ f4 ða1 Þ f ða Þ f ða Þ f ða Þ f ða Þ 1 2 2 2 3 2 4 2 T~ ¼ ··· · · · · · · · · · f1 ðan Þ f2 ðan Þ f3 ðan Þ f4 ðan Þ
(2) Selecting a preference function ðHðdÞÞ The evaluative difference of the ith candidate cell compared with the other candidate cells can be expressed as a n using 4 matrix d~i : d1 ðai ; a1 Þ d2 ðai ; a1 Þ d3 ðai ; a1 Þ d4 ðai ; a1 Þ d ða ; a Þ d ða ; a Þ d ða ; a Þ d ða ; a Þ 1 i 2 2 i 2 3 i 2 4 i 2 d~i ¼ ··· · · · · · · · · · d1 ðai ; an Þ d2 ðai ; an Þ d3 ðai ; an Þ d4 ðai ; an Þ where dðai ; an Þ ¼ f ðai Þ 2 f ðan Þ:
(3) Calculating the preference index (p) For each pair of alternatives, the preference index pðdÞ can be written as pðdÞ ¼
4 X
wj Hj ðdÞ
j¼1
where wj ; j ¼ 1; 2; 3; 4 are weights associated with each criterion. Because the priority for determining the optimal outlet from the four factors is h1 ¼ h2 . l1 . l2 : In this study, the weight of h1 ðw1 Þ is defined as 1, the weights of the other factors ðwj Þ can be written as w1 ¼ w2 ¼ 1 w3 ¼
1 1 or uh1 maxðl1 Þ uh2 maxðl1 Þ
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61 Table 1 The shape of the six possible types of generalized criteria Generalized criterion type Type I: Usual criterion ( 0 d¼0 HðdÞ ¼ 1 ldl . 0
Preference function ðHðdÞÞ
w4 ¼
53
w3 maxðl2 Þ
where uh1 is the variation in the outranking flow while changing h1 ; uh2 is the variation in the outranking flow while changing h2 ; maxðl1 Þ is the maximum of l1 ; maxðl2 Þ is the maximum of l2 : The weights can be ~ expressed as a 1 by 4 matrix w: ~ ¼ lw1 w2 w3 w4 l w
Type II: U-shape criterion ( 0 ldl # q HðdÞ ¼ 1 ldl . q
Type III:8V-shape criterion > < ldl ldl # p p HðdÞ ¼ > : 1 ldl . p
Type IV: Level criterion 8 0 ldl # q > > > <1 HðdÞ ¼ q , ldl # p > 2 > > : 1 ldl . p
Type V: V-shape criterion with indifference criterion 8 > 0 ldl # q > > > < ldl 2 q HðdÞ ¼ q , ldl # p > p2q > > > : 1 ldl . p
Type VI: Gaussian criterion HðdÞ ¼ 1 2 expð2d 2 =2s2 Þ
The preference index of the ith candidate cell compared with the other candidate cells can be written as ~ i ·w· ~ S~a pi ¼ H where S~a is a n by 1 scalar matrix 1 1 ðS~a ¼ Þ: ··· 1 (4) Calculating the outranking flow (f) According to the above preference index (p), the outranking flow for the ith candidate cell can be written as:
fi ¼
n X
pi
i¼1
(5) Generating a complete ranking From the outranking flow calculation, the candidate cell that has the smallest outranking flow will be the optimal outlet for the depression watershed. 2.3. Calculating the depressionless flow directions Two steps can be used to generate a new flow direction for the outlet and depression cells after determining the optimal outlet using the depression watershed method. (1) Changing the flow direction of neighbor cells to the outlet cell that will not flow into the outlet cell and reversing the flow directions of those cells located on the flow route from the outlet cell to the depression. (2) Assigning the depression flow direction to the outlet cell using the recursion algorithm.
54
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
3. Example The Shihmen Reservoir watershed in Northern Taiwan was used in this study (Fig. 3). The focus flat area (Fig. 4a – d) of lower Shihmen Reservoir watershed, contains major landform with reservoir encroachment, farmland and hardwood forest. A 16 by 20 line DEM data sample (Fig. 5) was chosen with average slope 2.97%, area 51.2 ha, no aspect area 56.1%. The geological data from Taiwan’s Central Geological Service is categorized Pleistocene Toen formation. It is composed of plateau gravel and clay with a thin layer of yellowish lateritic earth on top. It is also considered to be thin deposits on river terraces formed after the Tyuerki formation on the adjacent higher terraces. 3.1. Calculating depressionless flow directions The calculated depressless flow directions within sample area were shown in Fig. 6. The procedures for the depression watershed method can be divided into four portions as follows: (1) Calculating the incipient flow directions The non-depression flow directions can be calculated using elevation-differencing and surface-inclining methods. (2) Searching for depression watersheds
Fig. 4. (a) Elevation distribution of the study area. (b) Slope distribution of the study area. (c) Aspect distribution of the study area.
Fig. 3. Study area.
After calculating the incipient flow direction, searching for a depression watershed can be done by tracing the upstream cells that flow into the depression. Fig. 7 shows that there are four depression watersheds A, B, C and D in the DEM. According to the terrain surrounding the depression watersheds, the topography of the A and B depression watersheds are higher. The D depression watershed is lower. (3) Determining the optimal depression watershed outlet
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
55
Fig. 5. Original elevations of the sample DEM.
A breach is located at the bottom portion of the B and D depression watesheds. This breach is defined as the depression watershed outlet. The A and C depression watershed candidate cells are coded as shown in Fig. 7.
There are forty-five alternatives in the A depression watershed and twenty-one alternatives in the C depression. A pairwise comparison could provide a preference structure between these
Fig. 6. Illustration of depression watersheds.
56
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
Fig. 7. The code of the candidate cells in depression watersheds.
alternatives. The preference index with GUSSIAN function and criterion weights could then be calculated. According to the weight definitions, for h1 and h2 ; the weights can be written as: w1 ¼ w2 ¼ 1 In the l1 1 1 T~ ¼ 1 1
weight, uh1 can be obtained by substituting 2 3 1 1 1 1 1 1
into the PROMETHEE algorithm, and then maxðl1 Þ and maxðl2 Þ in the A and C depression watersheds can be derived as follows:
uh1 ¼ 1:2581 For A depression watershed: maxðl1 Þ ¼ 7:41 maxðl2 Þ ¼ 1:41 w3 ¼
1 ¼ 0:107 1:2581 £ 7:41
w4 ¼
0:107 ¼ 0:076 7:41
For C depression watershed: maxðl1 Þ ¼ 4:24 maxðl2 Þ ¼ 1:41 w3 ¼
1 ¼ 0:187 1:2581 £ 4:24
w4 ¼
0:187 ¼ 0:133 1:41
The criteria value, outranking flow and complete ranking for each candidate cell are shown in Tables 2a and 2b. In the A depression watershed, the candidate cell coded A7, with the smallest outranking flow, 2 80.719, is the optimal outlet cell. In the C depression watershed, the candidate cell coded C8, with the smallest outranking flow, 2 32.792, is the optimal outlet cell. According to the outranking flow analysis, the outlet for the A to D depression watersheds can be illustrated as shown in Fig. 8. (4) Flow direction assignment in a depression After determining the depression watershed outlet, the first procedure is assigning the flow direction from the outlet cell into a neighbor that does not flow into the outlet cell and has a lower elevation. The flow direction of cells that route from the outlet cell into the depression is then reversed. The second procedure is to assign
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
57
Table 2a The outranking flow and complete ranking of each candidate cells in A-depression watershed Alternatives
A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 A13 A14 A15 A16 A17 A18 A19 A20 A21 A22 A23 A24 A25 A26 A27 A28 A29 A30 A31 A32 A33 A34 A35 A36 A37 A38 A39 A40 A41 A42 A43 A44 A45
Criterion h1 ; w1 ¼ 1
h2 ; w2 ¼ 1
l1 ; w3 ¼ 0:107
l2 ; w4 ¼ 0:076
27 23 23 5 5 2 2 7 22 25 25 25 29 37 52 43 56 52 68 43 29 25 25 25 23 24 15 20 21 21 23 16 16 2 2 2 3 9 15 14 17 19 11 8 25
17 17 17 3 3 2 2 2 13 13 25 26 27 29 37 54 48 57 51 26 26 26 22 19 24 17 23 23 23 21 21 22 9 9 2 2 2 4 12 17 20 20 8 8 8
5.83 3.41 3.00 2.00 1.00 1.41 1.00 2.00 2.00 3.41 3.41 4.41 5.41 6.41 7.41 5.83 6.83 6.41 6.83 5.41 4.41 3.41 3.00 3.41 2.00 2.41 2.00 2.00 2.41 3.41 2.41 2.00 3.00 2.00 3.83 3.41 4.41 5.41 6.41 4.83 4.41 4.00 4.83 4.41 5.41
1.41 1.00 1.41 1.00 1.41 1.00 1.41 1.00 1.41 1.00 1.41 1.41 1.41 1.41 1.41 1.41 1.41 1.41 1.41 1.41 1.00 1.41 1.41 1.41 1.41 1.41 1.41 1.00 1.00 1.00 1.41 1.41 1.00 1.41 1.00 1.00 1.41 1.41 1.41 1.41 1.41 1.41 1.41 1.00 1.41
the flow direction of the depression cells into the outlet cell. Fig. 9 shows the results of the analyzed depressionless flow direction for the A to D depressions.
c
Rank
22.498 20.969 21.348 267.973 268.772 280.573 280.719 268.077 214.947 2.578 39.986 44.917 60.961 69.419 79.279 80.030 83.924 86.588 87.914 64.379 56.528 43.301 27.278 17.289 22.090 3.722 23.599 7.159 10.315 3.907 10.837 24.879 235.984 262.800 277.138 277.815 273.314 252.094 227.469 222.627 26.393 22.970 243.190 247.974 23.274
32 23 22 8 6 2 1 7 16 24 34 36 38 40 41 42 43 44 45 39 37 35 33 30 31 25 19 27 28 26 29 18 13 9 4 3 5 10 14 15 17 21 12 11 20
3.2. Watershed delineation application To verify the suitability of the proposed method, the depressionless flow direction was calculated using
58
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
Table 2b The outranking flow and complete ranking of each candidate cells in C-depression watershed Alternatives
C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 C13 C14 C15 C16 C17 C18 C19 C20 C21
Criterion h1 ; w1 ¼ 1
h2 ; w2 ¼ 1
l1 ; w3 ¼ 0:187
l2 ; w4 ¼ 0:133
22 24 13 6 2 2 1 1 5 9 15 20 20 15 2 1 2 3 13 21 23
21 15 20 8 4 4 1 1 2 2 4 9 23 9 3 3 3 1 1 5 16
4.24 3.83 2.83 2.41 2.00 1.00 1.41 1.00 1.00 2.00 2.41 3.41 3.00 3.41 2.41 1.41 1.00 2.00 2.41 2.41 3.83
1.41 1.41 1.41 1.41 1.00 1.41 1.00 1.41 1.00 1.41 1.41 1.41 1.41 1.41 1.41 1.00 1.41 1.00 1.41 1.41 1.41
Fig. 8. Illustration of the outlet of depression watersheds.
c
Rank
36.191 34.189 20.704 4.728 214.462 215.736 232.468 232.792 216.995 211.707 5.336 21.910 32.709 16.886 217.866 222.920 220.008 224.765 211.942 15.437 33.573
21 20 16 12 9 8 2 1 7 11 13 17 18 15 6 4 5 3 10 14 19
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
59
Fig. 9. Illustration of depressionless flow direction in a noised DEM.
for the Shihmen Reservoir using the depression watershed method, elevation-smoothing method (O’Callaghan and Mark, 1984), depression-filling method (Jenson and Domingue, 1988) and elevation-incrementing method (Martz and Garbrecht, 1998). Numerous large flat terrain areas that have more than one outflow point exist near the downstream area of the Shihmen Reservoir watershed. In these methods, the elevation-smoothing method removes only the shallow depressions. The deeper depressions remain. Therefore, these methods fail to fully accommodate flat areas. The elevationsmoothing curve for the Shihmen Reservoir watershed is shown in Fig. 10. The delineated watershed boundary using the depression watershed, elevationfilling and elevation-increasing methods are shown in Fig. 11. Fig. 11 shows that there is an obvious difference near the downstream area of the watershed. The watershed areas delineated using the depression watershed, depression-filling and elevation-incrementing methods were 75,633.60, 75,290.40 and 75,304.64 ha. To verify the suitability of the depression-processing methods, an electric map of the reservoir flood area, at a scale of 1/25,000, was
overlaid onto the watershed boundary. The overlay (Fig. 11) indicated that the result delineated using the depression watershed method was more realistic than the other methods.
Fig. 10. The convergent curve of applying elevation-smoothing method on raising depressions.
60
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
Fig. 11. Comparison of applying varied depression-processing methods on delineating the Shihmen Reservoir watershed.
4. Conclusions The proposed method in this study searches the best outlet for a depression without changing the grid elevation. It avoids the difficulty for determining outlet during increasing the depression elevation and obtaining an unrealistic watershed in a large flat area or several nearby depressions. Especially, when flat areas or depressions are located close to watershed boundary, the traditional methods may establish an incomplete watershed area. The new outlet algorithm is more complicated than usual. However, it improves the result obtained from automated watershed boundary delineation. The calculated flow directions in this study produce convolution flow patterns within each depression. A more rational flow pattern arrangement inside the depression need further compare and consult with the Jenson and Domingue (1988) or the Garbercht and Martz (1997) approach. The Shihmen Reservoir watershed in Northern Taiwan was delineated using several popular depression-processing techniques. The results indicated that the depression watershed method was more
realistic than the other methods. The algorithm proposed by this study produced quantitative analysis for depression outlet determination and can derive useful morphological information for hydrological applications.
References Brans, J.P., Mareschal, B., Vincke, P., 1984. PROMETHEE: A New Family of Outranking Methods in MCDM, IFORS 84, North Holland, pp. 477–490. Chang, K.T., 2002. Introduction to Geographic Information Systems, McGraw-Hill, New York, NY, p. 233. Collins, S.H., Moon, G.C., 1981. Algorithms for dense digital terrain models. Photogrammetric Engineering and Remote Sensing 47, 71–76. Doyle, F.J., 1978. Digital terrain models: an overview. Photogrammetric Engineering and Remote Sensing 44, 1481– 1485. Garbercht, J., Martz, L.W., 1997. Assignment of drainage direction over flat surfaces in raster digital elevation models. Journal of Hydrology 193, 204–213. Horn, B.K.P., 1981. Hill shading and the reflectance map. Proceedings of the IEEE 69, 14–47. Jenson, S.K., Trautwein, C.M., 1987. Methods and applications in surface depression analysis. Proceedings of Auto-Carto 8, 137 –144. Jenson, S.K., Domingue, J.O., 1988. Extracting topographic structure from digital elevation data for geographic information
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61 system analysis. Photogrammetric Engineering and Remote Sensing 54, 1593–1600. Marks, D., Dozier, J., Frew, J., 1984. Automated basin delineation from digital elevation sata. Geo-Processing 2, 299–311. Martz, L.W., Garbrecht, J., 1998. The treatment of flat areas and depressions in automated drainage analysis of raster digital elevation models. Hydrological Processes 12, 843–855.
61
O’Callaghan, J.F., Mark, D.M., 1984. The extraction of drainage networks from digital elevation data. Computer Graphics and Image Processing 28, 323–344. 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, 263– 293. US Geological Survey, 1987. Digital elevation models. US Geological Survey Data User’s Guide 5, p. 38.
Application of the PROMETHEE technique to determine depression outlet location and flow direction in DEM Tien-Yin Choua, Wen-Tzu Linb, Chao-Yuan Linc, Wen-Chieh Choud,*, Pi-Hui Huange a
Department of Land Management, Feng-Chia University, Taichung 407, Taiwan, ROC Graduate Institute of Environmental Planning and Design, Ming-Dao University, Changhua 523, Taiwan, ROC c Department of Soil and Water Conservation, National Chung-Hsing University, Taichung 402, Taiwan, ROC d Department of Civil Engineering, Chung-Hua University 707, Sec 2 Wu-Fu Rd, Hsinchu 300, Taiwan, ROC e Geographic Information Systems Research Center, Feng-Chia University, Taichung City 407, Taiwan
b
Received 27 August 2002; revised 22 September 2003; accepted 26 September 2003
Abstract With the fast growing progress of computer technologies, spatial information on watersheds such as flow direction, watershed boundaries and the drainage network can be automatically calculated or extracted from a digital elevation model (DEM). The stubborn problem that depressions exist in DEMs has been frequently encountered while extracting the spatial information of terrain. Several filling-up methods have been proposed for solving depressions. However, their suitability for large-scale flat areas is inadequate. This study proposes a depression watershed method coupled with the Preference Ranking Organization METHod for Enrichment Evaluations (PROMETHEEs) theory to determine the optimal outlet and calculate the flow direction in depressions. Three processing procedures are used to derive the depressionless flow direction: (1) calculating the incipient flow direction; (2) establishing the depression watershed by tracing the upstream drainage area and determining the depression outlet using PROMETHEE theory; (3) calculating the depressionless flow direction. The developed method was used to delineate the Shihmen Reservoir watershed located in Northern Taiwan. The results show that the depression watershed method can effectively solve the shortcomings such as depression outlet differentiating and looped flow direction between depressions. The suitability of the proposed approach was verified. q 2003 Elsevier B.V. All rights reserved. Keywords: Digital elevation model; Preference Ranking Organization METHod for Enrichment Evaluation; Hydrological parameter; Digital terrain analysis; TOPAZ; Depression-filling
1. Introduction Digital terrain model (DTM) represents the spatial distribution of terrain characteristics (Doyle, 1978). When there is only one property, elevation, the term * Corresponding author. Fax: þ 886-3-537-2188. E-mail address: [email protected] (W.-C. Chou). 0022-1694/$ - see front matter q 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2003.09.026
‘digital elevation model’ (DEM) is used (Collins and Moon, 1981). It can be used to derive a wealth of information about the morphology of a land surface (US Geological Survey, 1987). Because depressions sometimes exist in DEMs, a stubborn problem exists involving the automatic extraction of spatial information from the terrain. Some depressions are data errors introduced in
50
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
the surface generation process, while others represent real topographic features such as quarries or natural potholes (Jenson and Domingue, 1988). Consequently, a few depression-filling methods were presented to remove depressions from a DEM. O’Callaghan and Mark (1984) attempted to remove depressions by smoothing the DEM data. The smoothing approach removes shallow depressions, but deeper depressions remain (Jenson and Domingue, 1988). Marks et al. (1984) and Jenson and Trautwein (1987) presented an algorithm for filling depressions by increasing the values of cells in each depression to the value of the cell with the lowest value on the depression boundary. There are three situations in which flow direction assignment becomes problematic: (1) at cells where there are two or more possible flow directions; (2) in flat areas, where all of the cells have the same elevation; (3) at pits or closed depressions in the DEM (Tribe, 1992). Jenson and Domingue (1988) pointed out that depression-filling algorithms will change the natural terrain, enlarge the depression, loop the depression and/or produce an outflow point in the depression while processing a flat area. It is obvious from the calculated results from various depression-processing methods that watershed delineation, drainage network extraction and hydrologic event simulation accuracy will be affected. Recently, TOPAZ method proposed by Garbercht and Martz (1997) that moves generally away from higher and towards lower terrain in natural landscape drainage, and then repeatedly increases the depression elevation for drainage towards lower terrain. This approach can actually improve the shortcomings like depression outlet differentiating better than other depression-filling algorithms. However, the qualitative analysis for this method cannot be easily used to analyze the optimal outlet and flow direction for a large-scale flat area. Martz and Garbrecht (1998) proposed other algorithms that allow breaching of depression outlets and consider distribution of both higher and lower elevations in assigning flow directions on flat areas. Applying these algorithms can produce more realistic results in application. In this study, the watershed boundary delineated using the TOPAZ method was also compared to other methods. The objective of this study is to propose a new algorithm for determining the optimal outlet and flow
direction over large-scale flat area. The first step is to derive the incipient flow direction, non-depression flow direction, using an elevation-differentiating method coupled with a surface-inclining method. The second step is to calculate the optimal outlet and flow direction in depressions using a depression watershed method with the Preference Ranking Organization METHod for Enrichment Evaluations (PROMETHEEs) theory proposed by Brans et al. (1984). The third step is to generate the depressionless flow direction. To verify the suitability of proposed method, it was applied to delineate the watershed boundary of the Shihmen Reservoir and compared with other depression-processing methods.
2. Depression watershed method Several software systems have been developed using different depression-filling algorithms to calculate the flow direction, delineate the watershed and extract the drainage network. Methods such as the HYDROLOGY module in ArcView, the DWCON and TERRAIN ANALYSIS modules in the EASI/PACE (Chang, 2002) and TOPAZ landscape analysis tools have been used. As mentioned above, these depression-filling algorithms are not suitable for large-scale flat areas. The large-scale flat area are combining for several depressions nearby, the optimal outlet and flow direction are not easily determined in these areas. Exiting methods may result false simulation, Fig. 11 in the case study shows the unsuitable reason for some popular methods. For this reason, this study will propose an improved method to determine the optimal outlet and calculate the flow direction in depressions. Three procedures are used as follows. 2.1. Calculating the incipient flow direction A 3 by 3 moving window is used over a DEM to locate the flow direction for each cell. The elevation difference ðDhÞ from the center cell of the window to one of its eight neighbors is calculated as: ( Dh ¼
z5 2 zi ;
pffiffi ðz5 2 zi Þ= 2;
i ¼ 2; 4; 6; 8 i ¼ 1; 3; 7; 9
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
If Dh is less than zero, a zero flow direction is assigned to indicate that the flow direction is undefined. If Dh is greater than or equal to zero and occurs at one neighbor, the flow direction is assigned to that neighbor. If Dh is greater than zero and occurs at more than one neighbor, the flow direction is assigned logically according to the surface-incline method, which can be derived as follows. (1) Calculating the angle of inclination. If the coordinate of z5 is defined as (0, 0), the neighbor coordinates of z1 ; z2 ; z3 ; z4 ; z6 ; z7 ; z8 ; z9 are (2 1, 1), (0, 1), (1, 1), (2 1, 0), (1, 0), (2 1, 2 1), (0, 2 1) and (1, 2 1), respectively. The plane equation for a 3 £ 3 moving window can be expressed as z^ ¼ a þ bx þ cy: The angle of inclination ðuÞ for the plane is calculated using the least square method as u ¼ tan21 ðc=bÞ (Horn, 1981). (2) Determining flow direction of the center cell. The flow direction of the center cell is determined by the parameters b, c and u, shown in Fig. 1. If b is greater than zero, the elevation at the right neighbors is higher than that at the left neighbors. If c is greater than zero, the elevation at the top neighbors is higher than that at the bottom neighbors. In Fig. 1, the elevation of the colored cells is lower than the other cells, so the flow direction is assigned to the lowest elevation cell. If the colored cells have the same elevation, then the flow direction is determined by u.
Fig. 1. Illustration of determining the z5 flow direction.
51
2.2. Determining the optimal depression outlet using the depression watershed method After calculating the incipient flow direction, pits and/or flat areas coded as zero are necessary for advanced processing. Because drainage moves generally away from higher terrain towards lower terrain, the depression surfaces that include pits and/or flat areas will form a closed watershed, shown in Fig. 2. Four main factors (h1 ; h2 ; l1 and l2 ) can be used to determine the optimal depression watershed outlet. The h1 factor is the difference in elevation from the edge cells to the depression cell. The h2 factor is the difference in elevation from the edge cell to its neighbor cells with the lowest elevation, but not in the depression watershed. The l1 factor is the distance from the edge cells to the depression cell. The l2 factor is the distance from the edge cells to their neighbor cells that have the lowest elevation, but not in the depression watershed. The depression watershed outlet function (O) can be expressed as: O ¼ f ðh1 ; h2 ; l1 ; l2 Þ The priority for determining the optimal outlet from the four factors is h1 ¼ h2 . l1 . l2 : If the h1 and/or h2 of a cell on the edge of the depression watershed is the lowest and its l1 and/or l2 is the shortest, that cell will be the optimal outlet for the depression watershed. Multiple criteria decision-making methods are suitable for analyzing this problem because each cell around the edge of a depression watershed is likely to
Fig. 2. Illustration of depression watershed.
52
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
become the optimal candidate cell for the depression watershed outlet and be evaluated with these four factors. The PROMETHEE method presented by Brans et al. (1984) is an interactive multiple criteria decision-making technique designed to handle qualitative and discrete alternatives. In this study, the PROMETHEE method was chosen to quantify and determine the priority of the candidate cells for the depression watershed outlet. The steps applied to derive the optimal outlet using the PROMETHEE method are described as follows. (1) Establishing an alternatives and criterion matrix If the edge around a depression watershed has n candidate cells, there will be n alternatives ða1 ; a2 ; …; an Þ: Each candidate cell has four evaluative factors (h1 ; h2 ; l1 and l2 ). The evaluative function can be written as f1 ð·Þ; f2 ð·Þ; f3 ð·Þ and f4 ð·Þ: Both the alternatives and criterion can be expressed as a n by 4 ~ matrix T:
Brans et al. (1984) presented the shape of the six possible types of generalized criteria to assist the decision-maker with this selection, as shown in Table 1. He pointed out that the GAUSSIAN criterion was selected most by users for practical applications especially for continuing data. The major criteria in this study are elevation and distance, because of the criteria containing continuity, the GAUSSIAN criterion was chosen for evaluation. In the six type function, HðdÞ ¼ 1 2 expð2d2 =2s2 Þ; s is defined as the threshold value between the indifferent and strict preference areas. In real terrain, if the elevation difference in h1 or h2 is more than 1 m or the distance in l1 or l2 is more than one cell size, the effect on determining the depression outlet will be critical. Therefore, s is defined as 1 and the GAUSSIAN criterion function can be written as HðdÞ ¼ 1 2 expð2d2 =2Þ: The evaluative difference in the ith candidate cell compared with the other candidate ~ i: cells can be expressed as a n using the 4 matrix H 1 2 expð2 d1 ðai ; a1 Þ 1 2 expð2 d2 ðai ; a1 Þ Þ 1 2 expð2 d3 ðai ; a1 Þ Þ 1 2 expð2 d4 ðai ; a1 Þ Þ 2 2 2 2 d ða ; a Þ d ða ; a Þ d ða ; a Þ d ða ; a Þ 1 2 expð2 1 i 2 Þ 1 2 expð2 2 i 2 Þ 1 2 expð2 3 i 2 Þ 1 2 expð2 4 i 2 Þ ~i ¼ H 2 2 2 2 ··· ··· ··· ··· 1 2 expð2 d1 ðai ; an Þ Þ 1 2 expð2 d2 ðai ; an Þ Þ 1 2 expð2 d3 ðai ; an Þ Þ 1 2 expð2 d4 ðai ; an Þ Þ 2 2 2 2 f1 ða1 Þ f2 ða1 Þ f3 ða1 Þ f4 ða1 Þ f ða Þ f ða Þ f ða Þ f ða Þ 1 2 2 2 3 2 4 2 T~ ¼ ··· · · · · · · · · · f1 ðan Þ f2 ðan Þ f3 ðan Þ f4 ðan Þ
(2) Selecting a preference function ðHðdÞÞ The evaluative difference of the ith candidate cell compared with the other candidate cells can be expressed as a n using 4 matrix d~i : d1 ðai ; a1 Þ d2 ðai ; a1 Þ d3 ðai ; a1 Þ d4 ðai ; a1 Þ d ða ; a Þ d ða ; a Þ d ða ; a Þ d ða ; a Þ 1 i 2 2 i 2 3 i 2 4 i 2 d~i ¼ ··· · · · · · · · · · d1 ðai ; an Þ d2 ðai ; an Þ d3 ðai ; an Þ d4 ðai ; an Þ where dðai ; an Þ ¼ f ðai Þ 2 f ðan Þ:
(3) Calculating the preference index (p) For each pair of alternatives, the preference index pðdÞ can be written as pðdÞ ¼
4 X
wj Hj ðdÞ
j¼1
where wj ; j ¼ 1; 2; 3; 4 are weights associated with each criterion. Because the priority for determining the optimal outlet from the four factors is h1 ¼ h2 . l1 . l2 : In this study, the weight of h1 ðw1 Þ is defined as 1, the weights of the other factors ðwj Þ can be written as w1 ¼ w2 ¼ 1 w3 ¼
1 1 or uh1 maxðl1 Þ uh2 maxðl1 Þ
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61 Table 1 The shape of the six possible types of generalized criteria Generalized criterion type Type I: Usual criterion ( 0 d¼0 HðdÞ ¼ 1 ldl . 0
Preference function ðHðdÞÞ
w4 ¼
53
w3 maxðl2 Þ
where uh1 is the variation in the outranking flow while changing h1 ; uh2 is the variation in the outranking flow while changing h2 ; maxðl1 Þ is the maximum of l1 ; maxðl2 Þ is the maximum of l2 : The weights can be ~ expressed as a 1 by 4 matrix w: ~ ¼ lw1 w2 w3 w4 l w
Type II: U-shape criterion ( 0 ldl # q HðdÞ ¼ 1 ldl . q
Type III:8V-shape criterion > < ldl ldl # p p HðdÞ ¼ > : 1 ldl . p
Type IV: Level criterion 8 0 ldl # q > > > <1 HðdÞ ¼ q , ldl # p > 2 > > : 1 ldl . p
Type V: V-shape criterion with indifference criterion 8 > 0 ldl # q > > > < ldl 2 q HðdÞ ¼ q , ldl # p > p2q > > > : 1 ldl . p
Type VI: Gaussian criterion HðdÞ ¼ 1 2 expð2d 2 =2s2 Þ
The preference index of the ith candidate cell compared with the other candidate cells can be written as ~ i ·w· ~ S~a pi ¼ H where S~a is a n by 1 scalar matrix 1 1 ðS~a ¼ Þ: ··· 1 (4) Calculating the outranking flow (f) According to the above preference index (p), the outranking flow for the ith candidate cell can be written as:
fi ¼
n X
pi
i¼1
(5) Generating a complete ranking From the outranking flow calculation, the candidate cell that has the smallest outranking flow will be the optimal outlet for the depression watershed. 2.3. Calculating the depressionless flow directions Two steps can be used to generate a new flow direction for the outlet and depression cells after determining the optimal outlet using the depression watershed method. (1) Changing the flow direction of neighbor cells to the outlet cell that will not flow into the outlet cell and reversing the flow directions of those cells located on the flow route from the outlet cell to the depression. (2) Assigning the depression flow direction to the outlet cell using the recursion algorithm.
54
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
3. Example The Shihmen Reservoir watershed in Northern Taiwan was used in this study (Fig. 3). The focus flat area (Fig. 4a – d) of lower Shihmen Reservoir watershed, contains major landform with reservoir encroachment, farmland and hardwood forest. A 16 by 20 line DEM data sample (Fig. 5) was chosen with average slope 2.97%, area 51.2 ha, no aspect area 56.1%. The geological data from Taiwan’s Central Geological Service is categorized Pleistocene Toen formation. It is composed of plateau gravel and clay with a thin layer of yellowish lateritic earth on top. It is also considered to be thin deposits on river terraces formed after the Tyuerki formation on the adjacent higher terraces. 3.1. Calculating depressionless flow directions The calculated depressless flow directions within sample area were shown in Fig. 6. The procedures for the depression watershed method can be divided into four portions as follows: (1) Calculating the incipient flow directions The non-depression flow directions can be calculated using elevation-differencing and surface-inclining methods. (2) Searching for depression watersheds
Fig. 4. (a) Elevation distribution of the study area. (b) Slope distribution of the study area. (c) Aspect distribution of the study area.
Fig. 3. Study area.
After calculating the incipient flow direction, searching for a depression watershed can be done by tracing the upstream cells that flow into the depression. Fig. 7 shows that there are four depression watersheds A, B, C and D in the DEM. According to the terrain surrounding the depression watersheds, the topography of the A and B depression watersheds are higher. The D depression watershed is lower. (3) Determining the optimal depression watershed outlet
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
55
Fig. 5. Original elevations of the sample DEM.
A breach is located at the bottom portion of the B and D depression watesheds. This breach is defined as the depression watershed outlet. The A and C depression watershed candidate cells are coded as shown in Fig. 7.
There are forty-five alternatives in the A depression watershed and twenty-one alternatives in the C depression. A pairwise comparison could provide a preference structure between these
Fig. 6. Illustration of depression watersheds.
56
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
Fig. 7. The code of the candidate cells in depression watersheds.
alternatives. The preference index with GUSSIAN function and criterion weights could then be calculated. According to the weight definitions, for h1 and h2 ; the weights can be written as: w1 ¼ w2 ¼ 1 In the l1 1 1 T~ ¼ 1 1
weight, uh1 can be obtained by substituting 2 3 1 1 1 1 1 1
into the PROMETHEE algorithm, and then maxðl1 Þ and maxðl2 Þ in the A and C depression watersheds can be derived as follows:
uh1 ¼ 1:2581 For A depression watershed: maxðl1 Þ ¼ 7:41 maxðl2 Þ ¼ 1:41 w3 ¼
1 ¼ 0:107 1:2581 £ 7:41
w4 ¼
0:107 ¼ 0:076 7:41
For C depression watershed: maxðl1 Þ ¼ 4:24 maxðl2 Þ ¼ 1:41 w3 ¼
1 ¼ 0:187 1:2581 £ 4:24
w4 ¼
0:187 ¼ 0:133 1:41
The criteria value, outranking flow and complete ranking for each candidate cell are shown in Tables 2a and 2b. In the A depression watershed, the candidate cell coded A7, with the smallest outranking flow, 2 80.719, is the optimal outlet cell. In the C depression watershed, the candidate cell coded C8, with the smallest outranking flow, 2 32.792, is the optimal outlet cell. According to the outranking flow analysis, the outlet for the A to D depression watersheds can be illustrated as shown in Fig. 8. (4) Flow direction assignment in a depression After determining the depression watershed outlet, the first procedure is assigning the flow direction from the outlet cell into a neighbor that does not flow into the outlet cell and has a lower elevation. The flow direction of cells that route from the outlet cell into the depression is then reversed. The second procedure is to assign
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
57
Table 2a The outranking flow and complete ranking of each candidate cells in A-depression watershed Alternatives
A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 A13 A14 A15 A16 A17 A18 A19 A20 A21 A22 A23 A24 A25 A26 A27 A28 A29 A30 A31 A32 A33 A34 A35 A36 A37 A38 A39 A40 A41 A42 A43 A44 A45
Criterion h1 ; w1 ¼ 1
h2 ; w2 ¼ 1
l1 ; w3 ¼ 0:107
l2 ; w4 ¼ 0:076
27 23 23 5 5 2 2 7 22 25 25 25 29 37 52 43 56 52 68 43 29 25 25 25 23 24 15 20 21 21 23 16 16 2 2 2 3 9 15 14 17 19 11 8 25
17 17 17 3 3 2 2 2 13 13 25 26 27 29 37 54 48 57 51 26 26 26 22 19 24 17 23 23 23 21 21 22 9 9 2 2 2 4 12 17 20 20 8 8 8
5.83 3.41 3.00 2.00 1.00 1.41 1.00 2.00 2.00 3.41 3.41 4.41 5.41 6.41 7.41 5.83 6.83 6.41 6.83 5.41 4.41 3.41 3.00 3.41 2.00 2.41 2.00 2.00 2.41 3.41 2.41 2.00 3.00 2.00 3.83 3.41 4.41 5.41 6.41 4.83 4.41 4.00 4.83 4.41 5.41
1.41 1.00 1.41 1.00 1.41 1.00 1.41 1.00 1.41 1.00 1.41 1.41 1.41 1.41 1.41 1.41 1.41 1.41 1.41 1.41 1.00 1.41 1.41 1.41 1.41 1.41 1.41 1.00 1.00 1.00 1.41 1.41 1.00 1.41 1.00 1.00 1.41 1.41 1.41 1.41 1.41 1.41 1.41 1.00 1.41
the flow direction of the depression cells into the outlet cell. Fig. 9 shows the results of the analyzed depressionless flow direction for the A to D depressions.
c
Rank
22.498 20.969 21.348 267.973 268.772 280.573 280.719 268.077 214.947 2.578 39.986 44.917 60.961 69.419 79.279 80.030 83.924 86.588 87.914 64.379 56.528 43.301 27.278 17.289 22.090 3.722 23.599 7.159 10.315 3.907 10.837 24.879 235.984 262.800 277.138 277.815 273.314 252.094 227.469 222.627 26.393 22.970 243.190 247.974 23.274
32 23 22 8 6 2 1 7 16 24 34 36 38 40 41 42 43 44 45 39 37 35 33 30 31 25 19 27 28 26 29 18 13 9 4 3 5 10 14 15 17 21 12 11 20
3.2. Watershed delineation application To verify the suitability of the proposed method, the depressionless flow direction was calculated using
58
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
Table 2b The outranking flow and complete ranking of each candidate cells in C-depression watershed Alternatives
C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 C13 C14 C15 C16 C17 C18 C19 C20 C21
Criterion h1 ; w1 ¼ 1
h2 ; w2 ¼ 1
l1 ; w3 ¼ 0:187
l2 ; w4 ¼ 0:133
22 24 13 6 2 2 1 1 5 9 15 20 20 15 2 1 2 3 13 21 23
21 15 20 8 4 4 1 1 2 2 4 9 23 9 3 3 3 1 1 5 16
4.24 3.83 2.83 2.41 2.00 1.00 1.41 1.00 1.00 2.00 2.41 3.41 3.00 3.41 2.41 1.41 1.00 2.00 2.41 2.41 3.83
1.41 1.41 1.41 1.41 1.00 1.41 1.00 1.41 1.00 1.41 1.41 1.41 1.41 1.41 1.41 1.00 1.41 1.00 1.41 1.41 1.41
Fig. 8. Illustration of the outlet of depression watersheds.
c
Rank
36.191 34.189 20.704 4.728 214.462 215.736 232.468 232.792 216.995 211.707 5.336 21.910 32.709 16.886 217.866 222.920 220.008 224.765 211.942 15.437 33.573
21 20 16 12 9 8 2 1 7 11 13 17 18 15 6 4 5 3 10 14 19
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
59
Fig. 9. Illustration of depressionless flow direction in a noised DEM.
for the Shihmen Reservoir using the depression watershed method, elevation-smoothing method (O’Callaghan and Mark, 1984), depression-filling method (Jenson and Domingue, 1988) and elevation-incrementing method (Martz and Garbrecht, 1998). Numerous large flat terrain areas that have more than one outflow point exist near the downstream area of the Shihmen Reservoir watershed. In these methods, the elevation-smoothing method removes only the shallow depressions. The deeper depressions remain. Therefore, these methods fail to fully accommodate flat areas. The elevationsmoothing curve for the Shihmen Reservoir watershed is shown in Fig. 10. The delineated watershed boundary using the depression watershed, elevationfilling and elevation-increasing methods are shown in Fig. 11. Fig. 11 shows that there is an obvious difference near the downstream area of the watershed. The watershed areas delineated using the depression watershed, depression-filling and elevation-incrementing methods were 75,633.60, 75,290.40 and 75,304.64 ha. To verify the suitability of the depression-processing methods, an electric map of the reservoir flood area, at a scale of 1/25,000, was
overlaid onto the watershed boundary. The overlay (Fig. 11) indicated that the result delineated using the depression watershed method was more realistic than the other methods.
Fig. 10. The convergent curve of applying elevation-smoothing method on raising depressions.
60
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61
Fig. 11. Comparison of applying varied depression-processing methods on delineating the Shihmen Reservoir watershed.
4. Conclusions The proposed method in this study searches the best outlet for a depression without changing the grid elevation. It avoids the difficulty for determining outlet during increasing the depression elevation and obtaining an unrealistic watershed in a large flat area or several nearby depressions. Especially, when flat areas or depressions are located close to watershed boundary, the traditional methods may establish an incomplete watershed area. The new outlet algorithm is more complicated than usual. However, it improves the result obtained from automated watershed boundary delineation. The calculated flow directions in this study produce convolution flow patterns within each depression. A more rational flow pattern arrangement inside the depression need further compare and consult with the Jenson and Domingue (1988) or the Garbercht and Martz (1997) approach. The Shihmen Reservoir watershed in Northern Taiwan was delineated using several popular depression-processing techniques. The results indicated that the depression watershed method was more
realistic than the other methods. The algorithm proposed by this study produced quantitative analysis for depression outlet determination and can derive useful morphological information for hydrological applications.
References Brans, J.P., Mareschal, B., Vincke, P., 1984. PROMETHEE: A New Family of Outranking Methods in MCDM, IFORS 84, North Holland, pp. 477–490. Chang, K.T., 2002. Introduction to Geographic Information Systems, McGraw-Hill, New York, NY, p. 233. Collins, S.H., Moon, G.C., 1981. Algorithms for dense digital terrain models. Photogrammetric Engineering and Remote Sensing 47, 71–76. Doyle, F.J., 1978. Digital terrain models: an overview. Photogrammetric Engineering and Remote Sensing 44, 1481– 1485. Garbercht, J., Martz, L.W., 1997. Assignment of drainage direction over flat surfaces in raster digital elevation models. Journal of Hydrology 193, 204–213. Horn, B.K.P., 1981. Hill shading and the reflectance map. Proceedings of the IEEE 69, 14–47. Jenson, S.K., Trautwein, C.M., 1987. Methods and applications in surface depression analysis. Proceedings of Auto-Carto 8, 137 –144. Jenson, S.K., Domingue, J.O., 1988. Extracting topographic structure from digital elevation data for geographic information
T.-Y. Chou et al. / Journal of Hydrology 287 (2004) 49–61 system analysis. Photogrammetric Engineering and Remote Sensing 54, 1593–1600. Marks, D., Dozier, J., Frew, J., 1984. Automated basin delineation from digital elevation sata. Geo-Processing 2, 299–311. Martz, L.W., Garbrecht, J., 1998. The treatment of flat areas and depressions in automated drainage analysis of raster digital elevation models. Hydrological Processes 12, 843–855.
61
O’Callaghan, J.F., Mark, D.M., 1984. The extraction of drainage networks from digital elevation data. Computer Graphics and Image Processing 28, 323–344. 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, 263– 293. US Geological Survey, 1987. Digital elevation models. US Geological Survey Data User’s Guide 5, p. 38.