Ditch network extraction and hydrogeomorphological characterization using LiDAR-derived DTM in wetlands

Wetlands are often described as transitional ecosystems between aquatic and terrestrial ecosystems that perform a wide range of hydrological, biogeochemical, and ecological functions (Brinson 1993; Viaud et al. 2004; Maltby & Barker 2009; Montreuil et al. 2011; Maltby & Acreman 2011; Violle et al. 2011). The detailed and comprehensive mapping of the drainage network which includes streams and all ditches is necessary for a good assessment of wetland functions (Maltby & Acreman 2011; Stratford et al. 2011). Large ditches ensure the rapid entry or exit of water in wetlands. They participate in major transfers between rivers and wetlands and thus play an important role in flood control (Violle et al. 2006). Thus, detection of their presence allows evaluatation of water retention. Small ditches, whose density evolves according to human activities such as agriculture, urbanization, or infrastructure development, play a role in the biogeochemical dynamics (Acreman et al. 2011; Okruszko et al. 2011; Rouquette et al. 2011; Surridge et al. 2012). If the major rivers are usually already mapped in most parts of the world, this is not always the case for small streams, and even less for drainage ditches.

Most of drainage network inventories are made with aerial photography and validated by field observations and measures. They are often incomplete due to access problems, the time needed to conduct the surveys, and the multitude and small size of the elements to be mapped. For example, the French Ministry of Environment has decided to launch the mapping of all water bodies in France from aerial vertical photographs (Ministerial Decree No. 2009-1543 of 11 December 2009). This layer that is being developed by the French National Institute of Geographic and Forest Information (IGN) includes rivers and major drainage and irrigation channels. However, this spatial layer cannot be used for functional assessment of wetlands for two reasons: on the one hand, this map is incomplete because it omits the network of small ditches that form the majority of the drainage network; on the other hand, it lacks spatial accuracy, especially in closed landscapes such as woodland or scrubland. Thus, the use of other types of remote sensing data is envisaged to produce detailed maps of the drainage network of large geographic areas.

Besides aerial photography, a wide range of remotely sensed data can be used to map the drainage network, in both fields of passive or active remote sensing. Very high spatial resolution (VHSR) satellite optical data like Quickbird images can easily discriminate water bodies in the infrared wavelengths (Clément et al. 2008) with better results in winter than in summer (Bailly et al. 2006). VHSR optical remote sensing data have several limitations for the mapping of the drainage network: (1) the spatial resolution of these images remains insufficient to identify small ditches and streams; moreover, the mixed spectral signatures of shade and water cause many classification errors (Sawaya et al. 2003); (2) detection of small dry ditches and/or ditches under tree cover is difficult or impossible; and (3) characterization of ditches and streams is limited to their width.

Instruments such as airborne lasers, Light Detection and Ranging (LiDAR), provide innovative contributions to the detection and mapping of wetlands. LiDAR data have been used to characterize the ground microtopography with centimeter accuracy, including ditches and streams, even under tree cover (James et al. 2007) or herbaceous vegetation (Hopkinson et al. 2005). These active imaging sensors are used to measure the height of the canopy with respect to the ground surface. The LiDAR sensor sends a pulse in the near infrared: the first return only records the position of the first object (top of canopy, roof, etc.) while the last return indicates the height of the last object, which corresponds, according to the data collection configuration and canopy density, to soil or vegetation (Harding et al. 2001). Other returns also record the position of intermediate objects like branches or shrubs. If the height of the ground surface can be directly deduced from the points of the first return, a classification of pulses from the past returns is necessary to discriminate the points characterizing the lower intermediate layers of vegetation to those associated with soils (Axelsson 1999).

DTMs derived from LiDAR data are most often constructed in a given spatial resolution, according to the LIDAR data precision (point density). Different interpolation methods of ground points are used to create the DTMs derived from LiDAR data: Triangulated Irregular Network (TIN) (Hunt & James 2010; Lane & D'Amico 2010; Werbrouck et al. 2011), Nearest Neighbor (NN) (Bater & Coops 2009), kriging (Töyrä & Pietroniro 2005), spline (Cavalli et al. 2008; Maxa & Bolstad 2009), Inverse Distance Weighting (IDW) (Rosso et al. 2006; Liu et al. 2007). A study carried out on an area covered with conifers located on Vancouver Island in Canada showed that the NN interpolation method applied to LiDAR data acquired with a point density of 0.7 pt/m² provides better results compared to the IDW and spline methods (Bater & Coops 2009). Another study conducted on boreal wetlands using LiDAR data with a point density of 0.5 pt/m² highlighted that the kriging method is more efficient than IDW and TIN techniques (Töyrä et al. 2003). According to Lloyd & Atkinson (2002), the higher the point density decreases, the better the kriging method performs compared to the IDW method. Thus, it appears difficult to determine which interpolation method is most relevant from these experiments, especially since LiDAR data and landscape characteristics are different from one study to another (Fisher & Tate, 2006; Liu 2008).

Two main approaches are used to spatially distribute the drainage network from LiDAR data. The first approach is based on stream accumulation algorithms and highlights the potential drainage network. Studies show the benefits of this approach in natural areas to extract mountain torrents (Vianello et al. 2009) or small streams in forest areas (Hunt & James 2010). Although this approach allows a detailed mapping of the drainage network, DTM corrections using a computer-assisted photo-interpretation approach are often necessary to eliminate man-made features (e.g., bridges) that interrupt the continuity of the network (Liu & Zhang 2011). Moreover, this approach detects only the natural network. The artificial drainage network, which is not dependent on natural conditions of water circulation, is not identified. Finally, this approach ignores the spatial variability of the substratum impermeability. The second approach reconstructs the drainage network using form criteria derived from LiDAR data, either indirectly from the DTM (Cazorzi et al. 2013; Sofia et al. 2013) or directly from the point cloud (Bailly et al. 2008). This approach allows accurate detection of the actual drainage network, including small ditches and the historical drainage network (Werbrouck et al. 2011). Direct extraction from the point cloud is interesting because it offers opportunities to characterize landscape features in 3D. However, these approaches need to define beforehand areas where the potential hydrological saturation of soils caused by topography is high. Moreover, this approach was applied in open habitats, pastures, and vineyards, and has still to be evaluated in more forested landscapes.

Most studies are focused on the mapping of the drainage network without channel characterization (Vianello et al. 2009; Liu & Zhang 2011). However, information such as the width or depth of streams and ditches are needed for the functional assessment of wetlands, not only for hydrological functions, but also for biogeochemical and ecological ones.

Unlike pixel-based techniques which only use the layer pixel values, the object-oriented approach allows extraction of objects of interest from images using different parameters and classification of them based on different criteria, such as size, position, or shape (Blaschke et al. 2000). Moreover, each entity or object can easily be handled in vector systems, which facilitates mapping of ditch features and production of statistical indicators. While this approach has proved to be an effective means to extract geometric networks, such as roads from satellite images (Song & Civco 2004; Tian & Chen 2007), it has not been tested, to our knowledge, to extract and characterize drainage networks from remotely sensed images.

In this context, the objectives of this study are to: (1) extract the drainage network from LiDAR data estimating the required accuracy of data in various hydrological and vegetation environments; and (2) characterize the network elements in order to achieve a wetland functional assessment.

For these purposes, several DTMs have been derived from LiDAR data with different point densities and different interpolation methods. Then, the drainage network has been extracted from the most accurate DTM and characterized using an object-based image analysis.

The study was carried out on the Pleine-Fougères Long Term Ecological Research (LTER) network site (http://osur.univ-rennes1.fr/zoneatelier-armorique) located in north-eastern Brittany near the Mont-Saint-Michel Bay, France (Figure 1). Researchers of different disciplines dealing with environmental issues including earth, life, and geographical sciences examined the effects of landscape and climate change on biodiversity and abiotic fluxes in this agricultural site (Baudry et al. 2011). This site contains a broad river floodplain which is an environmentally protected area. In the wetland area, the research focused mainly on the variation in flooding and its impacts on the wetland functions (Marechal et al. 2012).

For this study, three test sites of 25 hectares each were selected in the wetland area to estimate the impact of the presence of water and wooded hedgerows on the accuracy of DTMs derived from the LiDAR data. Both landscape structures and water management are different in these grassland areas: in the Boucey marsh, ditches are often associated with hedges whereas Aucey and Sougéal marshes are open ones. All over the study site, ditches constitute a fine-mesh network managed for water evacuation (network density: 190 m/ha). One part of the Sougéal marsh is artificially flooded during 5 to 7 months of the year (it was flooded during LiDAR data acquisition) which is not the case for the two other marshes as they are only occasionally flooded.

LiDAR data were acquired with a point density of approximately 4 points per square meter on 5 April 2009 during leaf off period (Table 1). The original data files consisted of irregularly distributed laser points (X-, Y-, Z-coordinates) that were classified in ground and vegetation points with TerraScan software. After having visually checked them, we randomly subsampled these data with the CloudCompare software to provide additional files with a point density of approximately 1, 2, and 3 points per square meter. Field campaigns were performed from March to December 2009. Two topographic profiles were recorded with a total station theodolite (Trimble M3) and differential GPS equipment (Trimble GeoXH). One cross-section of 100 m in length is located in the Sougéal marsh, on a flat area (difference in height <0.30 m) while the other one, located in the Boucey marsh, covers a ditch section of 4 m length with a deeper slope (difference in height >0.50 m). All field measurements were integrated into a geographic information system (GIS) database.

Table 2 presents the classification scheme used to characterize the drainage network of the study area from the LiDAR data. This typology is based on morphological and phyto-sociological characteristics of ditches (Lanoë 2008). The drainage network of Aucey and Boucey marshes contains different types of ditches, ranging from small drainage ditches to large channels.

LiDAR data were acquired with a high density of ground points (4 pts/m²). We sought to determine the minimum ground-point density necessary to map the drainage network exhaustively. For this, we sampled the LiDAR point cloud data and observed the impact of the decrease in the density of points on the quality of the detection and characterization of the drainage network. We chose to process the LiDAR data not directly from the point cloud but from the DTM derived from the point cloud. Indeed, methods currently used for extracting the drainage network from raster data that are based on threshold values of the curvature (Tarolli & Dalla Fontana 2009; Pirotti & Tarolli 2010; Sofia et al. 2013), seem to be more appropriate to an operational approach compared to 3D analysis methods that are still exploratory (Bailly et al. 2006). Several interpolation algorithms are used to calculate DTM from LiDAR data. We chose to test four of them to compare their potential for the characterization of the drainage network according to the density of points.

The irregularly dispersed laser points assumed to correspond to the ground were used to derive high-accuracy DTMs with four interpolation methods, NN, IDW, spline, and ordinary kriging. The set of parameters used for these calculations is shown in Tables 3 and 4. For the kriging method, a variogram was calculated for each test area. The variogram analysis showed that each site presents spatial distributions of different elevation values that justify the use of two statistical models (exponential and power). The output grid cell datasets were fixed at 0.50 m. Calculations were run for the four point densities. Mean elevation differences between the field measurements and the LiDAR DTMs were measured by the root mean square error (RMSE).

In a second step, the drainage network mapping was performed using an object-oriented approach with eCognition software from the DTM that presents the lowest RMSE value.

First, the linear elements were highlighted with a linear filter. The convolution filters like Sobel directional or Laplacian filters (Haralick et al. 1987) are commonly used to detect linear features, but without taking into consideration their environment in the image. Steger (1998) proposed a method to extract linear features and their width using a geometric approach that considers the maximum width of linear features to retrieve and contrasts between the lines and their environment (Wessel & Wiedemann 2003). In comparison with convolution filters such as Sobel directional or Laplacian filters, this method has the advantage of taking into account parameters describing the shape but also the context of linear elements. All the parameters and configurations that have been selected in this study are presented in Table 5. The first configuration was used to detect small ditches, while the second was applied to identify large channels. One can observe that the two configurations are similar, except for the parameter ‘line width’. The first parameter is ‘line direction’. An iterative approach has been used to define the orientation from 0 to 180° with a range of 5° to the horizontal axis of the image. This amplitude allows the inclusion of all possible orientations. The second parameter is the ‘minimum length of the line’. It was set at 20 m to exclude all linear objects below this value (mostly noise, such as plowing a furrow in an agricultural plot). The third parameter is the ‘approximate width of the line’: it was set at 1 m for small ditches and 5 m for large ditches. The fourth parameter is the ‘border width’ which is the interface zone or transition between the linear and areal neighbors. This width was set at 2 m, which corresponds approximately to the transition between the ditch and its environment. The fifth parameter is the ‘maximum similarity of line and border’. This parameter defines the minimum difference between the pixel values of the linear and nonlinear neighboring pixels. Values range from 0 (low similarity) to 1 (strong similarity). The value of 0.85 was chosen as some ditches are shallow and have altitudes similar to their environment. The sixth parameter is the ‘minimum variance of the pixels of the line or the border’ which was set at 0.00 because some ditches may show homogeneous altitudes. Finally, the seventh parameter ‘minimum average difference of the pixels of the line over the edge’ can discriminate linear features higher than their neighbors (e.g., a tarred road) or lower than those in their neighborhood (e.g., a ditch). Bright lines are detected with values greater than 0.00, whereas dark lines are extracted with values less than 0.00. Both dark and bright lines could be detected when the value is 0.00. The value was set here at −0.02 because the ditch elevation values are often lower than those of their neighbors.

Second, ditches were extracted from the filtered image in two steps: first, the image was segmented into a series of segments which were composed of a cluster of contiguous and homogeneous pixels using different parameters (Baatz & Schäpe 2000). In our case, the parameter values are 10, 0.1, 0.9, 0.5, and 0.5 for scale, shape, color, compactness, and smoothness, respectively. Then, the segments with similar orientations (±40°) were aggregated to obtain one object per ditch.

Once identified, the drainage network was characterized as a function of depth and width of the ditch into four classes according to the phyto-morphological types defined before. The depth of each ditch was estimated from the DTM by subtracting the minimum altitude of the ditch from its maximum altitude.

This approach was applied on two close areas (test areas) of 25 hectares each (Figure 1) to evaluate the impact of wooded vegetation on DTMs' accuracy: in the first one, located in the Boucey marsh, ditches are often bordered by tree-lines, while there is no wooded hedge in the Aucey marsh. Classification results were validated with a drainage network field map. Results are expressed as percentages of the total length of well-detected network.

RMSE values of the DTMs calculated for the two cross-sections were compared for the four interpolation methods and the four ground-point densities (Table 6,). On the Sougéal marsh, all DTMs have the same global RMSE value (0.12 m) whatever the ground-point density or the interpolation method. On flat herbaceous meadows, RMSE values decrease by less than 0.05 m. In contrast, on flooded areas, where no LiDAR returns occur because of the presence of open water, RMSE values increase up to 0.30 m. On the flow ditch located at the Boucey marsh, RMSE values range from 0.12 to 0.27 m. As expected, the best results were obtained with a higher ground-point density. RMSE values significantly increased with a ground-point density of 1 point per m². DTMs derived using NN or the spline interpolation methods (RMSE: 0.12 m) have a better accuracy than DTMs derived using the IDW (RMSE: 0.15 m) and the kriging methods (RMSE: 0.19 m). These results confirm a previous study showing the superiority of the NN interpolation method compared to the IDW for DTMs derived from high point density datasets (Bater & Coops 2009). The comparison between the DTMs' profiles derived from LiDAR data and the field measurements shows that all DTMs underestimate the ditch depth (Figure 2). This could be explained by the presence of open water and dense herbaceous vegetation within the ditch. Figure 2 also highlights the fact that the DTMs obtained by the NN interpolation method with a ground-point density of 4 points per m² best fit with field measurements. This DTM was therefore selected for the detection of the drainage network.

In addition to the RMSE, which expresses the mean errors, we calculated other statistical indicators to better assess the impact of point density and interpolation method on the accuracy of the microtopography characterization. Figure 3 shows minimum, maximum, median, first, and third quartile errors of the DTMs derived from LiDAR data for the two cross-sections. For the Sougéal marsh cross-section, errors are the same whatever the interpolation method and the point density. For the Boucey marsh cross-section, the median error changes according to the point density. However, the decrease in point density leads to an increase of the maximum error. The largest error is observed for the kriging interpolation method, whatever the point density.

The detection results of the drainage network according to the ground-point density and the interpolation method are similar for both marshes, ranging from 40.5 to 57.9% and from 40.1 to 60.6% on the Boucey marsh and the Aucey marsh, respectively (Figure 4). Overestimation rates are low, between 0.9 and 3.6%, which allows a rapid correction by photo-interpretation. Differences in detection accuracy according to the ground-point density are not very prominent, although a break is observed between 2 and 1 pt/m². A decrease of 75% of the number of LiDAR ground points generates between −11 and −18% of well-detected elements according to the interpolation method considered. The differences between the interpolation methods are increasingly significant as the ground-point density decreases. At a density of 4 pts/m², the maximum difference is 2.4% of the linear network. This difference reaches 8.0% for a density of 1 pt/m². This observation is more pronounced on the Boucey marsh, where having the same mean point density, the number of ground points is lower than on the Aucey marsh due to the presence of trees. The NN and kriging interpolation methods provide the best results, while the spline method is the worst one.

The results of the classification performed with the object-oriented approach show significant accuracy differences between classes (Table 7). Channels and evacuation ditches are better identified than drainage and flow ditches, because the latter are less deep and wider. On the Boucey marsh, accuracy levels are lower because of the presence of hedgerows along the ditches that leads to no LiDAR returns. This is particularly true for drainage ditches, which are less deep and have a higher area covered by trees. The results also point out a very small over-detection rate for all classes. Drainage network maps are presented in Figure 5. Almost all ditches that were not detected by the object-oriented approach have been successfully mapped by photo-interpretation techniques.

Automatic processing of the DTM derived from LiDAR data on Aucey and Boucey marshes allows us to characterize 54.8 and 63.8%, respectively, of the total length of ditches previously detected, which represent 33.0 and 37.0% of the total length of the drainage network.

Characterization of ditches is determined by their width and depth. The ditch depth is calculated from the difference between the point (pixel) altitude value and the highest and lowest point values. The accuracy of the ditch depth depends on the accuracy of its boundary. An approximate ditch boundary generates errors for width and depth estimation (Figure 6). The presence of water and/or dense vegetation in the ditch leads to an underestimation of the depth. An overestimation of the ditch width will result in an overestimation of its depth. Conversely, an underestimation of the ditch width will lead to an underestimation of its depth. Errors can be easily corrected by visual photo-interpretation to produce a detailed and complete drainage network map that is necessary to the functional assessment of wetlands.

Results show that the LiDAR data can be used to characterize wetland ditch morphology. A ground-point density of between 4 and 2 pts/m² can reproduce the shape of the ditch with an average error lesser than 0.14 m. With a density lesser than 2 pts/m², the average error increases significantly to 0.27 m and the ditch depth is underestimated. The results are similar regardless of the interpolation method used. The decrease in ground-point density was conducted by a random sub-sampling. A different distribution of points could be obtained directly by decreasing the average mean point density during the LiDAR data acquisition, which could lead to slightly different results. Moreover, a ground-point density of 4 m or higher could provide additional very detailed information on the morphometric features of ditches.

LiDAR data were processed to detect the drainage network with an automatic approach, both in open and closed landscapes. Then, an object-oriented image analysis was used for ditch characterization. The results demonstrate that this object-based method has the advantage of classifying each ditch individually and is particularly well suited to the study context characterized by regular drainage network geometry. However, the approach used to extract ditches could be further improved. Indeed, the parameters of the linear filter applied in this study were determined through a time-consuming process of manual tuning. In addition, a long computation time involves processing data tile by tile allowing the efficient processing of large images. The identification of the drainage network could be performed using low-pass filter (Cazorzi et al. 2013) or curvature (Sofia et al. 2013) techniques, which are simpler, faster, and capable of handling larger images, compared to the linear filter method.

This automatic approach has been applied to the wetlands of the Pleine-Fougères LTER site (13 km²). A photo-interpretation step is still necessary to obtain a valid and complete drainage network of the entire wetland. Three limiting factors were identified for mapping the drainage network: presence of water even during periods of low water levels, vegetation density, and the accuracy of automatic delineation of ditches.

As part of this 2D mapping process, the decrease of the ground-point density leads to a relative decrease of the detection accuracy of the drainage network. More significant differences could be observed from 3D point cloud data. These 3D data could be used for the characterization of drainage network connectivity/hydraulic circulation in wetlands and, therefore, for wetland functional assessment. Moreover, morphometric features of ditches could be characterized directly from the cloud points, as was demonstrated by Lague et al. (2013) in investigating a 3D river bed from terrestrial LiDAR data.

In conclusion, results obtained in this study show that the drainage network map quality significantly depends primarily on the LIDAR data precision (point density) and, to a lesser extent, on the interpolation method used to create the DTM.

We derived several DTMs from LiDAR data with different point densities, ranging from 4 to 1 point/m² with four interpolation methods. Based on our experiments, we found that differences between LIDAR DTMs and field measurements vary from 0.12 to 0.27 m according to the LiDAR point density and the interpolation method. A minimum precision of 2 points per m² is required to properly map the ditch network in wetlands and thus to evaluate their hydrological functions. The NN interpolation method provides best results and in the fastest computation time, global computing time being less than 10 min for a 25-hectares area including DTM generation and automatic drainage extraction. This approach can detect more than half of the drainage network with a low overestimation rate, both in closed and open landscapes. The second phase characterization of the river system highlights that quality of the geometric characterization of objects is very sensitive to their spatial accuracy of the previously defined boundaries. The drainage network map can be corrected and supplemented by photo-interpretation. Two limitation factors have been identified for the mapping of drainage networks: the presence of open water and of vegetation, even during the leave-off period.

From these results, the mapping of the drainage network is currently being carried out on a 650 km² watershed in the Dordogne region, France, in order to evaluate wetland hydrological functions over this area that are characterized by low water levels in summer and to assess the reproducibility of the approach in different hydrological and landscape conditions.

We gratefully acknowledge Elven Lanoé (INRA) and Aurélien Bellanger (Communauté de communes de la Baie du Mont-Saint-Michel) for field support. The authors also thank Dr Olivier Dauteuil (Geosciences Rennes UMR 6118) for pre-processing GPS data.