r.watershed - GRASS GIS manual (2024)

NAME

r.watershed - Calculates hydrological parameters and RUSLE factors.

KEYWORDS

raster, hydrology, watershed, accumulation, drainage, stream network, stream power index, topographic index

SYNOPSIS

r.watershed

r.watershed --help

r.watershed [-s4mab] elevation=name [depression=name] [flow=name] [disturbed_land=name] [blocking=name] [retention=name] [threshold=integer] [max_slope_length=float] [accumulation=name] [tci=name] [spi=name] [drainage=name] [basin=name] [stream=name] [half_basin=name] [length_slope=name] [slope_steepness=name] [convergence=integer] [memory=memory in MB] [--overwrite] [--help] [--verbose] [--quiet] [--ui]

DESCRIPTION

r.watershed generates a set of maps indicating: 1) flowaccumulation, drainage direction, the location of streams andwatershed basins, and 2) the LS and S factors of the Revised UniversalSoil Loss Equation (RUSLE).

NOTES

Without flag -m set, the entire analysis is run in memorymaintained by the operating system. This can be limiting, but is veryfast. Setting this flag causes the program to manage memory on diskwhich allows very large maps to be processed but is slower.

Flag -s force the module to use single flow direction (SFD, D8)instead of multiple flow direction (MFD). MFD is enabled by default.

By -4 flag the user allow only horizontal and vertical flow ofwater. Stream and slope lengths are approximately the same as outputsfrom default surface flow (allows horizontal, vertical, and diagonalflow of water). This flag will also make the drainage basins lookmore hom*ogeneous.

When -a flag is specified the module will use positive flowaccumulation even for likely underestimates. When this flag is notset, cells with a flow accumulation value that is likely to be anunderestimate are converted to the negative. See below for a detaileddescription of flow accumulation output.

Option convergence specifies convergence factor for MFD. Lowervalues result in higher divergence, flow is more widelydistributed. Higher values result in higher convergence, flow is lesswidely distributed, becoming more similar to SFD.

Option elevation specifies the elevation data on which entireanalysis is based. NULL (nodata) cells are ignored, zero and negativevalues are valid elevation data. Gaps in the elevation map that arelocated within the area of interest must be filled beforehand,e.g. with r.fillnulls, toavoid distortions. The elevation map need not be sink-filled becausethe module uses a least-cost algorithm.

Option depression specifies the optional map of actualdepressions or sinkholes in the landscape that are largeenough to slow and store surface runoff from a storm event. All cellsthat are not NULL and not zero indicate depressions. Water will flowinto but not out of depressions.

Raster flow map specifies amount of overland flow per cell.This map indicates the amount of overland flow units that each cellwill contribute to the watershed basin model. Overland flow unitsrepresent the amount of overland flow each cell contributes to surfaceflow. If omitted, a value of one (1) is assumed.

Raster retention map specifies correction factors per cell forflow distribution. This map indicates the percentage of overland flowunits leaving each cell. Values should be between zero and 100. Ifomitted, a value of 100 is assumed.

Input Raster map or value containing the percent of disturbed land(i.e., croplands, and construction sites) where the raster or inputvalue of 17 equals 17%. If no map or value isgiven, r.watershed assumes no disturbed land. This input isused for the RUSLE calculations.

Option blocking specifies terrain that will block overlandsurface flow. Blocking cells and streams stop the slope length for theRUSLE. All cells that are not NULL and not zero indicate blockingterrain.

Option threshold specifies the minimum size of an exteriorwatershed basin in cells, if no flow map is input, or overland flowunits when a flow map is given. Warning: low threshold values willdramactically increase run time and generate difficult to read basinand half_basin results. This parameter also controls the level ofdetail in the stream segments map.

Value given by max_slope_length option indicates the maximumlength of overland surface flow in meters. If overland flow travelsgreater than the maximum length, the program assumes the maximumlength (it assumes that landscape characteristics not discernible inthe digital elevation model exist that maximize the slope length).This input is used for the RUSLE calculations and is a sensitiveparameter.

Output accumulation map contains the absolute value of theamount of overland flow that traverses each cell. This value will bethe number of upland cells plus one if no overland flow map is given.If the overland flow map is given, the value will be in overland flowunits. Negative numbers indicate that those cells possibly have surfacerunoff from outside of the current geographic region. Thus, any cellswith negative values cannot have their surface runoff and sedimentationyields calculated accurately.

Output tci raster map contains topographic index TCI,computed asln(α / tan(β)) where α is the cumulativeupslope area draining through a point per unit contour length andtan(β) is the local slope angle. The TCI reflects thetendency of water to accumulate at any point in the catchment and thetendency for gravitational forces to move that water downslope (Quinnet al. 1991). This value will be negative if α /tan(β) < 1.

Output spi raster map contains stream power index SPI,computed asα * tan(β) where α is the cumulativeupslope area draining through a point per unit contour length andtan(β) is the local slope angle. The SPI reflects thepower of water flow at any point in the catchment and thetendency for gravitational forces to move that water downslope (Mooreet al. 1991). This value will be negative if α < 0,i.e. for cells with possible surface runoff from outside of the currentgeographic region.

Output drainage raster map contains drainage direction.Provides the "aspect" for each cell measured CCW from East.

Figure: Drainage is 8 directions numbered counter-clockwisestarting from 1 in north-east direction(source)

Multiplying positive values by 45 will give the direction in degreesthat the surface runoff will travel from that cell. The value 0(zero) indicates that the cell is a depression area (defined by thedepression input map). Negative values indicate that surface runoffis leaving the boundaries of the current geographic region. Theabsolute value of these negative cells indicates the direction offlow. For MFD, drainage indicates the direction of maximum flow.

The output basin map contains unique label for each watershedbasin. Each basin will be given a unique positive even integer. Areasalong edges may not be large enough to create an exterior watershedbasin. NULL values indicate that the cell is not part of a completewatershed basin in the current geographic region.

The output stream contains stream segments. Values correspondto the watershed basin values. Can be vectorized after thinning(r.thin) withr.to.vect.

The output half_basin raster map stores each half-basin isgiven a unique value. Watershed basins are divided into left and rightsides. The right-hand side cell of the watershed basin (lookingupstream) are given even values corresponding to the values in basin.The left-hand side cells of the watershed basin are given odd valueswhich are one less than the value of the watershed basin.

The output length_slope raster map stores slope length andsteepness (LS) factor for the Revised Universal Soil Loss Equation(RUSLE). Equations taken from Revised Universal Soil LossEquation for Western Rangelands (Weltz et al. 1987). Since the LSfactor is a small number (usually less than one), the GRASS output mapis of type DCELL.

The output slope_steepness raster map stores slope steepness(S) factor for the Universal Soil Loss Equation (RUSLE). Equationstaken from article entitledRevised Slope Steepness Factor for the Universal SoilLoss Equation (McCool et al. 1987). Since the S factor is a smallnumber (usually less than one), the GRASS output map is of type DCELL.

A^T least-cost search algorithm

r.watershed uses an A^T least-cost search algorithm(see REFERENCES section) designed to minimize the impact of DEM dataerrors. Comparedto r.terraflow, this algorithmprovides more accurate results in areas of low slope as well as DEMsconstructed with techniques that mistake canopy tops as the groundelevation. Kinner et al. (2005), for example, used SRTM and IFSAR DEMsto compare r.watershedagainst r.terraflow results inPanama. r.terraflow was unableto replicate stream locations in the larger valleyswhile r.watershed performed much better. Thus, if forestcanopy exists in valleys, SRTM, IFSAR, and similar data products willcause major errors in r.terraflow stream output. Undersimilar conditions, r.watershed will generatebetter stream and half_basin results. If watersheddivides contain flat to low slope, r.watershed will generatebetter basin resultsthan r.terraflow. (r.terraflowuses the same type of algorithm as ESRI's ArcGIS watershed softwarewhich fails under these conditions.) Also, if watershed dividescontain forest canopy mixed with uncanopied areas using SRTM, IFSAR,and similar data products, r.watershed will generate betterbasin resultsthan r.terraflow. Thealgorithm produces results similar to those obtained when runningr.cost andr.drain on every cell on the raster map.

Multiple flow direction (MFD)

r.watershed offers two methods to calculate surface flow:single flow direction (SFD, D8) and multiple flow direction(MFD). With MFD, water flow is distributed to all neighbouring cellswith lower elevation, using slope towards neighbouring cells as aweighing factor for proportional distribution. The A^Tleast-cost path is always included. As a result, depressions andobstacles are traversed with a graceful flow convergence before theoverflow. The convergence factor causes flow accumulation to convergemore strongly with higher values. The supported range is 1 to 10,recommended is a convergence factor of 5 (Holmgren, 1994). If manysmall sliver basins are created with MFD, setting the convergencefactor to a higher value can reduce the amount of small sliver basins.

In-memory mode and disk swap mode

There are two versions of this program: ram and seg.ram is used by default, seg can be used by settingthe -m flag.

The ram version requires a maximum of 31 MB of RAM for 1million cells. Together with the amount of system memory (RAM)available, this value can be used to estimate whether the currentregion can be processed with the ram version.

The ram version uses virtual memory managed by the operatingsystem to store all the data structures and is faster thanthe seg version; seg uses the GRASS segmentationlibrary which manages data in disk files. seg uses only asmuch system memory (RAM) as specified with the memory option,allowing other processes to operate on the same system, even when thecurrent geographic region is huge.

Due to memory requirements of both programs, it is quite easy to runout of memory when working with huge map regions. If the ramversion runs out of memory and the resolution size of the currentgeographic region cannot be increased, either more memory needs to beadded to the computer, or the swap space size needs to beincreased. If seg runs out of memory, additional disk spaceneeds to be freed up for the program to run.The r.terraflow module wasspecifically designed with huge regions in mind and may be useful hereas an alternative, although disk space requirementsof r.terraflow are several times higher than of seg.

Large regions with many cells

The upper limit of the ram version is 2 billion(2³¹ - 1) cells, whereas the upper limit forthe seg version is 9 billion-billion (2⁶³ - 1 =9.223372e+18)cells.
In some situations, the region size (number of cells) maybe too large for the amount of time or memoryavailable. Running r.watershed may then require use of acoarser resolution. To make the results more closely resemble thefiner terrain data, create a map layer containing the lowest elevationvalues at the coarser resolution. This is done by: 1) Setting thecurrent geographic region equal to the elevation map layerwith g.region, and 2) Usethe r.neighbors orr.resamp.stats command tofind the lowest value for an area equal in size to the desiredresolution. For example, if the resolution of the elevation data is 30meters and the resolution of the geographic regionfor r.watershed will be 90 meters: use the minimum functionfor a 3 by 3 neighborhood. After changing to the resolution atwhich r.watershed will be run, r.watershed should berun using the values from the neighborhood output map layerthat represents the minimum elevation within the region of the coarsercell.

Basin threshold

The minimum size of drainage basins, defined by the thresholdparameter, is only relevant for those watersheds with a single streamhaving at least the threshold of cells flowing into it.(These watersheds are called exterior basins.) Interior drainagebasins contain stream segments below multiple tributaries. Interiordrainage basins can be of any size because the length of an interiorstream segment is determined by the distance between the tributariesflowing into it.

MASK and no data

The r.watershed program does not require the user to have thecurrent geographic region filled with elevation values. Areas withoutelevation data (masked or NULL cells) are ignored. It is NOT necessaryto create a raster map (or raster reclassification)named

for NULL cells. Areas without elevation data willbe treated as if they are off the edge of the region. Such areas willreduce the memory necessary to run the program. Masking outunimportant areas can significantly reduce processing time if thewatersheds of interest occupy a small percentage of the overall area.

Gaps (NULL cells) in the elevation map that are located within thearea of interest will heavily influence the analysis: water will flowinto but not out of these gaps. These gaps must be filled beforehand,e.g. with r.fillnulls.

Zero (0) and negative values will be treated as elevation data (notno_data).

Further processing of output layers

Problem areas, i.e. those parts of a basin with a likely underestimateof flow accumulation, can be easily identified with e.g.

 r.mapcalc "problems = if(flow_acc < 0, basin, null())"

If the region of interest contains such problem areas, and this is notdesired, the computational region must be expanded until the catchmentarea for the region of interest is completely included.

To isolate an individual river network using the output of thismodule, a number of approaches may be considered.

1. Use a resample of the basins catchment raster map as a MASK.
   The equivalent vector map method is similar using v.select or v.overlay.
2. Use the r.cost module with a point in the river as a starting point.
3. Use the v.net.iso module with a node in the river as a starting point.

All individual river networks in the stream segments output can beidentified through their ultimate outlet points. These points are allcells in the stream segments output with negative drainage direction.These points can be used as start pointsfor r.water.outlet orv.net.iso.

To create river mile segmentation from a vectorized streamsmap, trythe v.net.isoor v.lrs.segmentmodules.

The stream segments output can be easily vectorized after thinningwithr.thin. Each stream segment in thevector map will have the value of the associated basin. To isolatesubbasins and streams for a larger basin, a MASK for the larger basincan be created withr.water.outlet. The streamsegments output serves as a guide where to place the outlet point usedas input to r.water.outlet.The basin threshold must have been sufficiently small to isolate astream network and subbasins within the larger basin.

Given that the drainage is 8 directions numberedcounter-clockwise starting from 1 in north-east direction,multiplying the outputby 45 (by 45. to get a double precision floating point rastermap in r.mapcalc) givesthe directions in degrees. For most applications, zeroswhich indicate depressions specified by depressionand negative values which indicate runoff leaving the regionshould be replaced by NULL (null() inr.mapcalc).The following command performs these replacements:

r.mapcalc "drainage_degrees = if(drainage > 0, 45. * drainage, null())"

Alternatively, the user can use the -a flag or later the

function inr.mapcalc if the runoff is leavingthe region.

EXAMPLES

These examples use the Spearfish sample dataset.

Convert r.watershed streams map output to a vector map

If you want a detailed stream network, set the threshold option smallto create lots of catchment basins, as only one stream is presentedper catchment. The

flag preserves the catchmentID as the vector category number.

 r.watershed elev=elevation.dem stream=rwater.stream r.to.vect -v in=rwater.stream out=rwater_stream

Set a different color table for the accumulation map:

 MAP=rwater.accum r.watershed elev=elevation.dem accum=$MAP eval `r.univar -g "$MAP"` stddev_x_2=`echo $stddev | awk '{print $1 * 2}'` stddev_div_2=`echo $stddev | awk '{print $1 / 2}'` r.colors $MAP col=rules << EOF 0% red -$stddev_x_2 red -$stddev yellow -$stddev_div_2 cyan -$mean_of_abs blue 0 white $mean_of_abs blue $stddev_div_2 cyan $stddev yellow $stddev_x_2 red 100% red EOF

Create a more detailed stream map using the accumulation map andconvert it to a vector output map. The accumulation cut-off, andtherefore fractal dimension, is arbitrary; in this example we use themap's mean number of upstream catchment cells (calculated in the aboveexample by r.univar) as thecut-off value. This only works with SFD, not with MFD.

 r.watershed elev=elevation.dem accum=rwater.accum r.mapcalc 'MASK = if(!isnull(elevation.dem))' r.mapcalc "rwater.course = \ if( abs(rwater.accum) > $mean_of_abs, \ abs(rwater.accum), \ null() )" r.colors -g rwater.course col=bcyr g.remove -f type=raster name=MASK # Thinning is required before converting raster lines to vector r.thin in=rwater.course out=rwater.course.Thin r.colors -gn rwater.course.Thin color=grey r.to.vect in=rwater.course.Thin out=rwater_course type=line v.db.dropcolumn map=rwater_course column=label

Create watershed basins map and convert to a vector polygon map

 r.watershed elev=elevation.dem basin=rwater.basin thresh=15000 r.to.vect -s in=rwater.basin out=rwater_basins type=area v.db.dropcolumn map=rwater_basins column=label v.db.renamecolumn map=rwater_basins column=value,catchment

Display output in a nice way

 r.relief map=elevation.dem d.shade shade=elevation.dem.shade color=rwater.basin bright=40 d.vect rwater_course color=orange

REFERENCES

• Ehlschlaeger C. (1989). Using the A^T Search Algorithmto Develop Hydrologic Models from Digital Elevation Data,Proceedings of International Geographic Information Systems (IGIS)Symposium '89, pp 275-281 (Baltimore, MD, 18-19 March 1989).
  URL: http://chuck.ehlschlaeger.info/older/IGIS/paper.html
• Holmgren P. (1994). Multiple flow direction algorithms for runoffmodelling in grid based elevation models: An empirical evaluation.Hydrological Processes Vol 8(4), 327-334.
  DOI: 10.1002/hyp.3360080405
• Kinner D., Mitasova H., Harmon R., Toma L., Stallard R. (2005).GIS-based Stream Network Analysis for The Chagres River Basin,Republic of Panama. The Rio Chagres: A Multidisciplinary Profile ofa Tropical Watershed, R. Harmon (Ed.), Springer/Kluwer, p.83-95.
  URL: http://www4.ncsu.edu/~hmitaso/measwork/panama/panama.html
• McCool et al. (1987). Revised Slope Steepness Factor for the UniversalSoil Loss Equation, Transactions of the ASAE Vol 30(5).
• Metz M., Mitasova H., Harmon R. (2011). Efficient extraction ofdrainage networks from massive, radar-based elevation models with leastcost path search, Hydrol. Earth Syst. Sci. Vol 15, 667-678.
  DOI: 10.5194/hess-15-667-2011
• Moore I.D., Grayson R.B., Ladson A.R. (1991). Digital terrainmodelling: a review of hydrogical, geomorphological, and biologicalapplications, Hydrological Processes, Vol 5(1), 3-30
  DOI: 10.1002/hyp.3360050103
• Quinn P., K. Beven K., Chevallier P., Planchon O. (1991). Theprediction of hillslope flow paths for distributed hydrological modellingusing Digital Elevation Models, Hydrological Processes Vol 5(1),p.59-79.
  DOI: 10.1002/hyp.3360050106
• Weltz M. A., Renard K.G., Simanton J. R. (1987). Revised Universal SoilLoss Equation for Western Rangelands, U.S.A./Mexico Symposium ofStrategies for Classification and Management of Native Vegetation forFood Production In Arid Zones (Tucson, AZ, 12-16 Oct. 1987).

SEE ALSO

g.region,r.cost,r.drain,r.fillnulls,r.flow,r.mask,r.neighbors,r.param.scale,r.resamp.interp,r.terraflow,r.topidx,r.water.outlet,r.stream.extract

AUTHORS

Original version:Charles Ehlschlaeger, U.S. Army Construction Engineering Research Laboratory
Faster sorting algorithm and MFD support:Markus Metz <markus.metz.giswork at gmail.com>
Retention for flow distribution by Andreas Gericke (IGB Berlin)

SOURCE CODE

Available at: r.watershed source code (history)

Latest change: Thursday Jan 26 14:10:26 2023 in commit: cdd84c130cea04b204479e2efdc75c742efc4843

Main index |Raster index |Topics index |Keywords index |Graphical index |Full index

© 2003-2024GRASS Development Team,GRASS GIS 8.3.3dev Reference Manual