Accurate stream extraction from large , radar-based elevation models

Accurate stream extraction from large, radar-based elevation models M. Metz, H. Mitasova, and R. S. Harmon Institute of Experimental Ecology, University of Ulm, Ulm, Germany Department of Marine, Earth and Atmospheric Sciences, North Carolina State University, Raleigh, North Carolina, USA Environmental Sciences Division, Army Research Office, US Army Research Laboratory, Durham, North Carolina, USA

often low resolution coverage.These topographic data sets are increasingly used to improve mapping of geomorphic and hydrologic features in remote, hard to reach areas and at regional to global scales (e.g., Kinner et al., 2005;Lehner and D öll 2004;World Wildlife Fund, 2009).
Flow routing based extraction of hydrologic features from DEM poses several challenges: (a) elevation surfaces derived from remote sensing data such as SRTM or IFSARE include tree canopy that often requires flow routing through large depressions caused by gaps in vegetation or noise in elevation data; (b) depending on the size of the study region and resolution, data sets can be massive and require extensive processing time.Significant effort has been devoted to the development of new flow tracing and watershed analysis algorithms that support efficient processing of large DEMs and address the issue of routing through complex nested depressions (e.g., Arge et al., 2003;Danner et al., 2007).
The most widespread method for handling depressions is sink filling up to the level of the sink spill point combined with routing through the resulting flat area (Jenson and Domingue, 1988;Fairfield and Leymarie, 1991;Planchon and Darboux 2001;Wang and Liu 2006).The complementary method is known as carving or breaching (Rieger, 1998;Martz and Garbrecht, 1998) where a channel is carved out of each sink, breaking through the (artificial) obstacle.Both principles can be combined in an impact reduction approach (IRA; Lindsay and Creed, 2005) that determines for each sink the method that causes the least impact on the source dataset.All these approaches alter the elevation data in order to ensure full drainage, assuming that sinks are artefacts created either by too low (artificial pits) or too high (artificial drainage blocks) elevation values.An alternative to the modification of elevation data is to determine the least-cost drainage paths (LCP) through the unaltered terrain and out of sinks.In the implementation used here, the least-cost path search starts with potential outlet points and proceeds along the least steep uphill slope.If a depression is encountered, it follows the steepest downhill slope to the bottom of a depression and then proceeds again along the least steep uphill slope.At each step during the search, only neighbouring Introduction

Conclusions References
Tables Figures

Back Close
Full grid cells that are not yet in the search list are added to the list and drainage direction for these neighbouring grid cells is set towards the current grid cell.When applied to IFSARE derived DEMs, this approach has provided more robust flow routing through large, nested depressions with fewer artefacts than the traditional sink filling (Kinner et al., 2005).However, the implementation as a module r.watershed (Ehlschlaeger, 1989) in GRASS GIS (Neteler and Mitasova, 2008) was not designed for handling large data sets as are currently available (Arge et al., 2003).To make the method applicable to the new generation of watershed analysis tasks, redesign of the implementation for handling of massive DEMs and addition of more flow routing options with the potential to improve accuracy was necessary.Here we present the improvements and evaluation of their efficiency and accuracy when applied to radar-based DEMs.
The LCP algorithm (A* Search, Hart et al., 1968;module r.watershed in GRASS GIS, Ehlschlaeger, 1989) used for flow routing, flow accumulation, and watershed analysis was redesigned to increase the processing speed and decrease memory consumption.Its performance was then compared with the original implementation, an I/O efficient method (Arge et al., 2003), and Terrain Analysis System (TAS, Lindsay and Creed 2005).Impact of mapping technology (IFSARE, SRTM) and resampling on the extracted stream networks was also analysed.

Improved least-cost path search algorithm
The least-cost path search starts with potential outlet grid points (grid cells on map boundaries or cells with at least one neighbour with unknown elevation, e.g.masked ocean) that are kept in an open list sorted by costs.In this implementation, costs are measured as elevation and order of addition to the list (grid cells added earlier have higher precedence in case of equal elevation).The grid point with the smallest cost is extracted from the open list, put into a closed list and all neighbours of that grid point Introduction

Conclusions References
Tables Figures

Back Close
Full The search component of the original implementation was modified by augmenting it with a minimum heap used as priority queue (Atkinson et al., 1986) for both the search and the flow accumulation components have further enhanced the module's capacity to process larger DEMs.For massive datasets that can not be processed with the amount of memory available, the module can optionally use external memory with intermediate data stored on disk.These changes, aimed at improving the computational performance, do not affect the quality of the results.
Additionally, a multiple flow direction (MFD) algorithm, based on Holmgren (1994), was implemented using the path determined by the A* Search with an option to control the strength of flow convergence.For grid-based elevation models, an MFD approach conforms better to theoretical surface flow dispersion than single flow direction in which only one neighbouring cell can receive surface flow from the current cell (Freeman, 1991, Quinn et al., 1991;Holmgren, 1994).Initial D8 flow directions as determined by the LCP search are used to route flow out of depressions when all unprocessed neighbouring grid cells have higher elevation than the current grid cell.

Data sources and pre-processing
Countrywide elevation coverage of Panama was available as a 3 arc second resolution SRTM DEM.SRTM version 2 was selected for this study as the most reliable in terms of accuracy and minimal artefacts, after evaluating properties of the currently available SRTM products (v1, v2, v3, v4.1).A recent IFSARE survey has provided new, Introduction

Conclusions References
Tables Figures

Back Close
Full more detailed information about the topography in Central Panama at 10 m resolution (Fig. 1).The elevation data were not bare earth representations and over large regions include elevation surface defined by triple canopy tropical forest environment with tree heights of more than 30 m above the bare earth surface.SRTM v2 tiles covering all of Panama were combined and gaps in the dataset filled using the regularized spline with tension (RST) interpolation method (GRASS module r.fillnulls; Mitasova and Mitas, 1993).The seamless SRTM coverage was then reprojected with bicubic interpolation from geographic to a UTM zone 17N coordinate system with 90 m resolution to keep reprojection modifications to a minimum.The reprojected SRTM DEM was then resampled to 30 m resolution using the RST method for additional testing and to evaluate impact of resampling on stream extraction accuracy.
The original IFSARE data were collected at 2.5 m resolution and processed by standard procedures into 10 m resolution DEM (Kinner et al., 2005).The IFSARE DEM was converted to 30 m resolution using bicubic resampling in order to have DEMs from two sources with different level of detail but identical resolution.

Hydrological analysis and accuracy evaluation
Several tests were performed to evaluate the flow routing algorithms.First, we evaluated the performance of computing flow accumulation maps in terms of speed and capability to handle large DEMs for old and new implementations of the LCP algorithm (r.watershed) and an I/O efficient algorithm (r.terraflow).Second, we evaluated the accuracy of extracted streams for algorithms with different handling of depressions.Three approaches were tested: (i) sink filling as implemented in r.terraflow within GRASS GIS (Arge et al., 2003), (ii) the impact reduction approach (IRA) as implemented in TAS GIS (Lindsay and Creed, 2005), and (iii) the LCP algorithm as implemented in the enhanced version of r.watershed in GRASS GIS.Finally, we evaluated the impact of DEM resolution and level of detail on the accuracy of extracted streams for all tested algorithms (r.watershed, r.terraflow, and TAS).Introduction

Conclusions References
Tables Figures

Back Close
Full The tests were performed for the Central Panama region using both the IFSARE and SRTM data at the 10 m, 30 m and 90 m resolutions by computing flow accumulation using MFD routing.
Stream networks were then extracted from all MFD flow accumulation maps using a minimum upstream catchment area of 100 000 m 2 as a threshold for channel initiation.
A new channel was only initiated if it was not located within a wider stream tube which is common for MFD flow accumulation, especially in flood plains.Downstream channel tracing followed the steepest slope and maximum flow accumulation.To egress from depressions with the LCP method, streams followed the maximum flow accumulation which was in these cases identical to SFD flow accumulation.The streams extracted from MFD flow accumulation were consequently one grid cell wide and could easily be vectorized for further analysis.
Stream segments digitized from Landsat imagery (TM 5, year 2000, provided by the United States Geological Survey -USGS) and field measurements were used for accuracy assessment.Georeferenced stream data (GPS points) were collected in the field using a Corvallis Microtechnology MARCH v3.7 GPS unit with ∼2 m positional accuracy during the years 2002 to 2007 at sites in the Chagres river watershed (Fig. 1) and during 2005-2009 at lower reaches of most major rivers across Panama.Assessment of stream extraction accuracy was then performed by calculating distances between the stream segments derived from the IFSARE and SRTM DEMs and (a) 338 on-ground GPS measurements plus (b) 995 points along rivers digitized from Landsat imagery (Fig. 1).The GPS measurements included a larger proportion of points acquired in mountainous regions along smaller rivers, while most of the points digitized from Landsat imagery were located along larger rivers in flood plain and coastal plain landscapes.We did not use the SRTM water body database because this cannot be regarded as an independent data source since it was used in the creation of the SRTM DSM v2 and because even the larger rivers visible on Landsat imagery were still too narrow to be captured in the SRTM water body database.Introduction

Conclusions References
Tables Figures

Back Close
Full

Performance comparison
To illustrate the improvement in computational performance of the new implementation of the LCP method, we compared the processing times needed by the new and old version of r.watershed for flow accumulation computation at 90 m resolution for all of Panama and for the Central Panama (Table 1).The new version was 350 times faster than the old version for the Central Panama area represented by a relatively small DEM with 2 million grid cells at 90 m resolution.The improvement was even more dramatic for the countrywide coverage with 27 million grid cells: the new version was about 1750 times faster.Processing the larger regions used in this study was impractical with the old version because it could easily take days (Arge et al., 2003) whereas the new version required only minutes to execute.
As long as all intermediate data fit into memory, the processing time of the new r.watershed in all-in-memory mode is shorter than that of the two other programs tested.If data do not fit into memory, r.watershed can use segmented processing, with intermediate data being stored on disk.This leads to a longer processing time than for the all-in-memory mode.An overview of grid sizes and processing times is presented in Table 1.Observed processing times in the segmented mode were still shorter than for r.terraflow which uses theoretically I/O efficient algorithms specifically designed for large datasets (Arge et al., 2003).The size of intermediate data created by r.watershed in the segmented mode is less than 20% of the size of intermediate data created by r.terraflow.This can be of advantage if intermediate data are just a bit too large to fit into memory.In theory, r.terraflow with its advanced disk I/O algorithms should be faster than r.watershed for regions much larger than the ones processed here.
It was not possible to process the IFSARE DEM at 10 m resolution and the SRTM DEM at 30 m resolution for all of Panama by TAS GIS because memory requirements exceeded available system memory.TAS GIS is a 32 bit application for Microsoft Introduction

Conclusions References
Tables Figures

Back Close
Full Windows and can only utilize as much memory as a 32 bit application can manage.
The test system provided as much memory as TAS GIS could theoretically utilize, and TAS GIS, contrary to r.terraflow and r.watershed, keeps all intermediate data in memory, hence the out-of-memory error was regarded as a limitation of this application and not as a limitation of the test system.Therefore the SRTM DSM at 30 m resolution was processed in TAS GIS only for Central Panama as covered by the IFSARE DSM and results of the IRA were not available for IFSARE at 10 m resolution.At the time of this writing IRA was not available for WhiteboxGAT (Lindsay, 2009) that has replaced TAS.

Comparison between different methods of sink treatment
The median distances of reference points to extracted streams for each sink treatment method and each processed DEM are given in Table 2. Generally, the location of extracted streams improved with increasing detail available in the DEM.While the stream locations for the SRTM 90 m DEM were about 100 m away from reference points, the median distances improved to down to 34.4 m for the IFSARE 10 m DEM.Streams extracted from LCP flow accumulation were closer to the GPS field points (T-tests, all p<0.0001) for all DEMs at all tested resolutions (IFSARE 10 m, 30 m, and SRTM 30 m, 90 m) when compared with the sink filling method (Fig. 2).The streams extracted from LCP were also closer to the points digitized from Landsat for all DEMs, however, the distance difference was not significant for SRTM 90 m (T-test, p=0.23).
When compared with the IRA method (Fig. 3), streams extracted from LCP accu-

Comparison between different levels of detail at constant resolution
IFSARE 30 m and SRTM 30 m were interpolated to the same resolution from the source DEMs with vastly different levels of detail captured by their original resolutions of 10 m and 90 m respectively.We investigated how much the different levels of detail influenced the accuracy of extracted streams for all three sink treatment methods (Fig. 4).
The distance between the reference points and streams extracted from IFSARE 30 m was subtracted from the distance between the reference points and streams extracted from SRTM 30 m to measure the difference in accuracy of streams extracted from these two DEMs.For larger rivers in flat areas (Landsat points), the accuracy of stream locations was considerably higher for IFSARE 30 m than for SRTM 30 m (T-test, all p< 0.001; on average 35 m closer) for all methods.For smaller rivers in mountainous regions, only IRA showed a significant improvement in accuracy from SRTM 30 m to IF-SARE 30 m (T-test, p=0.001).However, as shown in the previous section for IFSARE 30 m, IRA stream locations in mountainous regions were significantly less accurate than LCP stream locations.

Effects of RST resampling on stream extraction accuracy
In order to assess the effect of SRTM resampling from 90 m to 30 m resolution by the RST method, the streams extracted from SRTM 90 m were compared with streams extracted from the resampled SRTM 30 m by calculating the difference in stream distances to reference points (Fig. 5).Corresponding median distances of reference points to stream locations are given in Table 2.The accuracy of stream locations was improved by resampling for larger rivers in flat landscape for all sink treatment methods (T-tests, all p<0.001), however, for GPS points in mountainous regions only the results obtained by the least-cost path search were improved (T-test, p<0.001).Introduction

Conclusions References
Tables Figures

Back Close
Full

Discussion
Sink (depression) removal, or hydrological conditioning, is a pre-processing step required by a number of applications for hydrological analysis and modelling.Handling of depressions during flow routing is an integral component of most flow accumulation and watershed analysis tools.It is safe to assume that no DEM, including the bare earth models, is a perfect representation of the terrain, and because every DEM is a model and approximation of the real terrain, it will contain a certain amount of random noise and errors.Nevertheless, actual terrain includes many true depressions, especially in riffle-and-pool portions of mountainous rivers, and on floodplains and coastal plains.Therefore, sink filling is not always an appropriate approach for DEM improvement because the measured elevation data are replaced with new elevation values based on the rather artificial condition that the terrain has no depressions.Although many hydrologic models require depression-less DEMs, removal of sinks can alter the true hydrological state that includes standing water in depressions.Moreover, in IFSARE and SRTM DEMs that include vegetation cover, many of the observed depressions are caused by gaps in canopy, or by vegetation partially obstructing narrow and deep valleys.The filling of such features requires altering elevations by several meters over areas of hundreds of m 2 leading to large flat areas and an artificial stream geometry that does not represent the ground truth.Figure 6 presents an example of the amount of modification that can result from such filling.Preserving most of the drainage directions of the original DEM is preferable for hydrological conditioning in this case.Our results suggest that the adverse effects of sink filling become more pronounced for higher resolution DEMs (see Fig. 7 for streams extracted from IFSARE 10 m).By contrast, at lower resolutions such as SRTM 90 m, the coarsest resolution tested in this study, there was often no significant difference in stream locations between the compared methods, particularly for Landsat reference points.Landsat points where predominantly located in relatively flat landscape where both the horizontal and the vertical resolution of SRTM 90 m were often not sufficient to locate rivers Introduction

Conclusions References
Tables Figures

Back Close
Full  Full  Full Mission (SRTM), Interferometric Synthetic Aperture Radar for Elevation (IFSARE) surveys and, since 2009, the ASTER Global Digital Elevation Model (ASTER GDEM) provide elevation data in regions that have had only limited, Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | that are neither in the open or closed list are added to the open list.Thus the search proceeds upstream along the least steep slope, in each step assigning flow directions to the neighbours of the current focus cell.The search terminates when all grid points are processed.The LCP search provides two results: (a) flow direction for each cell in standard D8 manner, and (b) the order in which cells must be processed for flow accumulation.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | mulation were closer to the GPS field points for IFSARE 30 m (T-test, p=0.018) and SRTM 30 m (T-test, p=0.002), but for the SRTM 90 m the difference was not significant (T-test, p=0.08).For the points digitized from Landsat, streams extracted from LCP accumulation were not closer than IRA streams for any of the tested DEMs (T-tests, all p>0.Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Table 1 .
Processing time for the different DEMs in minutes.

Table 2 .
Median distances of reference points to nearest extracted stream in meters.