xjtulyc/hydrology-pysheds

Hydrological Analysis with pysheds

One-line summary: Derive watersheds, stream networks, and flow accumulation from Digital Elevation Models (DEMs) using pysheds, rasterio, and geopandas.

---

When to Use This Skill

• When delineating watershed boundaries from a DEM for a pour point
• When extracting stream networks from topographic data
• When computing flow direction and accumulation grids
• When estimating runoff and drainage basin area
• When conditioning DEMs to remove pits and flat areas
• When analyzing upstream contributing area for flood modeling

Trigger keywords: watershed delineation, DEM, flow direction, flow accumulation, stream network, pysheds, hydrological analysis, basin area, pour point

---

Background & Key Concepts

DEM Processing Pipeline

$$ \text{Raw DEM} \xrightarrow{\text{pit-fill}} \text{Conditioned DEM} \xrightarrow{\text{flow direction}} \text{D8 grid} \xrightarrow{\text{accumulation}} \text{Flow acc.} \xrightarrow{\text{threshold}} \text{Stream network} $$

D8 Flow Direction

Each cell flows to one of 8 neighbors based on steepest descent. Flow direction encoding: N=64, NE=128, E=1, SE=2, S=4, SW=8, W=16, NW=32.

Flow Accumulation

$A_{ij}$ = number of upstream cells draining through cell $(i,j)$. Cells with high accumulation are stream channels.

Watershed Delineation

Starting from a pour point (outlet), trace upstream through the flow direction grid to find all contributing cells.

SCS Curve Number Method

For runoff estimation:

$$ Q = \frac{(P - I_a)^2}{P - I_a + S}, \quad S = \frac{25400}{CN} - 254, \quad I_a = 0.2S $$

where $P$ is precipitation depth (mm), $CN$ is the SCS curve number.

---

Environment Setup

Install Dependencies

pip install pysheds>=0.3 rasterio>=1.3 geopandas>=0.14 \
            matplotlib>=3.7 numpy>=1.24 scipy>=1.11

Download a Test DEM

# Download SRTM DEM tile (90m resolution) for a test region
# Using elevation package or manual download from USGS
pip install elevation
eio clip -o dem_test.tif --bounds -105.5 39.5 -104.5 40.5  # Colorado Front Range

Verify Installation

from pysheds.grid import Grid
import numpy as np
import matplotlib.pyplot as plt

print("pysheds installation OK")
# Create a synthetic DEM for testing
dem_arr = np.random.rand(50, 50) * 100
print(f"Test DEM shape: {dem_arr.shape}")

---

Core Workflow

Step 1: DEM Loading and Conditioning

from pysheds.grid import Grid
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors

def create_synthetic_dem(nrows=200, ncols=200, seed=42):
    """
    Create a synthetic DEM with a valley and drainage divide.
    Replace with: grid = Grid.from_raster('your_dem.tif')
    """
    rng = np.random.default_rng(seed)
    x = np.linspace(0, 2*np.pi, ncols)
    y = np.linspace(0, 2*np.pi, nrows)
    X, Y = np.meshgrid(x, y)

    # Base topography: slope + valley
    dem = 1000 - 200 * Y/Y.max() + 100 * np.sin(X) * np.cos(Y/2)
    # Add noise
    dem += 5 * rng.standard_normal((nrows, ncols))
    return dem

def load_and_condition_dem(dem_path=None):
    """
    Load a DEM and condition it for hydrological analysis.

    Parameters
    ----------
    dem_path : str or None
        Path to GeoTIFF DEM. If None, uses synthetic DEM.

    Returns
    -------
    grid : pysheds Grid
    dem : ndarray
    flooded_dem : ndarray (pit-filled)
    """
    if dem_path:
        grid = Grid.from_raster(dem_path)
        dem = grid.read_raster(dem_path)
    else:
        # Synthetic DEM wrapped in pysheds Grid
        dem_arr = create_synthetic_dem()
        # For real workflows, use Grid.from_raster(path)
        # Here we demonstrate with numpy arrays only
        print("Using synthetic DEM (replace with real DEM path)")
        return None, dem_arr, None

    # Pit filling (remove depressions that block flow)
    pit_filled = grid.fill_pits(dem)
    depressions_filled = grid.fill_depressions(pit_filled)
    # Resolve flats (areas with same elevation)
    inflated = grid.resolve_flats(depressions_filled)

    print(f"DEM shape: {dem.shape}")
    print(f"Elevation range: {dem.min():.1f} – {dem.max():.1f} m")
    print(f"Pit-filled cells: {(inflated != dem).sum():,}")

    # Visualize elevation
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    im0 = axes[0].imshow(dem, cmap='terrain')
    plt.colorbar(im0, ax=axes[0], label='Elevation (m)')
    axes[0].set_title("Original DEM")

    im1 = axes[1].imshow(inflated - dem, cmap='Blues')
    plt.colorbar(im1, ax=axes[1], label='Δ Elevation (m)')
    axes[1].set_title("Filled Depressions (difference)")

    plt.tight_layout()
    plt.savefig("dem_conditioning.png", dpi=150)
    plt.show()

    return grid, dem, inflated

# For demonstration with synthetic data
_, dem_raw, _ = load_and_condition_dem()
print(f"Synthetic DEM stats: mean={dem_raw.mean():.1f}, std={dem_raw.std():.1f}")

Step 2: Flow Direction and Accumulation

from pysheds.grid import Grid
import numpy as np
import matplotlib.pyplot as plt

# With a real DEM file:
# grid = Grid.from_raster('dem.tif')
# dem = grid.read_raster('dem.tif')
# pit_filled = grid.fill_pits(dem)
# depressions = grid.fill_depressions(pit_filled)
# inflated = grid.resolve_flats(depressions)

# Compute flow direction (D8 algorithm)
# fdir = grid.flowdir(inflated)
# print(f"Flow direction computed: unique values={np.unique(fdir)}")

# Compute flow accumulation
# acc = grid.accumulation(fdir)
# print(f"Max accumulation: {acc.max():,} cells")

# Stream network (threshold accumulation)
# stream_threshold = int(acc.max() * 0.001)  # 0.1% of max
# streams = acc > stream_threshold
# print(f"Stream cells: {streams.sum():,}")

# For demonstration without DEM file, show the logic:
print("Flow direction and accumulation workflow:")
print("1. grid.flowdir(conditioned_dem) → D8 direction grid")
print("2. grid.accumulation(fdir) → upstream area in cells")
print("3. acc > threshold → binary stream network")
print("4. grid.snap_to_mask(streams, xy) → snap pour point to stream")

# Synthetic example showing accumulation concept
np.random.seed(42)
n = 50
acc_synthetic = np.zeros((n, n))
for i in range(n):
    for j in range(n):
        acc_synthetic[i, j] = (i + 1) * (j + 1)  # simplified

fig, ax = plt.subplots(figsize=(6, 5))
im = ax.imshow(np.log1p(acc_synthetic), cmap='Blues')
plt.colorbar(im, ax=ax, label='log(accumulation + 1)')
ax.set_title("Synthetic Flow Accumulation")
plt.tight_layout()
plt.savefig("flow_accumulation.png", dpi=150)
plt.show()

Step 3: Watershed Delineation

from pysheds.grid import Grid
import geopandas as gpd
from shapely.geometry import shape
import numpy as np
import matplotlib.pyplot as plt

# With a real DEM:
# grid = Grid.from_raster('dem.tif')
# ...
# fdir = grid.flowdir(inflated)
# acc = grid.accumulation(fdir)
#
# # Snap pour point to nearest stream
# x_outlet, y_outlet = -105.0, 40.0  # longitude, latitude
# x_snap, y_snap = grid.snap_to_mask(acc > threshold, (x_outlet, y_outlet))
#
# # Delineate watershed
# catch = grid.catchment(x=x_snap, y=y_snap, fdir=fdir,
#                        xytype='coordinate')
# catch_view = grid.view(catch)
#
# # Convert to polygon
# shapes = grid.polygonize(catch_view)
# for geom, val in shapes:
#     if val:
#         watershed_polygon = shape(geom)
#         break
#
# # Area calculation
# gdf = gpd.GeoDataFrame({'geometry': [watershed_polygon]}, crs=grid.crs)
# gdf_proj = gdf.to_crs('EPSG:32614')  # UTM zone 14N
# area_km2 = gdf_proj.geometry.area.values[0] / 1e6
# print(f"Watershed area: {area_km2:.1f} km²")

# Demonstration output
print("Watershed delineation workflow complete")
print("Key outputs:")
print("  - watershed polygon (GeoDataFrame)")
print("  - drainage area in km²")
print("  - stream network as LineString features")

# SCS Curve Number runoff calculation
def scs_runoff(P_mm, CN):
    """
    Estimate direct runoff depth using SCS Curve Number method.

    Parameters
    ----------
    P_mm : float or ndarray
        Total precipitation depth (mm)
    CN : float
        SCS Curve Number (0-100)

    Returns
    -------
    Q_mm : float or ndarray
        Direct runoff depth (mm)
    """
    S = (25400 / CN) - 254   # potential maximum retention
    Ia = 0.2 * S              # initial abstraction
    P = np.asarray(P_mm)
    Q = np.where(P > Ia, (P - Ia)**2 / (P - Ia + S), 0.0)
    return Q

# Example: 25mm storm on different land cover
storm_mm = 25.0
land_covers = {
    "Row crops (poor)":  86,
    "Row crops (good)":  75,
    "Meadow (good)":     58,
    "Forest (good)":     55,
    "Impervious":        98,
}

print("\nSCS Runoff Estimation (25mm storm):")
print(f"{'Land Cover':<25} {'CN':>4} {'Runoff (mm)':>12} {'Runoff/P (%)':>13}")
print("-" * 56)
for lc, cn in land_covers.items():
    Q = scs_runoff(storm_mm, cn)
    print(f"{lc:<25} {cn:>4} {Q:>12.1f} {100*Q/storm_mm:>13.1f}")

---

Advanced Usage

Automated Multi-Watershed Analysis

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def estimate_flood_frequency(area_km2, mean_annual_precip_mm, CN):
    """
    Simplified regional flood frequency using SCS method.
    Returns peak discharge estimates (m³/s) for return periods.
    """
    return_periods = [2, 5, 10, 25, 50, 100]
    # Regional rainfall frequency (simplified, use NOAA Atlas 14 for real values)
    rp_factors = {2: 1.0, 5: 1.4, 10: 1.7, 25: 2.1, 50: 2.4, 100: 2.8}

    base_storm = mean_annual_precip_mm / 12  # monthly average
    results = []
    for rp in return_periods:
        P = base_storm * rp_factors[rp]
        Q_mm = scs_runoff(P, CN)
        # Rational method approximation: Q = CIA/360 (m³/s)
        runoff_vol_m3 = Q_mm / 1000 * area_km2 * 1e6
        Tc_h = 0.0078 * (area_km2**0.5 / 0.3)**0.77  # time of concentration (h)
        Q_peak = runoff_vol_m3 / (Tc_h * 3600)
        results.append({"return_period": rp, "Q_peak_m3s": Q_peak})

    return pd.DataFrame(results)

flood_table = estimate_flood_frequency(area_km2=150, mean_annual_precip_mm=600, CN=72)
print(flood_table.round(2))

fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogy(flood_table["return_period"], flood_table["Q_peak_m3s"], 'bs-', linewidth=2)
ax.set_xlabel("Return Period (years)")
ax.set_ylabel("Peak Discharge (m³/s)")
ax.set_title("Flood Frequency Curve")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("flood_frequency.png", dpi=150)
plt.show()

---

Troubleshooting

Error: ValueError: No data in raster

Cause: DEM has NoData values throughout or incorrect file path.

Fix:

import rasterio
with rasterio.open("dem.tif") as src:
    print(f"Nodata value: {src.nodata}")
    print(f"Data range: {src.read(1).min()} – {src.read(1).max()}")
    print(f"CRS: {src.crs}")

Issue: All cells flow in same direction (flat DEM)

Cause: DEM has no topographic relief (e.g., coastal plain or incorrect data).

Fix:

# After fill_depressions, always resolve_flats before flowdir
inflated = grid.resolve_flats(depressions)
fdir = grid.flowdir(inflated)

Version Compatibility

| Package | Tested versions | Known issues | |:--------|:----------------|:-------------| | pysheds | 0.3.5 | API changed in 0.3 (Grid class) | | rasterio | 1.3, 1.4 | None |

---

External Resources

Official Documentation

• pysheds documentation
• USGS SRTM DEM data

Key Papers

• Bartos, M. et al. (2021). pysheds: An open-source Python library for watershed delineation. JOSS.

---

Examples

Example 1: Extract Stream Network from SRTM DEM

# =============================================
# Stream network extraction workflow
# =============================================
# NOTE: Requires actual DEM file (download from USGS Earth Explorer)
# This shows the complete workflow assuming dem.tif exists

try:
    from pysheds.grid import Grid
    import numpy as np
    import matplotlib.pyplot as plt
    import geopandas as gpd

    # Step 1: Load
    grid = Grid.from_raster('dem.tif')
    dem = grid.read_raster('dem.tif')
    print(f"DEM loaded: {dem.shape}, range {dem.min():.0f}–{dem.max():.0f} m")

    # Step 2: Condition
    pit_filled = grid.fill_pits(dem)
    dep_filled = grid.fill_depressions(pit_filled)
    inflated   = grid.resolve_flats(dep_filled)

    # Step 3: Flow analysis
    fdir = grid.flowdir(inflated)
    acc  = grid.accumulation(fdir)

    # Step 4: Extract streams at multiple thresholds
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    for ax, thresh_pct in zip(axes, [0.1, 0.5, 2.0]):
        thresh = int(acc.max() * thresh_pct / 100)
        streams = acc > thresh
        ax.imshow(np.log1p(acc), cmap='Blues', alpha=0.5)
        ax.imshow(np.ma.masked_where(~streams, streams), cmap='Reds', alpha=0.8)
        ax.set_title(f"Threshold: {thresh_pct}% ({thresh:,} cells)")
        ax.axis('off')

    plt.suptitle("Stream Networks at Different Thresholds")
    plt.tight_layout()
    plt.savefig("stream_networks.png", dpi=150)
    plt.show()

except FileNotFoundError:
    print("DEM file not found. Download from USGS Earth Explorer (https://earthexplorer.usgs.gov/)")
    print("Then run: eio clip -o dem.tif --bounds <lon_min> <lat_min> <lon_max> <lat_max>")

Interpreting these results: Lower thresholds yield denser stream networks; higher thresholds show only major rivers. Use Strahler stream order to classify streams.

---

Last updated: 2026-03-17 | Maintainer: @xjtulyc Issues: GitHub Issues