Get drainage basin of a polygon (with QGIS, GRASS or other FOSS)

Question

Using a polygon shapefile and a raster DEM, I'd like to find a way to get the drainage basin of a polygon area — that is, the area of all cells which drain into or through the cells contained in the polygon. This is essentially the function provided by r.water.outlet in GRASS except that it only looks at the drainage basin of a single cell.

I figure there must be some way to do this, if only by summing together the drainage basins of all cells within the polygon. However this would require a means of identifying and iterating through all cells contained in a polygon and running r.water.outlet with each cell as a unique input. Ideally I would be able to automate this process with python scripting. Ultimately, I'd like to run it for a large number of polygons in a shapefile layer.

Is this possible?

EDIT: This work relates to semi-automated assessment of sinkholes, karst features and groundwater susceptibility with GRASS GIS & other FOSS. To illustrate what I am hoping to accomplish here, I've taken some screenshots.

A slopeshade raster overlaid on a map of low-lying areas and depressions (shown in blue, derived with r.neighbors). You can see a number of sinkholes here.

Given a vector layer containing polygon features representing the outlines of sinkholes, I'd like to calculate the geometry of the basin that drains into the sinkhole, which will generally be a very small fragment of the larger watershed basin.

Here is a map showing drainage basins for each sinkhole, calculated manually with r.water.outlet

Here you can see a split basin. The sinkhole contains two cells which are separate drainage endpoints in the DEM. I want to get the 'complete' drainage basin for a sinkhole polygon feature, that is the composite of the drainage basin for all points within the sinkhole polygon.

Answer 1

If you need only the total drainage area (not the actual geography of the area) then you can do this with flow accumulation. Run r.watershed with your DEM as input, and define a name for the flow accumulation output raster. Each pixel contains the upstream area (actually number of cells) that flows into that pixel. Next you need to extract the boundary of your polygon.

Now run v.rast.stats with these parameters: the polygon boundary as the input map= parameter, the flow accumulation as raster= and set method=sum. This module calculates univariate statistics from the raster for each category in the vector map, and adds new attribute columns to the vector. By summing flow accumlation, each line segment in the boundary will have a new column with the total flow accum crossing the boundary. You'll then need to multiply by the cell resolution in order to get total area in the units you prefer.

Edit...

From your description, I understand that your polygon can overlay parts of an unknown number of basins. The polygon could be cutting across watershed divides. So to get the actual geometry of this collection of parts of basins, I don't see any alternative to collecting all the basins that flow into the polygon, and merging them together.

Here's a possible way to go about this:

You'll need a DEM that covers all the possible catchments that intersect your polygon. The GRASS addon r.stream.basins can take a point map as drainage points and calculate the basin for each outlet point. When you run r.watershed give an output map name for the streams, and flow direction ("drainage" parameter). Now run r.thin on the streams raster, and convert to vector (r.to.vect ... type=line). Next intersect the streams and boundary line vectors. The way to do this is a bit strange: patch the two line vectors together (v.patch ) and then v.clean ... error=intersect_pts. The error option creates a point vector map with all the intersection points. Now run r.stream.basins using the flow direction map as input and the intersection points as points= parameter. That should give you in the output all basins that drain into the original polygon.