pysheds turns a digital elevation model into a delineated catchment in Python

Every step of the delineation is a call on the same grid object

The problem pysheds addresses is the middle of the hydrology pipeline. You have a digital elevation model, you want the boundary of the catchment that drains to a known point, and the tools that do this for you usually hand back a polygon without telling you what they did to the elevation data on the way. pysheds is a library for the version of that job where the intermediate grids matter.

The whole API hangs on a single object, Grid, and the steps read as a sequence of transformations on it. You build a grid from a raster file, then ask that same object to fill pits, fill depressions, resolve flats, compute flow direction, accumulate flow, and finally delineate a catchment. Because the grid is the same object throughout, the rasters stay aligned in extent and cell size, and there is no projection dance between steps.

The intended users are hydrologists, geomorphologists and GIS-adjacent Python developers working with digital elevation models. The repository carries a list of example notebooks for specific tasks, including crs_conversion, no_metadata, flowdir_vs_projection, laplacian, routing_modes, snap_to_mask, terrain, rfsm, hand, shapefiles, performance and quickstart, which is a fair picture of the questions the library is built to answer. The documentation lives at mdbartos.github.io/pysheds rather than in the README, which is itself a sign of how much there is to learn before the first run.

Conditioning the DEM comes before any flow direction is computed

The order of operations in the tutorial is not stylistic. Before any flow direction can be derived, the elevation surface has to be made hydrologically sensible, and pysheds exposes each conditioning step as its own call.

Pits are holes in the surface that trap water, so the first call fills them:

# Condition DEM
# ----------------------
# Fill pits in DEM
pit_filled_dem = grid.fill_pits(dem)

# Fill depressions in DEM
flooded_dem = grid.fill_depressions(pit_filled_dem)
    
# Resolve flats in DEM
inflated_dem = grid.resolve_flats(flooded_dem)

Three calls, each taking the previous result. Pits are local depressions, depressions are larger flat basins that need a consistent outlet, and flats are areas of equal elevation where a flow direction cannot be inferred from the surface at all. The last step is the one people skip by accident, and skipping it tends to produce a flow direction grid full of arbitrary values in flat terrain.

Skipping the first two is worse. An unconditioned DEM gives you a flow network with dead ends and spurious sinks, and a catchment boundary derived from that network will be wrong in ways that look plausible on a map. The cost is that you are working on a modified elevation surface, so any statistic you later compute on inflated_dem describes the conditioned raster rather than your original data. The library gives you both values, and which one you analyse is a decision you have to make.

dirmap is the contract that makes two flow direction grids comparable

The D8 method assigns each cell one of eight directions. The encoding of those eight directions is not standardised, so pysheds makes you pass the mapping explicitly as the dirmap argument, and every downstream call takes the same tuple. The tutorial uses:

# Determine D8 flow directions from DEM
# ----------------------
# Specify directional mapping
dirmap = (64, 128, 1, 2, 4, 8, 16, 32)
    
# Compute flow directions
# -------------------------------------
fdir = grid.flowdir(inflated_dem, dirmap=dirmap)

Those eight values are bit flags, and the tuple fixes which bit means which neighbour. This is the single most important detail in the library for anyone who will read results with another tool. A flow direction raster written with one encoding and read with another is silently wrong, and the failure does not look like an error, it looks like a plausible map of drainage. By making the mapping an explicit argument, pysheds forces the decision to be visible in the code rather than buried in a default, and it lets you keep several encodings side by side in one session.

The same convention carries into the accumulation step:

# Calculate flow accumulation
# --------------------------
acc = grid.accumulation(fdir, dirmap=dirmap)

Accumulation counts upstream cells for every cell in the network, and it is the intermediate that decides the pour point, the catchment extent and the stream network. The tutorial plots it with a logarithmic colour scale, which is the standard way to look at it, since the distribution of accumulation values spans orders of magnitude.

Pour points snap to a high accumulation cell before the catchment is drawn

Delineating a catchment means specifying a point and asking which cells drain to it. The point you have is rarely the point the raster actually wants, because the outlet in the flow network is a specific high accumulation cell, and handing that library a coordinate that is not on the network gives you either nothing or a surprising basin.

The tutorial fixes this with a snapping step:

# Delineate a catchment
# ---------------------
# Specify pour point
x, y = -97.294, 32.737

# Snap pour point to high accumulation cell
x_snap, y_snap = grid.snap_to_mask(acc > 1000, (x, y))

# Delineate the catchment
catch = grid.catchment(x=x_snap, y=y_snap, fdir=fdir, dirmap=dirmap, 
                       xytype='coordinate')

The threshold in acc > 1000 is the knob. Raise it and the snap moves to a cell further downstream with more upstream area, which can cross a drainage divide. Lower it and the outlet lands close to the coordinate you supplied. There is a dedicated example notebook, snap_to_mask.ipynb, because this is the step where results diverge from what a user expected, and the parameter is worth tuning against your own data rather than copying.

Once the catchment mask exists, the grid can be clipped to it, which limits everything downstream to the basin and cuts memory use for later steps:

# Crop and plot the catchment
# ---------------------------
# Clip the bounding box to the catchment
grid.clip_to(catch)
clipped_catch = grid.view(catch)

Note that clip_to changes the extent of the grid object itself, so anything you compute afterwards is computed on the clipped extent. If you need the full raster later, keep the original grid.

River networks come out of the accumulation raster, not out of the catchment

The stream network is derived from flow accumulation with a threshold, and the threshold is expressed in upstream cells, so its meaning depends on your raster resolution and on the flow direction encoding you chose:

# Extract river network
# ---------------------
branches = grid.extract_river_network(fdir, acc > 50, dirmap=dirmap)

The result comes back as a collection of branches whose features carry geometry coordinates, which the tutorial plots directly as line segments. The threshold of 50 upstream cells is a starting point, not a recommendation. Change the DEM resolution and the same number means a different physical distance; a coarser grid will treat a short reach as a channel, a finer one will miss it. Deciding that threshold against your own area of interest is the main judgement call in this step, and the fact that it is one argument means you can sweep it in a loop without touching the rest of the pipeline.

There is a related distance measurement in the library:

# Calculate distance to outlet from each cell
# -------------------------------------------
dist = grid.distance_to_outlet(x=x_snap, y=y_snap, fdir=fdir, dirmap=dirmap,
                               xytype='coordinate')

That call measures, for every cell, how far it is from the outlet along the flow paths, and the tutorial labels the resulting grid in cells rather than metres, which is a reminder that the output inherits the raster's cell units.

Together, extraction and distance cover the two things people actually want from a DEM: where the water goes, and how long it takes to get there in network terms. What pysheds does not do is produce a hydrologically correct stream order, and if that is your deliverable you will be adding it yourself.

rasterio, numba and pyproj decide what you can feed the library

setup.py lists the install requirements, and they describe the shape of the library better than the README does. The required set is affine, geojson, looseversion, numba, numpy, pandas, pyproj, rasterio>=1, scikit-image and scipy. Two of those deserve attention.

rasterio is the geospatial raster stack, and it is what makes Grid.from_raster and read_raster possible. The stated floor of rasterio>=1 matters when you are assembling an environment, because older rasterio and newer rasterio do not always agree on dataset handling, and pinning the wrong side of that boundary produces failures that look like pysheds bugs. The examples directory links a terrain raster, an impervious area archive and soil polygons, which suggests the library reads GeoTIFF input and expects a real geospatial raster with a CRS rather than a bare array of heights.

pyproj handles coordinates, and that is where most of the sharp edges are. The example notebooks include crs_conversion, no_metadata and flowdir_vs_projection, and their names describe the failure modes: rasters with no CRS attached, and the question of whether flow direction computed in a projected CRS differs from one computed in geographic coordinates. When you pass a pour point with xytype='coordinate', you are handing pysheds map coordinates, so the CRS of the raster and the CRS of your point have to agree. The example notebooks exist because this is where people get stuck, and the tutorial itself uses a longitude and latitude pair, so a geographic input is a valid path but not the only one.

numba is the third heavy dependency, used to compile the numeric kernels that do the fast work. It also means the first call on a new function carries compilation overhead, and that a mismatched numba or numpy combination in an existing environment can fail at import time rather than at the call you expected.

Release 0.5 shipped in 2025, and the code has moved since

The repository is not archived and the last push was on 2026-10-02. The published releases move slowly: 0.3.5 on 2023-05-16, 0.4 on 2024-05-09 and 0.5 on 2025-08-14, and setup.py still declares version 0.5. That pattern suggests a project with a stable public API and a low release cadence, where day to day changes land on master rather than on a tagged release. If you pin your environment, you are pinning to a release that is more than a year old, and if you track master you are tracking a moving target.

The version declared in setup.py is also behind the classifier list, which runs from Python 3.6 to 3.12, so the metadata describes a wider range than a fresh install necessarily targets. There is a pyproject.toml alongside setup.py that configures pytest rather than the build, with the test paths set to the tests directory, and the repository also carries tox.toml, a CI directory, a .travis.yml file and a CHANGELOG.md. The presence of both tox and CI configuration suggests the project treats compatibility across environments as a real concern, which matters for a library that must sit next to rasterio, numba and pyproj in someone else's environment.

The licence is GPL-3.0, and setup.py classifies it as GNU General Public License v3 or later. For internal analysis and research use this is not a constraint. If you intend to embed pysheds inside a proprietary product or a service you distribute, the copyleft terms of the GPL are the thing to examine, and the answer belongs with your legal team rather than with this article.

Where the library stops, and what you do instead

Three boundaries are worth being explicit about. First, this is a raster hydrology library, not a GIS. It does not manage vector topology, and it does not build a geodatabase. The optional recipes extra exists precisely because that work sits outside the library: it pulls in geopandas, ipython, matplotlib and seaborn, and the example notebooks for shapefiles and river networks do the vector and plotting work in the surrounding Python stack. Bringing your own geopandas workflow is the intended arrangement, not a workaround.

Second, the alternative to a dedicated GIS application is a different shape of work, not just a different shape of interface. A desktop or server GIS that offers catchment delineation as a tool does the whole chain for you in one operation, with a coverage format, a topological output and a data model you can rely on. The pysheds approach inverts that: you get the conditioned DEM, the direction grid, the accumulation raster and the mask, and you decide what a basin means. For a batch pipeline, a parameter sweep over thresholds, or a bespoke hydrological metric, that inversion is the point. For a one-off boundary that a colleague will open in a desktop GIS and edit, the coverage workflow is less work.

Third, the library is bound to a single flow method. The tutorial is D8, one of eight directions per cell, and that choice is baked into the direction encoding, the accumulation and the catchment walk. Other methods, such as multiple flow direction or a routing scheme that accounts for slope in the steepest descent choice, are a different implementation of the same idea, and the example notebooks named routing_modes and laplacian show that the surrounding project has been used to explore them. If your results have to be comparable to published work that used a specific routing method, check which one produced those numbers before assuming D8 is interchangeable.

Editorial conclusion

pysheds fits researchers and engineers who need reproducible catchment boundaries and stream networks from a DEM and want to inspect or modify each intermediate grid, and who are comfortable managing rasterio, numba and pyproj in the same environment. It is the wrong choice if you need a GUI, multi format coverage, or a tool that handles vector topology for you, because the library is deliberately narrow and leaves vector work to geopandas or shapely in your own code. Before you depend on it, check the CRS handling in examples/crs_conversion.ipynb against your own data, confirm the rasterio version in your environment satisfies the rasterio>=1 requirement, and read the GPL-3.0 terms if you plan to ship it inside a larger product.

Frequently asked questions

Official sources

1. Issues
2. License: GPL-3.0
3. pysheds/pysheds on GitHub
4. README
5. Releases