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:
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
r.to.vect in=wshed_patch out=wshed_patch feature=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=,`

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"

# 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