Handling Data Gaps in Sparse Bathymetry Interpolation #

Every gridding pass through the DEM interpolation stage faces the same silent hazard: the survey never covered the whole area. Widely spaced sounding lines, acoustic shadow zones behind a ridge, and no-data holidays where the swath dropped out all leave wide holes in the point set. An unconstrained interpolator does not know these holes are unsurveyed — it simply extrapolates, filling a kilometre-wide gap with a smooth but entirely fictional dome, or leaving a hard NaN void with a cliff around its rim. Both outcomes are dangerous: the false dome looks like a real shoal to a downstream navigation-surface consumer, and the void breaks any raster-algebra step that assumes continuity. The fix is to stop treating interpolation and extrapolation as the same operation. This page builds an explicit distance-to-sounding mask, caps interpolation to a defensible maximum gap radius, and records a per-cell data-density flag so every value carries its own confidence.

Why the False Dome Appears #

Inverse-distance weighting and most radial interpolators compute an output value as a weighted average of nearby samples. Inside a well-covered area that average is a genuine interpolation, bounded by real observations on all sides. Over a wide gap, however, the nearest samples all sit on one rim of the hole, so the estimator quietly becomes an extrapolation — it projects the trend of the rim inward with no data to constrain it. Because IDW weights are normalised, the interpolated surface always lands somewhere between the observed rim depths, and the smooth blending reads as a plausible bump. The estimator has no notion of “I am now too far from any sounding to answer.”

Formally, for a query node at distance $d_k$ from each of $N$ contributing soundings, IDW returns $\hat{z} = \sum_k w_k z_k / \sum_k w_k$ with $w_k = d_k^{-p}$. Nothing in that expression references how far the closest sounding actually is, so the estimate is defined even when every $d_k$ is enormous. The reproduction below grids a sparse point set with a deliberate rectangular hole and watches the dome form:

import numpy as np
from scipy.interpolate import griddata

rng = np.random.default_rng(42)
# Sparse soundings on a gently deepening slope, then punch a wide hole.
pts = rng.uniform(0, 1000, size=(400, 2))
depth = -8.0 - 0.01 * pts[:, 0]  # ~ -8 m to -18 m west-to-east
keep = ~((pts[:, 0] > 400) & (pts[:, 0] < 700) &
         (pts[:, 1] > 300) & (pts[:, 1] < 700))
pts, depth = pts[keep], depth[keep]  # 300 m x 400 m unsurveyed gap

gx, gy = np.meshgrid(np.linspace(0, 1000, 200), np.linspace(0, 1000, 200))
grid = griddata(pts, depth, (gx, gy), method="linear")
centre = grid[100, 110]           # inside the gap
edge = -8.0 - 0.01 * 550          # what the real slope would predict
print(round(centre - edge, 2))    # a non-zero, unconstrained fabrication

The value returned inside the hole is a pure invention of the triangulation, not a measurement. Choosing a smarter estimator does not help: the kriging vs IDW comparison shows kriging produces an equally confident-looking mean over the same gap, merely with a large kriging variance nobody downstream ever reads. The durable solution is a geometric one that any interpolator can share.

Step-by-Step Fix #

1. Build a distance-to-nearest-sounding surface #

Construct a scipy.spatial.cKDTree over the sounding positions and query it once per grid node. The returned nearest-neighbour distance is the key discriminator between interpolation (short distance, real neighbours on multiple sides) and extrapolation (long distance, deep inside a gap):

import logging
import numpy as np
from scipy.spatial import cKDTree

logger = logging.getLogger(__name__)


def distance_to_data(sounding_xy: np.ndarray,
                     grid_x: np.ndarray,
                     grid_y: np.ndarray) -> np.ndarray:
    """Distance from every grid node to its nearest real sounding.

    Args:
        sounding_xy: (N, 2) projected sounding coordinates in metres.
        grid_x: 2-D array of grid node X coordinates (from meshgrid).
        grid_y: 2-D array of grid node Y coordinates, same shape as grid_x.

    Returns:
        2-D float array of nearest-sounding distances, metres.

    Raises:
        ValueError: if inputs are empty or shapes disagree.
    """
    if sounding_xy.size == 0:
        raise ValueError("no soundings supplied; cannot build a gap mask")
    if grid_x.shape != grid_y.shape:
        raise ValueError("grid_x and grid_y must share a shape")

    tree = cKDTree(sounding_xy)
    query_pts = np.column_stack([grid_x.ravel(), grid_y.ravel()])
    dist, _ = tree.query(query_pts, k=1, workers=-1)
    dist_grid = dist.reshape(grid_x.shape)
    logger.info("distance-to-data: median=%.1f m, p99=%.1f m",
                float(np.median(dist_grid)), float(np.percentile(dist_grid, 99)))
    return dist_grid

2. Cap interpolation to a maximum gap radius #

Choose a maximum trusted search radius from the survey geometry — a common, defensible rule is a small multiple of the nominal along-track or line-spacing interval, so a cell may be interpolated across a normal inter-line gap but not across a whole missing survey block. Every cell whose distance-to-data exceeds that radius is demoted to no-data regardless of what the interpolator returned:

import logging
import numpy as np

logger = logging.getLogger(__name__)


def cap_to_gap_radius(interpolated: np.ndarray,
                      dist_to_data: np.ndarray,
                      max_gap_m: float,
                      nodata: float = np.nan) -> np.ndarray:
    """Mask interpolated cells that lie beyond the trusted search radius.

    Args:
        interpolated: 2-D grid returned by any interpolator.
        dist_to_data: matching distance-to-nearest-sounding grid, metres.
        max_gap_m: largest distance-to-data still treated as interpolation.
        nodata: sentinel written where the gap exceeds max_gap_m.

    Returns:
        Copy of ``interpolated`` with over-gap cells set to ``nodata``.

    Raises:
        ValueError: if max_gap_m is not positive or shapes disagree.
    """
    if max_gap_m <= 0:
        raise ValueError("max_gap_m must be > 0")
    if interpolated.shape != dist_to_data.shape:
        raise ValueError("interpolated and dist_to_data shapes differ")

    over_gap = dist_to_data > max_gap_m
    result = np.where(over_gap, nodata, interpolated).astype(np.float32)
    logger.info("masked %d/%d cells beyond %.0f m gap radius",
                int(over_gap.sum()), over_gap.size, max_gap_m)
    return result

Small holes inside the covered area — single beam dropouts, a few rejected pings — are legitimately fillable and need not be masked. If you want to close only those, dilate the valid region by the same max_gap_m in pixels using scipy.ndimage.binary_dilation and fill only where the dilated mask and the interpolation agree; anything the dilation cannot reach stays a true void that the surface smoothing stage must respect rather than average across.

3. Record a per-cell data-density flag #

Do not discard the information about why a cell is what it is. Emit a companion integer band classifying every cell as measured, interpolated-within-radius, or masked gap. Downstream navigation-surface and habitat-model consumers can then filter by confidence instead of trusting all cells equally:

import numpy as np

MEASURED, INTERPOLATED, GAP = np.uint8(1), np.uint8(2), np.uint8(0)


def data_density_flag(dist_to_data: np.ndarray,
                      cell_size_m: float,
                      max_gap_m: float) -> np.ndarray:
    """Classify each cell as measured, interpolated, or gap.

    A cell within roughly one cell-size of a sounding is treated as measured;
    within max_gap_m it is interpolated; beyond it is an unsurveyed gap.
    """
    flag = np.full(dist_to_data.shape, GAP, dtype=np.uint8)
    flag[dist_to_data <= max_gap_m] = INTERPOLATED
    flag[dist_to_data <= cell_size_m] = MEASURED
    return flag

Persist this flag as a second raster band or a sidecar so it travels with the DEM; align its CRS with the survey following the CRS alignment workflow so the geometry of the mask matches the geometry of the depths exactly.

Verification and Acceptance Test #

The gate is simple and absolute: after masking, no finite interpolated cell may lie farther from data than the maximum gap radius. Assert it and log the masked-cell count so a regression that widens the trusted radius cannot ship silently:

import logging
import numpy as np

logger = logging.getLogger(__name__)


def assert_gap_constraint(capped: np.ndarray,
                          dist_to_data: np.ndarray,
                          max_gap_m: float) -> dict:
    """Fail if any finite cell exceeds the permitted gap distance."""
    finite = np.isfinite(capped)
    violations = int(np.sum(finite & (dist_to_data > max_gap_m)))
    if violations:
        raise AssertionError(
            f"{violations} finite cells exceed {max_gap_m} m gap radius")
    masked = int(np.sum(~finite))
    logger.info("gap gate OK: %d cells masked, 0 violations", masked)
    return {"masked_cells": masked, "violations": violations}

Wire it into CI as a hard exit code so the build fails on any leaked extrapolation:

python -c "from gate import assert_gap_constraint as g; import numpy as np, sys; \
d = np.load('dist.npy'); z = np.load('capped.npy'); \
sys.exit(0 if g(z, d, 250.0)['violations'] == 0 else 1)"

Edge Cases and Gotchas #

• Anisotropic line spacing. Survey lines are usually far closer along-track than across-track, so a single isotropic max_gap_m either over-masks along dense lines or leaks across sparse ones. Scale the cKDTree query into a stretched coordinate frame (divide the across-track axis by the anisotropy ratio) before computing distances, so the trusted radius becomes an ellipse aligned to the survey geometry rather than a circle.
• Coastline and land masks. A gap that coincides with land is not an unsurveyed hole — it is genuinely dry. Rasterise the coastline with geopandas/rasterio and stamp land cells to a distinct sentinel before the gap gate, or the mask will “reclaim” intertidal cells that no sounding will ever cover. Re-check nearshore values after any vertical shift from the tidal datum transformation stage.
• Carry an uncertainty surface alongside. The distance-to-data grid is itself a cheap proxy for interpolation uncertainty: pair it with the kriging variance where available, and propagate both into the total vertical uncertainty so that IHO S-44 order classification reflects the true, gap-aware confidence rather than an even blanket over measured and invented cells alike.

Related #

• Kriging vs IDW for Bathymetry Interpolation — how each estimator behaves over the same gap, and why kriging variance alone does not save you.
• Preserving Seafloor Features While Smoothing — the downstream stage that must respect the voids this mask preserves.
• DEM Interpolation Techniques for Seafloor Mapping — parent workflow: estimator selection, grid resolution, and search-neighbourhood design.

Up: DEM Interpolation Techniques for Seafloor Mapping