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.


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.


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.

  • Thanks for the suggestion Micha. I did not think to use the flow accumulation map, and this will be a useful way to get the area values if no other solutions come up. However I am hoping to also get the geometry of the basin. I am working on a GIS-based sinkhole identification and groundwater susceptibility analysis, and knowing the geometry of the drainage basin would allow for some interesting possibilities. – corvus Jan 10 '17 at 21:27
  • See edit above. – Micha Jan 11 '17 at 14:25
  • Micha, thanks for the follow-up. I've tried to go through these steps but have become a bit perplexed. Could you clarify what you mean when you say "intersect the streams and boundary line vectors"? I'm supposing you mean the polygon boundaries of the sinks for which I want to get drainage basins? If this is the case, I'm wondering, what would happen if the polygon did not intersect with a stream feature? (These are mostly quite small polygons, representing hydrologically-isolated areas in a massive LiDAR DEM dataset and generally not near streams. I've edited my OP with some clarification) – corvus Jan 15 '17 at 16:58
  • I was a bit brief in the answer. But I think You understood. You extract the boundary of your polygon with v.extract ... type=boundary, then patch that boundary with the streams, as above. If there are no intersection points of streams with the boundary, then there's no flow into that polygon. Wasn't that the original question? – Micha Jan 16 '17 at 7:16
  • 1
    What about r.stream.basins using the boundary as point inputs? That should capture all flow into this boundary. – mankoff Feb 19 '20 at 15:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.