Optimal Pit Removal
in Watershed and Stream Network Delineation
Annjanette Dodd
Background
Grid digital elevation models (DEM's) are data structures in which the topographic elevation of each pixel is stored in a matrix node. Grid DEM's are readily available from the Internet (1:24,000 scale from http://edcnts12.cr.usgs.gov/ned/ or 1:250,000 scale from http://edc.usgs.gov/nsdi/gendem.htm). Two primary hydrologic quantities derived from DEM's are slopes and flow directions. TARDEM is a suite of programs that can be utilized for the analysis of digital elevation data to determine hydrologic quantities of watersheds (for details go to http://www.engineering.usu.edu/dtarb/). The steps involved are: (1) Pit filling corrections, (2) Computation of slopes and flow directions; (3) Computation of contributing area and specific catchment area and (4) Channel network extraction and computation of other quantities.
Here I will focus on the first and second steps only. The 'flood' routine in TARDEM involves a process of filling pits in the elevation data. Unnatural pits found in digital elevation data can affect the resulting flow directions determined from grid DEM's. Pits, rare in natural topography, are grid elements or sets of grid elements surrounded by higher terrain that, in terms of the DEM, do not drain (Figure 1). These arise from errors that result from the discrete nature of DEM's.
Figure 1: Example of an unnatural pit found in grid DEM's.
The 'flood' routine in TARDEM is an attempt to fill these unnatural pits by raising the elevation of each pit grid cell within the DEM to the elevation of the lowest pour point on the perimeter of the pit (Figure 2). This is not an optimal method for determining the elevation of each pit grid cell.
Figure 2: Example of pit removal using 'Flood' from TARDEM.
The 'D8' routine in TARDEM takes pit filled elevation data and creates flow directions using the eightdirection pour point model called the “D8 method”. The direction of flow is from a pixel, determined by the steepest downward slope, to one of its eight neighbors, either adjacent or diagonally. The slope is computed as the difference in elevation between two adjacent cells divided by the distance between them. The direction is identified as one of eight flow directions starting due east and moving counter clockwise (Figure 3). (In flat areas flow directions are assigned away from higher ground and towards lower ground using the method of Garbrecht and Martz).
Figure
3: Eight flow directions
utilized in the D8 method.
To improve the accuracy of flow direction along existing channel networks, the TARDEM includes the routine 'Streamtogrid' which formats EPA river reach files into grid DEM format and outputs a grid of flow directions along the rivers to enforce drainage along delineated river channels.
Statement of the
Problem
In a large DEM grid, multiple pits can exist (Figure 4). To better delineate flow paths, optimal pit removal (i.e., minimization of deviations from original elevations) should be performed.
Figure 4: Example of multiple pits in a large DEM. (Lecture presented by David Tarboton)
The purpose of this project is to create a program that minimizes the overall changes in elevation subject to grid drainage (determined by direction of steepest descent). This will be accomplished by performing the following on a grid DEM with original elevations (e_{i}):
1) Fill the pits in the original elevation data using the program 'flood'
2) Calculate the flow directions using the D8 method
The result is a grid of flow directions utilizing the pit filled elevation data (fd)
(This is a first estimate of flow directions in the channel network)
3) Set up objective function and constraint set of the system and optimize
The constraint set utilizes the flow directions (fd) and original elevations (e_{i})
The result is a new elevation grid with the set of elevations (e_{n})
4) Calculate the flow directions using the D8 method on the new elevations
5) Repeat steps 3 and 4 until a "no change" criterion between old and new elevations is met
The TARDEM software can be utilized to perform all of the above steps except the optimization. The rest of this paper is devoted to this optimization problem.
The Optimization
For a DEM grid with n rows and m columns, the problem can be written as a linearly constrained quadratic programming problem
(1)
where x_{i} is the deviation from the original elevation at node i and ci represents the downstream node that node i drains into. (There will actually be less than mn decision variables and mn1 constraints when empty cells of the DEM exist.) The new elevation for each node is determined by
. (2)
The constraint set, g, represents the imposed flow conditions determined using the D8 method on the old elevations. The last node mn will not have a constraint because it will not have a pour point within the grid being analyzed. The constraint set given in (1) can be written as
(3)
where b_{i} = e_{ci}  e_{i}. It is important to note that each constraint function is dependent only on the node it drains into. Written in matrix form, with a coefficient matrix A and a righthand side as a column vector representing the elevation differences of upstream and downstream pour nodes, the constraint set will have the form
.
(4)
Here A will be a (mn1) x (mn) sparse matrix with a coefficient of –1 in column i, a coefficient of 1 in the column represented by the pour node ci, and zeros everywhere else. The location of the pour node column (5) is only one of node i’s eight neighbors as determined from the D8 method.
Figure
5: Pour node identification
The coefficient matrix will be extremely large for large DEM's and will reduce computational efficiency. Compressing the system is necessary for efficient optimization of large DEM's. Separability of the problem in (1) can be exploited to reduce the size of the entire system.
Example 1: Setting up the constraint set.
The 4 x 3 DEM grid (Figure 5a) will have the flow direction grid (Figure 5b) that results from the flooded grid (Figure 5c) (where the superscripts represent the node number i). For the grid given in Figure 5a, the constraint set is summarized in Table 1 and an example of the coefficient matrix is given in Figure 6.
Figure 5: Example DEM grid pit filled using 'Flood' with flow directions determined using D8.
Table 1: Constraint set for Example 1
i 
x_{ci} 
x_{i} 
b_{i} 
1 
x_{5}x_{1} 
8 
2 
x_{5}x_{2} 
7 
3 
x_{6}x_{3} 
3 
4 
x_{5}x_{4} 
5 
5 
x_{9}x_{5} 
2 
6 
x_{9}x_{6} 
1 
7 
x_{11}x_{7} 
3 
8 
x_{12}x_{8} 
3 
9 
x_{12}x_{9} 
1 
10 
x_{11}x_{10} 
1 
11 
x_{12}x_{11} 
1 

1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
1 
1 
0 
0 
0 
1 
0 
0 
0 
0 
0 
0 
0 
2 
0 
1 
0 
0 
1 
0 
0 
0 
0 
0 
0 
0 
3 
0 
0 
1 
0 
0 
1 
0 
0 
0 
0 
0 
0 
4 
0 
0 
0 
1 
1 
0 
0 
0 
0 
0 
0 
0 
5 
0 
0 
0 
0 
1 
0 
0 
0 
1 
0 
0 
0 
6 
0 
0 
0 
0 
0 
1 
0 
0 
1 
0 
0 
0 
7 
0 
0 
0 
0 
0 
0 
1 
0 
0 
0 
1 
0 
8 
0 
0 
0 
0 
0 
0 
0 
1 
0 
0 
0 
1 
9 
0 
0 
0 
0 
0 
0 
0 
0 
1 
0 
0 
1 
10 
0 
0 
0 
0 
0 
0 
0 
0 
0 
1 
1 
0 
11 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
1 
1 
Figure 6: Constraint coefficient matrix, A, for Example 1.
Feasibility and Solution
The programming problem given in (1) has a constraint set that is closed, bounded, and nonempty. The objective function is realvalued and continuous. The Weierstrass and Localglobal theorems guarantee that a solution exists for (1). Since the constraints are linear, the feasible region is a convex set. The Hessian of the objective function (matrix of second partial derivatives) is 2*I (I is the identity matrix) which is positive semidefinite, thus the objective function is convex. The intersection of two convex sets is a convex set. Therefore, any local solution is a global solution to the programming problem (1). (Willis and Finney, 1998)
The development of the optimal solution begins by introducing a Lagrange multiplier, y_{i}, for each constraint. There are mn1 constraints, thus there will be mn1 Lagrange multipliers. The Lagrangian for (1) is
.
(5)
Program (1) is a separable problem because each function can be written as a sum of functions, each depending on only one variable. Due to separability of (1), the optimization problem decomposes into mn separate singlevariable minimization problems or
.
(6)
Equation (6) is called the dual function and is maximized with respect to y_{i}³0 for all i. The last summation is over the Lagrange multipliers for which x_{cj} = x_{i}, i.e., for which x_{i} is a downstream pour node.
The first order necessary conditions for a local maximum are
(7)
Using the first necessary condition in (7), the optimal decision variables are given by
(8)
The second necessary condition in (7) can be used to estimate y_{i}.
If a node is not a downstream node for which any other nodes flow into, then (8) reduces to
.
(9)
Thus, if x_{i}<0 then y_{i} must equal zero to satisfy the nonnegativity constraint for the Lagrange multipliers. If a node is a downstream flow node and all y_{j} for this node are zero (i.e., no other nodes pour into it) and x_{i}<0, then y_{i}=0. For these nodes, the second necessary condition is automatically satisfied. The remaining n_{sub} decision variables that are not predetermined to be zero are used to determine the optimum for the system that satisfies (7). This results in a system that is much reduced in size.
The original optimization problem in (1) is reduced to a subproblem that involves optimizing a system of n_{sub} decision variables. For the n_{sub} nodes that are left, the second necessary condition can be used to estimate y_{i} or
(10)
where, from (8),
(11)
Using an initial guess of y_{0}=0, results, from (8), in x_{0}=0. From (10), with x_{0}=0, the next value for y_{i} is b_{i}. Thus, for the next set of y_{i}³0, b_{i} must be nonnegative. If b_{i}£0, then, as discussed above, x_{i}£0 and the next values for both x_{i} and y_{i}_{ }are zero. Otherwise, y_{i}=b_{i} and the next x_{i} can be evaluated from (8). For each case that x_{i}=0, the updated value will remain zero from step to step. Thus, the problem becomes one of solving the original quadratic problem given in (1) with n_{sub} decision variables and the corresponding constraints. For the remaining y_{i}’s to be a solution to original quadratic programming problem, then the following equality must hold (Belegundu & Chandrupatla)
.
(12)
where the corresponding x_{i}’s are given by (8). Therefore, utilizing separability, Lagrange multipliers, and principals of a quadratic programming problem, the original optimization problem reduces to the solution of a set of n_{sub} linear equations that are easily solved.
Example 2: The system given in Example 1 is optimized using the method discussed above. For the DEM given in Example 1 all of the x_{i}'s and y_{i}'s are predetermined to be zero except for x_{5}, x_{9}, and x_{12} (Table 2). Thus, the resulting system to be solved is
Table 2: Evaluation of Lagrange Multipliers in Example 2
i 
x_{i} 
b_{i} 
1 
y_{1}/2 ̃0 
8 
2 
y_{2}/2 ̃0 
7 
3 
y_{3}/2 ̃0 
3 
4 
y_{4}/2 ̃0 
5 
5 
(y_{5}y_{1}y_{2}y_{4})/2 ̃ y_{5}/2 
2 
6 
(y_{6}y_{3})/2 ̃ y_{6}/2 ̃0 
1 
7 
y_{7}/2 ̃0 
3 
8 
y_{8}/2 ̃0 
3 
9 
(y_{9}y_{5}y_{6})/2 ̃ (y_{9}y_{5})/2 
1 
10 
y_{10}/2 ̃0 
1 
11 
(y_{11}y_{7}y_{10})/2 ̃ y_{11}/2̃0 
1 
12 
(y_{8}y_{9}y_{11})/2 ̃ y_{9}/2 

Thus, the Lagrange multipliers can be determined from substitution into the two remaining constraint equations given in (12), or
_{}̃_{}̃_{}̃_{}̃_{}
The objective function evaluated for this solution is z=2. The objective value from 'flood' was z=4. It is important to note that if the objective function was set up to minimize the sum of the absolute values of the perturbations, both objective values would be the same. Thus, for the linear programming problem, as given in this example, the solution would not be unique.
Figure 7:
Optimized DEM from Example 2.
Figure 8: Profile of original
pit and filled pit using optimization and 'Flood'
Discussion
The purpose of presenting the theory discussed here was an attempt to reduce the size of the original optimization problem for computational purposes. Something is missing in the assumptions and derivations made above and further development of this theory must occur before a full scale computer program is developed to optimize the pit filling routine in large DEM's.
Example 3: Testing the above theory on a slightly larger grid with multiple pits.
30 
29 
28 
29 

30 
29 
28 
29 

6 
6 
6 
7 

8 
7 
6 
3 
27 
22 
26 
28 

27 
24 
26 
28 

6 
11 
11 
12 

5 
2 
2 
2 
27 
26 
24 
26 

27 
26 
24 
26 

14 
14 
15 
16 

6 
5 
3 
3 
25 
21 
21 
23 

25 
23 
23 
23 

14 
15 
19 
20 

4 
0 
2 
0 
24 
24 
23 
23 

24 
24 
23 
23 

22 
23 
24 
24 

2 
3 
3 
3 
23 
22 
21 
20 

23 
22 
21 
20 

22 
23 
24 


1 
1 
1 

a b c d
Figure 9: a) Original Grid, b) Pit Filled Grid, c) Node number that cell pours into (ci), and d) grid representing b_{i}'s.
For this 4 x 6 DEM x_{i}'s and y_{i}'s are predetermined to be zero except for x_{6}, x_{11}, x_{15}, x_{19}, and x_{24 }(Table 3). Thus, the resulting system to be solved is
Table 3: Evaluation of Lagrange Multipliers in Example 3.
i 
x_{i} 
b_{i} 
15,810,13, 17,18,21 
y_{i}/2 ̃0 
negative 
6 
(y_{6}y_{1}y_{2}y_{3}y_{5})/2 ̃ y_{6}/2 
2 
7 
(y_{7}y_{4})/2 ̃ y_{7}/2 ̃0 
2 
11 
(y_{11}y_{6}y_{7})/2 ̃ (y_{11}y_{6})/2 
3 
12 
(y_{12}y_{8})/2 ̃ y_{12}/2 ̃0 
3 
14 
(y_{14}y_{9}y_{10}y_{13})/2 ̃ y_{14}/2 ̃0 
0 
15 
(y_{15}y_{11}y_{14})/2 ̃ (y_{15}y_{11})/2 
2 
16 
(y_{16}y_{12})/2 ̃ y_{16}/2 ̃0 
0 
19 
(y_{19}y_{15})/2 ̃ (y_{19}y_{15})/2 
3 
20 
(y_{20}y_{16})/2 ̃ y_{20}/2 ̃0 
3 
22 
(y_{22}y_{17}y_{21})/2 ̃ y_{22}/2 ̃0 
1 
23 
(y_{23}y_{18}y_{22})/2 ̃ y_{23}/2 ̃0 
1 
24 
(y_{19}y_{20}y_{23})/2 ̃ y_{19}/2 
1 
Thus, the Lagrange multipliers and optimal perturbations can be determined from substitution into the two remaining constraint equations given in (12), and are
i 
y_{i} 
x_{i} 
6 
0 
0 
11 
4 
2 
15 
2 
1 
19 
4 
1 
24 

2 
The optimized grid that results is
30 
29 
28 
29 
27 
22 
26 
28 
27 
26 
26 
26 
25 
21 
20 
23 
24 
24 
24 
23 
23 
22 
21 
18 
Figure 10: Optimized grid from Example 3
Obviously, the solution that is obtained in Example 3 does not satisfy the original flow directions as given in the constraint set and the pits are not removed. In this example, not only are they not removed, the pits are not improved or are made worse (Figure 10). This can only be a direct result of the assumptions made in eliminating Lagrange multipliers. Based on the theoretical analysis discussed above, it seems that the nodes that should have been perturbed are nodes 6, 11,14,15, 19 and 24. If node 14 is added into the constraint set, the perturbation at this node remains zero (because b_{14} = 0) and the system that results is inconsistent. This result indicates that a special constraint or condition is necessary when b_{i}=0.
The results of Example 3 make it evident that further analysis of the theory and system composition is required to fully understand the interaction between the Lagrange multipliers and the constraints.
References
Belegundu, A. and Chandrupatla, T. (1999). Optimization Concepts and Applications in Engineering. Prentice Hall, New Jersey.
Tarboton, D. http://www.engineering.usu.edu/dtarb/
Willis, R. and Finney, B. (1998). Environmental Systems Engineering. Dept. of Environmental Resources Engineering, Humboldt State University, California.