An Automated Processing Algorithm for Flat Areas Resulting from DEM Filling and Interpolation

Abstract

Correction of digital elevation models (DEMs) for flat areas is a critical process for hydrological analyses and modeling, such as the determination of flow directions and accumulations, and the delineation of drainage networks and sub-basins. In this study, a new algorithm is proposed for flat correction/removal. It uses the puddle delineation (PD) program to identify depressions (including their centers and overflow/spilling thresholds), compute topographic characteristics, and further fill the depressions. Three different levels of elevation increments are used for flat correction. The first and second level of increments create flows toward the thresholds and centers of the filled depressions or flats, while the third level of small random increments is introduced to cope with multiple threshold conditions. A set of artificial surfaces and two real-world landscapes were selected to test the new algorithm. The results showed that the proposed method was not limited by the shapes, the number of thresholds, and the surrounding topographic conditions of flat areas. Compared with the traditional methods, the new algorithm simplified the flat correction procedure and reduced the final elevation increments by 5.71–33.33%. This can be used to effectively remove/correct topographic flats and create flat-free DEMs.

1. Introduction

Digital elevation models (DEMs) have been widely used for hydrological analyses and modeling. Various hydrotopographic characteristics (e.g., depressions, slopes, flow directions and accumulations, stream channels, and basin boundaries) can be extracted from DEMs [1,2,3,4,5,6]. Due to the limitations in vertical resolution and accuracy, DEMs are usually reconstructed by using different interpolation methods before hydrological processing. However, because of various flaws and inherent shortcomings in the original elevation data and interpolation methods, numerous small depressions and flats (i.e., artifacts) are generated after interpolation in many cases [6,7,8]. Thus, flat areas are generally present in DEMs [9,10]. In addition, to eliminate the stagnation of water flowing from depressions to stream networks, surface depressions need to be fully filled before calculating flow directions and accumulations. However, filling depressions (i.e., increasing the elevations of all depressional cells) significantly increases the number of flat areas [9,10,11]. Due to the same elevation within a flat, it is difficult to determine flow directions [12]. Even though the D8 method is the most popular one for calculating flow directions and accumulations, it is incapable of dealing with flat areas [13]. The information on flow directions is essential to hydrological analyses. Except for slope, all other hydrologic parameters can be computed only after flow directions are determined [14]. Additionally, too many “insignificant” small flat areas will affect the efficiency of surface delineation [15] and even cause crashes of the delineation programs. Thus, development of efficient and effective methods for flat correction/removal is of theoretical and practical importance.

To address this issue, different methods have been proposed and used. Jenson and Domingue developed a method (JD method) with neighborhood techniques to assign flow directions for DEM cells [11], but the method created unrealistic parallel drainage lines [12]. Tribe [12] introduced a method involving a straight line between a valley entrance and outlet (s-o module), while it altered the elevations of cells outside of flat areas if they had irregular shapes. Jana et al. [16] used variable increments to correct the elevations of flat areas. Zhang et al. [17] overlaid a digital channel network on flat areas, and then used an interpolation method to make water flow toward channels. However, the digital channel network may not always be available. Pan et al. [14] developed a method that implemented linear interpolation between low-elevation outlet and inflow cells. The method required an iterative solution for flat areas of irregular shapes. Zhang et al. [18] proposed a water depth gradient burning method to extract watershed characteristics in flat regions. Su et al. [19,20] proposed the chain code matrix (CCM) and the distance transform (DT) to efficiently assign flow directions over a flat surface. Zhang et al. [21] proposed the maze-weight (MW) method, which was combined with non-flat flow accumulation values, to determine flow directions. Among the state-of-the-art methods, TOPAZ (topographic parameterization), an automated digital landscape analysis tool, has been extensively used in hydrological studies [22,23,24,25,26,27,28]. Although Barnes et al. [29] made some improvement with FlatMask by applying different weight gradients in the delineation processes, their approach tended to create more new small flat areas. Such small flats may seriously influence the efficiency of their surface delineation method [2,15,30], in which the search process starts from sinks and flat areas. Additionally, for multiple overflow thresholds of depressions, TOPAZ may create multiple flow direction cells under certain circumstances.

In this study, a new algorithm is proposed to correct/remove flat areas by adding small elevation increments for all flat cells. Three different levels of elevation increments are used to create artificial slopes and make water flow toward overflow thresholds. After the flat processing, no multi-direction cells exist for both single-threshold and multiple-threshold depressions. The new algorithm is tested by using a set of artificial surfaces and two real-world surfaces.

2. Methodology

2.1. Flat Area Correcting Procedure

Any depression has at least one overflow/spilling threshold. When a depression is fully filled, a flat is formed, and the water in the depression overflows through its threshold(s). In the new flat-correction algorithm, small elevation increments are added to all cells in flat areas to create gradients toward their thresholds. The flat-correcting procedure includes three main steps and one optional step. In step 1, gradients toward the threshold of a flat are created by adding different increment numbers. In step 2, additional gradients, which are away from the flat boundary cells, are created by adding different increment numbers for the flat cells. In step 3, final elevation increments, which are the product of the increment numbers (summation of the increment numbers added in steps 1 and 2) and an increment unit, are added to the original elevations of the flat cells. Step 4 is an optional one, which is used to correct the elevations of new small flats generated in the process and handle depressions with multiple thresholds.

To better describe the new algorithm, a square surface with 7 × 7 cells was selected as an example (Figure 1a). After implementing the depression filling process, a flat area (ranging from B3 to F6) with one threshold at D2 (Figure 1b) was generated on the surface. Except for cell D1 (downstream cell of D2 or outlet), the elevations of all surrounding cells are higher than the elevation of the flat area (Figure 1b). The flat correction steps for this surface are detailed as follows:

2.1.1. Step 1: Making Flow toward Thresholds

In step 1 of the flat-correcting procedure, gradients are created to make water flow toward the threshold of any depression by adding different elevation increment numbers. In this study, the increment numbers represent how many increment units are added to the flat area, as shown in Figure 2a. In the beginning of this step, the algorithm searches the threshold (D2) and marks the cells (C3–E3) that are next to the threshold in the flat area (Figure 2a). Then, a DO loop is performed to add increment numbers for the flat cells. In the first loop, if a cell is next to the threshold, no increment number is added; otherwise, an increment number of 2 is added to the flat cells. As shown in Figure 2b, the increment numbers of cells B3, F3, and B4–F6 become 2 after the first loop. Meanwhile, the number of cells with the added increment numbers in this loop is counted as N, which is used to determine whether or not the loop continues. The cells, which are next to the cells without the added increment numbers (e.g., C3–E3) in this loop, are treated as unadded increment number cells (e.g., B3, F3, and B4–F4) for the following loops. In the second loop, an increment number of 2 is added to the cells that have not been treated as unadded cells. As a result, the increment numbers of cells B5–F6 become 4 in the second loop (Figure 2c). Such a loop continues until N = 0. The final increment numbers of all cells from the first step are shown in Figure 2d and their flow directions are displayed in Figure 2e.

2.1.2. Step 2: Making Flow toward Centers

Unfortunately, step 1 may create many unrealistic “parallel flows” in the flat area, as mentioned in [11,12]. To solve this problem, increment gradients away from the boundaries of the flat will be created in the second step. The main theory is similar to the one proposed by Garbrecht and Martz [31], but some implementation details are different. In this study, the flat cells next to the surrounding outside cells are defined as boundary cells of the flat area. A DO loop is also performed in step 2 of the algorithm to add increment numbers for some flat cells. In the first loop, an increment number of 1 is added to all boundary cells. As a result, the increment numbers of cells B3–B6, C3–E3, C6–E6, and F3–F6 become 1 after the first loop (Figure 2f). These cells are marked as added cells for the following loops. In the second loop, an increment number of 1 is added to the marked cells and the cells adjacent to the marked cells (i.e., cells C4–E5) (Figure 2g). The loop continues until the increment numbers of all cells in the flat area (except the threshold cell) have been added at least once. This step ensures that water flows from all boundary cells toward the inside of the flat area, as shown in Figure 2h.

2.1.3. Step 3: Creating Artificial Slopes

In this step, the calculated increment numbers in steps 1 and 2 are added together (Figure 2i). Figure 2j shows the final flow directions from step 3. Since different levels of increment numbers (2 and 1) are used in the first two steps, the gradients, which are created by different levels of increment numbers, will not cancel each other when combining them. Thus, the current algorithm avoids the problem in the digital elevation drainage network model (DEDNM) [31] that the increment numbers in a cell may cancel each other out. In the new algorithm, a small increment of 1.0 × 10⁻⁵ m is selected as the increment unit, which is smaller than the increment unit (2.0 × 10⁻⁵ m) used in DEDNM. Compared with the vertical resolution of the original DEM, this increment unit is very small. Hence, the final increments, which are calculated by the product of the increment unit and the increment number, will not have any significant influences on the actual elevations of all DEM cells. The final elevation of each cell in the new modified DEM equals the summation of its original elevation and the final increment.

2.1.4. Step 4: Correcting New Flats and Multiple Flow Direction Cells

Although the results from step 3 are sufficient to determine flow directions for most cases, the new surface may still have some small flats such as cells C4–E4, and B6–F6 in Figure 2i. The new flat areas can cause serious problems in some algorithms which are very sensitive to flat areas [2,15]. Additionally, some cells in flat areas of the new surface with multiple overflow thresholds may have multiple flow directions. To deal with these special issues, an optional step is added in this new algorithm. To clearly describe the methodology, a surface with a four-threshold depression/flat was created (Figure 3).

The procedures for dealing with flat areas with multiple thresholds are similar to those for normal flat areas with a single threshold. In step 1, the algorithm searches simultaneously from different thresholds, and then creates gradients toward each threshold (Figure 3a). In step 2, the same procedures are implemented (Figure 3b). However, the situations become complex in step 3 for a case of multiple thresholds. As shown in Figure 3c, flow directions can be easily identified for most cells, but some cells (e.g., D4, D5, and E5) have multiple flow directions. To create a new surface without flat areas and multi-direction cells, random increment numbers ranging from 0.001 to 0.290 are assigned to the new flats, as justified below.

In Figure 3c, the smallest difference between the increment numbers is 1 resulted from step 3 of the new method. For example, cell C3 flows toward cell B3. The slopes from C3 to B3 and B2 are 1/L and 0.7071/L, respectively (L denotes the cell size). In order to make water flow from C3 to B3, the slope (1/L) between C3 and B3 must be greater than 0.7071/L. This also implies that the highest variation for cell B3 should be smaller than 0.2929 (i.e., 1 − 0.7071). To avoid creating the same increment number between neighbor cells, random increment numbers ranging from 0.001 to 0.290 are added to the increment numbers calculated in step 3. The final increment numbers are shown in Figure 3d. Thus, the final elevation increments, which are the product of the new increment number and the increment unit for each cell, are added to the original elevations for all DEM cells. The flow directions are shown in Figure 3e.

2.2. Introduction to DEDNM

DEDNM is the core component of TOPAZ, a DEM-based topographic analysis tool [31]. In the DEDNM module, two gradient increments are applied to create artificial slopes: one gradient toward the lower terrain, the other gradient away from the higher terrain [31]. The main difference between the new method and DEDNM is that DEDNM uses the same level of increment numbers (i.e., 1) in the first and second steps and 2.0 × 10⁻⁵ as an increment unit in the third step. More details can be found in [31]. However, the DEDNM module may sometimes suffer from an exceptional situation, in which the combined gradients will create new flat areas with no drainage direction. In this case, more loops are needed to correct the new flats by adding a half increment number (i.e., 0.5) to the combined increment numbers in step 3. As pointed out by Garbrecht and Martz [31], adding this half increment number may even encompass more adjacent cells and, hence, more loops are involved.

2.3. Testing of the New Flat Correction Algorithm

To test the new flat correction algorithm, a set of artificial surfaces with different flat shapes and different thresholds (both number and distribution) (Figure 4) were created in this study. Flat (a) represents a normal single-threshold topographic condition. Flats (b) and (c) consist of linearly-distributed flat cells. In reality, depressions and flats may have different shapes and relationships with their surrounding areas, as well as dissimilar connectivity via single or multiple thresholds. Flat (d) has a special spiral shape, while flat (e) is a combination of flats (d) and (a) with a large flat area. Note that many existing methods may not be able to efficiently handle flat (d) due to its complex spiral shape. Flat (f) has an irregular shape with two adjacent thresholds. Flat (g) is featured with a rectangular shape and two separated thresholds, while flat (h) has an irregular shape with multiple thresholds.

In addition, a saddle-shaped surface (Figure 5), which is similar to the one used in [31], was used to test the new algorithm and compare with the JD method in ArcGIS (ESRI, Redlands, CA, USA) [11,32] and DEDNM. The surface contains a flat area generated from depression filling. The flat area is surrounded by higher terrain, while the lower elevation points are located at the top and bottom boundaries (Figure 5).

To further test the new flat correction algorithm, a real-world surface in North Dakota, U.S. was selected. It is located in the northwest of Jamestown (Figure 6a). Its area is 2.5515 km² (1890 m × 1350 m) and the elevations of the area range from 549.81 to 589.36 m. Based on the 30-m DEM of the surface, the puddle delineation (PD) program [15] was used to identify all depressions and flats. The delineated depressions (including their centers and thresholds) are shown in Figure 6b. The grey cells represent depressions or flats (Figure 6b), and the green cells represent thresholds, through which water flows out of puddles (i.e., spilling). More details on the PD and its applications can be found in [15]. Due to the noises originated from the measuring techniques and interpolations, there are many small variations (<0.001 m) in the original DEM, which result in numerous tiny depressions and flats. These depressions and flats have no hydrologic significance and, in particular, affect the efficiency of surface delineation. Thus, the original elevation values were rounded to three decimal digits in this study. Any depressions with an area of smaller than 0.009 km² (i.e., 10 cells) were regarded as small “insignificant” depressions, which were first filled by using the PD program [15]. The resulting flats were then corrected by utilizing the new algorithm. Based on the original PD program [15,33], a new Windows-based program was developed to facilitate the aforementioned data processing, including delineation of depressions, filling of the identified depressions, and correction of the resulting flats.

3. Results and Analysis

3.1. Correction of Simple Flats

To demonstrate the capability of the new algorithm, both DEDNM and the new program were used to process the flat area with one threshold shown in Figure 4a and the results were compared. When combining the gradients from steps 1 and 2 in the DEDNM module, an exceptional situation occurred. Cells C4–E4, and D5 (Figure 7a) had the same gradients with opposite directions (thus, they canceled out), which led to no flow direction for these cells. As mentioned in [28], 0.5 was subsequently added to the increment numbers for cells C4–E4 and D5 (Figure 7b). However, as shown in Figure 7b, cell D5 still had no flow direction. Then, the module added another increment of 0.5 for cell D5. The increment number of cell D5 became 4 (Figure 7c). Unfortunately, the new increment number of cell D5 resulted in no flow direction for cell D6 (Figure 7c). More changes were needed until all cells had flow directions (Figure 8a). By the end, the product of the final increment number and the increment unit (2 × 10⁻⁵) was added to the original elevation for each DEM cell in the flat area. The new surface generated from this process is shown in Figure 8b. In contrast, the new algorithm with various levels of increment numbers in steps 1 and 2 did not have this problem. It was able to identify the flow direction for each cell without the requirement for any additional steps for the case of one threshold as shown in Figure 8d. Thus, the new algorithm simplified the procedure and showed an improved efficiency in this aspect.

As shown in Figure 8c,f, the flow directions determined by DEDNM and the new algorithm are toward the threshold of the flat area. Most flow directions from the two methods/programs are the same, except the flow directions of cells C6 and E6. The reason is that the increment number of 0.5 was added to cell D5 several times in DEDNM due to the issue of no flow direction. This operation in DEDNM changed the flow directions of cells C6 and E6 (Figure 7a and Figure 8c), while the new method did not have this issue. Since the new algorithm avoided the issue of no drainage directions, 1.0 × 10⁻⁵ was used as the increment unit, instead of 2.0 × 10⁻⁵ in DEDNM. Note that in Figure 8b,e, 0.02 and 0.01 were, respectively, used in DEDNM and the new algorithm, instead of 2.0 × 10⁻⁵ and 1.0×10⁻⁵ for real situations. The sum of increments in the new method is 0.0023, which is smaller than 0.00312 in DEDNM. The largest elevation increment in the new method is 0.00013 m, which is also smaller than 0.00016 m in DEDNM. The elevation increment is far smaller than the vertical resolution of the DEM, ensuring that the new elevations are lower than those of their surrounding cells and no flow directions are altered for the cells outside the flat.

3.2. Correction of Complex Flats

The seven artificial flat areas shown in Figure 4b–h were used to test the new method and further compare it with DEDNM. Figure 9 shows the final increment numbers of all cells in the flat areas from the two methods/programs. For flat areas (b) and (c), the new method and DEDNM had no significant difference, except the differences in the specific final increment numbers (Figure 9(1,2)). Both showed good performances. In both cases, all flow directions were toward the thresholds of the surfaces on the left end. The results demonstrated that like DEDNM, the new method was able to correct small flats of a linear shape.

For shape (d) (spiral flat area), DEDNM and the new method identified the flow directions toward the threshold correctly. As shown in Figure 9(3), some adjacent cells have the same increment numbers, resulting in new flats. The new algorithm can correct these flats by implementing the optional step (step 4), but DEDNM does not have this capability.

Shape (e), with a large rectangular flat area (Figure 9(4)), is more complex than shape (d). In DEDNM, the shaded cells showed no flow directions after combining the increment numbers calculated in steps 1 and 2. Hence, more computation loops were implemented by adding an increment of 0.5. Our test indicates that this case is very common in DEDNM, especially for the flats with a square area. The larger the square area is, the more repeating loops are needed (Figure 10). When many large flat areas are involved in this process, low efficiency can be a serious concern for DEDNM. However, the new method does not have this issue.

The results for surface (f), an irregular area with two adjacent thresholds (Figure 9(5)), again demonstrate the capability of the new flat correction algorithm. Similar to shape (e), the new method does not require more loops to correct the flats in step 3. Its final increments are smaller than those in DEDNM, as detailed in

Table 1. Final total elevation increments for eight surfaces (a–h) calculated by DEDNM and the new algorithm.

Surface/Shape	DEDNM	New Algorithm	Relative Difference (%)
(a)	3.12 × 10−3	2.30 × 10−3	26.28
(b)	6.00 × 10−5	4.00 × 10−5	33.33
(c)	3.00 × 10−4	2.50 × 10−4	16.67
(d)	3.50 × 10−3	3.30 × 10−3	5.71
(e)	9.84 × 10−3	8.63 × 10−3	12.30
(f)	3.29 × 10−3	2.41 × 10−3	26.75
(g)	1.20 × 10−3	8.46 × 10−4	29.50
(h)	1.14 × 10−3	8.04 × 10−4	29.52

.

For flat area (g) (Figure 9(6)), DEDNM suffered from another problem when handling the flat area with multiple thresholds. As shown in Figure 9(6A), no matter how many increment numbers were added to the two pink cells, they always had multiple flow directions. Additionally, a number of small flat areas were generated in DEDNM, such as the flats with the same increment number of 4 at the upper-right corner in Figure 9(6A). The new DEM contained many small, separate flat areas. Both multiple flow directions and new flat areas were considered in the new algorithm. By adding small random increment numbers to the increment numbers from step 3, the new method avoided creating new small flats and multiple flow direction cells (Figure 9(6B)). As shown in Table 1, the final increment from the new method for the flat area was still smaller than that from DEDNM, while it was sufficient to correct all flat areas with multiple thresholds.

Based on the increment numbers calculated by the new method for flat area (h) with a complex shape and multiple thresholds (Figure 9(7)), flow directions of all cells are toward the nearest thresholds (Figure 9(7B)). The results again demonstrate that the new algorithm is able to deal with very complex topographic surfaces. In addition, all the increments from the new method are smaller than those from DEDNM. The new method effectively reduces the increments by 5.71–33.33% (Table 1), which implies that the new algorithm has less influence on the actual elevations than DEDNM.

The results from the JD method in ArcGIS (ESRI, Redlands, CA, USA), DEDNM, and the new algorithm for the saddle-shaped surface are shown in Figure 11. The flow directions calculated by the JD method show a “parallel flow” pattern toward the three lower elevation locations (outlets) within the three drainage areas (Figure 11a). With the JD method in ArcGIS (ESRI, Redlands, CA, USA), water flows toward the lower elevation cells along the shortest flow paths. Better than the result of the JD method, the flow directions calculated by DEDNM display a convergent flow pattern (Figure 11b). However, multiple flow directions were generated in some cells (e.g., cell R11, Figure 11b). Additionally, DEDNM resulted in many step-type flats (Figure 11b) (note that the cells connected by a red dash line have the same elevation). Actually, elevation correction was performed only for the cells in the large flat area in DEDNM. Numerous small flats (flat steps) remained in the new DEM. These small flats can affect the efficiency of the delineation programs if they are sensitive to flat areas. The result from the new algorithm also shows a convergent flow pattern. The increments added in step 1 played a key role. Unlike DEDNM, more cells (e.g., cells B16, Q4, and L13) pointed toward lower-elevation cells (Figure 11c). In addition, the new algorithm eliminated the issues related to multiple flow directions and the step-type flats generated during the correcting process.

3.3. Correction of Real Depressions/Flats

As aforementioned, the original DEM of the real land surface in North Dakota was processed by keeping three decimal digits. According to the delineation results, this processed DEM showed more depressions and flats than the original DEM (Figure 6b and Figure 12a). The number of depressions and flats increased from 45 to 53 because rounding off the elevation decimals created more flats. Using the new algorithm, all small, insignificant flats were successfully corrected/removed (Figure 12b). The number of depressions and flats reduced significantly from 53 to 11, which demonstrated the effectiveness and efficiency of the new algorithm in flat correction.

4. Conclusions

Flow directions are usually calculated before hydrological analysis and modeling. However, correcting elevations for flat areas in raw DEMs or from the depression-filling process has been a challenge for the determination of flow directions and accumulations, as well as hydrotopographic analysis and applications. Particularly, some DEM correction algorithms (e.g., pit-filling methods) may create new unrealistic flats or artifacts. To solve these problems, a new algorithm, based on different levels of increment numbers, was proposed in this study. A set of artificial flat areas and two real land surfaces were used to test the new method. The results demonstrated that the new flat correction algorithm was not limited by the surrounding topography, the number of thresholds, and the shapes of flats. The new algorithm simplified the computational procedure of the traditional DEDNM by avoiding creating new flats and reduced the elevation increments by up to 33.33%. The new method also eliminated the problem associated with multiple flow directions and provided a new flat-free DEM. It avoided the “parallel flow” generated by the JD method and the step-type flats created by DEDNM. The real application showed that the new algorithm effectively corrected the elevations of small, insignificant depressions and flats (e.g., artifacts), and further improved the efficiency of surface delineation. It should be noted that since very small elevation increments are used in the new algorithm, the size of the processed/corrected DEM file can become larger. Additionally, more testing is needed after the new flat correction algorithm is incorporated into a hydrologic modeling system.