Why does a vertical offset between DEM tiles break watershed delineation?

A step of even a few centimetres along a tile boundary acts as a continuous scarp across the landscape. On steep terrain the real gradient dominates and nothing happens. On flat terrain the step is the steepest local gradient, so flow follows the seam instead of the valley, and the extracted network contains a suspiciously straight reach exactly along the tile edge.

Should DEM tiles be blended or hard-cut at the seam?

Blend when the two sources are of comparable quality and overlap by more than a few cells, because a feathered transition removes the step without inventing terrain. Hard-cut and prefer the better source when one is clearly superior — a 1 metre LiDAR tile beside a 30 metre global product should not be averaged with it, only trimmed to it. Estimate and remove the vertical offset first in either case.

How do I detect a vertical datum mismatch between elevation sources?

Compare the two sources over their overlap and look at the distribution of the difference. A near-constant offset of the order of tens of metres with a small spread is a datum mismatch — typically ellipsoidal versus orthometric heights, where the geoid separation is the offset. A difference that varies with land cover is a canopy effect, and one that varies with slope is a horizontal misregistration.

Merging Multi-Source DEM Tiles Without Seam Artifacts

A mosaic assembled from one acquisition is usually seamless. A mosaic assembled from several — different flights, different years, a LiDAR block dropped into a regional grid — carries steps at every boundary, and those steps become drainage divides. This guide covers detecting and removing them, as part of the SRTM and LiDAR data acquisition topic within hydrology data preparation and DEM processing.

Prerequisites

All tiles reprojected into one CRS and snapped to a common grid origin. Resampling during the merge itself hides the offsets this guide is about.
Overlapping data at the seams, ideally tens of cells wide. Butt-joined tiles with no overlap cannot have their offset estimated from the data.
rasterio, numpy, scipy.

Core Technique: Estimate, Remove, Feather

The three steps are always the same and are usually collapsed into one call to a merge function that does none of them.

Annotated Code Example

python

import logging

import numpy as np
import rasterio
from rasterio.merge import merge as rio_merge
from scipy.ndimage import distance_transform_edt

log = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")


def estimate_offset(a: np.ndarray, b: np.ndarray, nodata: float) -> dict:
    """
    Estimate the vertical offset between two overlapping surfaces.

    Uses the median rather than the mean because canopy, buildings and water
    produce a long one-sided tail that would drag a mean off the true offset.
    """
    m = (a != nodata) & (b != nodata) & np.isfinite(a) & np.isfinite(b)
    if m.sum() < 100:
        raise ValueError(f"only {int(m.sum())} overlapping cells — too few to "
                         "estimate an offset; the tiles may butt-join")
    d = (a[m] - b[m]).astype("float64")
    median = float(np.median(d))
    iqr = float(np.percentile(d, 75) - np.percentile(d, 25))

    log.info("Overlap: %d cells, median offset %.3f m, IQR %.3f m",
             int(m.sum()), median, iqr)
    # A tight distribution around a large offset is a datum difference;
    # a wide one is a genuine difference in what the two sensors measured.
    if abs(median) > 5.0 and iqr < 1.0:
        log.warning("A %.1f m offset with a %.2f m spread looks like a vertical "
                    "datum mismatch — check ellipsoidal vs orthometric heights",
                    median, iqr)
    elif iqr > 3.0:
        log.warning("Offset spread of %.2f m is large — the two sources may "
                    "differ in what they measure (canopy vs bare earth), and a "
                    "single constant offset will not reconcile them", iqr)
    return {"median": median, "iqr": iqr, "n_cells": int(m.sum())}


def feathered_merge(
    tile_paths: list[str],
    out_path: str,
    reference_index: int = 0,
    feather_cells: int = 40,
) -> dict:
    """
    Merge tiles with per-tile offset removal and a feathered transition.

    Parameters
    ----------
    tile_paths      : Tiles, already on a common CRS and grid.
    out_path        : Destination mosaic.
    reference_index : Which tile defines the vertical reference.
    feather_cells   : Width of the blend ramp, in cells.
    """
    srcs = [rasterio.open(p) for p in tile_paths]
    try:
        nodata = srcs[0].nodata if srcs[0].nodata is not None else -9999.0
        profile = srcs[0].profile.copy()

        # --- Pass 1: merge without adjustment, to find the overlaps. The
        # result is discarded; only its geometry is used. ---
        base, transform = rio_merge(srcs, nodata=nodata)
        shape = base.shape[1:]
        log.info("Mosaic extent: %d × %d cells from %d tiles",
                 shape[1], shape[0], len(srcs))

        accum = np.zeros(shape, dtype="float64")
        weight = np.zeros(shape, dtype="float64")
        offsets = {}

        ref_grid, _ = rio_merge([srcs[reference_index]], nodata=nodata,
                                bounds=rio_merge(srcs, nodata=nodata)[0].shape and None)
        # Re-read the reference on the mosaic grid so offsets are comparable.
        ref_grid, _ = rio_merge([srcs[reference_index]], nodata=nodata)

        for i, src in enumerate(srcs):
            grid, _ = rio_merge([src], nodata=nodata)
            arr = grid[0].astype("float64")
            valid = (arr != nodata) & np.isfinite(arr)

            if i != reference_index:
                # Offsets are estimated against the running mosaic, so tiles
                # chain correctly when no single tile overlaps all others.
                current = np.divide(accum, weight, out=np.full(shape, nodata),
                                    where=weight > 0)
                try:
                    off = estimate_offset(arr, current, nodata)
                    arr[valid] -= off["median"]
                    offsets[tile_paths[i]] = off["median"]
                    log.info("Tile %s adjusted by %+.3f m",
                             tile_paths[i], -off["median"])
                except ValueError as exc:
                    log.warning("Tile %s: %s — merged unadjusted",
                                tile_paths[i], exc)
                    offsets[tile_paths[i]] = 0.0
            else:
                offsets[tile_paths[i]] = 0.0

            # --- Feather weight: 0 at the tile edge rising to 1 inside, so
            # each tile contributes fully in its interior and tapers at its
            # margin. This is what removes the slope kink. ---
            dist = distance_transform_edt(valid)
            w = np.clip(dist / max(1, feather_cells), 0.0, 1.0)
            accum += np.where(valid, arr * w, 0.0)
            weight += np.where(valid, w, 0.0)

        mosaic = np.divide(accum, weight, out=np.full(shape, nodata, "float64"),
                           where=weight > 0).astype("float32")

        profile.update(height=shape[0], width=shape[1], transform=transform,
                       dtype="float32", nodata=nodata, compress="LZW",
                       tiled=True, blockxsize=512, blockysize=512)
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.write(mosaic, 1)
        log.info("Wrote mosaic: %s", out_path)
        return {"output": out_path, "offsets_m": offsets}
    finally:
        for s in srcs:
            s.close()


def seam_check(mosaic_path: str, seam_rows: list[int], window: int = 3) -> dict:
    """
    Measure the residual step across each known seam.

    A step under about a tenth of the DEM's vertical accuracy is invisible to
    the router; anything larger will steer flow on flat ground.
    """
    with rasterio.open(mosaic_path) as src:
        arr = src.read(1).astype("float64")
        nodata = src.nodata

    results = []
    for row in seam_rows:
        above = arr[row - window:row, :]
        below = arr[row:row + window, :]
        m = (above != nodata) & (below != nodata)
        if not m.any():
            continue
        step = float(np.nanmedian(below[m] - above[m]))
        results.append({"row": row, "residual_step_m": step})
        (log.info if abs(step) < 0.05 else log.error)(
            "Seam at row %d: residual step %.4f m", row, step)
    return {"seams": results}


# --- Example usage ---
# feathered_merge(["lidar_a.tif", "lidar_b.tif", "srtm_fill.tif"],
#                 "mosaic.tif", feather_cells=40)

Parameter Reference

Parameter	Typical	Effect
`feather_cells`	20–60	Too narrow leaves a visible kink; too wide blends real terrain across the seam
Offset statistic	Median	Mean is dragged by canopy and building returns
Reference tile	The best source	Everything else is adjusted toward it
Overlap required	≥ 100 cells	Below that, the offset estimate is noise
Acceptable residual step	< 0.05 m	Below the point where it can compete with real gradients

Worked Example: Reading the Offsets

A mosaic of four sources over one basin:

Tile	Source	Median offset	IQR	Treatment
A	1 m LiDAR, 2021	0.000 m	—	Reference
B	1 m LiDAR, 2016	+0.184 m	0.09 m	Constant offset removed
C	3 m IfSAR	−1.42 m	2.8 m	Offset removed, wide spread flagged
D	30 m global	+29.7 m	0.4 m	Datum mismatch — converted, not offset

Tile D is the instructive one. A 29.7 m offset with a 0.4 m spread is not a calibration difference; it is the geoid separation between ellipsoidal and orthometric heights at that latitude. Subtracting a constant would appear to work and would leave the tile in the wrong vertical datum, so any later comparison against surveyed elevations would be out by that amount. The right fix is a datum transformation, not an offset.

Whether a residual step matters at all depends on the terrain it sits in. The same 15 cm offset is invisible on a hillslope and decisive on a floodplain.

Gotchas and Edge Cases

rasterio.merge straight out of the box. It takes the first valid value at each cell and does nothing about offsets, so every seam becomes a step.
A mean offset instead of a median. Canopy and buildings drag it, typically by tens of centimetres on forested overlaps.
Feathering across a genuine cliff. If a real escarpment runs along the seam, feathering smooths it away. Check the seam against a hillshade before blending.
Butt-joined tiles. With no overlap there is nothing to estimate from. Either obtain overlapping data or accept an unadjusted mosaic and document it.
Merging before reprojecting. Resampling during the merge blurs the seam and makes the offset harder to see, without removing it.
Ignoring acquisition date. Two LiDAR flights five years apart differ genuinely where the land changed. The offset estimate treats that as noise, which is fine for a mosaic and wrong for change detection.
Not checking the seams afterwards. The residual step is one line of code and is the only proof the process worked.

SRTM & LiDAR Data Acquisition — the parent topic: sources, staging and the acquisition pipeline this feeds
SRTM vs LiDAR DEMs for Regional Watershed Modeling — why two sources differ systematically before any seam appears
Choosing a Projected CRS for Basin-Scale Hydrology — the common grid all tiles must reach before merging
Coordinate Reference System Alignment — grid snapping and the half-cell offsets that produce their own seams

Related topics