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.
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
Best Answer
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 inputmap=
parameter, the flow accumulation asraster=
and setmethod=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 runr.watershed
give an output map name for the streams, and flow direction ("drainage" parameter). Now runr.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 thenv.clean ... error=intersect_pts
. The error option creates a point vector map with all the intersection points. Now runr.stream.basins
using the flow direction map as input and the intersection points aspoints=
parameter. That should give you in the output all basins that drain into the original polygon.