(GRASS Raster Program)
r.terraflow takes as input a raster digital elevation model (DEM) and computes the flow direction raster and the flow accumulation raster, as well as the flooded elevation raster, sink-watershed raster (partition into watersheds around sinks) and tci (topographic convergence index) raster.
r.terraflow computes these rasters using well-known approaches, with the difference that its emphasis is on the computational complexity of the algorithms, rather than on modeling realistic flow. r.terraflow emerged from the necessity of having scalable software able to process efficiently very large terrains. It is based on theoretically optimal algorithms developed in the framework of I/O-efficient algorithms. r.terraflow was designed and optimized especially for massive grids and is able to process terrains which were impractical with similar functions existing in other GIS systems.
Flow directions are computed using either the MFD (Multiple Flow Direction) model or the SFD (Single Flow Direction, or D8) model, illustrated below. Both methods compute downslope flow directions by inspecting the 3-by-3 window around the current cell. The SFD method assigns a unique flow direction towards the steepest downslope neighbor. The MFD method assigns multiple flow directions towards all downslope neighbors.
Flow direction to steepest downslope neighbor (SFD). |
Flow direction to all downslope neighbors (MFD). |
---|
The SFD and the MFD method cannot compute flow directions for cells which have the same height as all their neighbors (flat areas) or cells which do not have downslope neighbors (one-cell pits).
In order to flood the terrain, r.terraflow identifies all sinks and partitions the terrain into sink-watersheds (a sink-watershed contains all the cells that flow into that sink), builds a graph representing the adjacency information of the sink-watersheds, and uses this sink-watershed graph to merge watersheds into each other along their lowest common boundary until all watersheds have a flow path outside the terrain. Flooding produces a sink-less terrain in which every cell has a downslope flow path leading outside the terrain and therefore every cell in the terrain can be assigned SFD/MFD flow directions as above.
Once flow directions are computed for every cell in the terrain, r.terraflow computes flow accumulation by routing water using the flow directions and keeping track of how much water flows through each cell. r.terraflow also computes the tci raster (topographic convergence index, defined as the logarithm of the ratio of flow accumulation and local slope).
For more details on the algorithms see [1,2,3].
r.terraflow [ -sq ] elev=name filled=name direction=name swatershed=name accumulation=name tci=name [d8cut=value] [memory=value] [STREAM_DIR=name] [stats=name]
Alternatively, the user can simply type r.terraflow on the command line and the program will ask for parameter values and flag settings interactively, using the standard GRASS parser interface.
One of the techniques used by r.terraflow is the space-time trade-off. In particular, in order to avoid searches, which are I/O-expensive, r.terraflow computes and works with an augmented elevation raster in which each cell stores relevant information about its 8 neighbors, in total up to 80B per cell. As a result r.terraflow works with intermediate temporary files that may be up to 80N bytes, where N is the number of cells (rows x columns) in the elevation raster (more precisely, 80K bytes, where K is the number of valid (not nodata) cells in the input elevation raster). All this intermediate temporary files are stored in the path specified by STREAM_DIR. Note: STREAM_DIR must contain enough free disk space in order to store up to 2 x 80N bytes.
The internal type used by r.terraflow to store elevations can be defined at compile-time. By default, r.terraflow is compiled to store elevations internally as floats. A version which is compiled to store elevations internally as shorts is available as r.terraflow.short. Other versions can be created by the user if needed.
r.terraflow.short uses less space (up to 60B per cell, up to 60N intermediate file) and therefore is more space and time efficient. r.terraflow is intended for use with floating point raster data (FCELL), and r.terraflow.short with integer raster data (CELL) in which the maximum elevation does not exceed the value of a short SHRT_MAX=32767 (this is not a constraint for any terrain data of the Earth, if elevation is stored in meters).
Both r.terraflow and r.terraflow.short work with input elevation rasters which can be either integer, floating point or double (CELL, FCELL, DCELL). If the input raster contains a value that exceeds the allowed internal range (short for r.terraflow.short, float for r.terraflow), the program exits with a warning message. Otherwise, if all values in the input elevation raster are in range, they will be converted (truncated) to the internal elevation type (short for r.terraflow.short, float for r.terraflow). In this case precision may be lost and artificial flat areas may be created.
For instance, if r.terraflow.short is used with floating point raster data (FCELL or DCELL), the values of the elevation will be truncated as shorts. This may create artificial flat areas, and the outpus of r.terraflow.short may be less realistic than those of r.terraflow on floating point raster data. The outputs of r.terraflow.short and r.terraflow are identical on integer raster data (CELL).