Watershed Modeling Python for hydrology automation

Best practices for filling sinks in high-resolution LiDAR data

The most reliable approach for filling sinks in sub-meter LiDAR DEMs combines a priority-flood algorithm with strict memory management, hydro-enforcement preprocessing, and post-fill flow-accumulation validation. This page is part of the DEM Pit Filling Algorithms guide within the broader Hydrology Data Preparation & DEM Processing domain. Avoid naive depression carving on high-resolution LiDAR; it creates artificial flat zones that distort overland flow routing. Instead, use a priority-flood implementation that preserves natural drainage gradients, process data in tiled or out-of-core chunks when rasters exceed available RAM, and always validate filled surfaces against known hydrography before deriving flow direction or accumulation grids.

Prerequisites

Beyond the baseline Python stack (rasterio, richdem, numpy) covered in the DEM Pit Filling Algorithms guide, this technique requires:

  • LiDAR DEM as a single-band float32 GeoTIFF in a projected CRS (UTM zone or state plane) with consistent nodata value — mixing float64 and int16 bands in the same pipeline causes silent precision loss.
  • A stream network vector layer (NHD or equivalent) co-registered to the DEM’s CRS for hydro-enforcement. See coordinate reference system alignment if your layers do not share a common projection.
  • Disk headroom of at least 3× the uncompressed DEM size for intermediate tiles and the final mosaic.
  • richdem ≥ 0.3.4 installed via pip install richdem (or the conda-forge build for systems without a C++ compiler).
  • System RAM ≥ 2× the size of a single processing tile (after overlap expansion).

Why standard depression filling fails at sub-meter resolution

High-resolution LiDAR (0.5–1 m cell size) captures micro-topography that introduces thousands of spurious depressions from vegetation returns, sensor noise, road crowns, and minor surface roughness. Traditional iterative carving artificially lowers terrain to force drainage, which flattens slopes, breaks hydraulic continuity, and produces unrealistic flow paths across what should be sloped hillsides.

The diagram below shows the algorithmic difference between carving (left branch) and priority-flood filling (right branch), and why the latter preserves terrain gradients:

Priority-flood vs carving: algorithm comparison for LiDAR sink removal Decision diagram comparing depression carving (left) against priority-flood filling (right), showing how carving lowers neighbour cells while priority-flood raises the pit cell to the minimum spill elevation. Sink Removal Strategy Comparison Pit cell detected (elevation below all neighbours) Carving Priority-flood Lower neighbour cells until outlet found Raise pit cell to minimum spill elevation Artificial flat zones & distorted slope gradients Terrain morphology preserved; gradients intact Flow direction & accumulation

Priority-flood algorithms process cells in order of elevation using a min-heap, propagating drainage outward from basin edges. This guarantees that only the minimum number of cells are modified to eliminate sinks, and that modified cells are raised — never lowered — preserving the original terrain morphology.

Core technique: priority-flood with hydro-enforcement

Step 1 — Hydro-enforce before filling

Never run a fill algorithm on raw LiDAR. Unenforced data treats culverts, bridges, road embankments, and engineered channels as solid barriers, producing hydrologically disconnected watersheds. Hydro-enforcement must happen first:

  • Vector-burn stream networks using depth-based incision: subtract 1–3 m from DEM cells that intersect the stream centerline, proportional to local channel width.
  • Drop culvert and bridge elevations to match surveyed invert levels. Even a 0.5 m mismatch will create artificial barriers at sub-meter resolution.
  • Mask known water bodies (lakes, reservoirs, tidal flats) to prevent the algorithm from filling legitimate open-water surfaces. Store masked cells as nodata before filling, then restore them afterward.

After burning, the fill step sees a DEM that already represents actual drainage pathways, so the pit-filling pass only corrects genuine micro-topographic artifacts.

Step 2 — Tiled windowed I/O

High-resolution LiDAR tiles routinely exceed 10–50 GB uncompressed. Loading entire rasters into memory causes swap thrashing or OOM kills. The diagram below shows the overlap-buffer strategy for memory-safe tiled processing:

Tiled overlap buffer strategy for large LiDAR DEM processing Diagram showing four DEM processing tiles arranged in a 2x2 grid, with the overlap buffer zone highlighted on each tile edge to prevent inter-tile sink formation. Full DEM extent Tile 1 chunk_px × chunk_px core Tile 2 chunk_px × chunk_px core Tile 3 Tile 4 overlap_m buffer (≥ 500 m) Overlap buffer zone Tile core (written to output)

Process in overlap-aware windows:

  • Read and write using rasterio windowed operations with a ≥ 500 m overlap buffer per tile edge.
  • The overlap buffer ensures that flow paths crossing tile boundaries have sufficient context for accurate D8 flow direction computation.
  • After all tiles are processed, run a single coarse-resolution fill pass on the stitched mosaic to resolve any remaining inter-tile sinks.

Step 3 — Post-fill validation

Compute D8 flow direction and accumulation on the filled DEM. Reject outputs that exhibit:

  • Unexplained flat zones larger than 3×3 cells outside known water bodies
  • Accumulation spikes or drop-offs at tile boundaries
  • Disconnected flow paths crossing mapped ridgelines
  • Drainage density deviation > 15% from pre-fill hydrography

Annotated code example

The following function implements all three steps with logging at every decision point:

python 
import logging
import os

import numpy as np
import rasterio
import richdem as rd
from rasterio.windows import Window

logger = logging.getLogger(__name__)


def fill_lidar_dem(
    input_tif: str,
    output_tif: str,
    chunk_px: int = 2048,
    overlap_m: float = 500.0,
) -> None:
    """Fill sinks in a high-resolution LiDAR DEM using priority-flood (richdem).

    Processes the raster in overlap-aware tiles to keep peak RAM usage bounded.
    A caller should run hydro-enforcement (stream burning, culvert drops) on
    ``input_tif`` before passing it here.

    Args:
        input_tif:  Path to a single-band float32 GeoTIFF (projected CRS).
        output_tif: Destination path for the filled DEM.
        chunk_px:   Tile dimension in pixels for windowed I/O (default 2048).
        overlap_m:  Buffer in map units added to each tile edge (default 500 m).
    """
    if not os.path.exists(input_tif):
        raise FileNotFoundError(f"Input DEM not found: {input_tif}")

    logger.info("Opening source DEM: %s", input_tif)

    with rasterio.open(input_tif) as src:
        if src.count != 1:
            raise ValueError(
                f"Expected single-band raster; got {src.count} bands."
            )

        # Compute overlap in pixels from the map-unit buffer.
        # src.res[0] is the pixel width in map units (e.g. metres for UTM).
        overlap_px = max(1, int(overlap_m / src.res[0]))
        nodata_val = src.nodata if src.nodata is not None else -9999.0
        logger.info(
            "DEM size: %d × %d px | resolution: %.3f m | overlap buffer: %d px",
            src.width,
            src.height,
            src.res[0],
            overlap_px,
        )

        # Copy source profile; force float32 + LZW + cloud-optimised tiling.
        profile = src.profile.copy()
        profile.update(
            dtype="float32",
            compress="lzw",
            nodata=nodata_val,
            tiled=True,
            blockxsize=512,
            blockysize=512,
        )

        tile_count = 0
        with rasterio.open(output_tif, "w", **profile) as dst:
            for row_off in range(0, src.height, chunk_px):
                for col_off in range(0, src.width, chunk_px):
                    # Expand the read window by the overlap buffer on all sides,
                    # clamped to the raster extent.
                    win_row = max(0, row_off - overlap_px)
                    win_col = max(0, col_off - overlap_px)
                    win_h = min(
                        chunk_px + 2 * overlap_px, src.height - win_row
                    )
                    win_w = min(
                        chunk_px + 2 * overlap_px, src.width - win_col
                    )

                    window = Window(win_col, win_row, win_w, win_h)
                    data = src.read(1, window=window).astype(np.float32)

                    # Identify nodata cells before filling so we can restore them.
                    if src.nodata is not None:
                        nodata_mask = data == src.nodata
                    else:
                        nodata_mask = np.zeros(data.shape, dtype=bool)

                    # richdem requires a sentinel nodata value on the array.
                    data[nodata_mask] = -9999.0

                    # Priority-flood fill via richdem.
                    # FillDepressions raises pit cells to their spill elevation —
                    # it never lowers neighbours, so slope gradients are preserved.
                    rd_arr = rd.rdarray(data, no_data=-9999.0)
                    filled = rd.FillDepressions(rd_arr, in_place=False)
                    filled_np = np.array(filled, dtype=np.float32)

                    # Restore nodata where the source had it.
                    filled_np[nodata_mask] = nodata_val

                    # Strip the overlap padding before writing so tile seams
                    # do not double-write or introduce boundary errors.
                    core_row = overlap_px if row_off > 0 else 0
                    core_col = overlap_px if col_off > 0 else 0
                    core_h = min(chunk_px, filled_np.shape[0] - core_row)
                    core_w = min(chunk_px, filled_np.shape[1] - core_col)

                    core = filled_np[
                        core_row : core_row + core_h,
                        core_col : core_col + core_w,
                    ]
                    out_window = Window(col_off, row_off, core_w, core_h)
                    dst.write(core, 1, window=out_window)

                    tile_count += 1
                    logger.debug(
                        "Wrote tile %d (row=%d, col=%d, size=%d×%d px)",
                        tile_count,
                        row_off,
                        col_off,
                        core_h,
                        core_w,
                    )

    logger.info(
        "Filled DEM written to %s (%d tiles processed)", output_tif, tile_count
    )

Parameter reference

Parameter Type Typical value Effect on hydrology
chunk_px int 1024–4096 Larger chunks reduce tile-seam artifacts but increase peak RAM; 2048 is safe up to ~16 GB datasets on a 32 GB machine.
overlap_m float 500–1000 m Controls how much terrain context surrounds each tile; too small and inter-tile sinks form; too large and RAM spikes.
rd.FillDepressions(in_place=False) bool False for production True modifies the input array in memory — convenient for exploration but unsafe if the same array is reused.
nodata_val float −9999.0 Must match the value written to the output profile; mismatches produce false-positive sink detections at nodata boundaries.
LZW compression str "lzw" Lossless, fast decode; use "deflate" with predictor=3 for float rasters if further size reduction is needed.
blockxsize / blockysize int 512 Matches typical COG tile size; aligns reads to cloud storage block boundaries (S3, GCS).

Worked example: output interpretation

A correctly filled 0.5 m LiDAR DEM will show the following characteristics when you compute D8 flow direction and accumulation afterward:

  • Zero isolated flat cells outside mapped water bodies. Any flat zone is a sign that the overlap buffer was too small (inter-tile artifact) or that hydro-enforcement was skipped (culvert barrier).
  • Smooth accumulation gradients up tributary networks. Sudden high-accumulation spikes mid-slope indicate a residual sink that was filled incorrectly to a higher elevation than its spill point.
  • Drainage density within 10–15% of NHD reference values for the same watershed. Large deviations (> 20%) suggest the stream-burning depth was too shallow, causing channels to remain blocked.
  • Tile boundary inspection: compute a binary difference raster (filled – unfilled) and check that the fill magnitude is consistent across tile boundaries. A sharp step in fill depth at a tile edge is a definitive sign of an insufficient overlap buffer.

When resampling DEMs before or after filling, always fill at the native resolution first — coarsening before fill merges micro-depressions into real basins and produces an over-filled DEM.

Gotchas and edge cases

  • Float precision and artificial terraces. Casting the DEM to int16 before filling truncates sub-centimeter gradients and creates artificial step-like terracing across slopes. Always process LiDAR as float32 minimum; float64 if your richdem build supports it and the dataset is small enough.
  • Nodata at tile edges. If the source LiDAR has irregular coverage (flightline gaps, coastal clipping), nodata cells inside the overlap buffer can be misinterpreted as basin outlets, directing flow toward them. Apply a pre-fill nodata edge fill using scipy.ndimage.distance_transform_edt to propagate the nearest valid elevation into gap regions before processing.
  • Inter-tile spill resolution. The tiled fill corrects intra-tile depressions but cannot resolve sinks whose spill point lies in an adjacent tile’s non-overlapping core. After mosaicking, run one coarse-resolution priority-flood pass (e.g., at 10 m resampled from the 0.5 m filled tiles) to catch these residual cases, then use the coarse fill only to flag problem tiles for a targeted re-run rather than blending coarse values into the high-res output.
  • Memory-mapped arrays for very large tiles. If a single tile with overlap exceeds available RAM, replace the src.read(...) call with a np.memmap-backed array written to a scratch directory or a temp path, then fill in-place. This trades I/O throughput for peak RAM headroom.
  • Hydro-enforcement depth calibration. A stream-burning depth that is too deep (> 5 m) can over-incise the channel, causing adjacent hillslope cells to drain directly to the burned channel rather than following natural overland flow paths. Calibrate burn depth against field measurements or channel cross-sections from USGS 3DEP data wherever available.
  • Vegetation-return artifacts in dense forest. Even ground-classified LiDAR returns in dense canopy can be 0.2–0.4 m higher than bare earth. If your LiDAR product is not a bare-earth DEM (last-return or ground-filtered), apply a morphological open filter before filling to suppress isolated high-return clusters that would otherwise create large artificial basins.

Frequently asked questions

Why does naive depression carving fail on sub-meter LiDAR?

Carving artificially lowers terrain to force drainage, which flattens slopes and breaks hydraulic continuity. At sub-meter resolution thousands of micro-depressions exist, so carving produces large, unrealistic flat zones that distort overland flow routing. Priority-flood filling raises pit cells to their minimum spill elevation instead, leaving surrounding terrain untouched.

How large an overlap buffer should I use between tiles?

At minimum 500 m of overlap (in map units). This ensures that flow paths crossing tile boundaries have enough context cells for accurate direction and accumulation computation. For very flat terrain — agricultural lowlands, coastal plains — increase to 1000 m to prevent inter-tile sinks from forming at seams.

Do I need a final mosaic-level fill pass after tiling?

Yes, for the highest accuracy. The per-tile fill resolves intra-tile depressions, but inter-tile spill points that span the buffer zone can still create residual sinks at the mosaic level. A single coarse-resolution pass on the stitched output — at 10 m resampled from the native 0.5 m — catches these without altering the high-resolution values.