Given a terrain, how to draw the stream flow path?

Question

Assuming I have a terrain, as usual the terrain has ridges, creeks and all the characteristics that you can find on a real life map. Water flows from the top of the mountain into lower area, the path that water flows is termed stream flow path.

The terrain is given in terms of triangular irregular network ( TIN), which each point p(x,y) has a z value. How to use this information to construct the stream flow path? What is the physics behind this?

From what I know, steepest descent method can be used to solve this problem. I am thinking about writing my own stream flow algorithm, so I am interested in the theoretical background rather than using the existing tools.

Answer 1

There are different possible implementations, but most procedures will start from a grid and not from a TIN.

The simplest one is most likely the D8 procedure: you calculate the direction where water would be flowing. There are 8 possibilities, the 8 cells that are next to a central grid cell. You can first calculate these directions, than how cells are connected, and finally you can draw the lines). An easy implementation is found in SAGA, it almost reads as pseudocode: http://saga-gis.svn.sourceforge.net/viewvc/saga-gis/trunk/saga-gis/src/modules_terrain_analysis/terrain_analysis/ta_channels/D8_Flow_Analysis.cpp?revision=911&view=markup

Although very easy, this is not very realistic: you will not have a stream starting in every cell. More advanced algorithms usually first close the pits (especially if you have a detailed DEM), then calculate the catchment area per cell, that is the number of cells that contribute water to a particular cell, and then use a threshold to determine whether a stream is present.

SAGA GIS implements a lot of these catchment area methods, you can find a description of them in this manual http://sourceforge.net/settings/mirror_choices?projectname=saga-gis&filename=SAGA%20-%20Documentation/SAGA%20Documents/SagaManual.pdf

It was written for an older version of SAGA GIS, but the description of the algorithms is still rather accurate, and I'll copy it here for quick reference (this is around page 120), since it is open source, you can check the implementation details by looking at the code.

• Deterministic 8 (D8): The classical. Flow goes from the center of a cell to the center of one(and only one) of the surrounding cells. Flow directions are, therefore, restricted to multiples of 45o , which is the main reason for most of the drawbacks of the method. (O’Callaghan & Mark 1984).
• Rho8: Same as above but with an stochastic component that should improve it. The flow direction is determined by a random argument that is dependent on the difference between aspect and the direction of the two adjacent neighbor cells. Not very useful. . . (Fairfield & Leymarie 1991).
• Deterministic infinity (D∞): Flow goes from one cell to two contiguous surrounding cells, thus considering a bidimensional flow and overcoming the drawbacks of the D8 method. (Tarboton 1998).
• Braunschweiger Digitales Reliefmodell: Another Multiple Flow Direction algo- rithm. The flow is split between the surrounding cell whose orientation is nearest to the aspect of the center cell and its two adjacent cells. (Bauer, Bork & Rohdenburg 1985).
• FD8 (found in SAGA simply as Multiple Flow Direction): A D8–derived bidimensional flow routing algorithm. (Quinn et al 1991).
• Kinematic Routing Algorithm (KRA). An unidimensional flow tracing algorithm. Flow behaves like a ball rolling down the DEM, without restricting its position to the center of cells. (Lea 1992).
• Digital Elevation Model Network (DEMON): The most complex one. A bidimensional flow tracing algorithm. Rather time consuming. (Costa-Cabral & Burgess 1994).

Even more models have been added recently:

• Triangular Multiple Flow Direction- Seibert, J. / McGlynn, B. (2007): 'A new triangular multiple flow direction algorithm for computing upslope areas from gridded digital elevation models', Water Resources Research, Vol. 43, W04501. This may be interesting for you because it might also work directly on a TIN
• The Mass-Flux Method (MFM) for the DEM based calculation of flow accumulation as proposed by Gruber and Peckham (2008). Gruber, S., Peckham, S. (2008): Land-Surface Parameters and Objects in Hydrology. In: Hengl, T. and Reuter, H.I. [Eds.]: Geomorphometry: Concepts, Software, Applications. Developments in Soil Science, Elsevier, Bd.33, S.293-308.
• The side algorithm: http://watershed.montana.edu/Hydrology/Home_files/2010WR009296.pdf and his code is on his website as well: http://thomasgrabs.com/side-algorithm/

Answer 2

An original approach is the one proposed in this paper:

Fisher, P., J. Wood, and T. Cheng (2004). Where is Helvellyn? Fuzziness of multiscale landscape morphometry. Transactions of the Institute of British Geographers 29, 106-128.

It proposes a method based on fuzzy and multi scale representation. I am not sure but this method may be the one implemented in LandSerf.

Answer 3

If you have access to Spatial Analyst in ArcGIS, then you have a series of tools to compute stream paths. A full workflow is provided in the ESRI reference, but the typical workflow includes:

1. Convert your TIN to an elevation raster.
2. Calculate the Flow Direction.
3. Fill small Sinks.
4. Calculate Flow Accumulation
5. Using a threshold, choose only the cells with a given amount of flow.
6. Use the Stream to Feature tool to export the streams to a vector shapefile.

Of course, there are numerous academic papers describing different methods, but this method is easy for all who have access to Spatial Analyst.

Answer 4

In grid-based digital elevation models, reliable determinations of the lines of slope are provided by the D8-LTD method:

Orlandini, S., and G. Moretti (2009), Determination of surface flow paths from gridded elevation data, Water Resour. Res., 45(3), W03417, doi: 10.1029/2008WR007099.

Orlandini, S., G. Moretti, M. Franchini, B. Aldighieri, and B. Testa (2003), Path-based methods for the determination of nondispersive drainage directions in grid-based digital elevation models, Water Resour. Res., 39(6), 1144, doi: 10.1029/2002WR001639.

In contour-based digital elevation models, lines of slope can be determined automatically by solving complex topographic structures using the (complex) model described in the following paper:

Moretti, G., and S. Orlandini (2008), Automatic delineation of drainage basins from contour elevation data using skeleton construction techniques, Water Resour. Res., 44(5), W05403, doi: 10.1029/2007WR006309.

Answer 5

Seems that this will be quite the job to write a tool from scratch. ESRI has been at it for decades and they still don't have right.

AutoCAD (Civil 3D) can this this using a TIN. I am not aware of what is going on behind the scenes there but in ArcGIS identifying stream networks is handled via raster analysis.

In a nutshell an input DEM raster (where each cell has X,Y,Z values) is used as input and an algorithm calculates quoting "accumulated flow (as the accumulated weight) of all cells flowing into each downslope cell in the input raster." The product is a raster where each cell has a flow accumulation value. To identify the stream network you then isolate the cells of high flow which are the areas of "concentrated flow". There are other considerations such as optional weight factor, hydrologicaly correct input DEM, etc.

I will just toss in a few ideas: In term of the "mechanics" of such algorithm I suppose it could be quite straight forward; recursively and for each cell, determine the location and elevation of all surrounding cells and based on its elevation add up the number of cells flowing into it. As for TIN you could probably construct a line from two points on each triangle (highest and lowest vertex) then join all these together into a network.