Watershed Analysis with GRASS

Arcview users, needing to delineate watersheds and stream networks, would do well to choose the extension called  “Arc Hydro” (requires, at least the Spatial Analyst extension from ESRI). 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. 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.

#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, 10×10 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.

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

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 “errors” and can be saved to a new point vector. So by merging the roads and streams vectors, we have lines (streams) crossing lines (roads) with no node at that location. Then we run v.clean to get all those 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

Once we have the crossing points in a file, we simply run r.water.outlet in a loop to create the watersheds for each cross point.  The module r.water.outlet sets the raster value to 1 in the watershed and 0 elsewhere. We’ll use r.null to force the raster to ‘null’ outside of the watershed. THen we convert each raster to a vector (polygon).

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.null -q wshed_$i setnull=0 \
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"

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"
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

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>