Delineate Watersheds in PyQGIS
Every point in a landscape drains somewhere. A watershed — or catchment, or drainage basin — is all the land that drains to a given outlet: a gauging station, a culvert, a lake inlet, a bridge. Delineating it from an elevation model is a standard hydrological workflow: fill spurious depressions, compute the direction water flows from each cell, accumulate flow downstream, extract streams, and trace upstream from an outlet. QGIS runs all of it through the GRASS algorithms, and PyQGIS makes it repeatable for any number of outlets.
This recipe belongs to Terrain & Interpolation Analysis. It conditions a DEM, computes flow accumulation and drainage direction with r.watershed, extracts a stream network with a chosen threshold, snaps outlets onto streams, delineates basins with r.water.outlet, polygonizes them and checks their areas.
Prerequisites
- QGIS 3.34 LTR or newer, or the QGIS 4 series, with the GRASS Processing provider (included in the standard Windows and macOS installers; a separate package on some Linux distributions).
- A DEM in a projected CRS with metre units, covering the whole area that drains to your outlets — a catchment cut off by the DEM edge will be truncated.
- Outlet points, such as gauging stations, as a point layer.
Find the GRASS algorithms
The GRASS provider's algorithm ids have a version-dependent prefix. A small helper finds them by name so scripts work on every release.
import processing
from qgis.core import QgsApplication
def grass(name):
for alg in QgsApplication.processingRegistry().algorithms():
if alg.id().split(":")[-1] == name and alg.provider().id().startswith("grass"):
return alg.id()
raise RuntimeError(f"GRASS algorithm {name} not available - is the GRASS provider enabled?")
R_WATERSHED = grass("r.watershed")
R_OUTLET = grass("r.water.outlet")
R_FILL = grass("r.fill.dir")
print(R_WATERSHED, R_OUTLET, R_FILL)
Breakdown: Matching on the part after the colon finds grass:r.watershed or grass7:r.watershed alike. Raising a clear error when GRASS is missing saves confusion later — the most common reason these workflows fail on a new machine is a QGIS installed without the GRASS provider. processing.algorithmHelp(R_WATERSHED) lists parameter names, which follow GRASS's own option names.
Condition the DEM
Real DEMs contain small artificial depressions — from interpolation, from bridges and culverts that the DEM does not represent — where water would pool instead of flowing on. Filling them makes every cell drain to the edge or an outlet.
from qgis.core import QgsRasterLayer
dem = "/data/terrain/dem_10m.tif"
filled = processing.run(R_FILL, {
"input": dem,
"format": 0,
"output": "/data/hydro/dem_filled.tif",
"direction": "/data/hydro/fill_direction.tif",
"areas": "/data/hydro/problem_areas.tif",
})["output"]
print("filled DEM written:", filled)
Breakdown: r.fill.dir raises depression cells to their pour point so flow can continue, and writes the direction raster and a raster of problem areas it could not resolve. Filling is standard for small sinks but wrong for real depressions such as quarries and closed lake basins; check the problem areas and the difference between filled and original DEM where it is large. r.watershed itself handles depressions with a least-cost search and often works on unfilled DEMs, so compare results with and without filling for your terrain.
Compute accumulation and drainage
r.watershed computes, for every cell, the direction water leaves it and how many cells drain through it — flow accumulation, the basis for streams and catchments.
ws = processing.run(R_WATERSHED, {
"elevation": filled,
"threshold": 10000, # cells; 10,000 × 100 m² = 1 km²
"-s": True, # single flow direction (D8)
"-a": True, # positive accumulation only
"accumulation": "/data/hydro/accumulation.tif",
"drainage": "/data/hydro/drainage.tif",
"stream": "/data/hydro/streams.tif",
"basin": "/data/hydro/subbasins.tif",
})
print({k: v for k, v in ws.items() if isinstance(v, str)})
Breakdown: threshold is the minimum drainage area, in cells, for a stream and for subbasin generation; at 10 m resolution, 10,000 cells is 1 km². Single flow direction (-s) sends all flow to the steepest neighbour, which gives clean, connected streams; the default multiple-flow algorithm spreads flow and suits hillslope hydrology better. The drainage raster encodes direction and is what r.water.outlet needs. Subbasins delineated at the threshold are a useful by-product: every stream segment's own catchment.
Choose a stream threshold
The threshold decides how dense the stream network is. Too small and every gully is a stream; too large and real streams are missed. Comparing against a mapped river network is the best calibration.
for km2 in (0.25, 1, 4):
cells = int(km2 * 1e6 / 100) # 10 m cells
processing.run(R_WATERSHED, {
"elevation": filled, "threshold": cells, "-s": True, "-a": True,
"stream": f"/data/hydro/streams_{km2}km2.tif"})
print(f"threshold {km2} km² = {cells} cells")
Breakdown: Running a few thresholds and overlaying each stream raster on mapped rivers shows which matches the official network's density; that threshold is then the right one for this landscape. Steep, wet terrain supports channels at small drainage areas; flat or dry terrain needs larger ones. Convert the chosen stream raster to lines with r.to.vect or r.thin plus polygonizing for mapping.
Snap outlets to the stream network
An outlet point that lies a few metres off the computed stream sits on a hillslope, and its "catchment" is a tiny patch. Moving each outlet to the highest-accumulation cell nearby fixes it.
import numpy as np
from osgeo import gdal
from qgis.core import QgsVectorLayer, QgsPointXY
acc_ds = gdal.Open("/data/hydro/accumulation.tif")
acc = acc_ds.ReadAsArray()
gt = acc_ds.GetGeoTransform()
def snap(x, y, radius_m=50):
col, row = int((x - gt[0]) / gt[1]), int((y - gt[3]) / gt[5])
r = int(radius_m / gt[1])
window = acc[max(row - r, 0):row + r + 1, max(col - r, 0):col + r + 1]
i, j = np.unravel_index(np.argmax(window), window.shape)
row2, col2 = max(row - r, 0) + i, max(col - r, 0) + j
return gt[0] + (col2 + 0.5) * gt[1], gt[3] + (row2 + 0.5) * gt[5]
gauges = QgsVectorLayer("/data/hydro/gauges.gpkg", "gauges", "ogr")
snapped = {f["gauge_id"]: snap(*f.geometry().asPoint()) for f in gauges.getFeatures()}
print(list(snapped.items())[:3])
Breakdown: Converting the outlet's coordinates to a row and column, searching a small window and taking the cell with the largest accumulation moves the outlet onto the main channel. The radius should be small — tens of metres — or an outlet near a confluence may jump to the larger river. Returning cell-centre coordinates avoids edge ambiguity. Check snapped outlets against the stream raster on the map before delineating.
Delineate and polygonize basins
r.water.outlet traces every cell that drains to an outlet. Running it per outlet and polygonizing gives one catchment polygon per gauge.
results = []
for gid, (x, y) in snapped.items():
basin = processing.run(R_OUTLET, {
"input": "/data/hydro/drainage.tif",
"coordinates": f"{x},{y}",
"output": f"/data/hydro/basin_{gid}.tif"})["output"]
poly = processing.run("gdal:polygonize", {
"INPUT": basin, "BAND": 1, "FIELD": "basin", "EIGHT_CONNECTEDNESS": False,
"OUTPUT": "memory:"})["OUTPUT"]
for f in poly.getFeatures():
if f["basin"] == 1:
results.append((gid, f.geometry().area() / 1e6))
for gid, km2 in results:
print(f"{gid}: {km2:,.1f} km²")
Breakdown: r.water.outlet writes 1 for cells draining to the outlet and NoData elsewhere; polygonizing and keeping the polygon with value 1 gives the catchment outline. The area in square kilometres is the first check: gauging stations publish their catchment areas, and a delineated area within a few percent confirms the workflow; a large difference points to a misplaced outlet, a DEM edge cutting the basin, or flat areas where flow direction is ambiguous. The polygonize recipe covers cleaning staircase outlines.
Summarise each catchment
A catchment polygon becomes useful when it carries numbers: mean slope, land-cover shares, rainfall, the length of streams inside it. Zonal statistics over the basins add them in one step per raster.
basins = QgsVectorLayer("/data/hydro/catchments.gpkg", "catchments", "ogr")
for raster, prefix in (("/data/terrain/slope_deg.tif", "slope_"),
("/data/climate/rain_mm.tif", "rain_")):
processing.run("native:zonalstatisticsfb", {
"INPUT": basins, "INPUT_RASTER": raster, "RASTER_BAND": 1,
"COLUMN_PREFIX": prefix, "STATISTICS": [2, 4, 6], # mean, std dev, max
"OUTPUT": "memory:"})
Breakdown: The feature-based zonal statistics algorithm adds the chosen statistics as fields with the given prefix — mean, standard deviation and maximum here, by their positions in the algorithm's statistics list. Mean slope and rainfall per catchment are classic inputs for runoff estimates and for comparing gauged with ungauged basins. Writing the outputs to a GeoPackage alongside the catchment polygons gives one table describing every basin; zonal statistics covers statistics codes and performance.
QGIS version compatibility
r.watershed, r.water.outlet and r.fill.dir are available through the GRASS provider on QGIS 3.34 LTR, 3.40 LTR and QGIS 4 when GRASS is installed; the provider prefix varies, hence the lookup helper. GRASS option names such as elevation, threshold and drainage follow GRASS 8.
Troubleshooting
- The catchment is tiny. The outlet is off the stream; snap it to high accumulation.
- The catchment is cut by a straight line. The DEM does not cover the whole basin; extend it upstream.
- Streams appear in flat floodplains as parallel lines. Flow direction is ambiguous on flat areas; burn in mapped streams or use a coarser threshold.
- GRASS algorithms are missing. The GRASS provider is not installed or not enabled.
Conclusion
Find GRASS algorithms by name, fill spurious sinks while checking real depressions, compute accumulation and drainage with r.watershed, calibrate the stream threshold against mapped rivers, snap outlets onto the computed streams, trace and polygonize basins with r.water.outlet, and check areas against published catchment sizes.
Frequently Asked Questions
Can I do this without GRASS? SAGA and WhiteboxTools offer equivalent tools through their Processing providers; the steps are the same.
What resolution should the DEM have? Fine enough to resolve the channels that matter — 10–25 m for regional catchments, 1–2 m for urban drainage.
How do I handle culverts under roads? Burn them into the DEM by lowering cells along the culvert line before computing flow.
Can I get the stream network as lines?
Yes — thin the stream raster and convert to vector with r.thin and r.to.vect.