A probabilistic approach to river network detection in digital elevation models

Catena 87 (2011) 341–350

Contents lists available at ScienceDirect

Catena j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / c a t e n a

A probabilistic approach to river network detection in digital elevation models Laura Poggio ⁎, Pierre Soille The James Hutton Institute, Aberdeen, UK Joint Research Centre, European Commission, Ispra, Italy

a r t i c l e

i n f o

Article history: Received 16 June 2010 Received in revised form 1 July 2011 Accepted 3 July 2011 Keywords: River networks position DEM Probabilistic approach Channel probability map

a b s t r a c t River networks are often derived from digital terrain models and are affected by uncertainty and errors of the corresponding elevation data. The analysis of the spatial distribution of the errors provides information on the confidence level of the derived networks. However indications on the most probable river network as a whole are missing. This study proposes a method to indicate which is the river network maximising the sum of the probability values along the network, given a map reporting the likelihood that a cell belongs to the network itself. The method is considering the inverse of the channel probability map as pseudo-DEM from which drainage networks are derived. A reference network is used to assess the spatial match of the extracted river networks using the Euclidean distance as simple comparison parameter. The network extracted from the inverse of the channel probability map is the closest to the reference. The use of a probabilistic approach to error modelling significantly increases the values of channel probability for extracted river networks and the spatial match with a ground reference dataset. © 2011 Elsevier B.V. All rights reserved.

1. Introduction The availability of digital elevation data and advances in computer processing power have contributed to the increasing interest in conducting terrain and hydrologic analysis using digital elevation models (DEMs) on local and global scales (Tarboton, 2003). Morphological features such as slope, aspect, flow length, contributing areas, drainage divides, and channel networks can be extracted from DEMs. For the delineation of watershed and drainage networks, a large number of techniques and algorithms, such as 3-by-3 moving windows calculation for detecting surface specific points, topographic features tracing, and accumulative flow analysis have been developed (e.g. Garbrecht and Martz, 2000; Lin et al., 2006). The accuracy of the results for features derived from elevation data depends on numerous factors, including the underlying terrain, the algorithms employed, the resolution and the uncertainty of the data. DEM datasets contain errors originating from a variety of sources: sampling, measurement, and interpolation (Burrough and McDonell, 2000; Pike, 2002; Wise, 1998). The modelling of digital terrain data introduces increasing uncertainty in the elevation values due to data manipulation. Even small elevation errors can greatly affect derived features (Wilson and Gallant, 2000). Therefore, it is necessary to take into account the uncertainty in the analysis based on DEM information (Darnell et al., 2008; Fisher, 1998; Lindsay, 2006).

⁎ Corresponding author. E-mail address: [email protected] (L. Poggio). 0341-8162/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.catena.2011.07.001

Monte Carlo simulation assumes that the considered map is only one realisation of a host of potential realisations. Each cell therefore can be represented by a probability distribution function (PDF) and each cell has a known mean and variance. A value is drawn from the PDF for each cell. This process repeated many times generates a set of realisation maps. A stochastic Monte Carlo simulation gives a set of potential realisations of the DEM and of the derived features (Hunter and Goodchild, 1997; Lindsay, 2006; Wechsler, 2007). Such a modelling approach can give more information on the spatial uncertainty of the original datasets and on the confidence level of the results for derived features (Darnell et al., 2008; Oksanen and Sarjakoski, 2006; Temme et al., 2008; Wechsler, 2007). The statistical analysis of the set of realisations provides information on the probability that a cell belongs to a river network. However this analysis does not indicate the actual location of the most probable river network, i.e. the one that maximises the probabilities obtained through the modelling of the probability distribution functions on the dataset realisations. The underlying assumptions of the Monte Carlo simulation procedure as applied to DEM uncertainty assessment are as follows: a) DEM error exists and constitutes uncertainty that is propagated with manipulation of the elevation data; b) the nature and extent of these errors are unknown; c) DEM error can be represented by a distribution of DEM realisations; and d) the true elevation lies somewhere within this distribution. This study is proposing a method to identify the river network that maximises the probability values obtained through an error simulation of digital terrain models. The method is considering the notchannel probability map as pseudo-DEM from which drainage

342

L. Poggio, P. Soille / Catena 87 (2011) 341–350

networks are derived. The elevation errors computed as the difference between the four considered DEMs and a reference one were used as input for Monte Carlo simulations to obtain a probability estimation of the position of the river networks. Finally the results were compared with an independent river network used as reference. The paper is organised as follows: Section 2 describes the test area and the datasets used; Section 3 holds the methodological part with the description of the DEM processing and the various methods used in the study, such as pit removal and flow-routing; Section 4 reports about the extraction of the most probable river network and the main findings; and Section 5 address the comparison of the results obtained in the previous section with a reference river network. The conclusions are given in Section 6.

Table 1 Summary statistics of test area.

2. Test area and datasets

et al., 2003), iii) coregistered, and, if necessary, iv) resampled to a common resolution of 30 m using cubic convolution interpolation (Keys, 1981; Park and Schowengerdt, 1983). The EuroDEM, (2008) was used as reference elevation dataset (EuroGeographics, 2008; Hovenbitzer, 2008). A reference river network (DTK in the following text) was obtained from the Digital Topographic Maps (Digitalen Topographischen Karten, DTK) provided by the Web Map Service (WMS) of the Rheinland-Pfalz region at 1:5000 scale. Recently, the availability of LIDAR (Light Detection And Ranging) and other very high resolution data has increased and tests were performed to measure the suitability of LIDAR derived DEM for river networks extraction (Anders et al., 2009; Passalacqua et al., 2010; Pirotti and Tarolli, 2010). Very high resolution LIDAR data resolve all fine structures of the surface of the Earth so that LIDAR elevation data need to be preprocessed to remove objects above the ground such as buildings and vegetation. Although this process can be partly automated (Chen et al., 2007; Zhang et al., 2003), user interaction is often required to remove remaining objects. Furthermore, remaining anthropic modifications of the landscape like road ditching and buried drainage systems need to be taken into account for successfully deriving drainage networks from LIDAR data (Murphy et al., 2008). Unfortunately, information addressing this need is seldom available and cannot be mapped in the field for large areas. LIDAR data can

Area (km2) Elevation (m)

Slope (degrees)

Drainage density Stream frequency Number of streams

The area used for the test is part of the Rhine basin with a size of approximately 120 km 2. It presents a variety of land uses and morphological characteristics (Fig. 1). The main morphological and hydrological features are summarised in Table 1. The DEMs used in this study were chosen with different original cell resolution, origins and error structure, and for their availability at European continental scale (Table 2 for a summary of the main characteristics): 1. GDEM: Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) dataset version 1.0 with an average of 9 stereopairs used on the test area; 2. SRTM: The Shuttle Radar Topographic Mission (SRTM) original digital elevation data produced by NASA (Rodríguez et al., 2006) were further processed to fill in no-data voids (Jarvis et al., 2006); 3. IMAP_DSM: NEXTMap® Digital Surface Model (DSM) commercial products of Intermap Technologies Inc. 4. IMAP_DTM: NEXTMap® Digital Terrain Model (DTM) commercial products of Intermap Technologies Inc. The datasets were pre-processed as follow: i) mosaicked for the available scenes, ii) projected to the ETRS89-LAEA projection (Annoni

s

0

5

120 142 698 375 39 0 11 1.8 2.9 350

b

10 Kilometers Dhro

m

re

Nim

Prü

Sû

a

Maximum Minimum Average Maximum Minimum Average

Mose

n

Ky

ll

lle

r Felle

DE

r

Al

we

Ru

tb ac

h

B ac h

Syre

LU

River network

Test area

Prims

NAMEDRIVERS

River network Ga nd

Sa

er

rre

Test area

FR Fig. 1. Test area. DE for Germany, FR for France and LU for Luxembourg. Overall river network derived from CCM2 (Catchment characterisation and modelling; (Vogt et al., 2007)).

L. Poggio, P. Soille / Catena 87 (2011) 341–350

3. DEM processing

Table 2 Summary characteristics of the considered datasets. Dataset name

Acquisition method

Resolution Source

SRTM

Synthetic aperture radar interferometry Stereopairs of optical images

3 arc sec

GDEM

1 arc sec

IMAP DSM and DTM

Airborne interferometric 0.3 arc sec synthetic aperture radar EuroDEM (Reference Mosaicked from different 1 arc sec elevation dataset) sources

343

CIAT (Jarvis et al., 2006) ERSDAC (Earth Remote Sensing Data Analysis Centre) Intermap Technologies Inc. Eurogeographic (Hovenbitzer, 2008)

provide detailed information on the morphology of the area. However their global coverage is still limited and they cannot be used for large scale applications.

The selected datasets were processed in order to obtain the channel probability maps. The processing flow (Fig. 2a) is composed of three main steps detailed below: i) probabilistic error simulations, ii) removal of spurious pits, and iii) extraction of connected stream networks from the datasets listed in Section 4.

3.1. DEM probabilistic error simulations The probability distribution was analysed using Monte Carlo approach with consideration of spatial variability of the errors (Temme et al., 2008). The procedure followed can be summarised in three main steps (Fig. 2b): (a) calculation of elevation errors, (b) modelling of elevation errors, and (c) simulations.

(a) DEM processing (the steps in the grey area are detailed below)

(b) DEM simulations

Fig. 2. Flow charts of the DEM processing (a) and simulations (b) methods, to be applied to each of the DEM used in the study (2).

L. Poggio, P. Soille / Catena 87 (2011) 341–350

DEM′ = DEM + kriged error

100

80

60

40

20

1000

2000

3000

4000

Distance (m)

(b) SRTM

150

100

50

The procedure was repeated 100 times. The number of realisations had to be sufficiently large to obtain stable results depending on how accurate the results of the uncertainty analysis should be. Indeed, the accuracy of the method is proportional to the square root of the number of runs (Wechsler and Kroll, 2006). In order to assess if 100 realisations were sufficient, a high number of random DEMs was calculated for each dataset in order to calculate the variance explained by the model, and thus the confidence level of the results. A random Gaussian noise was generated with mean, με, and standard deviation, σε, calculated for the errors model of each dataset. It was then added to the original data. The process was repeated N = 500 times and the standard deviation was calculated across all the N generated grids. The standard deviation obtained for all N random DEMs was considered as the total variation. Then, the standard deviations of N − 1, N − 2, … random DEMs were calculated in order to determine the minimum number k of iterations such that σk ≥ tσN where t is a given threshold value. Fig. 4 shows the results for two of the considered datasets. The values of explained variance for 100 iterations are always above 90% (solid thick vertical line in the figure), i.e., σ100 ≥ 0.9σN = 500 for all DEMs. In Fig. 4 the other vertical lines

Table 3 Error models analysis: parameters.

με σε Range (m) Nugget Partial Sill Model

(a) GDEM

semivariance

Initial experiments for the application of the procedure to river network extraction are presented and discussed in Poggio and Soille (2008) following the method described by Temme et al. (2008) for topographic and geomorphological parameters. The error propagation approach for assessing the uncertainty of stream networks derived from elevation data is also applied in Hengl et al. (2010). The errors of each of the DEMs considered were calculated as the difference between the elevation data and a set of points randomly extracted from the EuroDEM reference DEM. The spatial variability of the elevation errors was assessed and a variogram was derived assuming that the modelling of the DEM error was possible with a second-order stationary Gaussian random field (Oksanen and Sarjakoski, 2005; Temme et al., 2008). For each dataset a variogram was derived and defined in terms of lag distance, range, nugget, and partial sill (Table 3.). The exponential and spherical spatial autocorrelation models (Cressie, 1993) were selected to represent the correlation structure of the DEM error, as they have been found to be realistic in earlier works (Holmes et al., 2000). The geoR package of the R software (http://www.r-project.org/) was used for fitting the variograms (Ribeiro and Diggle, 2001). The fitted variograms (Fig. 3) for all the considered DEMs show a short range and a nugget that is significantly smaller than the sill, confirming the spatial autocorrelation in the datasets (Cambardella et al., 1994). The points were interpolated using the derived variogram parameters within a Gaussian simulation approach, as implemented in the GSTAT package (Pebesma, 2004) of the R software. The simulated DEM values used in further analyses were calculated as the sum of each original DEM with the simulated interpolated elevation errors (kriged error in the following text) obtained with the conditional Gaussian simulations Temme et al., 2008:

semivariance

344

GDEM

SRTM

IMAP_DSM

IMAP_DTM

6.64 30.90 3863 47.4 58.9 Spherical

− 0.03 25.75 2226 59.3 124.3 Spherical

2.02 19.87 1967 54.6 121.7 Spherical

3.12 18.79 2145 46.5 57.5 Spherical

1000

2000

3000

4000

Distance (m) Fig. 3. Error models analysis: variograms.

represent the number of iteration needed to explain 95% of the variance for each dataset. The selected number of iterations (100) corresponds to more than 90% explained variance for all datasets considered.

3.2. Pit removal All hydrologic models rely on some form of overland flow simulation to define drainage courses and watershed structure (Garbrecht and Martz, 2000). To create a fully connected and fully labelled drainage network and watershed partition, water outflow at every grid cell of the DEM needs to be routed to an outlet on the border of the DEM. Nevertheless, the frequent presence of surface depressions in the DEM prevents simulated water flow from draining into outlets, resulting in disconnected stream-flow patterns and spurious interior sub-watersheds pouring into these depressions. Due to the undesirable results, surface depressions in DEMs are treated as nuisance features in hydrologic modelling. The common practice is to locate and remove surface depressions in the DEM at the very first step of hydrologic analysis. The removal of pits can be achieved by a variety of methods (Reuter et al., 2008) that can be grouped into three

345

20

40

60

80

90 95 100

L. Poggio, P. Soille / Catena 87 (2011) 341–350

0

SRTM 0

100

200

GDEM SRTM

GDEM 300

400

500

Fig. 4. Total variance for N = 500 iterations. For all considered DEMs, 100 iterations explain more than 90% of the variance obtained for 500 iterations and considered as the total variation.

main categories depending on how the original elevation values are modified: 1. incremental methods suppress pits by incrementing their elevation value until their pour point is reached. Examples of incremental methods are the depression filling (Jenson and Domingue, 1988) and morphological pit filling in Soille and Ansoult(1990); 2. decremental methods whereby values along a path starting from the bottom of the pit and reaching a pixel of lower elevation value are decremented by setting their elevation value to that of the bottom of the pit. The two main decremental methods are known as the breaching (Martz and Garbrecht, 1998) and carving (Soille et al., 2003) methods; 3. hybrid methods combine incremental and decremental methods. An optimal combination of pit filling and carving is proposed in Soille(2004b), whilst a combination of depression filling and breaching is put forward in Lindsay and Creed(2005). To study the impact of pit removal method when searching for the most probable river network, a method was selected from each category. The selected methods rely on concepts of morphological image analysis (Soille, 2004c). Indeed, morphological image analysis offers a common framework for the development of effective and efficient implementations coping with the presence of nested pits and natural depressions (Soille, 2007). The methods used are briefly described hereafter: 1. Pit filling (F1). Spurious pits are removed by considering a marker function flagging all non-spurious pits (typically those connected to the border of the considered DEM) and then performing a morphological reconstruction by erosion (Soille, 2004c) of the input DEM from this marker function. This technique originally proposed in Soille and Ansoult(1990) is very effective and fast for it removes nested pits without any extra computation and can be implemented with efficient algorithms (Soille and Gratin, 1994); 2. Carving (F2). The carving method (Soille et al., 2003) relies on a flooding simulation. The pits are not filled, but the terrain is carved to make pits flowing further down, i.e. carving decreases the elevation of pixels occurring along a path linking the bottom of the pit to the nearest pixel of lower elevation. An efficient implementation of carving is detailed in Soille(2004a).

3. Optimal hybrid (F3). The optimal approach combines morphological pit filling and carving (Soille, 2004b) in order to minimise the differences between the original DEM and the elaborated one. In this combined approach pits are filled up to a certain level and then carving proceeds from that level. The level is set so as to either minimise (i) the sum of the height differences between the input and the output depressionless DEM or (ii) the number of modified pixels. In this study the first option was used. Modelling at landscape scale does not always require pit removal (Hancock, 2008; Temme et al., 2006). However, the removal of spurious pits is especially necessary in this case, as the DEMs were generated by Gaussian simulation that introduced noise in the data. 3.3. River network extraction The numerous methods that have been developed for the calculations of flow directions can be broadly divided in single (non dispersive) and multiple (dispersive) neighbour flow (Gruber and Peckham, 2008). Examples of single neighbour flow algorithms are (a) deterministic D8 (O'Callaghan and Mark, 1984), (b) Rho8 (Fairfield and Leymarie, 1991), a stochastic extension of the D8, (c) D8-LTD (Orlandini et al., 2003), or (d) a global search method that maximises the use of all information stored in the dataset to define the flow path (Paik, 2008). Examples of multiple neighbour flow algorithms are (a) MFD approaches that handle divergent flow and divide the flow coming out of a single cell to all lower neighbours (e.g.: Freeman, 1991; Quinn et al., 1991; Holmgren, 1994) and (b) D∞ in which an infinite number of flow directions is allowed as the flow direction angle is continuously defined (Tarboton, 1997). Some of the dispersive algorithms can improve the estimation of flow accumulation and are well suited for estimating the break in slope of the log–log relationship between contributing drainage area and slope (Vogt et al., 2003). However multi-flow directions cannot define specific flow paths and thus are not suitable for processes where channel definition is important (Orlandini et al., 2003; Paik, 2008; Poggio and Soille, 2009). McMaster et al.(2002) showed that the multiple flow direction methods do not provide a significant improvement in network detection in steep terrains with well-defined channels. Submeter landscape features are one of the main reasons for dispersive

346

L. Poggio, P. Soille / Catena 87 (2011) 341–350

flow (Wilson et al., 2007) and they can direct the run-off in cells not on the steepest descent path (Endreny and Wood, 2003). However the medium resolution datasets used here cannot resolve such microtopography features. Moussa (2009) concluded that for certain applications, such as stream networks extraction, the use of nondispersive methods may appear preferable as more consistent with the physical definition of upstream drainage area. For all these reasons and given the rather steep terrains of the test area, the D8 method was used to calculate flow direction and flow accumulation. The river network was then defined as all cells with a flow accumulation value higher than 100 pixels. This threshold was set according to the morphological characteristics of the considered region (Colombo et al., 2007; Vogt et al., 2003). The river networks obtained for all iterations were added to each other to calculate in how many realisations a cell was part of the river network. The value of each cell of this channel probability map corresponds to the probability that a river network passed through the cell.

4. Most probable river network In order to identify the most probable network, the following datasets were considered: i) the inverse of the channel probability map (i.e. not-channel probability, inv in the following text), calculated on a cell by cell basis as: not−channel probability = 1−channel probability value

This surface represents the probability of a grid cell not being part of a channel (Fig. 5a). The not-channel probability map was considered as a pseudo-DEM that was subsequently treated as if it were a DEM. That is, flow paths through the areas with lowest not-channel probability maps were

Network High=1 Low=0

(a) Extracted river network plotted over the probability map. Here an example of IMAP_DTM and Inv_F3

(b) GDEM Elev

(c) GDEM Itermax

(d) GDEM Mean

(f) IMAP _DTM Elev (g) IMAP _DTM Iter- (h) IMAP _DTM max

ð1Þ

Mean

(e) GDEM Inv _F3

(i) IMAP _DTM Inv _F3

Fig. 5. Maps of the most probable river networks. The black square in (a) delineates the area detailed in following figures.

L. Poggio, P. Soille / Catena 87 (2011) 341–350

calculated to obtain a network maximising the channel probabilities according to the simulations of the data set errors. The idea of applying hydrological concepts to non-elevation data was also considered in Soille and Grazzini (2007) to detect arborescent structures in optical image data. However, in this latter case, the pseudo-DEM was calculated by computing a cost function on the input optical image. ii) the average (mean in the following text) of all the 100 simulated DEMs obtained with the model realisations. These datasets are considered to represent a summary of the modelled data. iii) the river network extracted from the original dataset without error modelling (elev in the following text); iv) the river network with the highest overall probability amongst all the drainage networks that were extracted in the 100 realisations (itermax in the following text).

Table 4 Overall probability values of the river networks normalised by the number of pixels (0–1 range).

DTK InvF3 Mean IterMax Elev

GDEM

SRTM

IMAP_DSM

IMAP_DTM

0.30 0.26 0.19 0.17 0.18

0.28 0.24 0.21 0.16 0.17

0.36 0.32 0.29 0.24 0.33

0.40 0.36 0.33 0.28 0.29

method provided results very similar to the ones obtained for the network extracted from the mean dataset. The river networks extracted from the different options tested were plotted on the channel probability distribution map (Fig. 5). The networks extracted from the original DEM and from itermax datasets show numerous artefacts in the junctions, where multiple parallel lines are derived. Similar problems are visible also in the main large segments. The networks do not always lie in the areas with higher probability values. In the networks derived from the mean dataset the features lie closer to the regions with higher values in the channel probability distribution map. However, in case of segments of lower Strahler order, the uncertainty is higher. The networks derived from the inverse of the channel probability maps and filled with optimal hybrid method follow closely the areas with higher probability values and present less parallel lines. 5. Comparison to reference network The overall probability values of each extracted network were compared with the overall probability values of the DTK reference network (Table 4). For all four considered DEMs, the overall probability value of the reference network is higher than that of extracted networks. This suggests that the probability maps are consistent with the reference network (Fig. 7). Drainage networks are difficult objects to compare quantitatively, and there is no standard method to assess the quality of extracted networks (Molloy and Stepinski, 2007). In order to assess in a more quantitative way the spatial match of the extracted networks with the selected reference, a distance map from the DTK reference network was produced using the algorithm proposed by Saito and Toriwaki

F1: Pit filling F2: Carving F3: Optimal hybrid

SRTM GDEM IMAP_DTM IMAP_DSM

0.35 0.30 0.25 0.20 0.15

Probabilities%

0.40

0.45

0.50

In order to remove the spurious pits the three methods previously described (Section 3.2) were used for the inverse of the channel probability map (i) above and the mean dataset (ii). The overall probability of the river network was calculated by summing the values of the channel probability map for all the cells belonging to the network considered. The most probable network was defined as the one maximising the overall probability value. The results of the overall probability for the considered datasets are presented in Fig. 6 where the values are normalised by the number of pixels in the river network. The overall probability of the networks extracted from the DEMs after the error modelling (i.e. mean and inv datasets) is significantly higher than the values obtained for the networks extracted from the original dataset (elev) or for the network that maximised the overall probability amongst the different realisations (itermax). The steepest gain is obtained in case of IMAP_DTM and IMAP_DSM, due to the originally higher resolution and thus to the rougher surface in the original dataset. The network extracted from the inverse of the channel probability map for river networks has the highest overall probability values, with significant differences due to the filling algorithms used. By construction the pseudo-DEM is particularly cluttered with spurious pits and the impact of pit removal methods is thus more evident. The use of the optimal hybrid algorithm (F3) obtained the highest values, whilst the use of the pit filling (F1)

347

elev

itermax

mean

inv

Fig. 6. Overall river networks probability of the considered datasets.

348

L. Poggio, P. Soille / Catena 87 (2011) 341–350

(a) Distance Values SRTM GDEM IMAP_DTM IMAP_DSM

80

100

120

140

160

mean distance of river networks from the reference (m)

F1: Pit filling F2: Carving F3: Optimal hybrid

elev

itermax

aver

inv

80

90

100

110

120

SRTM GDEM IMAP_DTM IMAP_DSM

70

Mean distance of river networks from the reference (m)

(b) Inverse of the channel probability maps

inv_F1

inv_F2

inv_F3

Fig. 7. Mean distance of river networks from the reference.

(1994). Fig. 7 presents the distances calculated for the extracted networks, where lower values indicate that the network is closer to the reference source. The pattern shown is similar to the one found for the probability values, with mean and inv significantly closer to the reference network. The itermax and the elev networks had significantly higher values than the networks extracted from the notchannel probability map. The inv_F3 network had the lowest values, thus being the closest to the reference network as shown in Poggio

Table 5 Mean distance values (m) between extracted networks and reference DTK network.

InvF1 InvF2 InvF3 Mean IterMax Elev

GDEM

SRTM

IMAP_DSM

IMAP_DTM

105 98 95 111 143 144

97 90 85 104 146 132

108 102 100 115 140 142

84 76 73 91 128 126

and Soille(2009) where only non-simulated datasets were considered (Table 5). 6. Conclusions The spatial positioning and the geomorphological properties of hydrographic features are important in many environmental and management studies. Digital elevation data are commonly used to model drainage divides, hydrographic channel networks and river attributes such as contributing drainage area, channel sinuosity or tributary junctions. The uncertainty and the errors of derived features are often neglected by users, because of complexity issue and computational loads. The proposed probabilistic approach provided a significant increase in the accuracy of the spatial positioning of the river networks extracted. Furthermore, this study presents an alternative method to extract river networks considering the notchannel probability map as a pseudo-DEM from which a flow map can be calculated through the areas of lowest not-channel values, thus obtaining a network maximising the probabilities of having a river,

L. Poggio, P. Soille / Catena 87 (2011) 341–350

according to the simulations of DEM errors. The channel networks derived from the channel probability maps: i. reached a significant higher probability than the best river network out of the 100 realisations; ii. had a lower variability in lower Strahler orders; iii. presented an improved spatial matching of features, being closer to the reference network used. The influence of the filling method is manifested when considering the inverse channel probability map with more homogeneous morphology especially in flat regions. The datasets treated with the optimal hybrid method provided the most probable river networks that were closer to the reference network. The analysis of the distribution of the probabilities gives information about the uncertainty of derived channels and thus about the confidence of the results. The derived most probable network gives an indication on the location where there are the highest probabilities to find the features on the terrain. The method presented is useful to derive hydrographic features that minimise the uncertainty due to errors in the elevation data source. Acknowledgements The authors are grateful to Intermap Technologies Inc for providing sample of the NEXTMap ® products for an evaluation with scientific and non-commercial purposes. References Anders, N.S., Seijmonsbergen, A.C., Bouten, W., 2009. Modelling channel incision and alpine hillslope development using laser altimetry data. Geomorphology (ISSN: 0169-555X) 1130 (1–2, Sp. Iss. SI), 35–46. doi:10.1016/j.geomorph.2009.03.022. DEC 1. Annoni, A., Luzet, C., Gubler, E., Ihde, J. (Eds.), 2003. Map projections for Europe, volume EUR 2012 0EN. European Commission, DG Joint Research Centre. http://www.ec-gis. org/sdi/publist/pdfs/annoni-etal2003eur.pdf. URL. Burrough, P.A., McDonell, R.A., 2000. Principles of Geographical Information System. Oxford University Press. Cambardella, C.A., Moorman, T.B., Novak, J.M., Parkin, T.B., Karlen, D.L., Turco, R.F., Konopka, A.E., 1994. Field-scale variability of soil properties in Central Iowa soils. Soil Scientific Society American Journal 58, 1501–1511. Chen, Q., Gong, P., Baldocchi, D., Xie, G., 2007. Filtering airborne laser scanning data with morphological methods. Photogrammetric Engineering and Remote Sensing 730 (2), 175–185. Colombo, R., Vogt, J.V., Soille, P., Paracchini, M.L., de Jager, A., 2007. Deriving river networks and catchments at the European scale from medium resolution digital elevation data. Catena 70, 296–305. doi:10.1016/j.catena.2006.10.001. Cressie, N.A.C., 1993. Statistics for Spatial Data. Wiley, New York. Darnell, A.R., Tate, N.J., Brunsdon, C., 2008. Improving user assessment of error implications in digital elevation models. Computers, Environment and Urban Systems 32, 268–277. doi:10.1016/j.compenvurbsys.2008.02.003. Endreny, T.A., Wood, E.F., 2003. Maximizing spatial congruence of observed and DEMdelineated overland flow networks. International Journal of Geographical Information Science 17, 699–713. doi:10.1080/1365881031000135483. EuroDEM, 2008. Sample from EuroDEM (c) EuroGeographics and Bundesamt fuer Kartographie und Geodaesie. Technical report, EuroGeographics. http://www. eurogeographics.org/eng/EuroDEM.asp. URL. EuroGeographics, 2008. EuroDEM user guide. Technical Report, EuroGeographics. April. Fairfield, J., Leymarie, P., 1991. Drainage networks from grid digital elevations models. Water Resources Research 270 (5), 709–717. Fisher, P.F., 1998. Improved modeling of elevation error with geostatistics. GeoInformatica 2, 215–233. Freeman, T.G., 1991. Calculating catchment area with divergent flow based on a regular grid. Computer and Geosciences 170 (3), 413–422. Garbrecht, J., Martz, L., 2000. Digital elevation model issues in water resources modeling. In: Maidment, D., Djokic, D. (Eds.), Hydrologic and Hydraulic Modeling Support with Geographic Information Systems. Environmental Systems Research Institute Press, Redlands. Gruber, S., Peckham, S., 2008. Land-surface parameters and objects in hydrology. In: Hengl, T., Reuter, H. (Eds.), Geomorphometry: Concepts, Software, and Applications, Volume 33 of Developments in Soil Science, Chapter 5. Elsevier, pp. 121–140. Hancock, G.R., 2008. The impact of depression removal on catchment geomorphology, soil erosion and landscape evolution. Earth Surface Processes and Landforms (ISSN: 0197-9337) 330 (3), 459–474. doi:10.1002/esp.1598. Hengl, T., Heuvelink, G.B.M., van Loon, E.E., 2010. On the uncertainty of stream networks derived from elevation data: the error propagation approach [in review]. Hydrol. Earth Syst. Sci. Discuss. 7, 767–799.

349

Holmes, K.W., Chadwick, O.A., Kyriakidis, P.C., 2000. Error in a USGS 30-meter digital elevation model and its impact on terrain modeling. Journal of Hydrology 233, 154–173. Holmgren, P., 1994. Multiple flow direction algorithms for runoff modelling in grid based elevation models: an empirical evaluation. Hydrological Processes 8, 327–334. doi:10.1002/hyp.3360080405. Hovenbitzer, M., 2008. The European DEM (EuroDEM). The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XXXVII. ISPRS, pp. 1853–1856. ISPRS Congress Beijing 2008. Hunter, G., Goodchild, M., 1997. Modeling the uncertainty of slope and aspect estimates derived from spatial databases. Geographical Analysis 29, 35–49. Jarvis, A., Reuter, H.I., Nelson, A., Guevara, E., 2006. Hole-filled seamless SRTM data V3. Technical Report, International Centre for Tropical Agriculture (CIAT). Jenson, S.K., Domingue, J.O., 1988. Extracting topographic structure from digital elevation data for geographic information system analysis. Photogrammetric Engineering and Remote Sensing 540 (11), 1593–1600. Keys, R.G., 1981. Cubic convolution interpolation for digital image processing. IEEE Transactions in Acoustics, Speech and Signal Processing 290 (6), 1153–1160. Lin, W.-T., Chou, W.-C., Lin, C.-Y., Huang, P.H., Tsai, J.-S., 2006. Automated suitable drainage network extraction from digital elevation models in Taiwans upstream watersheds. Hydrological Processes 20 (289306). doi:10.1002/hyp.5911. Lindsay, J., 2006. Sensitivity of channel mapping techniques to uncertainty in digital elevation data. International Journal of Geographical Information Science 20 (6), 669–692. doi:10.1080/13658810600661433. Lindsay, J.B., Creed, I.F., 2005. Removal of artifact depressions from digital elevation models: towards a minimum impact approach. Hydrological Processes 19, 3113–3126. doi:10.1002/hyp.5835. Martz, L.W., Garbrecht, J., 1998. The treatment of flat areas and depressions in automated drainage analysis of raster digital elevation models. Hydrological Processes 12, 843–855. McMaster, K.J., et al., 2002. Effects of digital elevation model resolution on derived stream network positions. Water Resources Research (ISSN: 0043-1397) 380 (4), 13. doi:10.1029/2000WR000150. Molloy, I., Stepinski, T.F., 2007. Automatic mapping of valley networks on Mars. Computers & Geosciences (ISSN: 0098-3004) 330 (6), 728–738. doi:10.1016/j. cageo.2006.09.009. Moussa, R., 2009. Definition of new equivalent indices of Horton–Strahler ratios for the derivation of the geomorphological instantaneous unit hydrograph. Water Resources Research 45, W09406. doi:10.1029/2008WR007330. Murphy, P.N.C., Ogilvie, J., Meng, F.-R., Arp, P., 2008. Stream network modelling using lidar and photogrammetric digital elevation models: a comparison and field verification. Hydrological Processes (ISSN: 1099-1085) 220 (12), 1747–1754http:// dx.doi.org/10.1002/hyp.6770. doi:10.1002/hyp.6770 URL. O'Callaghan, J.F., Mark, D.M., 1984. The extraction of drainage networks from digital elevation data. Computer Vision Graphics Image Processing 28 (3), 323–344. Oksanen, J., Sarjakoski, T., 2005. Error propagation of DEM-based surface derivatives. Computers and Geosciences 310 (8), 1015–1027. doi:10.1016/j.cageo.2005.02.014. Oksanen, J., Sarjakoski, T., 2006. Uncovering the statistical and spatial characteristics of fine toposcale DEM error. International Journal of Geographical Information Science 200 (4), 345–369. doi:10.1080/13658810500433891. Orlandini, S., Moretti, G., Franchini, M., Aldighieri, B., Testa, B., 2003. Path-based methods for the determination of nondispersive drainage directions in grid-based digital elevation models. Water Resources Research 390 (6), 1144. Paik, K., 2008. Global search algorithm for nondispersive flow path extraction. Journal of Geophysical Research 113, F04001. doi:10.1029/2007JF000964. Park, S.K., Schowengerdt, R.A., 1983. Image reconstruction by parametric cubic convolution. Computer Vision, Graphics, and Image Processing 23, 258–272. Passalacqua, P., Trung, T.-D., Foufoula-Georgiou, E., Sapiro, G., Dietrich, W.E., 2010. A geometric framework for channel network extraction from lidar: nonlinear diffusion and geodesic paths. Journal of Geophysical Research—Earth Surface (ISSN: 0148-0227) 115. doi:10.1029/2009JF001254 JAN 7. Pebesma, Edzer J., 2004. Multivariable geostatistics in S: the gstat package. Computers & Geosciences 30, 683–691. Pike, R.J., 2002. A bibliography of terrain modeling (geomorphometry), the quantitative representation of topography — supplement 4.0. Technical Report OPEN-FILE REPORT 02–465, U.S. Department of the Interior — U.S. Geological Survey. Pirotti, F., Tarolli, P., 2010. Suitability of LiDAR point density and derived landform curvature maps for channel network extraction. Hydrological Processes (ISSN: 0885-6087) 240 (9), 1187–1197. doi:10.1002/hyp.7582. APR 30. Poggio, L., Soille, P., 2008. Quality assessment of hydro-geomorphological features derived from digital terrain models. Number EUR23489EN. European Commission, DG Joint Research Centre. http://mdigest.jrc.ec.europa.eu/soille/poggio-soille2008eur.pdf. Poggio, L., Soille, P., 2009. Influence of spurious pit removal methods from SRTM on river network positioning. In: Purves, R., Gruber, S., Straumann, R., Hengl, T. (Eds.), Proceedings of Geomorphometry 2009. University of Zurich, Zurich. http:// geomorphometry.org/poggiosoille2009geomorphometry. Quinn, P., Beven, K., Chevallier, P., Planchon, O., 1991. The prediction of hillslope flow paths for distributed hydrological modeling using digital terrain models. Hydrological Processes 5, 59–80. Reuter, H., Hengl, T., Gessler, P., Soille, P., 2008. Preparation of DEMs for geomorphometric analysis. In: Hengl, T., Reuter, H. (Eds.), Geomorphometry: Concepts, Software, and Applications, Volume 33 of Developments in Soil Science, Chapter 4. Elsevier, pp. 87–120. Ribeiro, P.J., Diggle, P.J., 2001. geoR: a package for geostatistical analysis. R-NEWS 10 (2), 14–18. http://CRAN.R-project.org/doc/Rnews/. Rodríguez, E., Morris, C., Belz, J., 2006. A global assessment of the SRTM performance. Photogrammetric Engineering and Remote Sensing 720 (3), 249–260.

350

L. Poggio, P. Soille / Catena 87 (2011) 341–350

Saito, T., Toriwaki, J.-I., 1994. New algorithms for euclidean distance transformation of an n-dimensional digitized picture with applications. Pattern Recognition 270 (11), 1551–1565. doi:10.1016/0031-3203(94)90133-3. Soille, P., 2004a. Morphological carving. Pattern Recognition Letters 25, 543–550. doi:10.1016/j.patrec.2003.12.007. Soille, P., 2004b. Optimal removal of spurious pits in grid digital elevation models. Water Resources Research 40 (12), W12509. doi:10.1029/2004WR003060. Soille, P., 2004c. Morphological Image Analysis: Principles and Applications. SpringerVerlag, Berlin and New York. corrected 2nd printing of the 2nd edition edition. Soille, P., 2007. From mathematical morphology to morphological terrain features. In: Peckham, R., Jordan, G. (Eds.), Digital Terrain Modelling. Springer, Berlin, pp. 45–66. doi:10.1007/978-3-540-36731-4-2. Soille, P., Ansoult, M., 1990. Automated basin delineation from digital elevation models using mathematical morphology. Signal Processing 20, 171–182. doi:10.1016/ 0165-1684(90)90127-K. Soille, P., Gratin, C., 1994. An efficient algorithm for drainage networks extraction on DEMs. Journal of Visual Communication and Image Representation 50 (2), 181–189. doi:10.1006/jvci.1994.1017. Soille, P., Grazzini, J., 2007. Extraction of river networks from satellite images by combining mathematical morphology and hydrology. Lecture Notes in Computer Science 4673, 636–644. doi:10.1007/978-3-540-74272-2_79. Soille, P., Vogt, J.V., Colombo, R., 2003. Carving and adaptative drainage enforcement of grid digital elevation models. Water Resources Research 39 (12), 1366. doi:10.1029/2002WR001879. Tarboton, D.G., 1997. A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research 33, 309–319. Tarboton, D.G., 2003. Terrain analysis using digital elevation models in hydrology. 23rd ESRI International Users Conference, San Diego, California. July 7–11.

Temme, A.J.A.M., Schoorl, J.M., Veldkamp, A., 2006. Algorithm for dealing with depressions in dynamic landscape evolution models. Computers & Geosciences (ISSN: 0098-3004) 320 (4), 452–461. doi:10.1016/j.cageo.2005.08.001. Temme, A.J.A.M., Heuvelink, G.B.M., Schoorl, J.M., Claessens, L., 2008. Geostatistical simulation and error propagation in geomorphometry. In: Hengl, T., Reuter, H. (Eds.), Geomorphometry: Concepts, Software, and Applications, Volume 33 of Developments in Soil Science, Chapter 5. Elsevier, pp. 121–140. Vogt, J.V., Colombo, R., Bertolo, F., 2003. Deriving drainage networks and catchment boundaries: a new method combining digital elevation data and environmental characteristics. Geomorphology 53, 281–298. doi:10.1016/S0169-555X(02)00319-7. Vogt, J., Soille, P., de Jager, A., Rimavičiute, E., Mehl, S., Foisneau, S., Bódis, K., Dusart, J., Paracchini, M., Haastrup, P., Bamps, C., 2007. A pan-European river and catchment database. Number EUR 22920 EN. European Publications Office. doi:10.2788/35907. Wechsler, S.P., 2007. Uncertainties associated with digital elevation models for hydrologic applications: a review. Hydrology and Earth System Sciences 110 (4), 1481–1500. www.hydrol-earth-syst-sci.net/11/1481/2007/. Wechsler, S.P., Kroll, C.N., 2006. Quantifying DEM uncertainty and its effect on topographic parameters. Photogrammetric Engineering & Remote Sensing 720 (9), 1081–1090. Wilson, J.P., Gallant, J.C., 2000. Terrain Analysis — Principles and Applications. Wiley, New York. Wilson, J.P., Lam, C.S., Deng, Y., 2007. Comparison of the performance of flow-routing algorithms used in GIS-based hydrologic analysis. Hydrological Processes 21, 1026–1044. doi:10.1002/hyp.6277. Wise, S., 1998. The effect of GIS interpolation errors on the use of digital elevation models in geomorphology. In: Lane, S.N., Richards, K.S., Chandler, J.H. (Eds.), Landform Monitoring, Modelling and Analysis. John Wiley and Sons, pp. 139–164. Zhang, K., Chen, S.-C., Whitman, D., Shyu, M.-L., Yan, J., Zhang, C., 2003. A progressive morphological filter for removing nonground measurements form airborne LIDAR data 410 (4), 872–882. doi:10.1109/TGRS.2003.810682.