Resampling DEMs Without Losing Hydrologic Connectivity

Grid coarsening is unavoidable in regional and continental hydrology workflows — 1 m LiDAR tiles become computationally prohibitive at basin scale, and model inputs often demand a single consistent resolution. But standard GIS resampling methods destroy the drainage skeleton, creating artificial sinks, severed channels, and false divides that cause routing algorithms to fail silently. This page covers the flow-accumulation-weighted approach that prevents those failures. It sits within Spatial Resolution Tradeoffs under the broader Hydrology Data Preparation & DEM Processing domain.

Prerequisites

This technique builds on the Spatial Resolution Tradeoffs baseline. You additionally need:

• A flow accumulation raster computed on the fine-resolution DEM — generated by DEM Pit Filling Algorithms followed by a D8 routing step.
• rasterio >= 1.3, numpy >= 1.24, scipy >= 1.10, and richdem >= 0.3 (or WhiteboxTools for the re-conditioning step).
• Input DEM in a projected CRS with square pixels — geographic coordinates (degrees) produce incorrect scale factors. See Coordinate Reference System Alignment if you need to reproject first.
• The target resolution should ideally be an integer multiple of the source. Non-integer ratios require a different code path (described in the Gotchas section below).

Why Standard Resampling Breaks Drainage Networks

Conventional GIS resampling treats elevation as a scalar field. Hydrologic routing depends on gradient continuity and monotonic descent — properties that are invisible to generic interpolators:

• Bilinear and cubic interpolation smooth sharp channel incisions, artificially raising thalwegs. The result is false divides that split watersheds and truncate stream networks even when the output looks visually clean.
• Nearest-neighbor preserves pixel values but introduces stair-step artifacts that fragment flow paths, generate micro-sinks, and cause D8 flow direction algorithms to stall at artificial flat zones.
• Simple statistical averaging dilutes channel depth relative to hillslope cells, flattening the hydraulic gradient because channel cells are vastly outnumbered by hillslope cells in any block.

The root cause is that all three methods weight fine-resolution cells equally. A 30-metre coarse cell might aggregate 9 fine cells, of which 8 represent hillslope and only 1 is the actual thalweg. Equal weighting means the thalweg’s low elevation is averaged away. Flow-accumulation weighting inverts that ratio, giving the channel cell influence proportional to how much drainage area it carries.

The diagram below shows the contrast between equal-weight averaging and flow-accumulation-weighted averaging for a single 3 × 3 aggregation block. In the equal-weight case, the thalweg cell (elevation 48 m, accumulation 312 cells) contributes only 1/9 of the output; in the weighted case it contributes roughly 55% of the weight.

Core Technique: Flow-Accumulation-Weighted Aggregation

The technique has four logical components:

1. Log-transform weighting. Raw flow accumulation spans many orders of magnitude. A headwater cell may carry 1 contributing cell; a trunk-stream cell may carry millions. Using raw accumulation as a weight would cause trunk streams to dominate so completely that headwater channels are erased. log1p(accumulation) compresses the range while still assigning channels much higher weight than hillslopes.

2. Weighted block mean. For each coarse output cell, compute the weighted mean elevation of the contributing fine cells:

z_coarse = Σ(w_i · z_i) / Σ(w_i)

where w_i = log1p(FAC_i) and the sum runs over all fine cells that fall within the coarse cell’s footprint.

3. Strided downsampling. After computing the weighted moving average across the full fine grid, extract every n-th pixel (where n = target_res / source_res) to produce the coarse grid. This is computationally efficient because scipy.ndimage.uniform_filter runs the convolution on the already-weighted grid in a single pass.

4. Post-aggregation re-conditioning. Coarsening always introduces minor topological noise — flat zones appear where the weighted average of two adjacent blocks lands at the same elevation. These must be resolved by DEM pit filling before running any flow routing on the output.

Annotated Code Example

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

log = logging.getLogger(__name__)


def hydrologic_resample(
    dem_path: str | Path,
    flow_acc_path: str | Path,
    target_resolution: float,
    output_path: str | Path,
    nodata: float = -9999.0,
) -> None:
    """
    Resample a DEM using flow-accumulation weighting to preserve
    hydrologic connectivity across resolution boundaries.

    Parameters
    ----------
    dem_path : path to the fine-resolution DEM (GeoTIFF, projected CRS).
    flow_acc_path : path to the D8/D-Inf flow accumulation raster computed
        on the same fine DEM after sink filling.
    target_resolution : output pixel size in the DEM's CRS units (metres).
    output_path : destination GeoTIFF path.
    nodata : nodata sentinel written to the output raster.
    """
    dem_path = Path(dem_path)
    flow_acc_path = Path(flow_acc_path)
    output_path = Path(output_path)

    log.info("Reading fine DEM: %s", dem_path)
    with rasterio.open(dem_path) as src:
        dem = src.read(1).astype(np.float32)
        src_nodata = src.nodata if src.nodata is not None else nodata
        src_transform = src.transform
        src_crs = src.crs
        src_res = src.res[0]          # assumes square pixels
        src_height, src_width = src.height, src.width

    # Mask nodata as NaN for arithmetic; restored before writing
    dem[dem == src_nodata] = np.nan

    log.info("Reading flow accumulation: %s", flow_acc_path)
    with rasterio.open(flow_acc_path) as fac_src:
        fac = fac_src.read(1).astype(np.float32)
        fac_nodata = fac_src.nodata if fac_src.nodata is not None else nodata
        fac[fac == fac_nodata] = np.nan

    # --- Step 1: compute aggregation window size --------------------------
    scale_factor = target_resolution / src_res
    window_size = max(int(round(scale_factor)), 1)
    log.info(
        "Scale factor %.2f → window_size=%d (%.1f m → %.1f m)",
        scale_factor, window_size, src_res, target_resolution,
    )

    # --- Step 2: build log-transformed weights ----------------------------
    # log1p compresses the accumulation range so headwater channels retain
    # influence; hillslope cells (FAC≈0) receive weight≈0 but never exactly
    # 0 to avoid NaN-propagation in the weighted sum.
    weights = np.log1p(np.nan_to_num(fac, nan=0.0))

    # --- Step 3: weighted block aggregation via uniform_filter ------------
    # uniform_filter is equivalent to a box-blur (moving-window mean) and
    # runs in O(n) time regardless of window size.  Applying it to the
    # element-wise product dem*weights, then dividing by the smoothed
    # weights, gives the correct weighted mean per block.
    dem_filled = np.nan_to_num(dem, nan=0.0)
    dem_weighted = dem_filled * weights

    agg_dem = uniform_filter(dem_weighted, size=window_size, mode="nearest")
    agg_wts = uniform_filter(weights,      size=window_size, mode="nearest")

    # Guard against division by zero in completely un-accumulated areas
    agg_wts = np.where(agg_wts == 0.0, 1e-9, agg_wts)
    smoothed = agg_dem / agg_wts

    # Re-apply nodata mask before downsampling
    smoothed[np.isnan(dem)] = nodata

    # --- Step 4: strided downsample to target resolution -----------------
    # Slice every window_size-th pixel from the smoothed grid. Works only
    # when target_resolution is an integer multiple of src_res; see Gotchas
    # for the non-integer case.
    coarse = smoothed[::window_size, ::window_size]

    # Trim to the expected output dimensions in case rounding added a row
    exp_h = int(np.ceil(src_height / window_size))
    exp_w = int(np.ceil(src_width  / window_size))
    coarse = coarse[:exp_h, :exp_w]

    # --- Step 5: build updated geotransform ------------------------------
    new_transform = Affine(
        target_resolution,  0.0, src_transform.c,
        0.0, -target_resolution, src_transform.f,
    )

    # --- Step 6: write output --------------------------------------------
    output_path.parent.mkdir(parents=True, exist_ok=True)
    with rasterio.open(
        output_path, "w",
        driver="GTiff",
        height=coarse.shape[0],
        width=coarse.shape[1],
        count=1,
        dtype=coarse.dtype,
        crs=src_crs,
        transform=new_transform,
        nodata=nodata,
        compress="lzw",
    ) as dst:
        dst.write(coarse, 1)

    log.info(
        "Wrote coarse DEM %dx%d → %s (re-condition before routing)",
        coarse.shape[1], coarse.shape[0], output_path,
    )
    log.warning(
        "Post-processing required: fill depressions in %s before "
        "computing flow direction at %g m resolution.",
        output_path, target_resolution,
    )


# --- Usage example -------------------------------------------------------
# import logging
# logging.basicConfig(level=logging.INFO)
# hydrologic_resample(
#     "dem_10m.tif", "flow_acc_10m.tif", 30.0, "dem_30m_hydro.tif"
# )

Parameter Reference

Worked Example: 10 m LiDAR to 30 m Analysis Grid

A 10 m LiDAR-derived DEM covering a 200 km² catchment is being coarsened to a 30 m analysis grid for a HEC-HMS simulation. window_size is 3 (30 / 10), so each output cell aggregates a 3 × 3 block of fine cells.

Before resampling (standard bilinear): D8 flow direction routing on the 30 m bilinear output generates 47 spurious sinks and 3 broken flow paths that terminate mid-slope. Stream extraction at a 500-cell accumulation threshold yields a network 12% shorter than the reference NHD reach dataset.

After flow-weighted resampling and re-conditioning: Zero residual sinks after one richdem.FillDepressions() pass. Stream extraction at the same threshold produces a network within 2% of NHD length. The Hausdorff distance between the coarse and fine stream centrelines drops from 85 m to 18 m.

Reading the output: In a GIS viewer, overlay the coarse DEM’s hillshade against the fine-resolution stream network. Thalwegs should remain topographically low relative to their immediate neighbours — no ridge-like artefacts should cross the expected channel path. If you see any, the weight transform was too weak; switch from log1p to raw accumulation for a single trial run to verify the channel geometry exists in the source data.

Gotchas & Edge Cases

• Non-integer resolution ratios. The strided downsample only works when target_res is an integer multiple of src_res. For ratios like 1.5× or 2.5×, replace the strided slice with rasterio.warp.reproject(..., resampling=Resampling.average) applied to the weight-smoothed grid, then update the geotransform accordingly.

• Always re-condition before routing. The script performs weighted aggregation but does not include a full sink-filling routine. Pass the output through richdem.FillDepressions() or WhiteboxTools FillDepressions before computing D-Infinity routing or D8 flow direction on the coarse grid.

• Continental-scale catchments. When flow accumulation spans six or more orders of magnitude (e.g., sub-continental rivers), log1p may still allow trunk streams to dominate so strongly that second-order tributaries are erased. Try np.sqrt() or np.cbrt() as the weight transform and compare stream-network length statistics at both settings.

• Flat-area artifacts at resolution boundaries. Where the source DEM has large pre-existing flat areas (e.g., reservoir surfaces, coastal plains), the weighted average can produce extended flat zones in the coarse output. These require an enforcement pass — WhiteboxTools’ BreachDepressions or RichDEM’s ResolveFlats — before flow routing can proceed. See removing flat-area artifacts from flow direction grids.

• CRS mismatch between DEM and flow accumulation raster. If the flow accumulation raster was generated in a different CRS than the target DEM, the weight grid will not align correctly and the function will silently produce incorrect results. Verify that both inputs share the same CRS and transform before calling the function; see fixing CRS mismatches in watershed shapefiles for correction strategies.

Frequently Asked Questions

Why does bilinear resampling break stream networks?

Bilinear interpolation smooths elevation as a continuous scalar field, which raises thalwegs relative to surrounding terrain. The interpolated surface treats a channel incision the same as any other elevation gradient and averages it toward its neighbours. The result is false divides at the coarser resolution that split or truncate stream networks even when the output DEM looks visually clean in a hillshade.

Do I always need to re-fill sinks after resampling?

Yes. Even flow-weighted aggregation introduces minor topological noise at coarser resolutions — flat zones and micro-sinks appear wherever adjacent blocks happen to converge on the same weighted mean elevation. Always pass the output through a DEM pit filling step before computing new flow directions, regardless of how careful the aggregation was.

What log transform should I use for extreme accumulation ranges?

log1p works well for most catchments because it compresses the enormous FAC range while still giving channels far more influence than hillslopes. For continental rivers where accumulation spans many orders of magnitude, try np.sqrt() to reduce the dominance of very large trunk-stream cells and preserve headwater tributary detail. Compare stream-network length statistics from both transforms against a reference dataset such as the NHD to decide which is appropriate.

Can I apply this technique when the source DEM has already been resampled?

Yes, but the quality depends heavily on the quality of the flow accumulation raster. If the source DEM was previously downsampled with a generic method (bilinear, nearest), its flow accumulation will already reflect the distorted drainage network. Re-deriving the accumulation from a corrected source is always preferable. If you must work from an already-coarsened DEM, run best practices for filling sinks in high-resolution LiDAR data on it first, re-derive accumulation, then apply the weighted aggregation to reach the final target resolution.

Related Topics

• Spatial Resolution Tradeoffs — parent: choosing between 1 m, 10 m, and 30 m DEMs for different modeling goals
• DEM Pit Filling Algorithms — the mandatory re-conditioning step after any resampling operation
• Best Practices for Filling Sinks in High-Resolution LiDAR Data — filling strategies for the fine-resolution DEM before computing weights
• Removing Flat-Area Artifacts from Flow Direction Grids — handling flat zones introduced during aggregation
• Coordinate Reference System Alignment — ensuring DEM and FAC inputs share the same projected CRS before resampling