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.

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:

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.

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).