Watershed Modeling Python for hydrology automation

Resampling DEMs without losing hydrologic connectivity

To preserve drainage networks during grid coarsening, resampling DEMs without losing hydrologic connectivity requires replacing blind statistical interpolation with flow-weighted aggregation. Standard methods (nearest-neighbor, bilinear, cubic) optimize for visual smoothness or pixel-level accuracy, but they ignore directional derivatives. The correct workflow pre-conditions the fine-resolution DEM, applies flow-accumulation-weighted averaging during downsampling, and re-conditions the coarse output. This preserves channel thalwegs, prevents artificial sink creation, and guarantees continuous flow paths across resolution boundaries.

Why Standard Interpolation Breaks Drainage Networks

Conventional GIS resampling treats elevation as a scalar field. Hydrologic routing, however, depends on gradient continuity and monotonic descent. Applying standard algorithms to raw DEMs introduces systematic topological errors:

  • Bilinear/Cubic Interpolation: Smooths sharp incisions and artificially raises thalwegs. This creates false divides that split watersheds and truncate stream networks.
  • Nearest-Neighbor: Introduces stair-stepping artifacts that fragment flow paths, generate micro-sinks, and cause flow-routing algorithms to stall or diverge.
  • Statistical Averaging: Dilutes channel depth relative to hillslopes, flattening the hydraulic gradient and violating the assumption that water follows continuous downhill paths.

As documented in GDAL resampling specifications, these methods optimize for radiometric or geometric fidelity, not topological consistency. For hydrologic modeling, drainage preservation must be explicitly enforced.

The Hydrologically Correct Resampling Pipeline

A robust pipeline decouples statistical aggregation from hydraulic conditioning. Follow these four steps to maintain network integrity:

  1. Pre-condition the Source DEM: Fill sinks and compute D8 (or D∞) flow direction and flow accumulation on the fine-resolution grid. This establishes the true drainage skeleton before aggregation.
  2. Generate Spatial Weights: Normalize or log-transform the flow accumulation raster. High-accumulation cells (channels) receive higher weights, ensuring they dominate the aggregated elevation value.
  3. Apply Flow-Weighted Aggregation: Downsample using a weighted moving window or block reduction. The output elevation at each coarse cell is the weighted mean of contributing fine cells, preserving channel geometry.
  4. Re-condition the Coarse Output: Re-fill depressions and recompute flow direction on the new grid. Coarsening inevitably introduces minor topological noise; re-conditioning guarantees monotonic descent and continuous routing.

This methodology aligns with established Spatial Resolution Tradeoffs where grid coarsening demands explicit drainage preservation rather than blind averaging. For complete preprocessing workflows, see our guide to Hydrology Data Preparation & DEM Processing.

Working Python Implementation

The following script implements flow-weighted resampling using rasterio, numpy, and scipy.ndimage. It assumes you have already generated a flow accumulation raster (flow_acc.tif) using a hydrology library like WhiteboxTools or TauDEM.

python 
import rasterio
import numpy as np
from scipy.ndimage import uniform_filter
from rasterio.transform import Affine
from pathlib import Path

def hydrologic_resample(
    dem_path: str,
    flow_acc_path: str,
    target_resolution: float,
    output_path: str,
    nodata: float = -9999.0
) -> None:
    """
    Resample a DEM using flow-accumulation weighting to preserve hydrologic connectivity.
    Uses a uniform filter for weighted block aggregation.
    """
    # 1. Read source DEM and flow accumulation
    with rasterio.open(dem_path) as src:
        dem = src.read(1).astype(np.float32)
        dem[dem == nodata] = np.nan
        src_transform = src.transform
        src_crs = src.crs
        src_res = src.res[0]  # Assume square pixels
        src_height = src.height
        src_width = src.width
        
    with rasterio.open(flow_acc_path) as fac_src:
        fac = fac_src.read(1).astype(np.float32)
        fac[fac == nodata] = np.nan

    # 2. Calculate aggregation window size
    scale_factor = target_resolution / src_res
    window_size = max(int(np.round(scale_factor)), 1)

    # 3. Create weights from flow accumulation
    # Log-transform prevents extreme channel dominance while preserving hierarchy
    weights = np.log1p(np.nan_to_num(fac, nan=0.0))
    
    # 4. Weighted aggregation via uniform filter, then strided downsample to the
    # coarse grid. The uniform_filter smooths within each window; strided
    # indexing collapses the smoothed grid to the target resolution.
    dem_weighted = dem * weights
    agg_dem = uniform_filter(dem_weighted, size=window_size, mode='nearest')
    agg_weights = uniform_filter(weights, size=window_size, mode='nearest')
    
    # Avoid division by zero in hillslope areas
    agg_weights[agg_weights == 0] = 1e-9
    smoothed_dem = agg_dem / agg_weights
    
    # Restore nodata
    smoothed_dem[np.isnan(smoothed_dem)] = nodata

    # 5. Strided downsample to target resolution
    coarse_dem = smoothed_dem[::window_size, ::window_size]

    # 6. Update geotransform for new resolution
    new_transform = Affine(
        target_resolution, 0.0, src_transform.c,
        0.0, -target_resolution, src_transform.f
    )
    
    # Sanity check on output dimensions
    expected_height = int(np.ceil(src_height * src_res / target_resolution))
    expected_width = int(np.ceil(src_width * src_res / target_resolution))
    coarse_dem = coarse_dem[:expected_height, :expected_width]

    # 6. Write output
    with rasterio.open(
        output_path, 'w',
        driver='GTiff',
        height=coarse_dem.shape[0],
        width=coarse_dem.shape[1],
        count=1,
        dtype=coarse_dem.dtype,
        crs=src_crs,
        transform=new_transform,
        nodata=nodata,
        compress='lzw'
    ) as dst:
        dst.write(coarse_dem, 1)

# Usage example:
# hydrologic_resample("dem_10m.tif", "flow_acc_10m.tif", 30.0, "dem_30m_hydro.tif")

Important Production Notes

  • Re-conditioning Required: The script performs weighted aggregation but does not include a full sink-filling routine. For production, pass the output through richdem.fill_depressions() or WhiteboxTools’ FillDepressions before computing new flow directions.
  • Non-Integer Scaling: If your target resolution isn’t an exact multiple of the source, replace uniform_filter with scipy.ndimage.zoom or use rasterio’s reproject with a custom weighting kernel.
  • Channel Prioritization: The log1p transform balances channel preservation with hillslope representation. Adjust the exponent or use np.sqrt() if your catchment has extreme accumulation ranges.

Validation & Best Practices

After resampling, verify connectivity before routing water:

  1. Flow Path Continuity: Run a D8 flow direction algorithm and trace paths from headwaters to outlets. Paths should not terminate at artificial sinks.
  2. Stream Network Comparison: Extract streams (e.g., flow_acc > threshold) and compute the Jaccard index or Hausdorff distance against the fine-resolution network.
  3. Gradient Checks: Ensure no flat or reverse-gradient cells exist along primary thalwegs. Use WhiteboxTools hydrological analysis tools to validate monotonic descent.

When scaling DEMs for regional or continental modeling, always prioritize topological fidelity over pixel-level RMSE. Hydrologic routing algorithms fail predictably when drainage networks are fragmented, making flow-weighted aggregation a non-negotiable step in any multi-resolution workflow.