Coastal & Marine Spatial Analysis Pipelines

Handling VDatum Coverage Gaps Offshore

VDatum covers the US coastal zone. It does not cover the open ocean, it thins out across the shelf break, and its regional grids meet at seams where two models disagree by centimetres. A survey that runs from the beach to the shelf edge therefore has a separation value for part of its extent and no value for the rest, and the question of what to write in the rest is where most of the damage in a tidal-datum workflow happens. This page is that decision, within the tidal datum transformation workflow.

The default behaviours are all wrong in the same direction: they produce a complete-looking raster. A nodata fill of zero yields a chart that says the seabed is exactly at the tidal datum offshore. A nearest-neighbour fill yields a chart that says the offshore separation equals the last coastal value — and the coastal value is the least representative one available, because the separation gradient is steepest right where coverage ends.

Three Ways to Fill the Gap, Two of Them Fabrications A profile running left to right from the shoreline to beyond the shelf break. The true NAVD88-to-MLLW separation curve rises from about 0.3 metres nearshore to about 1.1 metres at the shelf break and then flattens. A vertical line marks where VDatum coverage ends, roughly two thirds of the way out. To the left of it the curve is drawn solid; to the right, three candidate continuations are drawn: a flat line at zero, labelled zero fill and marked as claiming no correction at all; a flat line held at the last covered value, labelled nearest fill and marked as freezing the steepest part of the gradient; and an absent line with hatching, labelled nodata and marked as the only one that makes no claim. A caption notes that the first two produce a complete-looking raster and the third does not, which is why the first two are chosen by default. distance offshore separation (m) VDatum coverage ends zero fill — no correction at all nearest fill — freezes the edge nodata the only option making no claim measured separation The first two produce a complete-looking raster and the third does not. That is why the first two get chosen.

Root cause: the gradient is steepest where the coverage stops

Coverage does not end at a random place. It ends where the hydrodynamic model that produced the separation stops being trustworthy, which is where the tidal regime changes — typically at the shelf break, where the separation surface is bending fastest. Extrapolating from the last covered cell therefore takes the least stable value in the grid and projects it indefinitely.

import numpy as np
import rasterio

with rasterio.open("vdatum_navd88_to_mllw.tif") as src:
    sep = src.read(1, masked=True)

# The last two covered columns, and the gradient between them.
covered = ~sep.mask
last = np.argmax(covered[::-1].any(axis=0)[::-1])
edge = sep[:, -last - 2:-last]
print(float(np.ma.mean(np.diff(edge, axis=1))))     # 0.0141 m per cell at the edge
print(float(np.ma.mean(np.diff(sep[:, :50], axis=1))))  # 0.0009 m per cell nearshore

A gradient fifteen times steeper at the boundary than in the interior is the signature to expect. Holding that last value constant over 40 km of shelf is not a small approximation; it is an error that grows in exactly the region where nobody has a reference to check it against.

Step-by-step fix with production code

1. Classify, do not fill

Replace the binary valid/invalid flag with a three-way classification. The middle class is the one that carries the information a consumer needs:

import logging
from enum import IntEnum

import numpy as np
from scipy import ndimage

logger = logging.getLogger(__name__)


class Coverage(IntEnum):
    UNCOVERED = 0        # no separation available; the datum shift was not applied
    EDGE = 1             # covered, but within `edge_cells` of the coverage boundary
    INTERIOR = 2         # covered, and far enough in to be trusted


def classify_coverage(sep: np.ma.MaskedArray, *, edge_cells: int = 5) -> np.ndarray:
    """Three-way coverage classification for a sampled separation grid."""
    if edge_cells < 1:
        raise ValueError("edge_cells must be >= 1")
    covered = ~np.ma.getmaskarray(sep)
    if not covered.any():
        raise ValueError("no VDatum coverage anywhere in this extent")

    # Distance, in cells, from each covered cell to the nearest uncovered one.
    dist = ndimage.distance_transform_edt(covered)
    out = np.full(covered.shape, Coverage.UNCOVERED, dtype="uint8")
    out[covered] = Coverage.EDGE
    out[covered & (dist > edge_cells)] = Coverage.INTERIOR

    frac = {c.name: float((out == c).mean()) for c in Coverage}
    logger.info("coverage: interior %.1f%%, edge %.1f%%, uncovered %.1f%%",
                frac["INTERIOR"] * 100, frac["EDGE"] * 100, frac["UNCOVERED"] * 100)
    return out

The edge class exists because coverage is not a step function in quality. The last few cells inside the boundary are computed from a model that is already extrapolating internally, and a consumer computing a dredge volume will want to exclude them even though a valid number is present.

2. Handle a region seam by blending, not by choosing

Where two VDatum regions overlap, both grids have values and they disagree. Choosing one produces a visible step along the seam; blending across the overlap produces a surface that is wrong by less than the disagreement everywhere and by nothing at either end:

def blend_overlap(a: np.ma.MaskedArray, b: np.ma.MaskedArray) -> np.ma.MaskedArray:
    """Linear blend of two separation grids across the region where both exist."""
    both = ~np.ma.getmaskarray(a) & ~np.ma.getmaskarray(b)
    if not both.any():
        return np.ma.where(np.ma.getmaskarray(a), b, a)

    disagreement = float(np.abs(a[both] - b[both]).max())
    logger.info("region overlap: %d cells, max disagreement %.3f m",
                int(both.sum()), disagreement)
    if disagreement > 0.15:
        raise ValueError(
            f"regions disagree by {disagreement:.3f} m in their overlap — that is "
            "too large to blend away; check that both grids use the same epoch")

    # Weight by distance into the overlap so each grid dominates on its own side.
    da = ndimage.distance_transform_edt(both | ~np.ma.getmaskarray(a))
    db = ndimage.distance_transform_edt(both | ~np.ma.getmaskarray(b))
    w = np.divide(da, da + db, out=np.full(da.shape, 0.5), where=(da + db) > 0)
    blended = np.ma.where(both, a * (1 - w) + b * w, np.ma.where(np.ma.getmaskarray(a), b, a))
    return blended

The disagreement check is the important line. A blend hides a small inconsistency, which is what it is for; a large one means the two grids were built against different epochs or different geoid models, and hiding that is exactly wrong.

A Seam You Can See, or a Blend You Cannot A profile crossing the boundary between two adjacent VDatum regions whose separation grids overlap for a band in the middle. In the upper trace, region A is used up to the boundary and region B afterwards; the two values differ by about six centimetres so the trace shows a vertical step at the seam. In the lower trace, the two grids are blended across the overlap with a distance weighting; the trace matches region A exactly at the left edge of the overlap, region B exactly at the right edge, and passes smoothly between them with no step. A note records that a blend is only legitimate when the disagreement is small, and that a large one indicates the grids were built against different epochs. overlap choose one 6 cm step blend exact at both edges, smooth between Blend only when the disagreement is small. A large one means different epochs, and hiding it is the wrong fix.

3. Decide what the product says where it cannot answer

A coverage gap forces a product decision that is not the pipeline’s to make silently. There are three defensible answers and the difference between them is what the deliverable claims, not what the code does.

Clip the product to coverage. The published raster ends where VDatum ends. This is the strongest statement — everything in the file is referenced to the stated datum — and the weakest deliverable, because a client who asked for a survey to the shelf edge receives one that stops short.

Publish the full extent with an explicit datum band. Every cell carries which vertical reference it is on: MLLW inside coverage, the source datum outside. The file is complete, and no cell claims a datum it is not on. This needs a consumer who reads the band, and most do not.

Publish two products. The MLLW product clipped to coverage, and the source-datum product over the full extent, cross-referenced in both manifests. Twice the storage, no ambiguity, and the only option where a careless consumer cannot go wrong.

The choice is the deliverable owner’s. What the pipeline owes them is that the choice is available — which means the uncovered region must survive the transformation stage as an identifiable class rather than being resolved into a number early. Filling at the point of sampling forecloses all three options before anyone is asked.

def apply_shift_with_coverage(
    raw: np.ma.MaskedArray, sep: np.ma.MaskedArray, cov: np.ndarray,
) -> tuple[np.ma.MaskedArray, np.ndarray]:
    """Shift only where a separation exists; leave the rest on the source datum."""
    shifted = raw.copy()
    covered = cov != Coverage.UNCOVERED
    shifted[covered] = raw[covered] - sep[covered]
    datum_band = np.where(covered, DatumCode.MLLW, DatumCode.SOURCE).astype("uint8")
    logger.info("shift applied to %d of %d cells", int(covered.sum()), covered.size)
    return shifted, datum_band
Three Honest Deliverables, One Coverage Gap Three options laid out side by side for the same survey. Clip to coverage produces a smaller raster in which every cell is on MLLW; the cost is that the client asked for the shelf edge and does not get it. Publish the full extent with a datum band produces a complete raster plus a per-cell code saying which vertical reference each cell is on; the cost is that it depends on the consumer reading the band. Publish two products gives an MLLW raster clipped to coverage and a source-datum raster over the full extent, cross-referenced in both manifests; the cost is twice the storage and the benefit is that a careless consumer cannot misuse either. A footer states that the pipeline's job is to keep all three available by not resolving the gap into a number at sampling time. clip to coverage all MLLW strongest claim smaller than requested full extent + datum band MLLW source complete, nothing lies needs a consumer who reads it two products MLLW, clipped source datum, full extent cross-referenced manifests nobody can misuse either The choice belongs to whoever owns the deliverable. What the pipeline owes them is that all three stay possible — which means not resolving the gap into a number at sampling time. Fill it early and the only remaining option is the one nobody chose.

Verification and acceptance test 

import pytest

MAX_UNCOVERED_FRACTION = 0.05      # this product's spec


def test_uncovered_fraction_is_inside_spec(sep_grid):
    cov = classify_coverage(sep_grid)
    uncovered = float((cov == Coverage.UNCOVERED).mean())
    assert uncovered <= MAX_UNCOVERED_FRACTION, (
        f"{uncovered:.1%} of the survey has no VDatum coverage, over the "
        f"{MAX_UNCOVERED_FRACTION:.0%} allowed — extend the model or reduce the extent")


def test_no_shift_was_applied_where_there_is_no_separation(shifted, cov, raw):
    uncovered = cov == Coverage.UNCOVERED
    assert np.allclose(shifted[uncovered], raw[uncovered], equal_nan=True), (
        "cells with no separation were modified — something filled the gap")

The second test is the one that catches a fill introduced three refactors later. It asserts the negative: where there was no answer, nothing changed.

Edge cases and gotchas

  • A masked array that loses its mask. np.ma masks do not survive every NumPy operation, and a silently unmasked array fills with whatever the underlying buffer holds — often zeros. Carry the coverage classification as its own integer band rather than relying on the mask to persist.
  • Coverage is not the same as the survey footprint. A cell can be inside VDatum coverage and outside the surveyed area, or the reverse. Two independent masks, intersected explicitly, is clearer than one combined flag that nobody can decompose later.
  • Publish the class, not just the mask. A consumer who receives only valid/invalid cannot exclude the edge band, and the edge band is where a dredge quantity acquires most of its error. The three-way classification costs one byte per cell.

Up: Tidal Datum Transformations in Python