Creating watersheds

Arcview users, needing to delineate watersheds and stream networks, choose the extension called "Arc Hydro" (requires at least Spatial Analyst). This extension introduces the concepts of "Batch Points" and "Adjoint Catchments". Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points. Here we'll demonstrate how to get similar results with GRASS GIS 6.

Creating watersheds with specific drainage outlets

As an example, we'll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We'll assume our starting data includes an elevation raster called "dem" and a line vector called "roads". We first create the regular hydrology layers.

Preparation

#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem

# threshold in map cells (try 10000 as a start). Note: DEM Sink-filling not needed, since GRASS GIS uses the A* least cost routing, see r.watershed manual:
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>

So far, pretty straightforward. There's abundant information on r.watershed. I'll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10x10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

Display first results

When the process finishes we'll have three new raster maps: the flow direction map, the streams and the catchments. Let's see what we've got so far:

#Convert the steams and catchments to vectors
# r.to.vect in=catch out=catchments feature=area        # Grass 6
r.to.vect in=catch out=catchments type=area             # Grass 7

# the stream raster usually requires thinning
r.thin in=str out=str_thin
# r.to.vect in=str_thin out=streams feature=line        # Grass 6
r.to.vect in=str_thin out=streams type=line             # Grass 7
r.colors dem col=elevation

# Make a hillshade raster for displaying "3D" 
# r.shaded.relief map=dem shade=dem_shade zmult=1.5     # Grass 6
r.relief input=dem output=dem_shade zscale=1.5          # Grass 7

# Now display layers
# d.mon x0          # Grass 6
d.mon start=wx0     # Grass 7
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2

Determine drainage points

Now we need to find all the points where streams cross roads. The v.overlay module does not deal with point vectors (hint: v.select does). Instead we use a trick in v.clean. When cleaning a line vector, all points where lines cross and no node exists are considered topological "errors" and can be saved to a new point vector. So by merging the roads and streams vectors, we create a vector with lines (streams) crossing other lines (roads) without a node. Then we run v.clean, and we get all those intersection points in a new vector.

# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points

#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12

# Save crossing points to a text file
# v.out.ascii in=cross_points out=cross_points.txt format=point fs=space                   # Grass 6
v.out.ascii in=cross_points out=cross_points.txt format=point layer=-1 separator=space     # Grass 7

Looping thru drainage outlet points

Once we have the crossing points in a file, we simply run r.water.outlet in a loop to create a watershed for each cross point. However the raster result of r.water.outlet has value '1' in each cell that is upstream of the drainage point, and '0' everywhere else. For our purposes, we want to patch the rasters together after running the loop, so we need to have null values outside of the watersheds, and each watershed must use a different value in the upstream cells for its drainage point. To achieve these results, we use the r.null module to set '0' value cells to null. Then, we take advantage of the r.reclass function to make a reclassed raster with different values for each watershed. Here's how it works for GRASS 6:

i=0   #an iterator to give consecutive names to the new watersheds
while read X Y; do \
    i=$(( ${i} + 1 ))  \
    r.water.outlet drain=fdir east=$X north=$Y basin=tmp_$i  \
    r.null tmp_$i setnull=0 \
    echo 1=$i | r.reclass in=tmp_$i out=tmp_reclass_$i
done < cross_points.txt
echo "Created $i watersheds"

And here's how it works for GRASS 7:

i=0   #an iterator to give consecutive names to the new watersheds
while read X Y; do
    i=$(( ${i} + 1 ))
    r.water.outlet --overwrite input=$drain coordinates=$X,$Y output=tmp_$i 
    r.null tmp_$i setnull=0
    echo 1=$i | r.reclass --overwrite in=tmp_$i out=tmp_reclass_$i rules=- 
done < cross_points.txt
echo "Created $i watersheds"

Combining watersheds into one patched vector

Next we patch together all the reclassed rasters (watersheds), convert to vector and clean the merged watersheds vector.

# r.patch in=`g.mlist rast pattern=tmp_reclass* separator=,` out=wshed_patch     # Grass 6
# r.to.vect in=wshed_patch out=wshed_patch feature=area                          # Grass 6
r.patch in=`g.list raster pattern=tmp_reclass* separator=,` out=wshed_patch      # Grass 7
r.to.vect in=wshed_patch out=wshed_patch type=area
# Use v.clean to remove tiny areas (that were a string of single cells in the raster)
v.clean wshed_patch out=wshed_final tool=rmarea thresh=150

Choose an appropriate threshold value based on your region resolution. With a region resolution of 10, each individual cell will be 100 sqm, so choosing 150 as the threshold for v.clean allows removing these small areas. Additional manual cleaning may be required.

Clean up tmp rasters:

# g.remove rast=`g.mlist pattern=tmp* sep=,`                          # Grass 6
g.remove -f type=raster name=`g.list raster pattern=tmp* sep=,  `     # Grass 7

Most likely we'll want to calculate the area for each watershed.

# Add a table with a column for area in sq.km.
# v.db.addcol map=wshed_final col="area_sqkm double"      # Grass 6
v.db.addcolumn map=wshed_final col="area_sqkm double"     # Grass 7

# Use unit=k(ilometers) to get area in sq. km.
v.to.db map=wshed_final option=area col=area_sqkm unit=k

And finally, we can view the catchments, and their area values (you may use the wxGUI):

d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange

Calculating total drainage area for each stream reach

Some other software, such as the ArcGIS extension "ArcHydro" can calculate the total upstream drainage area flowing through each stream reach. The same result is achieved in GRASS with a few hydrology modules and a database join. Begin by rerunning the r.watershed command as in the first section, but adding a flow accumulation raster as the output:

#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem
# Choose threshold in map cells (try 10000 as a start). Note: DEM Sink-filling not needed:
r.watershed elev=dem accumulation=facc thresh=<your-threshold>

Now create a stream vector layer (or use the vector layer created above) with the r.stream.extract module. Then add to the vector attribute table columns for end point coordinates of each section of the stream, and for the flow accumulation area:

r.stream.extract --o elev=dem accum=facc thresh=<your-threshold> stream_vector=streams
v.db.addcolumn streams col="end_x double, end_y double, accum_area double"
v.to.db streams option=end columns="end_x,end_y"

We now have a stream vector map with the X-Y coordinates of the end point of each section of the stream, and an empty column ready to accept the upstream drainage area. The next step involves exporting the end point coordinates to a point layer, and using the v.what.rast module to get the values at each point from the flow accumulation grid. Then the end_points vector will have an attribute column containing the total number of upstream pixels flowing thru each point.

v.db.select -c streams column="cat,end_x,end_y" | v.in.ascii --o input=- output=end_points x=2 y=3
v.db.addcolumn end_points column="accum_pixels double, accum_area double"
v.what.rast end_points raster=facc column=accum_pixels

The following queries the current region setting for the cell resolution. The total number of upstream cells (from the flow accumulation raster) must be multiplied by the cell size in order to get upstream area in map units. Finally an SQL expression updates the total upstream area to the streams attribute table, using the number of pixels from the end_points multiplied by the area of each pixel.

eval `g.region -g`
SQ_M=$(( ${ewres}*${nsres} ))
v.db.update end_points column=accum_area query_col="accum_pixels*${SQ_M}"
db.execute sql="UPDATE streams SET accum_area=(SELECT accum_area FROM end_points where cat=streams.cat)"

This leaves the accum_area column populated with values for total upstream drainage area for each stream reach.

References

• Jasiewicz, J., & Metz, M. (2011). A new GRASS GIS toolkit for Hortonian analysis of drainage networks. Computers & Geosciences, 37(8), 1162-1173. (DOI)
• Metz, M. (2015). A new least impact approach for hydrological conditioning of DEMs. Geomorphometry for Geosciences, 95-96. (PDF)
• Metz, M., Mitasova, H., & Harmon, R. S. (2011). Efficient extraction of drainage networks from massive, radar-based elevation models with least cost path search. Hydrol. Earth Syst. Sci., 15, 667-678, 2011 (DOI)
• Metz, M., Mitasova, H., & Harmon, R. S. (2010). Accurate stream extraction from large, radar-based elevation models. Hydrology & Earth System Sciences Discussions, 7(3). (DOI)