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).

<!-- Interactive mode not activated in GRASS 7.
<p>
<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>NOTES</h2>

Without flag <b>-m</b> 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.

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

<p>
By <b>-4</b> flag the user 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.

<p>
When <b>-a</b> flag is specified the module will 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.

<p>
Option <b>convergence</b> specifies 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.

<p>
Option <b>elevation</b> specifies the elevation data 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><a href="r.fillnulls.html">r.fillnulls</a></em>, to
avoid distortions.  The elevation map need not be sink-filled because
the module uses a least-cost algorithm.

<p>
Map 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.

<p>
Raster <b>flow</b> map specifies 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.

<p>
Input Raster map 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.

<p>
Option <b>blocking</b> specifies terrain that will block overland
surface flow.  errain 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.

<p>
Option <b>threshold</b> specifies 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 <b>stream</b> segments map.

<p>
Value given by <b>max_slope_length</b> option indicates 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.

<p>
Output <b>accumulation</b> map contains the absolute value of each
cell in this output map 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.

<p>
Output <b>tci</b> raster map contains topographic index TCI is
computed as
<tt>ln(&alpha; / tan(&beta;))</tt> where &alpha; is the cumulative
upslope area draining through a point per unit contour length and
<tt>tan(&beta;)</tt> 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 <tt>&alpha; /
tan(&beta;) &lt; 1</tt>.

<p>
Output <b>drainage</b> raster map contains drainage direction.
Provides the &quot;aspect&quot; 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.

<p>
The output <b>basin</b> map contains 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.

<p>
The output <b>stream</b> contains stream segments. Values correspond
to the watershed basin values.  Can be vectorized after thinning
(<em><a href="r.thin.html">r.thin</a></em>) with
<em><a href="r.to.vect.html">r.to.vect</a></em>.

<p>
The output <b>half_basin</b> raster map stores 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.

<p>
The output <b>length_slope</b> raster map stores 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.

<p>
The output <b>slope_steepness</b> raster map stores 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.

<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 REFERENCES section) designed to minimize the impact of DEM data
errors. Compared
to <em><a href="r.terraflow.html">r.terraflow</a></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><a href="r.terraflow.html">r.terraflow</a></em> results in
Panama. <em><a href="r.terraflow.html">r.terraflow</a></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><a href="r.terraflow.html">r.terraflow</a></em>. (<em><a href="r.terraflow.html">r.terraflow</a></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><a href="r.terraflow.html">r.terraflow</a></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 raster 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 <b>-m</b> flag.

<p>
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.

<p>
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 <b>memory</b> option,
allowing other processes to operate on the same system, even when the
current geographic region is huge.

<p>
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><a href="r.terraflow.html">r.terraflow</a></em> module was
specifically designed with huge regions in mind and may be useful here
as an alternative, although disk space requirements
of <em><a href="r.terraflow.html">r.terraflow</a></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><a href="g.region.html">g.region</a></em>, and 2) Use
the <em><a href="r.neighbors.html">r.neighbors</a></em> or
<em><a href="r.resamp.stats.html">r.resamp.stats</a></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 <b>neighborhood</b> 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 <b>threshold</b>
parameter, is only relevant for those watersheds with a single stream
having at least the <b>threshold</b> 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>


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><a href="r.fillnulls.html">r.fillnulls</a></em>.

<p>
Zero (0) and negative values will be treated as elevation data (not
no_data).

<h3>Further processing of output layers</h3>

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><a href="v.select.html">v.select</a></em> or
  <em><a href="v.overlay.html">v.overlay</a></em>.
<li>Use the <em><a href="r.cost.html">r.cost</a></em> module with a
  point in the river as a starting point.
<li>Use the <em><a href="v.net.iso.html">v.net.iso</a></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><a href="r.water.outlet.html">r.water.outlet</a></em> or
<em><a href="v.net.iso.html">v.net.iso<a/></em>.

<p>
To create <i>river mile</i> segmentation from a vectorized streams
map, try
the <em><a href="v.net.iso.html">v.net.iso<a/></em> or <em><a href="v.lrs.segment.html">v.lrs.segment</a></em>
modules.

<p>
The stream segments output can be easily vectorized after thinning
with
<em><a href="r.thin.html">r.thin</a></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><a href="r.water.outlet.html">r.water.outlet</a></em>. The stream
segments output serves as a guide where to place the outlet point used
as input to <em><a href="r.water.outlet.html">r.water.outlet</a></em>.
The basin threshold must have been sufficiently small to isolate a
stream network and subbasins within the larger basin.

<h2>EXAMPLES</h2>

These examples use the Spearfish sample dataset.

<h3>Convert <em>r.watershed</em> streams map output to a vector map</h3>

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 <tt>r.to.vect -v</tt> 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>

<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>

<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><a href="r.univar.html">r.univar</a></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 -f type=rast name=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 type=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? -->

<h3>Create watershed basins map and convert to a vector polygon map</h3>

<div class="code"><pre>
  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
</pre></div>

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

<h2>REFERENCES</h2>

<ul>
<li>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>

<li>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>

<li>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://www4.ncsu.edu/~hmitaso/measwork/panama/panama.html">
http://www4.ncsu.edu/~hmitaso/measwork/panama/panama.html</a>

<li>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).

<li>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>

<li>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>

<li>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).
</li>
</ul>

<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>,
<a href="r.stream.extract.html">r.stream.extract</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>