Geografic-Information-System/raster/r.watershed/front/r.watershed.html

<h2>DESCRIPTION</h2>

<em>r.watershed</em> generates a set of maps indicating:
1) flow accumulation, drainage direction, the location of streams and watershed basins, and
2) the LS and S factors of the Revised Universal Soil Loss Equation (RUSLE).

<p>
<!-- Interactive mode not activated in GRASS 7.
<em>r.watershed</em> can be run either interactively or non-interactively.
The interactive version of
<em>r.watershed</em> can prepare inputs to lumped-parameter hydrologic models.
After a verbose interactive session, <em>r.watershed</em> will query the user
for a number of
map layers.  Each map layer's values will be tabulated by watershed basin and sent
to an output file.  This output file is organized to ease data entry into a
lumped-parameter hydrologic model program.  The non-interactive version of
<em>r.watershed</em> cannot create this file.
-->

<h2>OPTIONS</h2>

<dl>
<dt><em>-m</em> 

<dd>Without this flag set, the entire analysis is run in memory
maintained by the operating system.  This can be limiting, but is
very fast.  Setting this flag causes the program to manage memory
on disk which allows very large maps to be processed but is slower.

<dt><em>-s</em> 

<dd>Use single flow direction (SFD, D8) instead of multiple flow direction (MFD).
MFD is enabled by default.

<dt><em>-4</em> 

<dd>Allow only horizontal and vertical flow of water.
Stream and slope lengths are approximately the same as outputs from default
surface flow (allows horizontal, vertical, and diagonal flow of water).
This flag will also make the drainage basins look more homogeneous.

<dt><em>-a</em> 

<dd>Use positive flow accumulation even for likely underestimates. When this
flag is not set, cells with a flow accumulation value that is likely to be
an underestimate are converted to the negative. See below for a detailed
description of flow accumulation output.

<dt><em>memory</em> 

<dd>Maximum amount of memory in MB to be used with -m set. More memory 
speeds up the processes.

<dt><em>convergence</em> 

<dd>Convergence factor for MFD. Lower values result in higher divergence,
flow is more widely distributed. Higher values result in higher convergence, 
flow is less widely distributed, becoming more similar to SFD. 

<dt><em>elevation</em> 

<dd>Input map: Elevation on which entire analysis is based. NULL (nodata) 
cells are ignored, zero and negative values are valid elevation data. 
Gaps in the elevation map that are located within the area of interest 
must be filled beforehand, e.g. with <em>r.fillnulls</em>, to avoid 
distortions.

<dt><em>depression</em> 

<dd>Input map:  Map layer of actual depressions or sinkholes in the
landscape that are large enough to slow and store surface runoff from 
a storm event.  All cells that are not NULL and not zero indicate 
depressions. Water will flow into but not out of depressions.

<dt><em>flow</em> 

<dd>Input map: amount of overland flow per cell.  This map indicates the
amount of overland flow units that each cell will contribute to the
watershed basin model.  Overland flow units represent the amount of
overland flow each cell contributes to surface flow.  If omitted, a
value of one (1) is assumed.

<dt><em>disturbed_land</em> 

<dd>Raster map input layer or value containing the percent of disturbed
land (i.e., croplands, and construction sites) where the raster or input
value of 17 equals 17%.  If no map or value is given, <em>r.watershed</em>
assumes no disturbed land.  This input is used for the RUSLE calculations.

<dt><em>blocking</em> 

<dd>Input map: terrain that will block overland surface flow.  Terrain
that will block overland surface flow and restart the slope length
for the RUSLE.  All cells that are not NULL and not zero indicate blocking 
terrain.

<dt><em>threshold</em> 

<dd>The minimum size of an exterior watershed basin in cells, if no flow
map is input, or overland flow units when a flow map is given.
Warning: low threshold values will dramactically increase run time and
generate difficult to read basin and half_basin results.
This parameter also controls the level of detail in the <em>stream</em>
segments map.

<dt><em>max_slope_length</em> 

<dd>Input value indicating the maximum length of overland surface flow
in meters.  If overland flow travels greater than the maximum length,
the program assumes the maximum length (it assumes that landscape
characteristics not discernible in the digital elevation model exist
that maximize the slope length).  This input is used for the RUSLE calculations
and is a sensitive parameter.

<dt><em>accumulation</em> 

<dd>Output map: The absolute value of each cell in this output map layer is
the amount of overland flow that traverses the cell. This value will be
the 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 flow units.
Negative numbers indicate that those cells possibly have surface runoff
from outside of the current geographic region. Thus, any cells with
negative values cannot have their surface runoff and sedimentation yields
calculated accurately.

<dt><em>tci</em> 

<dd>Output map: The topographic index TCI is computed as 
<em>ln(&alpha; / tan(&beta;))</em> where &alpha; a is the cumulative 
uplsope area draining through a point per unit contour length and 
tan(&beta;) is the local slope angle. The TCI reflects the tendency of 
water to accumulate at any point in the catchment and the tendency for 
gravitaional forces to move that water downslope (Quinn et al. 1991). 
This value will be negative if &alpha; / tan(&beta;) &lt; 1.

<dt><em>drainage</em> 

<dd>Output map: drainage direction.  Provides the "aspect" for each
cell measured CCW from East.  Multiplying positive values by 45 will give 
the direction in degrees that the surface runoff will travel from that 
cell.  The value 0 (zero) indicates that the cell is a depression area 
(defined by the depression input map).  Negative values indicate that
surface runoff is leaving the boundaries of the current geographic
region.  The absolute value of these negative cells indicates the
direction of flow.

<dt><em>basin</em> 

<dd>Output map: Unique label for each watershed basin.  Each basin will
be given a unique positive even integer.  Areas along edges may not
be large enough to create an exterior watershed basin.  0 values
indicate that the cell is not part of a complete watershed basin
in the current geographic region.

<dt><em>stream</em> 

<dd>Output map: stream segments.  Values correspond to the watershed
basin values.  Can be vectorized after thinning (<em>r.thin</em>) with 
<em>r.to.vect</em>.

<dt><em>half_basin</em> 

<dd>Output map: each half-basin is given a unique value.  Watershed
basins are divided into left and right sides.  The right-hand side
cell of the watershed basin (looking upstream) are given even values
corresponding to the values in basin.  The left-hand side
cells of the watershed basin are given odd values which are one less
than the value of the watershed basin.

<dt><em>length_slope</em> 

<dd>Output map: slope length and steepness (LS) factor for the Revised 
Universal Soil Loss Equation (RUSLE).  Equations taken from <em>Revised 
Universal Soil Loss Equation for Western Rangelands</em>
(Weltz et al. 1987). Since the LS factor is a small number (usually less 
than one), the GRASS output map is of type DCELL.

<dt><em>slope_steepness</em> 

<dd>Output map: slope steepness (S) factor for the Universal Soil
Loss Equation (RUSLE).  Equations taken from article entitled
<em>Revised Slope Steepness Factor for the Universal Soil
Loss Equation</em> (McCool et al. 1987).  Since the S factor is a small 
number (usually less than one), the GRASS output map is of type DCELL.
</dd>
</dl>


<h2>NOTES</h2>

<h3>A<sup>T</sup> least-cost search algorithm</h3>
<em>r.watershed</em> uses an A<sup>T</sup> least-cost search algorithm 
(see <a href="#references">REFERENCES</a> section) designed to minimize 
the impact of DEM data errors. Compared to <em>r.terraflow</em>, this 
algorithm provides more accurate results in areas of low slope as well 
as DEMs constructed with techniques that mistake canopy tops as the 
ground elevation. Kinner et al. (2005), for example, used SRTM and IFSAR 
DEMs to compare <em>r.watershed</em> against <em>r.terraflow</em> 
results in Panama. <em>r.terraflow</em> was unable to replicate stream 
locations in the larger valleys while <em>r.watershed</em> performed 
much better. Thus, if forest canopy exists in valleys, SRTM, IFSAR, and 
similar data products will cause major errors in <em>r.terraflow</em> 
stream output. Under similar conditions, <em>r.watershed</em> will 
generate better <b>stream</b> and <b>half_basin</b> results. If 
watershed divides contain flat to low slope, <em>r.watershed</em>
will generate better basin results than <em>r.terraflow</em>.
(<em>r.terraflow</em> uses the same type of algorithm as ESRI's ArcGIS
watershed software which fails under these conditions.) Also, if watershed
divides contain forest canopy mixed with uncanopied areas using SRTM, IFSAR,
and similar data products, <em>r.watershed</em> will generate better basin
results than <em>r.terraflow</em>.
The algorithm produces results similar to those obtained when running
<em><a href="r.cost.html">r.cost</a></em> and
<em><a href="r.drain.html">r.drain</a></em> on every cell on the map.

<h3>Multiple flow direction (MFD)</h3>

<em>r.watershed</em> 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 cells with lower 
elevation, using slope towards neighbouring cells as a weighing factor 
for proportional distribution. The A<sup>T</sup> least-cost path is 
always included. As a result, depressions and obstacles are traversed 
with a gracefull flow convergence before the overflow. The convergence 
factor causes flow accumulation to converge more strongly with higher 
values. The supported range is 1 to 10, recommended is a convergence 
factor of 5 (Holmgren, 1994). If many small sliver basins are created 
with MFD, setting the convergence factor to a higher value can reduce 
the amount of small sliver basins.

<h3>In-memory mode and disk swap mode</h3>
There are two versions of this program: <em>ram</em> and <em>seg</em>.
<em>ram</em> is used by default, <em>seg</em> can be used by setting 
the <em>-m</em> flag.
<br>
The <em>ram</em> version requires a maximum of 31 MB of RAM for 1 million 
cells. Together with the amount of system memory (RAM) available, this 
value can be used to estimate whether the current region can be 
processed with the <em>ram</em> version.
<br>
The <em>ram</em> version uses virtual memory managed by the operating
system to store all the data structures and is faster than the <em>seg</em>
version; <em>seg</em> uses the GRASS segmentation library which manages 
data in disk files. <em>seg</em> uses only as much system memory (RAM) as 
specified with the <em>memory</em> option, allowing other processes to 
operate on the same system, even when the current geographic region is huge.
<br>
Due to memory requirements of both programs, it is quite easy to run out of
memory when working with huge map regions. If the <em>ram</em> version runs
out of memory and the resolution size of the current geographic region
cannot be increased, either more memory needs to be added to the computer,
or the swap space size needs to be increased. If <em>seg</em> runs out of
memory, additional disk space needs to be freed up for the program to run.
The <em>r.terraflow</em> module was specifically designed with huge
regions in mind and may be useful here as an alternative, although disk
space requirements of <em>r.terraflow</em> are several times higher than
of <em>seg</em>.

<h3>Large regions with many cells</h3>
The upper limit of the <em>ram</em> version is 2 billion 
(2<sup>31</sup> - 1) cells, whereas the upper limit for the <em>seg</em> 
version is 9 billion billion (2<sup>63</sup> - 1) cells.<br>
In some situations, the region size (number of cells) may be too large for
the amount of time or memory available. Running <em>r.watershed</em> may
then require use of a coarser resolution. To make the results more closely
resemble the finer terrain data, create a map layer containing the
lowest elevation values at the coarser resolution. This is done by:
1) Setting the current geographic region equal to the elevation map
layer with <em>g.region</em>, and 2) Use the <em>r.neighbors</em> or
<em>r.resamp.stats</em> command to find the lowest value for an area
equal in size to the desired resolution. For example, if the resolution
of the elevation data is 30 meters and the resolution of the geographic
region for <em>r.watershed</em> will be 90 meters: use the minimum 
function for a 3 by 3 neighborhood. After changing to the resolution at
which <em>r.watershed</em> will be run, <em>r.watershed</em> should be run
using the values from the <em>neighborhood</em> output map layer that
represents the minimum elevation within the region of the coarser cell.

<h3>Basin threshold</h3>
The minimum size of drainage basins, defined by the <em>threshold</em>
parameter, is only relevant for those watersheds with a single stream
having at least the <em>threshold</em> of cells flowing into it.
(These watersheds are called exterior basins.)
Interior drainage basins contain stream segments below multiple tributaries.
Interior drainage basins can be of any size because the length of
an interior stream segment is determined by the distance between the
tributaries flowing into it.

<h3>MASK and no data</h3>
<p>
The <em>r.watershed</em> program does not require the user to have the
current geographic region filled with elevation values.  Areas without
elevation data (masked or NULL cells) are ignored. It is NOT necessary to
create a raster map (or raster reclassification) named <tt>MASK</tt> for 
NULL cells.  Areas without elevation data will be treated as if they are 
off the edge of the region. Such areas will reduce the memory necessary 
to run the program.  Masking out unimportant areas can significantly 
reduce processing time if the watersheds of interest occupy a small 
percentage of the overall area.
<p>
Gaps (NULL cells) in the elevation map that are located within the area 
of interest will heavily influence the analysis: water will 
flow into but not out of these gaps. These gaps must be filled beforehand, 
e.g. with <em>r.fillnulls</em>.
<p>
Zero (0) and negative values will be treated as elevation data (not no_data).

<h3>Further processing of output layers</h3>
<p>
Problem areas, i.e. those parts of a basin with a likely underestimate of
flow accumulation, can be easily identified with e.g.
<div class="code"><pre>
  r.mapcalc "problems = if(flow_acc &lt; 0, basin, null())"
</pre></div>
If the region of interest contains such problem areas, and this is not
desired, the computational region must be expanded until the catchment
area for the region of interest is completely included.

<p>
To isolate an individual river network using the output of this module,
a number of approaches may be considered.
<ol>
<li>Use a resample of the basins catchment raster map as a MASK.<br>
  The equivalent vector map method is similar using <em>v.select</em> or
  <em>v.overlay</em>.
<li>Use the <em>r.cost</em> module with a point in the river as a starting
  point.
<li>Use the <em>v.net.iso</em> module with a node in the river as a
  starting point.
</ol>

All individual river networks in the stream segments output can be
identified through their ultimate outlet points. These points are all
cells in the stream segments output with negative drainage direction.
These points can be used as start points for <em>r.water.outlet</em> or
<em>v.net.iso</em>.

<p>
To create <i>river mile</i> segmentation from a vectorized streams map,
try the <em>v.net.iso</em> or <em>v.lrs.segment</em> modules.
<p>
The stream segments output can be easily vectorized after thinning with 
<em>r.thin</em>. Each stream segment in the vector map will have the 
value of the associated basin. To isolate subbasins and streams for a 
larger basin, a MASK for the larger basin can be created with 
<em>r.water.outlet</em>. The stream segments output serves as a guide 
where to place the outlet point used as input to <em>r.water.outlet</em>. 
The basin threshold must have been sufficiently small to isolate a 
stream network and subbasins within the larger basin.

<h2>EXAMPLES</h2>
<i>These examples use the Spearfish sample dataset.</i>
<p>
Convert <em>r.watershed</em> streams map output to a vector layer.
<p>
If you want a detailed stream network, set the threshold option
small to create lots of catchment basins, as only one stream is
presented per catchment. The r.to.vect -v flag preserves the
catchment ID as the vector category number.

<div class="code"><pre>
  r.watershed elev=elevation.dem stream=rwater.stream
  r.to.vect -v in=rwater.stream out=rwater_stream
</pre></div>
<br>

<p>
Set a different color table for the accumulation map:
<div class="code"><pre>
  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 &lt;&lt; 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
</pre></div>
<br>


<p>
Create a more detailed stream map using the accumulation map and convert
it to a vector output map. The accumulation cut-off, and therefore fractal
dimension, is arbitrary; in this example we use the map's mean number of
upstream catchment cells (calculated in the above example by <em>r.univar</em>)
as the cut-off value. This only works with SFD, not with MFD.
<div class="code"><pre>
  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 MASK

  # <i>Thinning is required before converting raster lines to vector</i>
  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 feature=line
  v.db.dropcolumn map=rwater_course column=label
</pre></div>
<!-- can't set line attribute to catchment it is in as v.what.rast and 
  v.distance only work for point features. Could create endpoint node
  points map and upload to that ?? -->
<!-- Note value column containing accumulation cells in output vector
  may not necessarily reference the downstream end of the line! drop it? -->
<br>

<p>
Create watershed basins map and convert to a vector polygon map
<div class="code"><pre>
  r.watershed elev=elevation.dem basin=rwater.basin thresh=15000
  r.to.vect -s in=rwater.basin out=rwater_basins feature=area
  v.db.dropcolumn map=rwater_basins column=label
  v.db.renamecolumn map=rwater_basins column=value,catchment
</pre></div>
<br>

<p>
Display output in a nice way
<div class="code"><pre>
  r.shaded.relief map=elevation.dem
  d.shadedmap rel=elevation.dem.shade drape=rwater.basin bright=40
  d.vect rwater_course color=orange
</pre></div>
<br>

<a name="references"></a>
<h2>REFERENCES</h2>


Ehlschlaeger C. (1989). <i>Using the A<sup>T</sup> Search Algorithm
to Develop Hydrologic Models from Digital Elevation Data</i>,
<b>Proceedings of International Geographic Information Systems (IGIS)
Symposium '89</b>, pp 275-281 (Baltimore, MD, 18-19 March 1989).<br>
URL: <a href="http://chuck.ehlschlaeger.info/older/IGIS/paper.html">
http://chuck.ehlschlaeger.info/older/IGIS/paper.html</a>

<p>
Holmgren P. (1994). <i>Multiple flow direction algorithms for runoff 
modelling in grid based elevation models: An empirical evaluation.</i>
<b>Hydrological Processes</b> Vol 8(4), 327-334.<br>
DOI: <a href="http://dx.doi.org/10.1002/hyp.3360080405">10.1002/hyp.3360080405</a>

<p>
Kinner D., Mitasova H., Harmon R., Toma L., Stallard R. (2005).
<i>GIS-based Stream Network Analysis for The Chagres River Basin,
Republic of Panama</i>. <b>The Rio Chagres: A Multidisciplinary Profile of
a Tropical Watershed</b>, R. Harmon (Ed.), Springer/Kluwer, p.83-95.<br>
URL: <a href="http://skagit.meas.ncsu.edu/%7Ehelena/measwork/panama/panama.html">
http://skagit.meas.ncsu.edu/~helena/measwork/panama/panama.html</a>

<p>
McCool et al. (1987). <i>Revised Slope Steepness Factor for the Universal
Soil Loss Equation</i>, <b>Transactions of the ASAE</b> Vol 30(5).

<p>
Metz M., Mitasova H., Harmon R. (2011). <i>Efficient extraction of 
drainage networks from massive, radar-based elevation models with least 
cost path search</i>, <b>Hydrol. Earth Syst. Sci.</b> Vol 15, 667-678.<br>
DOI: <a href="http://dx.doi.org/10.5194/hess-15-667-2011">10.5194/hess-15-667-2011</a>

<p>
Quinn P., K. Beven K., Chevallier P., Planchon O. (1991). <i>The 
prediction of hillslope flow paths for distributed hydrological modelling 
using Digital Elevation Models</i>, <b>Hydrological Processes</b> Vol 5(1), 
p.59-79.<br>
DOI: <a href="http://dx.doi.org/10.1002/hyp.3360050106">10.1002/hyp.3360050106</a>

<p>
Weltz M. A., Renard K.G., Simanton J. R. (1987). <i>Revised Universal Soil
Loss Equation for Western Rangelands</i>, <b>U.S.A./Mexico Symposium of
Strategies for Classification and Management of Native Vegetation for
Food Production In Arid Zones</b> (Tucson, AZ, 12-16 Oct. 1987).

<a name="seealso"></a>
<h2>SEE ALSO</h2>

<em>
<a href="g.region.html">g.region</a>,
<a href="r.cost.html">r.cost</a>,
<a href="r.drain.html">r.drain</a>,
<a href="r.fillnulls.html">r.fillnulls</a>,
<a href="r.flow.html">r.flow</a>,
<!-- <a href="r.flowmd.html">r.flowmd</a>, -->
<a href="r.mask.html">r.mask</a>,
<a href="r.neighbors.html">r.neighbors</a>,
<a href="r.param.scale.html">r.param.scale</a>,
<a href="r.resamp.interp.html">r.resamp.interp</a>,
<a href="r.terraflow.html">r.terraflow</a>,
<a href="r.topidx.html">r.topidx</a>,
<a href="r.water.outlet.html">r.water.outlet</a>
</em>


<h2>AUTHORS</h2>

Original version:
Charles Ehlschlaeger, U.S. Army Construction Engineering Research Laboratory
<br>
Faster sorting algorithm and MFD support:
Markus Metz &lt;markus.metz.giswork at gmail.com&gt;

<p>
<i>Last changed: $Date$</i>