1r.watershed(1) GRASS GIS User's Manual r.watershed(1)
2
6 r.watershed - Calculates hydrological parameters and RUSLE factors.
7
9 raster, hydrology, watershed, accumulation, drainage, stream network,
10 stream power index, topographic index
11
13 r.watershed
14 r.watershed --help
15 r.watershed [-s4mab] elevation=name [depression=name] [flow=name]
16 [disturbed_land=name] [blocking=name] [retention=name] [thresh‐
17 old=integer] [max_slope_length=float] [accumulation=name]
18 [tci=name] [spi=name] [drainage=name] [basin=name]
19 [stream=name] [half_basin=name] [length_slope=name] [slope_steep‐
20 ness=name] [convergence=integer] [memory=memory in MB] [--over‐
21 write] [--help] [--verbose] [--quiet] [--ui]
22 Flags:
24 -s
25 SFD (D8) flow (default is MFD)
26 SFD: single flow direction, MFD: multiple flow direction
27 -4
29 Allow only horizontal and vertical flow of water
30 -m
32 Enable disk swap memory option: Operation is slow
33 Only needed if memory requirements exceed available RAM; see manual
34 on how to calculate memory requirements
35 -a
37 Use positive flow accumulation even for likely underestimates
38 See manual for a detailed description of flow accumulation output
39 -b
41 Beautify flat areas
42 Flow direction in flat areas is modified to look prettier
43 --overwrite
45 Allow output files to overwrite existing files
46 --help
48 Print usage summary
49 --verbose
51 Verbose module output
52 --quiet
54 Quiet module output
55 --ui
57 Force launching GUI dialog
58 Parameters:
60 elevation=name [required]
61 Name of input elevation raster map
62 depression=name
64 Name of input depressions raster map
65 All non-NULL and non-zero cells are considered as real depressions
66 flow=name
68 Name of input raster representing amount of overland flow per cell
69 disturbed_land=name
71 Name of input raster map percent of disturbed land
72 For USLE
73 blocking=name
75 Name of input raster map blocking overland surface flow
76 For USLE. All non-NULL and non-zero cells are considered as block‐
77 ing terrain.
78 retention=name
80 Name of input raster map with percentages for flow accumulation.
81 threshold=integer
83 Minimum size of exterior watershed basin
84 max_slope_length=float
86 Maximum length of surface flow in map units
87 For USLE
88 accumulation=name
90 Name for output accumulation raster map
91 Number of cells that drain through each cell
92 tci=name
94 Name for output topographic index ln(a / tan(b)) map
95 spi=name
97 Name for output stream power index a * tan(b)
98 Name for output raster map
99 drainage=name
101 Name for output drainage direction raster map
102 Directions numbered from 1 to 8
103 basin=name
105 Name for output basins raster map
106 stream=name
108 Name for output stream segments raster map
109 half_basin=name
111 Name for output half basins raster map
112 Each half-basin is given a unique value
113 length_slope=name
115 Name for output slope length raster map
116 Slope length and steepness (LS) factor for USLE
117 slope_steepness=name
119 Name for output slope steepness raster map
120 Slope steepness (S) factor for USLE
121 convergence=integer
123 Convergence factor for MFD (1-10)
124 1 = most diverging flow, 10 = most converging flow. Recommended: 5
125 Default: 5
126 memory=memory in MB
128 Maximum memory to be used (in MB)
129 Cache size for raster rows
130 Default: 300
131
133 r.watershed generates a set of maps indicating: 1) flow accumulation,
134 drainage direction, the location of streams and watershed basins, and
135 2) the LS and S factors of the Revised Universal Soil Loss Equation
136 (RUSLE).
137
139 Without flag -m set, the entire analysis is run in memory maintained by
140 the operating system. This can be limiting, but is very fast. Setting
141 this flag causes the program to manage memory on disk which allows very
142 large maps to be processed but is slower.
143 Flag -s force the module to use single flow direction (SFD, D8) instead
145 of multiple flow direction (MFD). MFD is enabled by default.
146 By -4 flag the user allow only horizontal and vertical flow of water.
148 Stream and slope lengths are approximately the same as outputs from
149 default surface flow (allows horizontal, vertical, and diagonal flow of
150 water). This flag will also make the drainage basins look more homoge‐
151 neous.
152 When -a flag is specified the module will use positive flow accumula‐
154 tion even for likely underestimates. When this flag is not set, cells
155 with a flow accumulation value that is likely to be an underestimate
156 are converted to the negative. See below for a detailed description of
157 flow accumulation output.
158 Option convergence specifies convergence factor for MFD. Lower values
160 result in higher divergence, flow is more widely distributed. Higher
161 values result in higher convergence, flow is less widely distributed,
162 becoming more similar to SFD.
163 Option elevation specifies the elevation data on which entire analysis
165 is based. NULL (nodata) cells are ignored, zero and negative values are
166 valid elevation data. Gaps in the elevation map that are located
167 within the area of interest must be filled beforehand, e.g. with
168 r.fillnulls, to avoid distortions. The elevation map need not be
169 sink-filled because the module uses a least-cost algorithm.
170 Option depression specifies the optional map of actual depressions or
172 sinkholes in the landscape that are large enough to slow and store sur‐
173 face runoff from a storm event. All cells that are not NULL and not
174 zero indicate depressions. Water will flow into but not out of depres‐
175 sions.
176 Raster flow map specifies amount of overland flow per cell. This map
178 indicates the amount of overland flow units that each cell will con‐
179 tribute to the watershed basin model. Overland flow units represent the
180 amount of overland flow each cell contributes to surface flow. If omit‐
181 ted, a value of one (1) is assumed.
182 Raster retention map specifies correction factors per cell for flow
184 distribution. This map indicates the percentage of overland flow units
185 leaving each cell. Values should be between zero and 100. If omitted, a
186 value of 100 is assumed.
187 Input Raster map or value containing the percent of disturbed land
189 (i.e., croplands, and construction sites) where the raster or input
190 value of 17 equals 17%. If no map or value is given, r.watershed
191 assumes no disturbed land. This input is used for the RUSLE calcula‐
192 tions.
193 Option blocking specifies terrain that will block overland surface
195 flow. Blocking cells and streams stop the slope length for the RUSLE.
196 All cells that are not NULL and not zero indicate blocking terrain.
197 Option threshold specifies the minimum size of an exterior watershed
199 basin in cells, if no flow map is input, or overland flow units when a
200 flow map is given. Warning: low threshold values will dramactically
201 increase run time and generate difficult to read basin and half_basin
202 results. This parameter also controls the level of detail in the
203 stream segments map.
204 Value given by max_slope_length option indicates the maximum length of
206 overland surface flow in meters. If overland flow travels greater than
207 the maximum length, the program assumes the maximum length (it assumes
208 that landscape characteristics not discernible in the digital elevation
209 model exist that maximize the slope length). This input is used for
210 the RUSLE calculations and is a sensitive parameter.
211 Output accumulation map contains the absolute value of the amount of
213 overland flow that traverses each cell. This value will be the number
214 of upland cells plus one if no overland flow map is given. If the
215 overland flow map is given, the value will be in overland flow units.
216 Negative numbers indicate that those cells possibly have surface runoff
217 from outside of the current geographic region. Thus, any cells with
218 negative values cannot have their surface runoff and sedimentation
219 yields calculated accurately.
220 Output tci raster map contains topographic index TCI, computed as ln(α
222 / tan(β)) where α is the cumulative upslope area draining through a
223 point per unit contour length and tan(β) is the local slope angle. The
224 TCI reflects the tendency of water to accumulate at any point in the
225 catchment and the tendency for gravitational forces to move that water
226 downslope (Quinn et al. 1991). This value will be negative if α /
227 tan(β) < 1.
228 Output spi raster map contains stream power index SPI, computed as α *
230 tan(β) where α is the cumulative upslope area draining through a
231 point per unit contour length and tan(β) is the local slope angle. The
232 SPI reflects the power of water flow at any point in the catchment and
233 the tendency for gravitational forces to move that water downslope
234 (Moore et al. 1991). This value will be negative if α < 0, i.e. for
235 cells with possible surface runoff from outside of the current geo‐
236 graphic region.
237 Output drainage raster map contains drainage direction. Provides the
239 "aspect" for each cell measured CCW from East.
240 Figure: Drainage is 8 directions numbered counter-clockwise starting
241 from 1 in north-east direction (source) Multiplying positive values by
242 45 will give the direction in degrees that the surface runoff will
243 travel from that cell. The value 0 (zero) indicates that the cell is a
244 depression area (defined by the depression input map). Negative values
245 indicate that surface runoff is leaving the boundaries of the current
246 geographic region. The absolute value of these negative cells indi‐
247 cates the direction of flow. For MFD, drainage indicates the direction
248 of maximum flow.
249 The output basin map contains unique label for each watershed basin.
251 Each basin will be given a unique positive even integer. Areas along
252 edges may not be large enough to create an exterior watershed basin.
253 NULL values indicate that the cell is not part of a complete watershed
254 basin in the current geographic region.
255 The output stream contains stream segments. Values correspond to the
257 watershed basin values. Can be vectorized after thinning (r.thin) with
258 r.to.vect.
259 The output half_basin raster map stores each half-basin is given a
261 unique value. Watershed basins are divided into left and right sides.
262 The right-hand side cell of the watershed basin (looking upstream) are
263 given even values corresponding to the values in basin. The left-hand
264 side cells of the watershed basin are given odd values which are one
265 less than the value of the watershed basin.
266 The output length_slope raster map stores slope length and steepness
268 (LS) factor for the Revised Universal Soil Loss Equation (RUSLE).
269 Equations taken from Revised Universal Soil Loss Equation for Western
270 Rangelands (Weltz et al. 1987). Since the LS factor is a small number
271 (usually less than one), the GRASS output map is of type DCELL.
272 The output slope_steepness raster map stores slope steepness (S) factor
274 for the Universal Soil Loss Equation (RUSLE). Equations taken from
275 article entitled Revised Slope Steepness Factor for the Universal Soil
276 Loss Equation (McCool et al. 1987). Since the S factor is a small num‐
277 ber (usually less than one), the GRASS output map is of type DCELL.
278 AT least-cost search algorithm
280 r.watershed uses an AT least-cost search algorithm (see REFERENCES sec‐
281 tion) designed to minimize the impact of DEM data errors. Compared to
282 r.terraflow, this algorithm provides more accurate results in areas of
283 low slope as well as DEMs constructed with techniques that mistake
284 canopy tops as the ground elevation. Kinner et al. (2005), for example,
285 used SRTM and IFSAR DEMs to compare r.watershed against r.terraflow
286 results in Panama. r.terraflow was unable to replicate stream locations
287 in the larger valleys while r.watershed performed much better. Thus, if
288 forest canopy exists in valleys, SRTM, IFSAR, and similar data products
289 will cause major errors in r.terraflow stream output. Under similar
290 conditions, r.watershed will generate better stream and half_basin
291 results. If watershed divides contain flat to low slope, r.watershed
292 will generate better basin results than r.terraflow. (r.terraflow uses
293 the same type of algorithm as ESRI’s ArcGIS watershed software which
294 fails under these conditions.) Also, if watershed divides contain for‐
295 est canopy mixed with uncanopied areas using SRTM, IFSAR, and similar
296 data products, r.watershed will generate better basin results than
297 r.terraflow. The algorithm produces results similar to those obtained
298 when running r.cost and r.drain on every cell on the raster map.
299 Multiple flow direction (MFD)
301 r.watershed offers two methods to calculate surface flow: single flow
302 direction (SFD, D8) and multiple flow direction (MFD). With MFD, water
303 flow is distributed to all neighbouring cells with lower elevation,
304 using slope towards neighbouring cells as a weighing factor for propor‐
305 tional distribution. The AT least-cost path is always included. As a
306 result, depressions and obstacles are traversed with a graceful flow
307 convergence before the overflow. The convergence factor causes flow
308 accumulation to converge more strongly with higher values. The sup‐
309 ported range is 1 to 10, recommended is a convergence factor of 5
310 (Holmgren, 1994). If many small sliver basins are created with MFD,
311 setting the convergence factor to a higher value can reduce the amount
312 of small sliver basins.
313 In-memory mode and disk swap mode
315 There are two versions of this program: ram and seg. ram is used by
316 default, seg can be used by setting the -m flag.
317 The ram version requires a maximum of 31 MB of RAM for 1 million cells.
319 Together with the amount of system memory (RAM) available, this value
320 can be used to estimate whether the current region can be processed
321 with the ram version.
322 The ram version uses virtual memory managed by the operating system to
324 store all the data structures and is faster than the seg version; seg
325 uses the GRASS segmentation library which manages data in disk files.
326 seg uses only as much system memory (RAM) as specified with the memory
327 option, allowing other processes to operate on the same system, even
328 when the current geographic region is huge.
329 Due to memory requirements of both programs, it is quite easy to run
331 out of memory when working with huge map regions. If the ram version
332 runs out of memory and the resolution size of the current geographic
333 region cannot be increased, either more memory needs to be added to the
334 computer, or the swap space size needs to be increased. If seg runs out
335 of memory, additional disk space needs to be freed up for the program
336 to run. The r.terraflow module was specifically designed with huge
337 regions in mind and may be useful here as an alternative, although disk
338 space requirements of r.terraflow are several times higher than of seg.
339 Large regions with many cells
341 The upper limit of the ram version is 2 billion (231 - 1) cells,
342 whereas the upper limit for the seg version is 9 billion-billion (263 -
343 1 = 9.223372e+18) cells.
344 In some situations, the region size (number of cells) may be too large
345 for the amount of time or memory available. Running r.watershed may
346 then require use of a coarser resolution. To make the results more
347 closely resemble the finer terrain data, create a map layer containing
348 the lowest elevation values at the coarser resolution. This is done by:
349 1) Setting the current geographic region equal to the elevation map
350 layer with g.region, and 2) Use the r.neighbors or r.resamp.stats com‐
351 mand to find the lowest value for an area equal in size to the desired
352 resolution. For example, if the resolution of the elevation data is 30
353 meters and the resolution of the geographic region for r.watershed will
354 be 90 meters: use the minimum function for a 3 by 3 neighborhood. After
355 changing to the resolution at which r.watershed will be run, r.water‐
356 shed should be run using the values from the neighborhood output map
357 layer that represents the minimum elevation within the region of the
358 coarser cell.
359 Basin threshold
361 The minimum size of drainage basins, defined by the threshold parame‐
362 ter, is only relevant for those watersheds with a single stream having
363 at least the threshold of cells flowing into it. (These watersheds are
364 called exterior basins.) Interior drainage basins contain stream seg‐
365 ments below multiple tributaries. Interior drainage basins can be of
366 any size because the length of an interior stream segment is determined
367 by the distance between the tributaries flowing into it.
368 MASK and no data
370 The r.watershed program does not require the user to have the current
371 geographic region filled with elevation values. Areas without eleva‐
372 tion data (masked or NULL cells) are ignored. It is NOT necessary to
373 create a raster map (or raster reclassification) named MASK for NULL
374 cells. Areas without elevation data will be treated as if they are off
375 the edge of the region. Such areas will reduce the memory necessary to
376 run the program. Masking out unimportant areas can significantly
377 reduce processing time if the watersheds of interest occupy a small
378 percentage of the overall area.
379 Gaps (NULL cells) in the elevation map that are located within the area
381 of interest will heavily influence the analysis: water will flow into
382 but not out of these gaps. These gaps must be filled beforehand, e.g.
383 with r.fillnulls.
384 Zero (0) and negative values will be treated as elevation data (not
386 no_data).
387 Further processing of output layers
389 Problem areas, i.e. those parts of a basin with a likely underestimate
390 of flow accumulation, can be easily identified with e.g.
391 r.mapcalc "problems = if(flow_acc < 0, basin, null())"
392 If the region of interest contains such problem areas, and this is not
393 desired, the computational region must be expanded until the catchment
394 area for the region of interest is completely included.
395 To isolate an individual river network using the output of this module,
397 a number of approaches may be considered.
398 1 Use a resample of the basins catchment raster map as a MASK.
400 The equivalent vector map method is similar using v.select or
401 v.overlay.
402 2 Use the r.cost module with a point in the river as a starting
404 point.
405 3 Use the v.net.iso module with a node in the river as a starting
407 point.
408 All individual river networks in the stream segments output can be
410 identified through their ultimate outlet points. These points are all
411 cells in the stream segments output with negative drainage direction.
412 These points can be used as start points for r.water.outlet or
413 v.net.iso.
414 To create river mile segmentation from a vectorized streams map, try
416 the v.net.iso or v.lrs.segment modules.
417 The stream segments output can be easily vectorized after thinning with
419 r.thin. Each stream segment in the vector map will have the value of
420 the associated basin. To isolate subbasins and streams for a larger
421 basin, a MASK for the larger basin can be created with r.water.outlet.
422 The stream segments output serves as a guide where to place the outlet
423 point used as input to r.water.outlet. The basin threshold must have
424 been sufficiently small to isolate a stream network and subbasins
425 within the larger basin.
426 Given that the drainage is 8 directions numbered counter-clockwise
428 starting from 1 in north-east direction, multiplying the output by 45
429 (by 45. to get a double precision floating point raster map in r.map‐
430 calc) gives the directions in degrees. For most applications, zeros
431 which indicate depressions specified by depression and negative values
432 which indicate runoff leaving the region should be replaced by NULL
433 (null() in r.mapcalc). The following command performs these replace‐
434 ments:
435 r.mapcalc "drainage_degrees = if(drainage > 0, 45. * drainage, null())"
436 Alternatively, the user can use the -a flag or later the abs() function
437 in r.mapcalc if the runoff is leaving the region.
438
440 These examples use the Spearfish sample dataset.
441 Convert r.watershed streams map output to a vector map
443 If you want a detailed stream network, set the threshold option small
444 to create lots of catchment basins, as only one stream is presented per
445 catchment. The r.to.vect -v flag preserves the catchment ID as the vec‐
446 tor category number.
447 r.watershed elev=elevation.dem stream=rwater.stream
448 r.to.vect -v in=rwater.stream out=rwater_stream
449 Set a different color table for the accumulation map:
451 MAP=rwater.accum
452 r.watershed elev=elevation.dem accum=$MAP
453 eval `r.univar -g "$MAP"`
454 stddev_x_2=`echo $stddev | awk ’{print $1 * 2}’`
455 stddev_div_2=`echo $stddev | awk ’{print $1 / 2}’`
456 r.colors $MAP col=rules << EOF
457 0% red
458 -$stddev_x_2 red
459 -$stddev yellow
460 -$stddev_div_2 cyan
461 -$mean_of_abs blue
462 0 white
463 $mean_of_abs blue
464 $stddev_div_2 cyan
465 $stddev yellow
466 $stddev_x_2 red
467 100% red
468 EOF
469 Create a more detailed stream map using the accumulation map and con‐
471 vert it to a vector output map. The accumulation cut-off, and therefore
472 fractal dimension, is arbitrary; in this example we use the map’s mean
473 number of upstream catchment cells (calculated in the above example by
474 r.univar) as the cut-off value. This only works with SFD, not with MFD.
475 r.watershed elev=elevation.dem accum=rwater.accum
476 r.mapcalc ’MASK = if(!isnull(elevation.dem))’
477 r.mapcalc "rwater.course = \
478 if( abs(rwater.accum) > $mean_of_abs, \
479 abs(rwater.accum), \
480 null() )"
481 r.colors -g rwater.course col=bcyr
482 g.remove -f type=raster name=MASK
483 # Thinning is required before converting raster lines to vector
484 r.thin in=rwater.course out=rwater.course.Thin
485 r.colors -gn rwater.course.Thin color=grey
486 r.to.vect in=rwater.course.Thin out=rwater_course type=line
487 v.db.dropcolumn map=rwater_course column=label
488 Create watershed basins map and convert to a vector polygon map
490 r.watershed elev=elevation.dem basin=rwater.basin thresh=15000
491 r.to.vect -s in=rwater.basin out=rwater_basins type=area
492 v.db.dropcolumn map=rwater_basins column=label
493 v.db.renamecolumn map=rwater_basins column=value,catchment
494 Display output in a nice way
496 r.relief map=elevation.dem
497 d.shade shade=elevation.dem.shade color=rwater.basin bright=40
498 d.vect rwater_course color=orange
499
501 · Ehlschlaeger C. (1989). Using the AT Search Algorithm to
502 Develop Hydrologic Models from Digital Elevation Data, Proceed‐
503 ings of International Geographic Information Systems (IGIS)
504 Symposium ’89, pp 275-281 (Baltimore, MD, 18-19 March 1989).
505 URL: http://chuck.ehlschlaeger.info/older/IGIS/paper.html
506 · Holmgren P. (1994). Multiple flow direction algorithms for
508 runoff modelling in grid based elevation models: An empirical
509 evaluation. Hydrological Processes Vol 8(4), 327-334.
510 DOI: 10.1002/hyp.3360080405
511 · Kinner D., Mitasova H., Harmon R., Toma L., Stallard R. (2005).
513 GIS-based Stream Network Analysis for The Chagres River Basin,
514 Republic of Panama. The Rio Chagres: A Multidisciplinary Pro‐
515 file of a Tropical Watershed, R. Harmon (Ed.), Springer/Kluwer,
516 p.83-95.
517 URL: http://www4.ncsu.edu/~hmitaso/measwork/panama/panama.html
518 · McCool et al. (1987). Revised Slope Steepness Factor for the
520 Universal Soil Loss Equation, Transactions of the ASAE Vol
521 30(5).
522 · Metz M., Mitasova H., Harmon R. (2011). Efficient extraction of
524 drainage networks from massive, radar-based elevation models
525 with least cost path search, Hydrol. Earth Syst. Sci. Vol 15,
526 667-678.
527 DOI: 10.5194/hess-15-667-2011
528 · Moore I.D., Grayson R.B., Ladson A.R. (1991). Digital terrain
530 modelling: a review of hydrogical, geomorphological, and bio‐
531 logical applications, Hydrological Processes, Vol 5(1), 3-30
532 DOI: 10.1002/hyp.3360050103
533 · Quinn P., K. Beven K., Chevallier P., Planchon O. (1991). The
535 prediction of hillslope flow paths for distributed hydrological
536 modelling using Digital Elevation Models, Hydrological Pro‐
537 cesses Vol 5(1), p.59-79.
538 DOI: 10.1002/hyp.3360050106
539 · Weltz M. A., Renard K.G., Simanton J. R. (1987). Revised Uni‐
541 versal Soil Loss Equation for Western Rangelands, U.S.A./Mexico
542 Symposium of Strategies for Classification and Management of
543 Native Vegetation for Food Production In Arid Zones (Tucson,
544 AZ, 12-16 Oct. 1987).
545
547 g.region, r.cost, r.drain, r.fillnulls, r.flow, r.mask, r.neighbors,
548 r.param.scale, r.resamp.interp, r.terraflow, r.topidx, r.water.outlet,
549 r.stream.extract
550
552 Original version: Charles Ehlschlaeger, U.S. Army Construction Engi‐
553 neering Research Laboratory
554 Faster sorting algorithm and MFD support: Markus Metz <markus.metz.gis‐
555 work at gmail.com>
556 Retention for flow distribution by Andreas Gericke (IGB Berlin)
557
559 Available at: r.watershed source code (history)
560 Main index | Raster index | Topics index | Keywords index | Graphical
562 index | Full index
563 © 2003-2020 GRASS Development Team, GRASS GIS 7.8.5 Reference Manual
565GRASS 7.8.5 r.watershed(1)