THALWEG DETECTION FOR RIVER NETWORK CARTOGRAPHY IN FOREST FROM HIGH-RESOLUTION LIDAR DATA

Abstract. This paper addresses the problem of extracting the drainage network in forested areas. A precise description of the drainage network including intermittent streams is required for the planning of logging operations and environmental conservation. LiDAR provides now high-resolution point clouds from which the terrain is modelled and the drainage extracted but it also brings some challenges for traditional approaches. First, the raster DTM is interpolated from LiDAR ground points and has to be split in tiles for processing, adding approximations. Second, drainage enforcement techniques alter the terrain and rely on parameters difficult to fix and limiting the optimisation of the process. In this context, we discuss a new approach aiming at: (1) Designing a data structure to model the terrain with a Triangulated Irregular Network in order to avoid interpolation. This data structure must enable the distribution of data and processes across several nodes in Big data architectures and eventually, the processing of complete watersheds with no tiling. (2) Modelling the river network through thalwegs and avoiding the filling and breaching operations. Thalweg detection is more robust, removing the need for filling and breaching. However, it yields a very dense network requiring a simplification step. Combining this model and the architecture will enable the design and modelling of a new tool for river network computation directly from LiDAR ground points. In this paper, we mainly discuss the second point and propose to model the drainage by a network of thalwegs computed from the terrain. Thalwegs are extracted from the surface network, a topological structure formed of peaks, pits and saddles as vertices and ridges and thalwegs as vertices. We present preliminary results comparing the thalweg network and the drainage network.



INTRODUCTION
In forest management, knowledge of the hydrology is required for planning logging operations and designing haul routes.Traditionally, hydrological maps were drawn from photointerpretation on aerial photographies and from digital terrain models.However, the resolution and the quality of these maps were low because most of the details is hidden by the canopy.Hence, a survey of the logging area is performed before the operations start to locate streams and wetlands.However, such surveys are long, fastidious and, depending on the time of year and weather conditions, may overlook some ephemeral or intermittent streams or some humid areas.Now, LiDAR provides datasets at a much higher resolution and can see the floor under the canopy.As a consequence, many mapping agencies or governmental organisations have launched programmes to collect high-density point clouds and build highresolution DTM from which the hydrology is computed.The principle of the approach consists in interpolating a grid DTMfrom the points on the ground, usually at a resolution around one metre.The river network is then extracted by traditional flow accumulation techniques.
However, LiDAR point clouds are massive, and a single watershed amounts to several terabytes of data.Current software and approaches are not designed to handle such volumes at once.The terrain is split into tiles but a watershed may be spread over several tiles.The DTM is interpolated from scattered ground points but their density varies according to the type of terrain, the thickness of the canopy and the quality of the classification.Furthermore, network computation methods often rely on preprocessing * Corresponding author techniques enforcing the drainage by altering the terrain.Their parameters may depend on more or less arbitrary factors and require validation by an expert.
In order to move towards a fully automated processing of LiDAR point clouds, we need to limit manual operations and avoid arbitrary parameters which affect the reliability of the result.We also need to develop data structures that optimise data storage and access and to preserve the quality of original measurements.The objective is to be able to work on complete watersheds without tiling and to avoid or linit interpolation and drainage enforcement operations that modify the terrain model.
In this paper, we discuss the need for a change of paradigm and propose (1) to model the terrain with a TIN directly from Li-DAR points and (2) to identify rivers from thalwegs in order to avoid limitations of flow accumulation approaches.The next section presents the current approaches with its limitations.Section 3 more specifically addresses thalweg computation and presents preliminary results comparing thalweg network and drainage network.The last section concludes and presents on-going work and future developments.

The processing chain
Current approaches rely on the extraction of the river network from a raster DTM.Since LiDAR points are irregularly distributed, the DTM is computed by interpolation.However, the amount of points is too large for current software to process them all at once and they have to be tiled.LiDAR point classification and DTM interpolation are then performed in each tiles separately.Computation of the river network is done in a separate stage for each watershed.In the case a watershed spreads over several tiles, these tiles have to be assembled, sometimes requiring some adjustment on the overlapping areas between tiles.If the watershed is too big to be processed at once, tiles are processed separately and sub-networks are connected afterward.
Flow accumulation methods simulate a water flow over a terrain.Rivers are defined by pixels whose flow accumulation is above a given threshold.While the principle is simple, the method requires some preprocessing to enforce the drainage.The processing chain is summarised in Figure 1.Errors and approximations mainly come from the steps in red on the figure where the terrain model is modified.The quality of the terrain model also depends on the quality of the ground point classification algorithm (in yellow).This step is still an open issue, especially in forest (Meng et al., 2010), but will not be addressed in this paper.

DTM construction
Traditionally, terrain analysis is conducted on raster DTM since terrain data were mainly obtained from aerial photography or satellite images.However, LiDAR ground points are not distributed regularly.Instead of working on an irregular distribution, the raster DTM is interpolated from ground points.There are two main reasons.First, triangulated irregular networks are more complicated to handle, they require an explicit data structure where rasters have an implicit data structure allowing for faster computation.Second, existing processes were designed to handle raster DTM and handling triangulated terrain would require to implement new tools for this.Since the point density is much higher than before, rasters produced from LiDAR points are of a much higher resolution.Rasters built for drainage extraction usually have a resolution of 1 metre while rasters used for road detection can be of a resolution of 0.25 metres (Kiss et al., 2016), compared to previous DTM at 10 or 30 metre resolution (James et al., 2007).The networks are therefore much more detailed than before but the quality of the DTM depends mainly of: • The thickness of the canopy; • The type of terrain; • The interpolation method.
Aerial LiDAR mesures returns from the ground or objects on the ground.Logically, in open areas, all points reach the ground while, if the canopy is thick, most points hit the canopy and few hit the ground.It is also more difficult in steep slopes to classify ground points, leading to more unclassified points.In previous works, (Sherba et al., 2014) obtained a ground point spacing of 0.91 metre in average with a LiDAR emitting 6 pulses per square metre.However, (White et al., 2010) conducted a survey at highdensity (12 pts/m 2 ) in forest with steep slope and a thick canopy and observed that 94% of the points did not reach the ground.Hence, the density of ground points varies considerably on a terrain.(Montealegre et al., 2015) showed that in such areas, the density is below 0.3 pts/m 2 .Hence, the raster DTM homogeneous resolution does not reflect the ground point distribution.In places where few LiDAR points hit the ground, the raster gives a false impression of precision since the DTM at 1 m was built in some areas from samples of points distant of several metres, sometimes more than 10 metres apart.The only piece of information about the quality given to the user is the point density over the whole terrain.
To a lesser extent, the DTM also depends on the interpolation method computing the raster from the ground points.There is no best method here and the quality depends also on the point density and the terrain.Common methods are inverse distance weight, kriging and triangulation followed by a linear interpolation.(Montealegre et al., 2015) obtained best results with the triangulation.(Stereńczak et al., 2016) also found that triangulation yields smaller errors but all methods can achieve similar errors.However, bias and RMSE depend on the slope and the undergrowth height, meaning that the quality of the DTM is not homogeneous.
As mentioned above, interpolation is conducted separately on each tile.When computing the drainage network of a given watershed, only the part of the terrain corresponding to this watershed is processed.If the watershed overlaps two or more tiles, the DTM is built by merging the tiles.Since elevations on overlapping areas do not match exactly (they were computed in two different interpolation processes), a resampling may be done on these areas.In the case the watershed covers too many tiles to be processed at once, the watershed must be split in several pieces which will be treated separately.This problem occurs for example for very large watersheds like the St-Lawrence river in Canada but also for its main tributaries.

River network construction
Extracting the river network is commonly done by computing the flow accumulation (O'Callaghan and Mark, 1984).It consists in simulating a water flow and counting the amount of water that flows through each pixel.In most approaches, water runs along the steepest slope chosen among the eight neighbours of a pixel but some approaches also consider multiple directions.Pixels belonging to a stream are then defined by a flow accumulation threshold.The classification is then vectorised to produce a stream network and compute the stream order.One difficulty is in fixing this threshold.A too high threshold may erase stream heads and miss some connections between streams while a low threshold will lead to too many streams and a difficulty to identify them.But the main issue is the spurious pits.They are points on the terrain lower than their neighbours.As such, there is no direction out of these pits and the flow is interrupted.
In order to correct this problem, several methods have been developed.A first method consists in integrating hydrographic data in the DTM construction process as in the ANUDEM algorithm Figure 2: Drainages obtained after filling (green) and after breaching (blue) on a one-metre resolution DTM (Hutchinson, 1989).However, it requires an a priori knowledge of the hydrography and existing maps (when they exist) lack of data and precision compared to LiDAR data.A second approach is the filtering of the terrain with a low-pass filter or a focal operation which assigns to a pixel the lowest elevation within a window.Filtering removes small pits by smoothing the terrain but it works only for pits smaller than the window.Nonetheless, it is often used to remove noise and enlarge existing channels before applying other techniques.
Three main drainage enforcement techniques are applied to remove spurious pits and avoid flow interruptions: burning, filling and breaching.Burning the DTM (Lindsay, 2016b) means that users carves the terrain in places where they want to direct the flow.It is often used to force the flow according to some hydrographic data.For example, streams are often diverted under forest tracks by culverts.When these culverts are known, the terrain is burnt at these location for streams to cross the tracks.This method requires data of a higher accuracy than the DTM (Lindsay and Dhun, 2015) and the procedure is arbitrary and done according to the user's knowledge.As such, it cannot be considered in an automatic process.
The most popular approach is the filling of depressions (Jenson andDomingue, 1988, Planchon andDarboux, 2002).The terrain inside the depression is raised to the level of the first outflow point, creating a flat area.The method is applied automatically on the terrain to fill either all the pits or only pits below a certain size.Hence it guarantees that streams are not interrupted but it introduces large deformations in the terrain with large flat areas and unnatural streams running in straight lines as shown in Figure 2.
Breaching lowers the terrain in order to create paths from outflow points to lower depressions.Two main approaches are considered: breaches can follow a natural path, connecting one point to a lower neighbour based on a least cost function (Soille et al., 2003) or they can follow a more direct route (Lindsay and Dhun, 2015).The latter leads to paths carved through eminences and less natural patterns but provides more relevant results when breaching across tracks.
Breaching has the benefit of limiting terrain deformations compared to filling.In Figure 2, a large depression has been filled and lead to straight streams where breaching carved the terrain in a more limited area.We also note that depression filling connect streams but breaching leaves them disconnected, showing that stream order and watershed can both be affected by the enforcement method.
Like filling, breaching does not necessarily require parameters but this may lead to connecting the wrong streams together.Often, breaching is constrained by a maximum length and carving depth but these parameters may depend on the type of terrain and the kind of breach to create.For example, in Figure 2, breaching was constrained by a maximum distance.By breaching over a longer distances, the two streams on the top right corner may have been connected.Furthermore, the result is not the same with different breaching techniques and some pits may remain afterward.Hence, filling may still be performed after breaching to remove any remaining pit (Lindsay, 2016a).

Application to large datasets
In order to extract river networks from large datasets, our objective is to provide an automated and robust method.The main limitation comes from the need for enforcement methods: the choice of a method and the setting of parameters depend on the type of terrain and the experience of the user.Drainage enforcement is still required because flow computation relies on the flow direction.The flow direction is defined by the terrain gradient computed for each pixel and enforcement helps in having consistent flow directions.
While breaching seems the most appropriate technique, it cannot guaranty that all streams were detected and it may add streams that do not exist.Enforcement methods also alter the terrain model in a significant way creating an elevation difference of several metres with the original raster DTM (up to 2 metres for the example of Figure 2) and stream location can vary of several metres with the method.
Hence, while LiDAR brings in data with a considerable precision, drainage enforcement approaches significantly modify the terrain and provide inaccurate results in comparison.Main river streams are usually detected correctly but small streams, often ephemeral or intermittent, cannot be mapped accurately.We consider that, because of the need for enforcement, flow accumulation methods lack of robustness at such a high resolution and may require a user's supervision and possibly intervention.Alternative approaches shall be explored in view of a full automation of the process.The points to take into account are: • Alterations on the terrain must be avoided to preserve terrain accuracy and limit deviations; • User-set parameters should be avoided to improve the robustness of the computation.
A first terrain approximation occurs when interpolating the raster DTM.Another issue with raster DTM is the large volume of data imposing a tiling.As mentioned above, in large watersheds, the river network has to be assembled from different tiles which is a source of inconsistencies.In order to avoid interpolation that brings inaccuracies in the DTM and facilitates partitioning in watersheds, a solution is to work on a TIN terrain built directly from LiDAR ground points.Drainage computation algorithms have also been designed for TIN (Freitas et al., 2016) but they follow a similar approach and require drainage enforcement (Silveira and van Oostrum, 2007) which brings significant alterations in some places.
Indeed, the drainage network is closely related to the topography: rivers run towards lower areas and watercourses follow valley thalwegs.Thalwegs on a DTM can be detected by morphometric classification but these approaches still rely on pixel classification.The result is scale-dependent and yields areal features that need to be vectorised.We propose instead to extract directly the set of thalwegs as part of the surface network and to associate it with the river network.In the next section, we present the surface network and compare it with the drainage network.

The surface network
A surface network (Wolf, 2004) is a tripartite graph connecting three independent sets of vertices, the pits, peaks and saddles (or passes) of a terrain by edges representing its ridges and thalwegs (Figure 3).A surface network can be computed on any 2D manifold and observes several topological constraints: 1.A ridge connects a saddle and a peak; 2. A thalweg connects a saddle and a pit; 3. The number of ridges and of thalwegs connected to a saddle are equal and cannot be lower than 2; 4. Vertices obey the Euler-Poincaré formula: Equation ( 1) takes into consideration the multiplicity of a saddle.It is given by the number of ridges (or thalwegs) connecting at this saddle minus one.For example, if three ridges and three thalwegs start at the same saddle, this saddle is of multiplicity 2. This equation applies to surfaces topologically equivalent to the surface of a sphere.Hence, a virtual pit below the terrain is added to close the DTM.In the case the DTM contains holes, a virtual pit needs to be added to each hole to close the surface (Cortés Murcia et al., 2016).
Surface networks can be computed on both raster and TIN.In the case of TIN, vertices are critical points of the triangulation: pits are points lower than their neighbours, peaks are higher and saddles correspond to a minimum in one direction and a maximum in another.On a raster, two approaches can be used: the raster can be triangulated by splitting each grid square in two triangles and the network is built on the subsequent triangulation (Takahashi et al., 1995) or vertices are calculated by bilinear or bicubic interpolation from the raster (Schneider, 2003).In the second case, they do not necessarily correspond to nodes on the raster grid.Such a definition of peaks, pits and saddles lacks of robustness in the case two neighbouring points have the same elevation.A solution proposed by (Takahashi et al., 1995) is to introduce a lexicographical order between the points: if two points have the same elevation, the one on the right is considered higher or, if both have the same x coordinate, the one with the highest y is higher.This simple rule is effective on natural terrains since completely flat areas hardly occur.
Computation of ridges and thalwegs comes in a second stage.A ridge is a line starting from a saddle and moving upward along the steepest slope to a peak.Respectively, a thalweg is computed by marching downward from a saddle to a pit.Whether a TIN or a raster, computation of the surface network yields a connected graph.Since a pit is located at the bottom of a depression and outflow points correspond to saddles, pits are usually connected by two or more thalwegs.That means that water courses associated to thalweg lines are not interrupted by pits and no drainage enforcement is performed in the process.Furthermore, since no vectorisation is required, no accumulation threshold is defined in the process.These aspects guarantee the robustness of the approach compared to flow accumulation.
However, surface networks face several limitations and cannot substitute yet the drainage network in hydrological applications.While pits and peaks are of importance in the characterisation of a terrain, most saddles have little significance (Brändli, 1996).
Ridges are seen as equivalent to water divides separating watersheds but they connect to thalwegs at saddle points while streams and drainage divides connect at points of confluence which are not modelled in the network.This is reflected by the fact that several thalwegs connecting to the same pit may partially overlap when coming close to the pit.
Hence, in order to remove overlaps which creates redundancy, we propose to insert confluences as nodes in the graph.Thalwegs are then segmented into streams where a stream is a line connecting a saddle and a pit (the stream is equal to the thalweg), a saddle and a confluence (in the upper part of the thalweg), two confluences (in its middle part) or a confluence and a pit (in its lower part).
In Figure 4, 3 thalwegs have been divided into 6 streams, three upper streams, two middle streams and one lower stream.At this stage, saddle points are kept to maintain the connection with ridges which can later be used to build the drainage divide.The surface network algorithm was implemented in Python using Numpy and the QGIS 2.18 API, using the same approach as (Cortés Murcia et al., 2016) but was applied on a raster DTM instead of a TIN for comparison with existing data.Since the objective of the project is to compute networks on TIN, the raster was triangulated.This approach is also more realistic for fluvially eroded terrains (Brändli, 1996) where there are few depressions.Point elevation is not modified by this operation.
The surface network was computed on the original DTM and on the filtered DTM used to compute the drainage.About twice as many features were found on the original DTM (Table 1).Indeed, filtering removed most noise on the terrain and surface network is very sensible to it.
The thalweg network is much denser than the drainage network (Figure 5).Since no thresholding or selection has been done, no detail has been removed while the drainage depends on an accumulation threshold.However, many thalwegs are not significant streams and the distribution of thalwegs is heterogeneous.In the upper part and lower parts of Figure 5, the two networks are quite similar.They correspond to smooth slopes where all streams run towards the centre of the terrain.The most noticeable difference is that thalwegs are initiated at a higher point than streams because stream heads have been removed by the thresholding.
Most thalwegs are concentrated at some places on the borders of the watershed and along the main stream in the centre.These areas correspond to relatively flat areas located on higher or lower grounds (Figure 6).As such, a lot of saddles and pits are found, generating many thalwegs because of the roughness of the terrain.Variations in elevation between adjacent saddles and pits are around 20 cm, the precision of the LiDAR.Much breaching is done in these areas to connect low points together but few streams are defined because the accumulation is too low.
The results were also compared with sample data collected on the ground.These data are point positions acquired when walking up the streams.Most computed streams fit with field data.The few points for which no stream was computed were located further up the streams and were matched by thalwegs (Figure 7).Hence thalwegs can represent relevant streams that were lost otherwise.These streams mainly correspond to ephemeral or intermittent streams whose flow accumulation is low or are located in humid zones where the position of the stream is difficult to establish.
Finally, the surface network algorithm failed to detect fifteen streams in box 3 of Figure 5, all flowing northward.These streams were detected on the original DTM but missed on the filtered DTM.Indeed, there were not enough saddle points to initiate the thalwegs.The filter created many flat areas of 3 × 3 pixels and, because of the orientation of the slope and the lexicographic rule chosen above, no saddle point was detected in this area.On the original terrain, this rule was relevant because adjacent points hardly have the same elevation but is not appropriate when the terrain contains many flat areas due to some filtering or enforcement techniques.

CONCLUSION AND PERSPECTIVES
Overall, the detection of thalwegs instead of streams can be an interesting solution in place of the stream approach because it does not require any drainage enforcement altering the terrain.But, at this stage, it detects too many thalwegs and a selection has to be done.Surface network simplication methods exist (Bremer et al., 2003, Danovaro et al., 2003, Rana and Morley, 2002) and take into accont topological (to maintain network consistency) and geometrical criteria such as the difference of elevation between vertices or change of slope along edges.The latter was used by (Cortés Murcia et al., 2016) to characterise submarine canyons.In comparison to flow accumulation or other indicators such as slope computed on the terrain, they are not pixel-based and less sensitive to noise.
The algorithm needs to be improved to handle flat or nearly flat areas.In the first case, some saddles may be missed meaning that some streams are missed.The definition of vertices needs to be modified to look on a larger window but also to handle water surfaces, such as a lake or a large stream whose two banks can be delineated, which we did not deal with so far.
In the second case, many details are detected that relate to noise rather than real streams.A simplification method should remove most if not all this noise but, especially in lower ground, these areas can be humid zones.These zones need also to be detected and preserved.Hence it is more relevant to delineate them and integrate them in the data structure instead of computing a stream whose location is fuzzy.
Finally, we mentioned that the large volume of data imposes tiling of the terrain.Current systems cannot handle such volumes and a new architecture needs to be defined.We suggested using a TIN data structure since the DTM can be built directly from LiDAR ground points, avoiding an interpolation step.TIN is also more adaptive and makes it easier to extract a part of a terrain, such as a watershed to work with it.However, TIN are heavier to manipulate and geospatial big data architectures must be considered to distribute the workload.Few solutions exist that guarrantee scalability.Some promising work has been done for vector data and topological data structure (Engélinus and Badard, 2018) but more specific architectures handling TIN need to be investigated.

Figure 1 :
Figure 1: River network extraction through DTM interpolation and flow accumulation

Figure 4 :
Figure 4: Segmentation of 3 thalwegs connected to a same pit into streams

Figure 5 :
Figure 5: DEM with drainage network (blue) superimposed on thalweg network (red) 3.2 Preliminary results and analysis The surface network was compared with a drainage network computed from a raster DTM interpolated from LiDAR data points shown in Figure 5.The one-metre resolution DTM was provided by the Ministry of Forests, Wildlife and Parks of Québec.It represents a watershed in the Montmorency forest in Québec, Canada of an area of 1.05 km 2 .The watershed sink is the rightmost point on the figure.Drainage computation followed the subsequent process:

Figure 7 :
Figure 7: Field points acquired along real streams with thalwegs (red) and streams (blue) from box 1 of Figure 5

Table 1 :
Kok et al., 2007) performed to minimise the number of pits and better align edges with thalwegs (deKok et al., 2007): each quadrangle formed by four grid points is split in two triangles by tracing the diagonal joining the lowest node.Number of vertices and edges in the drainage network.Saddles are counted with their multiplicity.