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. I'd like to demonstrate how to get similar results with GRASS GIS.

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
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
# the stream raster usually requires thinning
r.thin in=str out=str_thin
r.to.vect in=str_thin out=streams feature=line
r.colors dem col=elevation
# Make a hillshade raster for displaying "3D"
r.shaded.relief map=dem shade=dem_shade zmult=1.5
# Now display layers
d.mon start=x0
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. Instead we use a trick in v.clean. When cleaning a line vector, all points where lines intersect 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, then running v.clean, we get all the intersection points in a new layer.

# 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

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.

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=wshed_$i  \
    r.to.vect -s in=wshed_$i out=wshed_$i feature=area \
    g.remove rast=wshed_$i
done < cross_points.txt
echo "Created $i catchments"

Combining watersheds into one patched vector

Next we patch together all the area vectors and clean the merged watershed vector.

v.patch in=`g.mlist vect pattern=wshed_* separator=,` out=wshed_patch
v.clean wshed_patch out=wshed_clean tool=snap,rmdupl thresh=20

Choose an appropriate snap threshold value based on your region resolution. Additional manual cleaning may be required. After cleaning, reset the cat values (only for centroids):

v.category in=wshed_clean out=wshed_clean2 opt=del
v.category in=wshed_clean2 out=wshed_final type=centroid opt=add
g.remove vect=wshed_clean,wshed_clean2

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.addtable 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:

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

Creating watersheds

Contents