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