--- title: "Watershed delineation in Python with pysheds: the River Derwent above Ladybower" description: "Delineate a watershed from a DEM in Python with pysheds: fill depressions, D8 flow, accumulation, pour-point snapping, Strahler order and an NRFA check." url: https://www.technicalgis.com/tutorials/watershed-delineation-python/ author: TechnicalGIS Team publisher: TechnicalGIS (https://www.technicalgis.com/about/) published: 2026-10-17 updated: 2026-10-17 language: en-GB --- # Watershed delineation in Python with pysheds: the River Derwent above Ladybower > Delineate a watershed from a DEM in Python with pysheds: fill depressions, D8 flow, accumulation, pour-point snapping, Strahler order and an NRFA check. - Level: Intermediate - Tested with: Python 3.12 - Source code: [Download scripts](https://www.technicalgis.com/code/watershed-delineation-python/) ## Key takeaways 1. A catchment is five calls in pysheds: fill\_pits, fill\_depressions, resolve\_flats, flowdir and catchment. The hard parts are a good DEM and a pour point on the right stream cell. 2. From the free 30 m Copernicus GLO-30 DEM, the River Derwent above Yorkshire Bridge drains 126.38 km², 0.3% more than the official NRFA area of 126.0 km², with an IoU of 0.988 against the NRFA boundary. 3. Always snap the pour point to the flow-accumulation network. A gauge a single cell off the river can give a catchment of a fraction of a square kilometre instead of 126. 4. The stream network is a choice, not a measurement: channel heads at 0.5 km² give 109.5 km of streams and Strahler order 4; at 0.1 km² the same valley has order 5 and twice the drainage density. ## Before you start - Comfortable with NumPy and rasterio; the rasterio terrain analysis tutorial covers DEMs, hillshades and coordinate systems - Python 3.10+ with pysheds, numba, rasterio, geopandas, shapely, pyproj, pandas, matplotlib, pystac-client, planetary-computer and requests - An internet connection: the DEM comes from the Microsoft Planetary Computer and the reference boundaries from UKCEH and the Environment Agency, with no accounts or keys **Short answer:** to delineate a watershed in Python, load a DEM with `pysheds`, then **fill pits and depressions** and **resolve flats** so every cell can drain. Compute **D8 flow direction** and **flow accumulation**, and **snap your outlet** to the nearest cell with large accumulation. Then call `grid.catchment()` for that point and convert the mask to a polygon with `rasterio.features.shapes`. For the stream network, keep cells whose upstream area exceeds a threshold (here 0.5 km²) and run `grid.stream_order()` for Strahler order. Always check the result against an official boundary or a published catchment area if one exists. We delineate a real, well-documented catchment: the **River Derwent above Yorkshire Bridge** in the Peak District, just below Ladybower dam. This is the catchment of Howden, Derwent and Ladybower reservoirs, the same reservoirs we monitored from space in the [Sentinel-2 reservoir water quality tutorial](https://www.technicalgis.com/tutorials/reservoir-water-quality-sentinel-2/). The UK National River Flow Archive (NRFA) gauges the river there as station 28001 and publishes its catchment boundary, so we can measure how close a free 30 m DEM gets. You will learn how to: 1. get a gauge location and the official catchment from the **NRFA** web services; 2. download the **Copernicus GLO-30** DEM and reproject it to British National Grid; 3. **condition** the DEM: fill pits and depressions, and resolve flats; 4. compute **D8 flow direction** and **flow accumulation**; 5. **snap the pour point** and delineate the **catchment polygon** and its area; 6. extract a **stream network** with **Strahler order** and drainage density, and test the threshold; 7. compare the result with the **NRFA boundary** and the **Environment Agency’s Water Framework Directive** catchments, and export a GeoPackage. All outputs below are real results from [`watershed_delineation.py`](https://www.technicalgis.com/code/watershed-delineation-python/), run on 14 October 2026 with Python 3.12, pysheds 0.5, numpy 2.5.4, rasterio 1.5.2 and geopandas 1.2.0. ## Step 1: the gauge and the official answer Start with the answer you want to reproduce. The NRFA station-information API returns the gauge’s grid reference and catchment descriptors as JSON, and its GeoServer serves the official catchment boundary over WFS: `watershed_delineation.py` ```python STATION = 28001 # River Derwent at Yorkshire Bridge, below Ladybower dam r = requests.get("https://nrfaapps.ceh.ac.uk/nrfa/ws/station-info", params={"station": STATION, "format": "json-object", "fields": "all"}, timeout=60) info = r.json()["data"][0] gx, gy = info["easting"], info["northing"] r = requests.get("https://hydrologygeoserver.ceh.ac.uk/geoserver/ows", params={"service": "WFS", "version": "2.0.0", "request": "GetFeature", "typename": "NRFA:NRFA_catchments", "CQL_FILTER": f"ID={STATION}", "outputformat": "application/json", "srsname": "EPSG:27700"}, timeout=60) nrfa = gpd.read_file(io.StringIO(r.text)) # the script caches both responses in data/ ``` `Output — NRFA station 28001` ```text === NRFA station 28001: Derwent at Yorkshire Bridge === gauge at British National Grid 419810 E, 385090 N | catchment area 126.0 km² (IHDTM 126.71 km²) | altitude 160.0–631.2 m | drainage density 1.138 km/km² NRFA boundary polygon: 126.70 km² ``` The NRFA publishes two areas: 126.0 km² is the official catchment area, and 126.71 km² is the area derived from the 50 m Integrated Hydrological DTM (IHDTM), which is also the area of the boundary polygon. ## Step 2: download and reproject the DEM We use **Copernicus GLO-30**, a free global DEM with 1 arc-second (about 30 m) cells, from the Planetary Computer. Hydrology needs square cells in metres, so we reproject it to British National Grid (EPSG:27700) at 30 m for a 25 km × 21 km box around the catchment: `watershed_delineation.py` ```python RES = 30.0 X0, Y0, X1, Y1 = 402_000, 381_000, 427_000, 402_000 # 25 km x 21 km of the Dark Peak in BNG lon0, lat0, lon1, lat1 = Transformer.from_crs(27700, 4326, always_xy=True).transform_bounds(X0, Y0, X1, Y1) cat = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", modifier=planetary_computer.sign_inplace) items = list(cat.search(collections=["cop-dem-glo-30"], bbox=[lon0, lat0, lon1, lat1]).items()) for item in items: with rasterio.open(item.assets["data"].href) as src: reproject(rasterio.band(src, 1), part, dst_transform=from_origin(X0, Y1, RES, RES), dst_crs="EPSG:27700", dst_nodata=np.nan, resampling=Resampling.bilinear) grid = Grid.from_raster(str(dem_path)) dem = grid.read_raster(str(dem_path)) ``` `Output — DEM` ```text === DEM === 833 x 700 cells of 30 m (EPSG:27700) | elevation 129.7–634.5 m | missing cells 0 ``` Make the box comfortably larger than the catchment. If the divide touches the edge of the DEM, part of the catchment is cut off and nothing warns you. ![A hillshaded elevation map of the Dark Peak, coloured from blue in the valleys (about 150 m) through green and yellow to brown and white on the moorland tops (over 600 m). A dashed black line outlines the NRFA catchment, which covers the upper Derwent and Ashop valleys. Howden, Derwent Upper and Ladybower reservoirs are drawn in blue down the middle of the catchment, and a red triangle marks the Yorkshire Bridge gauge at its southern tip.](https://www.technicalgis.com/_astro/01-dem-pour-point.Cthuo8ww_Y2lsQ.webp) _The Copernicus GLO-30 DEM with the NRFA catchment boundary for station 28001 and the three reservoirs from the Environment Agency's WFD lake layer._ ## Step 3: condition the DEM Real DEMs contain closed hollows: single-cell pits from noise, and larger depressions from errors, bridges, dams and genuine features such as quarries. Water flowing into a hollow has nowhere to go, so the flow network breaks. pysheds fixes this in three steps: `watershed_delineation.py` ```python pit_filled = grid.fill_pits(dem) pits = np.asarray(pit_filled) > np.asarray(dem) # copy now: the next step can modify in place flooded = grid.fill_depressions(pit_filled) inflated = grid.resolve_flats(flooded) fill = np.asarray(flooded) - np.asarray(dem) ``` `Output — conditioning` ```text single-cell pits raised: 1,440 cells raised by depression filling (beyond the pits): 6,936 (1.2%) | max fill 60.2 m | volume 22.10 million m³ deepest fill at (415635, 381525), DEM 224.0 m raised to 284.2 m flat cells after filling: 14,869; resolve_flats adds a tiny gradient so they drain ``` Filling raised 8,376 cells. Most changes are under a few metres, but the deepest is 60 m at the bottom edge of the box, near the Hope cement works quarry, outside our catchment. Filling turns every depression into a flat lake surface, so `resolve_flats` then adds a tiny gradient across the 14,869 flat cells so that water drains towards the outlet instead of stopping. ![Left: a grey hillshade of the study area with scattered purple and red patches where cells were raised by filling, including a large patch along a valley in the north-west and dark blobs in the south. Right: a histogram of fill depth on a logarithmic axis from 0.001 to 100 m, with most cells raised by roughly 1 to 4 m and a separate spike near 0.05 m.](https://www.technicalgis.com/_astro/02-depression-filling.DsJjT_Oz_1sDz0l.webp) _Where and by how much the DEM was raised so that every cell drains. The reservoir outlines are drawn in blue: their surfaces are flat in the DEM and need resolve_flats._ ## Step 4: flow direction and flow accumulation **D8** sends all the water in each cell to its steepest downhill neighbour. pysheds encodes the eight directions as powers of two, and the `dirmap` argument says which code means which direction (here the ESRI convention). **Flow accumulation** counts how many cells drain through each cell; multiplied by the cell area, it gives the upstream area: `watershed_delineation.py` ```python DIRMAP = (64, 128, 1, 2, 4, 8, 16, 32) # N, NE, E, SE, S, SW, W, NW (ESRI convention) fdir = grid.flowdir(inflated, dirmap=DIRMAP) acc = grid.accumulation(fdir, dirmap=DIRMAP) acc_km2 = np.asarray(acc) * RES ** 2 / 1e6 ``` `Output — D8` ```text largest upstream area in the grid: 196.3 km² flow directions: E 11.9%, SE 12.3%, S 15.0%, SW 12.0%, W 11.1%, NW 9.6%, N 14.6%, NE 13.5% ``` The largest value, 196.3 km², is where the Derwent leaves the bottom of our box after the River Noe joins it. If you use a different `dirmap` in different calls, the results are garbage without any error, so define it once and pass it everywhere. ![Left: a close-up of the D8 flow direction grid around the gauge, with each cell coloured by direction and black arrows showing where water flows, converging on the river channel in the middle. Right: flow accumulation for the whole area over a hillshade, with small yellow channels on the hillsides growing into dark blue rivers in the valleys, and a red box marking the close-up.](https://www.technicalgis.com/_astro/03-flow-direction-accumulation.D5tthbMN_1wayil.webp) _Left: D8 flow directions in the red box around the gauge. Right: upstream area on a logarithmic scale; cells draining less than 0.05 km² are hidden._ ## Step 5: snap the pour point A gauge’s published coordinates rarely fall exactly on the DEM’s river cell. One cell to the side, and you get the catchment of a hillside a few hundred metres square. The standard fix is to **snap** the point to the nearest cell whose upstream area is above a threshold suited to the river: `watershed_delineation.py` ```python SNAP_KM2 = 50 # snap to a cell that drains at least 50 km² sx, sy = grid.snap_to_mask(acc > SNAP_KM2 * 1e6 / RES ** 2, (gx, gy)) ``` `Output — pour point` ```text gauge cell drains 126.37 km² before snapping snapped to (419790, 385080), 22 m from the gauge; that cell drains 126.4 km² ``` Here the gauge already sits on the right cell and snapping moves it only 22 m. That is not always the case. Use a threshold well below the expected catchment area but well above the tributaries next to the outlet. If it is too low, the point can snap to a side stream. ## Step 6: delineate the catchment With a conditioned flow-direction grid and a good pour point, the catchment is one call. We turn the mask into a polygon and compare it with the NRFA boundary: `watershed_delineation.py` ```python catch = grid.catchment(x=sx, y=sy, fdir=fdir, dirmap=DIRMAP, xytype="coordinate") cmask = np.asarray(catch).astype(bool) area = cmask.sum() * RES ** 2 / 1e6 polys = [shape(g) for g, v in features.shapes(cmask.astype("uint8"), mask=cmask, transform=grid.affine, connectivity=8)] ours = gpd.GeoDataFrame({"area_km2": [round(area, 2)]}, geometry=[shapely.union_all(polys)], crs=27700) ref = nrfa.geometry.iloc[0] inter, union = ours.geometry.iloc[0].intersection(ref).area, ours.geometry.iloc[0].union(ref).area ``` `Output — catchment` ```text delineated area 126.38 km² vs NRFA 126.0 km² (+0.3%); vs NRFA boundary polygon 126.70 km² IoU with the NRFA boundary 0.988 | only in ours 0.61 km² | only in NRFA 0.94 km² min 10% median 90% max GLO-30 (m) 162.9 285.0 422.8 525.8 629.1 NRFA (m) 160.0 279.3 424.7 528.7 631.2 mean elevation 417 m (NRFA ALTBAR 416 m) longest flow path 19.9 km; mean flow path 11.34 km (NRFA DPLBAR 11.21 km) ``` That is a very good match. The area is within 0.3% of the official figure and the IoU with the official boundary is 0.988. The elevation percentiles are within 6 m and the mean flow-path length (DPLBAR) within 1.2%. The disagreements are slivers along the divide, where the moorland tops are so flat that 30 m and 50 m DEMs put the watershed in slightly different places. ![Left: the delineated catchment in green over a hillshade, with thin red slivers where only the GLO-30 catchment extends and thin blue slivers where only the NRFA boundary extends, all along the edge. Right: the Environment Agency's river water-body catchments in different colours inside our black catchment outline, labelled Derwent from Source to Westend, Westend Catchment, Alport Catchment, Ashop from Source to Alport, Ashop from Alport to Derwent and Highshore Clough Catchment, each 99 to 100 percent inside. The large Derwent from Westend to Wye water body extends beyond the outlet to the south-east.](https://www.technicalgis.com/_astro/04-catchment-vs-nrfa-wfd.BjqFrV2f_28Hzl4.webp) _Left: our catchment against the NRFA boundary (IoU 0.988). Right: the Environment Agency WFD river water-body catchments and our boundary._ ### A second check: Environment Agency WFD catchments The Environment Agency divides England’s rivers into Water Framework Directive (WFD) water bodies, each with its own catchment, on the Catchment Data Explorer. Six of them lie upstream of Ladybower dam: `Output — EA WFD catchments` ```text id name kind area_km2 share_inside GB104028057960 Derwent from Source to Westend river 20.56 0.995 GB104028057950 Westend Catchment (trib of Derwent) river 14.92 0.999 GB104028057930 Ashop from Source to Alport river 27.82 0.990 GB104028057940 Alport Catchment (trib of Ashop) river 11.34 0.992 GB104028057910 Ashop from Alport to Derwent river 10.19 0.992 GB104028057900 Highshore Clough Catchment (trib of Derwent) river 9.84 0.994 GB104028057880 Derwent from Westend to Wye river 142.84 0.219 GB30432299 Howden Reservoir lake 0.62 1.000 GB30432359 Derwent Upper Reservoir lake 0.70 1.000 GB30432459 Ladybower Reservoir lake 1.97 1.000 river water bodies more than half inside: 6, 94.67 km² together; IoU with our catchment 0.740 ``` Every upstream water body is 99–100% inside our boundary. They add up to only 94.67 km², though, because the main-stem Derwent water body runs from Westend all the way down to the River Wye (142.84 km²). Only 22% of it is above the gauge, and it fills the green gap down the middle of the right-hand map. WFD water bodies are management units, not gauge catchments, so they make a good consistency check but not a replacement for the NRFA boundary. ## Step 7: stream network and Strahler order A stream network from a DEM is simply the cells whose upstream area exceeds a **channel initiation threshold**. pysheds then traces the branches as lines and computes Strahler order: `watershed_delineation.py` ```python STREAM_KM2 = 0.5 # a channel starts where 0.5 km² drains to a cell BOOL_VF = ViewFinder(affine=grid.affine, shape=grid.shape, crs=grid.crs, nodata=False) # for mask rasters streams = Raster((acc_km2 > STREAM_KM2) & cmask, viewfinder=BOOL_VF) # pysheds also labels the cell just past the outlet (order + 1), so keep only stream cells order = np.where(np.asarray(streams), np.asarray(grid.stream_order(fdir, streams, dirmap=DIRMAP)), 0) branches = grid.extract_river_network(fdir, streams, dirmap=DIRMAP) ``` Each branch takes the most common order of its cells, and the script builds a GeoDataFrame of lines with lengths: `Output — stream network` ```text channel heads at 0.5 km²: 127 segments, 109.5 km, highest Strahler order 4 segments length_km strahler 1 64 48.1 2 44 37.6 3 17 22.2 4 2 1.5 drainage density 0.867 km/km² (NRFA DRAINDENS 1.138) ``` The upper Derwent and the Ashop are both order 3. They join in Ladybower, only 1.5 km above the dam, which is why order 4 is so short. ![The delineated catchment outline over a hillshade, with streams drawn in four colours and widths: thin purple first-order headwaters, green second-order streams, thick orange third-order rivers along the Derwent and Ashop valleys, and a short dark red fourth-order reach from their confluence in Ladybower to the dam. Pale blue reservoirs lie under the main channels.](https://www.technicalgis.com/_astro/05-strahler-order.Cmxp2QwQ_15olce.webp) _Streams starting at 0.5 km² of upstream area, coloured by Strahler order, with the lengths of each order in the legend._ ### How much does the threshold matter? A lot. The NRFA’s drainage density of 1.138 km/km² comes from a 1:50,000 river map, not from a DEM, and nothing in the DEM tells you where a real channel starts. Re-running with different thresholds shows the effect: `Output — threshold sensitivity` ```text threshold_km2 drainage_density max_strahler 0.10 1.652 5 0.25 1.155 5 0.50 0.867 4 1.00 0.665 4 2.00 0.497 3 5.00 0.362 3 ``` A threshold of about 0.25 km² reproduces the NRFA drainage density in this valley. The highest Strahler order changes from 3 to 5 across the range, so **always report the threshold** with any stream order or drainage density. A better approach for a real study is to calibrate the threshold against mapped channel heads or a river layer such as OS Open Rivers. ![Left: drainage density falling from 1.65 to 0.36 km per square kilometre as the channel threshold rises from 0.1 to 5 square kilometres on a logarithmic axis, with each point labelled with its maximum Strahler order, and a dashed red line at the NRFA value of 1.138 crossing near 0.25. Right: a hypsometric curve of elevation against percentage of catchment area above it, falling from 630 m to 160 m, with red NRFA altitude percentiles sitting on the curve.](https://www.technicalgis.com/_astro/06-threshold-hypsometry.Bfqq8eLS_KW1aY.webp) _Left: drainage density and maximum Strahler order for six channel thresholds. Right: the catchment's hypsometric curve from GLO-30 with the NRFA altitude percentiles._ ## Step 8: export The script writes the catchment, streams, both reference layers and the pour points to one GeoPackage, plus the flow-accumulation raster in km²: `watershed_delineation.py` ```python ours.to_file(gpkg, layer="catchment") net.to_file(gpkg, layer="streams") nrfa.drop(columns=["id"], errors="ignore").to_file(gpkg, layer="nrfa_catchment") wfd.to_file(gpkg, layer="ea_wfd_waterbodies") ``` `Output — export` ```text wrote derwent_catchment.gpkg (layers: catchment, streams, nrfa_catchment, ea_wfd_waterbodies, pour_point) and flow_accumulation_km2.tif ``` The NRFA layer has both `ID` and `id` columns, and GeoPackage field names are case-insensitive, so writing it unchanged raises `FieldError: Error adding field 'ID'`. Dropping one of them fixes it. ## Pitfalls and limitations - **GLO-30 is a surface model.** Copernicus GLO-30 is a DSM: it includes trees and buildings. In the open moorland of the Dark Peak that hardly matters, but in wooded or urban catchments, blocks of forest and buildings act as dams. Use a bare-earth DTM where you can. For England, the Environment Agency’s 1 m lidar composite DTM is free from the [Defra Survey Data Download](https://www.technicalgis.com/data-sources/united-kingdom/#src-defra-survey-data-download). Resample it to 5–10 m first, or the flow grid becomes huge and dominated by walls and ditches. - **Filling is not always right.** Depression filling assumes every hollow is an error. Real sinks, such as quarries, karst and lakes with no outlet, get filled too, and in flat areas filled DEMs can send water along the wrong route. Breaching (carving through the obstruction), as in WhiteboxTools’ `BreachDepressionsLeastCost`, often gives more realistic paths, especially where roads and embankments cross valleys. - **Reservoirs and dams.** A dam is a wall across the valley in the DEM, and the reservoir surface is a big flat. Here the gauge is below the dam so it works, but if you put a pour point in the middle of a reservoir, check that the flats were resolved towards the real outlet. - **Natural versus real catchment.** The NRFA notes that the Derwent has diversions within the basin (Ashop), imports (water piped in from the River Noe, outside the catchment) and public water supply exports. A DEM gives the topographic catchment only. The water arriving at the gauge comes from a slightly different area, and the flow is dominated by Ladybower releases. - **D8 limitations.** D8 cannot split flow, so on smooth slopes it draws parallel straight lines and it overestimates convergence on ridges. That does not affect catchment boundaries much, but for wetness indices or erosion modelling, use D-infinity or MFD (`grid.flowdir(..., routing="dinf")` in pysheds). - **Resolution at the divide.** On flat moorland plateaus the exact watershed is uncertain by tens to hundreds of metres. The 0.61 km² and 0.94 km² slivers in Step 6 are that uncertainty, not errors in either source. - **pysheds and NumPy.** pysheds 0.5 needs the `np.in1d = np.isin` shim with NumPy 2.4 or later, and it uses numba, so the first run compiles functions and takes longer. ## Frequently asked questions ### How do you delineate a watershed in Python? Load a DEM with `pysheds.grid.Grid.from_raster`, fill pits and depressions and resolve flats, then compute D8 flow direction and flow accumulation. Snap your outlet to a high-accumulation cell with `grid.snap_to_mask`, and call `grid.catchment(x, y, fdir)` to get a boolean catchment mask. Convert the mask to a polygon with `rasterio.features.shapes` and save it with GeoPandas. Steps 3 to 6 above show the full code. ### What DEM is best for watershed delineation? Use the best bare-earth DTM you can get at a sensible resolution for the catchment size. For catchments of tens to hundreds of square kilometres, a 10–30 m DEM is enough: Copernicus GLO-30, as here, is free worldwide and matched the official boundary at an IoU of 0.988. For small catchments, urban drainage or flat land, use lidar DTMs such as the Environment Agency’s 1 m composite in England or USGS 3DEP in the US, resampled to a few metres. ### What is the difference between pysheds and WhiteboxTools? pysheds is a pure-Python library (NumPy and numba) with a compact API for flow direction, accumulation, catchments and stream networks, which makes it easy to use in a notebook. WhiteboxTools (and its `whitebox-workflows` Python package) is a large compiled toolbox with hundreds of terrain tools, including depression breaching, more flow algorithms, and tools for watershed sets, basins and stream ordering. It is faster on large rasters. Both give very similar catchment boundaries for a well-conditioned DEM. ### Why is my delineated watershed tiny? Almost always because the pour point is not on the river cell in the DEM. A point one cell to the side of the channel drains only its own hillside. Snap the point to the flow-accumulation network first, as in Step 5, with a threshold appropriate for the river size. Also check that you use the same `dirmap` for `flowdir`, `accumulation` and `catchment`. ### How do I choose the flow accumulation threshold for streams? There is no universal value. Common starting points are 0.1–1 km² of upstream area for a 30 m DEM. Calibrate it by comparing drainage density with a mapped river network (here, 0.25 km² reproduces the NRFA value), or by checking channel heads on aerial imagery. Report the threshold, because stream order and drainage density depend directly on it. ### Where can I get official catchment boundaries in the UK? The National River Flow Archive publishes catchment boundaries for its gauging stations through a WFS and on each station page. The Environment Agency’s Catchment Data Explorer publishes WFD river and lake water-body catchments for England under the Open Government Licence. For Scotland, SEPA publishes river basin and water-body catchments. For the rest of the world, HydroSHEDS and its HydroBASINS layer give consistent nested basins from a 90 m DEM. ## Next steps - Look inside the reservoirs this catchment feeds in [monitoring reservoir water quality with Sentinel-2](https://www.technicalgis.com/tutorials/reservoir-water-quality-sentinel-2/), and relate colour to the peat cover upstream. - Add slope, aspect and hillshade maps of the catchment with [DEM terrain analysis in rasterio](https://www.technicalgis.com/tutorials/raster-terrain-analysis-rasterio/), or print it in 3D with [DEM to STL in Python](https://www.technicalgis.com/tutorials/dem-to-stl-python/). - Swap GLO-30 for a lidar DTM. Our [building footprints from LiDAR tutorial](https://www.technicalgis.com/tutorials/building-footprints-lidar-python/) shows how to build one from a classified point cloud. - Delineate every NRFA gauge in a region in a loop, and compare all the areas with the published ones using [GeoPandas](https://www.technicalgis.com/tutorials/geopandas-getting-started/). - Put the catchment and stream network on an interactive web map with [Folium](https://www.technicalgis.com/tutorials/interactive-web-maps-folium/). ## Data and licences - **DEM**: Copernicus DEM GLO-30, © DLR e.V. 2010–2014 and © Airbus Defence and Space GmbH 2014–2018, provided under COPERNICUS by the European Union and ESA; all rights reserved, free to use under the Copernicus DEM licence. Accessed from the [Microsoft Planetary Computer](https://planetarycomputer.microsoft.com/dataset/cop-dem-glo-30) on 14 October 2026. - **Gauging station data and catchment boundary**: National River Flow Archive, UK Centre for Ecology & Hydrology, station 28001 Derwent at Yorkshire Bridge. NRFA data are available under the [NRFA terms and conditions](https://nrfa.ceh.ac.uk/data-terms-and-conditions). Accessed on 14 October 2026. - **WFD water bodies**: Environment Agency Catchment Data Explorer, © Environment Agency, licensed under the Open Government Licence v3.0. Accessed on 14 October 2026. - Code: [`watershed_delineation.py`](https://www.technicalgis.com/code/watershed-delineation-python/), tested with Python 3.12.0, pysheds 0.5, numba 0.68.0, numpy 2.5.4, rasterio 1.5.2, geopandas 1.2.0, shapely 2.2.0 and pyproj 3.8.0. ## Key terms - **Watershed (catchment)**: All the land that drains to one point on a river, the pour point. British English usually says catchment; in American English a watershed is the same thing, while in British usage watershed can also mean the dividing ridge. - **Pour point**: The outlet cell for which a catchment is delineated, often a gauging station, a dam or a confluence. - **D8**: A flow-direction method that sends all the water from each cell to the steepest of its eight neighbours. Simple and robust, but it cannot split flow, so it draws parallel lines on smooth slopes. - **Flow accumulation**: The number of cells (or the area) upstream of each cell. Rivers are the cells with large values. - **Depression filling**: Raising closed hollows in a DEM to the level of their lowest outlet so that every cell can drain to the edge. Most hollows are errors in the DEM, but some are real, such as lakes and quarries. - **Strahler order**: A stream ranking where headwater channels are order 1, and two channels of the same order join to make the next order up. A channel joined by a lower-order stream keeps its order. - **Drainage density**: Total stream length divided by catchment area, in km/km². It depends heavily on where you start drawing channels. ## Data for this tutorial - [Copernicus DEM (GLO-30 / GLO-90)](https://www.technicalgis.com/data-sources/#src-copernicus-dem-glo-30-glo-90) · The 30 m DEM used here - [Microsoft Planetary Computer](https://www.technicalgis.com/data-sources/#src-microsoft-planetary-computer) · GLO-30 as Cloud Optimised GeoTIFFs via STAC - [Defra Survey Data Download](https://www.technicalgis.com/data-sources/united-kingdom/#src-defra-survey-data-download) · Environment Agency 1 m lidar DTM for England - [HydroSHEDS](https://www.technicalgis.com/data-sources/#src-hydrosheds) · Ready-made global basins and river networks Licences and formats for each source are listed in the [GIS data sources catalogue](https://www.technicalgis.com/data-sources/). --- Source: https://www.technicalgis.com/tutorials/watershed-delineation-python/