Next Article in Journal
A 3D Digital Cadastre for New Zealand and the International Opportunity
Next Article in Special Issue
Roughness Spectra Derived from Multi-Scale LiDAR Point Clouds of a Gravel Surface: A Comparison and Sensitivity Analysis
Previous Article in Journal
Spatio-Temporal Series Remote Sensing Image Prediction Based on Multi-Dictionary Bayesian Fusion
Previous Article in Special Issue
Visualization of Features in 3D Terrain
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
School of Water Resources and Environment, China University of Geosciences (Beijing), Beijing 100083, China
2
Department of Civil and Environmental Engineering (Dept 2470), North Dakota State University, P.O. Box 6050, Fargo, ND 58108-6050, USA
*
Author to whom correspondence should be addressed.
ISPRS Int. J. Geo-Inf. 2017, 6(11), 376; https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi6110376
Submission received: 13 September 2017 / Revised: 19 November 2017 / Accepted: 20 November 2017 / Published: 21 November 2017
(This article belongs to the Special Issue Leading Progress in Digital Terrain Analysis and Modeling)

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−5 m is selected as the increment unit, which is smaller than the increment unit (2.0 × 10−5 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−5 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 km2 (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 km2 (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−5) 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−5 was used as the increment unit, instead of 2.0 × 10−5 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−5 and 1.0×10−5 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.
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.

Acknowledgments

This material is based upon work supported by the National Science Foundation under grant no. NSF EPSCoR Award IIA-1355466 and the Fundamental Research Funds for the Central Universities (no. 2652017180). The authors would like to thank Jianli Zhang for his contributions to the related work.

Author Contributions

Xuefeng Chu and Xingwei Liu proposed and designed the research. All authors made fairly equal contributions to the writing of the paper.

Conflicts of Interest

The author declare no conflict of interest.

References

  1. Nasab, M.T.; Zhang, J.; Chu, X. A new depression-dominated delineation (D-cubed) method for improved watershed modeling. Hydrol. Process. 2017. [Google Scholar] [CrossRef]
  2. Chu, X.; Yang, J.; Chi, Y.; Zhang, J. Dynamic puddle delineation and modeling of puddle-to-puddle filling-spilling-merging-splitting overland flow processes. Water Resour. Res. 2013, 49, 3825–3829. [Google Scholar] [CrossRef]
  3. Zhang, W.H.; Montgomery, D.R. Digital elevation model grid size, landscape representation, and hydrologic simulations. Water Resour. Res. 1994, 30, 1019–1028. [Google Scholar] [CrossRef]
  4. Jain, M.K.; Singh, V.P. DEM-based modelling of surface runoff using diffusion wave equation. J. Hydrol. 2005, 302, 107–126. [Google Scholar] [CrossRef]
  5. Fairfield, J.; Leymarie, P. Drainage networks from grid digital elevation models. Water Resour. Res. 1991, 27, 709–717. [Google Scholar] [CrossRef]
  6. Grimaldi, S.; Teles, V.; Bras, R.L. Sensitivity of a physically based method for terrain interpolation to initial conditions and its conditioning on stream location. Earth Surf. Process. Landf. 2004, 29, 587–597. [Google Scholar] [CrossRef]
  7. Kenny, F.; Matthews, B.; Todd, K. Routing overland flow through sinks and flats in interpolated raster terrain surfaces. Comput. Geosci. 2008, 34, 1417–1430. [Google Scholar] [CrossRef]
  8. Grimaldi, S.; Teles, V.; Bras, R.L. Preserving first and second moments of the slope area relationship during the interpolation of digital elevation models. Adv. Water Resour. 2005, 28, 583–588. [Google Scholar] [CrossRef]
  9. Grimaldi, S.; Nardi, F.; Di Benedetto, F.; Istanbulluoglu, E.; Bras, R.L. A physically-based method for removing pits in digital elevation models. Adv. Water Resour. 2007, 30, 2151–2158. [Google Scholar] [CrossRef]
  10. Nardi, F.; Grimaldi, S.; Santini, M.; Petroselli, A.; Ubertini, L. Hydrogeomorphic properties of simulated drainage patterns using digital elevation models: The flat area issue/propriétés hydro-géomorphologiques de réseaux de drainage simulés à partir de modèles numériques de terrain: La question des zones planes. Hydrol. Sci. J. 2008, 53, 1176–1193. [Google Scholar] [CrossRef]
  11. Jenson, S.K.; Domingue, J.O. Extracting topographic structure from digital elevation data for geographic information system analysis. Photogramm. Eng. Remote Sens. 1988, 54, 1593–1600. [Google Scholar]
  12. Tribe, A. Automated recognition of valley lines and drainage networks from grid digital elevation models: A review and a new method. J. Hydrol. 1992, 139, 263–293. [Google Scholar] [CrossRef]
  13. O’Callaghan, J.F.; Mark, D.M. The extraction of drainage networks from digital elevation data. Comput. Vis. Graph. Image Process. 1984, 28, 323–344. [Google Scholar] [CrossRef]
  14. Pan, F.; Stieglitz, M.; McKane, R.B. An algorithm for treating flat areas and depressions in digital elevation models using linear interpolation. Water Resour. Res. 2012, 48. [Google Scholar] [CrossRef]
  15. Chu, X.; Zhang, J.; Chi, Y.; Yang, J. An improved method for watershed delineation and computation of surface depression storage. In Proceedings of the Watershed Management 2010: Innovations in Watershed Management under Land Use and Climate Change, Madison, WI, USA, 23–27 August 2010; pp. 1113–1122. [Google Scholar]
  16. Jana, R.; Reshmidevi, T.V.; Arun, P.S.; Eldho, T.I. An enhanced technique in construction of the discrete drainage network from low-resolution spatial database. Comput. Geosci. 2007, 33, 717–727. [Google Scholar] [CrossRef]
  17. Zhang, H.; Huang, G. Building channel networks for flat regions in digital elevation models. Hydrol. Process. 2009, 23, 2879–2887. [Google Scholar] [CrossRef]
  18. Zhang, S.; Zhao, B.; Erdun, E. Watershed characteristics extraction and subsequent terrain analysis based on digital elevation model in flat region. J. Hydrol. Eng. 2014, 19. [Google Scholar] [CrossRef]
  19. Su, C.; Wang, X.Z.; Feng, C.J.; Huang, Z.C.; Zhang, X.C. An integrated algorithm for depression filling and assignment of drainage directions over flat surfaces in digital elevation models. Earth Sci. Inform. 2015, 8, 895–905. [Google Scholar] [CrossRef]
  20. Su, C.; Feng, C.J.; Wang, X.Z.; Huang, Z.C.; Zhang, X.C. An efficient algorithm for assignment of flow direction over flat surfaces in raster DEMs based on distance transform. Earth Sci. Inform. 2016, 9, 225–233. [Google Scholar] [CrossRef]
  21. Zhang, H.; Yao, Z.; Yang, Q.; Li, S.; Baartman, J.E.; Gai, L.; Yao, M.; Yang, X.; Ritsema, C.J.; Geissen, V. An integrated algorithm to evaluate flow direction and flow accumulation in flat regions of hydrologically corrected DEMs. Catena 2017, 151, 174–181. [Google Scholar] [CrossRef]
  22. Garbrecht, J.; Martz, L.W. Automated channel ordering and node indexing for raster channel networks. Comput. Geosci. 1997, 23, 961–966. [Google Scholar] [CrossRef]
  23. Kite, G. Modelling the Mekong: Hydrological simulation for environmental impact studies. J. Hydrol. 2001, 253, 1–13. [Google Scholar] [CrossRef]
  24. Alsdorf, D.E. Water storage of the central amazon floodplain measured with GIS and remote sensing imagery. Ann. Assoc. Am. Geogr. 2003, 93, 55–66. [Google Scholar] [CrossRef]
  25. White, A.B.; Kumar, P.; Saco, P.M.; Rhoads, B.L.; Yen, B.C. Hydrodynamic and geomorphologic dispersion: Scale effects in the Illinois River basin. J. Hydrol. 2004, 288, 237–257. [Google Scholar] [CrossRef]
  26. Clennon, J.A.; King, C.H.; Muchiri, E.M.; Kitron, U. Hydrological modelling of snail dispersal patterns in Msambweni, Kenya and potential resurgence of Schistosoma haematobium transmission. Parasitology 2007, 134, 683–693. [Google Scholar] [CrossRef] [PubMed]
  27. Phillips, J.D.; Slattery, M.C. Downstream trends in discharge, slope, and stream power in a tower coastal plain river. J. Hydrol. 2007, 334, 290–303. [Google Scholar] [CrossRef]
  28. Clarke, S.E.; Burnett, K.M.; Miller, D.J. Modeling streams and hydrogeomorphic attributes in Oregon from digital and field data. J. Am. Water Resour. Assoc. 2008, 44, 459–477. [Google Scholar] [CrossRef]
  29. Barnes, R.; Lehman, C.; Mulla, D. An efficient assignment of drainage direction over flat surfaces in raster digital elevation models. Comput. Geosci. 2014, 62, 128–135. [Google Scholar] [CrossRef]
  30. Yang, J.; Chu, X.F. A new modeling approach for simulating microtopography-dominated, discontinuous overland flow on infiltrating surfaces. Adv. Water Resour. 2015, 78, 80–93. [Google Scholar] [CrossRef]
  31. Garbrecht, J.; Martz, L.W. The assignment of drainage direction over flat surfaces in raster digital elevation models. J. Hydrol. 1997, 193, 204–213. [Google Scholar] [CrossRef]
  32. Ormsby, T.; Napoleon, E.J.; Burke, R.; Groessl, C.; Bowden, L. Getting to Know Arcgis Desktop; Esri Press: Redlands, CA, USA, 2010; pp. 527–536. [Google Scholar]
  33. Chu, X. Delineation of pothole-dominated wetlands and modeling of their threshold behaviors. J. Hydrol. Eng. 2017, 22. [Google Scholar] [CrossRef]
Figure 1. (a) Original DEM; and (b) DEM after depression filling.
Figure 1. (a) Original DEM; and (b) DEM after depression filling.
Ijgi 06 00376 g001
Figure 2. Increment numbers for steps 1–3 in the new algorithm (arrows represent the corrected flow directions across the flat area): (a) increment numbers at the beginning in step 1; (bd) increment numbers after the first, second, and third loops respectively in step 1; (e) flow directions as a result of step 1; (f-g) increment numbers after the first and second loops respectively in step 2; (h) flow directions as a result of step 2; (i) final increment numbers after combining the gradients from step 1 and step 2; and (j) final flow directions of the corrected flat area.
Figure 2. Increment numbers for steps 1–3 in the new algorithm (arrows represent the corrected flow directions across the flat area): (a) increment numbers at the beginning in step 1; (bd) increment numbers after the first, second, and third loops respectively in step 1; (e) flow directions as a result of step 1; (f-g) increment numbers after the first and second loops respectively in step 2; (h) flow directions as a result of step 2; (i) final increment numbers after combining the gradients from step 1 and step 2; and (j) final flow directions of the corrected flat area.
Ijgi 06 00376 g002
Figure 3. Procedures for flat areas with multiple thresholds: (a) increment numbers in step 1; (b) increment numbers in step 2; (c) combined increment numbers from steps 1 and 2; (d) final increment numbers after adding random increments; and (e) final flow directions of the corrected DEM.
Figure 3. Procedures for flat areas with multiple thresholds: (a) increment numbers in step 1; (b) increment numbers in step 2; (c) combined increment numbers from steps 1 and 2; (d) final increment numbers after adding random increments; and (e) final flow directions of the corrected DEM.
Ijgi 06 00376 g003
Figure 4. Topographic characteristics of different testing surfaces: (a) a normal single-threshold flat area; (b-c) linearly distributed flat areas; (d) spiral-shaped flat area; (e) combination of flats (d) and (a); (f) flat area with two adjacent thresholds; (g) rectangular flat area with two separated thresholds; and (h) irregular flat area with multiple thresholds.
Figure 4. Topographic characteristics of different testing surfaces: (a) a normal single-threshold flat area; (b-c) linearly distributed flat areas; (d) spiral-shaped flat area; (e) combination of flats (d) and (a); (f) flat area with two adjacent thresholds; (g) rectangular flat area with two separated thresholds; and (h) irregular flat area with multiple thresholds.
Ijgi 06 00376 g004
Figure 5. Saddle-shaped surface surrounded by higher terrain.
Figure 5. Saddle-shaped surface surrounded by higher terrain.
Ijgi 06 00376 g005
Figure 6. (a) Location of the selected surface; (b) delineated depressions/flats based on the original DEM.
Figure 6. (a) Location of the selected surface; (b) delineated depressions/flats based on the original DEM.
Ijgi 06 00376 g006
Figure 7. Step 3 of DEDNM for flat correction.
Figure 7. Step 3 of DEDNM for flat correction.
Ijgi 06 00376 g007
Figure 8. Final increment numbers, new elevations, and flow directions for the normal flat determined by DEDNM (ac) and the new algorithm (df).
Figure 8. Final increment numbers, new elevations, and flow directions for the normal flat determined by DEDNM (ac) and the new algorithm (df).
Ijgi 06 00376 g008
Figure 9. Final increment numbers from two methods (A: DEDNM; B: new algorithm) for seven different flat surfaces (bh) (note that M represents the repeating times to correct no drainage cells in DEDNM).
Figure 9. Final increment numbers from two methods (A: DEDNM; B: new algorithm) for seven different flat surfaces (bh) (note that M represents the repeating times to correct no drainage cells in DEDNM).
Ijgi 06 00376 g009
Figure 10. Application of DEDNM to flats with a square area (N = length of the square area; M = repeating times for processing no drainage cells).
Figure 10. Application of DEDNM to flats with a square area (N = length of the square area; M = repeating times for processing no drainage cells).
Ijgi 06 00376 g010
Figure 11. Flow directions after flat correction determined by (a) JD method in ArcGIS (ESRI, Redlands, CA, USA); (b) DEDNM, and (c) the new algorithm (the arrows indicate flow directions; the cells with a red dash line have the same elevation; and the dark solid lines show drainage boundaries).
Figure 11. Flow directions after flat correction determined by (a) JD method in ArcGIS (ESRI, Redlands, CA, USA); (b) DEDNM, and (c) the new algorithm (the arrows indicate flow directions; the cells with a red dash line have the same elevation; and the dark solid lines show drainage boundaries).
Ijgi 06 00376 g011
Figure 12. (a) Delineated depressions and flats based on the DEM with an accuracy of three decimal digits; and (b) delineated depressions and flats using the new DEM after correcting small depressions/flats.
Figure 12. (a) Delineated depressions and flats based on the DEM with an accuracy of three decimal digits; and (b) delineated depressions and flats using the new DEM after correcting small depressions/flats.
Ijgi 06 00376 g012
Table 1. Final total elevation increments for eight surfaces (a–h) calculated by DEDNM and the new algorithm.
Table 1. Final total elevation increments for eight surfaces (a–h) calculated by DEDNM and the new algorithm.
Surface/ShapeDEDNMNew AlgorithmRelative Difference (%)
(a)3.12 × 10−32.30 × 10−326.28
(b)6.00 × 10−54.00 × 10−533.33
(c)3.00 × 10−42.50 × 10−416.67
(d)3.50 × 10−33.30 × 10−35.71
(e)9.84 × 10−38.63 × 10−312.30
(f)3.29 × 10−32.41 × 10−326.75
(g)1.20 × 10−38.46 × 10−429.50
(h)1.14 × 10−38.04 × 10−429.52

Share and Cite

MDPI and ACS Style

Liu, X.; Wang, N.; Shao, J.; Chu, X. An Automated Processing Algorithm for Flat Areas Resulting from DEM Filling and Interpolation. ISPRS Int. J. Geo-Inf. 2017, 6, 376. https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi6110376

AMA Style

Liu X, Wang N, Shao J, Chu X. An Automated Processing Algorithm for Flat Areas Resulting from DEM Filling and Interpolation. ISPRS International Journal of Geo-Information. 2017; 6(11):376. https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi6110376

Chicago/Turabian Style

Liu, Xingwei, Ning Wang, Jingli Shao, and Xuefeng Chu. 2017. "An Automated Processing Algorithm for Flat Areas Resulting from DEM Filling and Interpolation" ISPRS International Journal of Geo-Information 6, no. 11: 376. https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi6110376

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop