Coastal & Marine Spatial Analysis Pipelines

Converting NAVD88 to MLLW with VDatum

Airborne lidar, RTK topo-bathy, and merged coastal DEMs almost always arrive on NAVD88 — an orthometric, gravity-referenced datum — while nautical charts, dredging pay-volume calculations, and shoreline baselines are computed against Mean Lower Low Water (MLLW), a tidal datum. Converting between them is not a scalar tide correction: the separation between NAVD88 and MLLW is a surface that varies spatially, sometimes by more than a metre across a single estuary. This page is the point-and-grid conversion recipe within the broader tidal datum transformation workflow; it complements Applying MLLW to Coastal Survey Data, which handles the raster-to-raster case, by focusing on obtaining and sampling the NOAA VDatum separation surface per point with correct sign and coverage handling.

Why a Single NAVD88-to-MLLW Offset Fails Across a Survey A schematic coastal cross-section runs left (open shelf) to right (estuary head). A horizontal line marks NAVD88. A second line marks MLLW; it starts close to NAVD88 offshore and tilts increasingly below it toward the estuary head, so the vertical separation sep(x,y) grows from about 0.4 metres to about 1.4 metres. A dashed line shows a single constant 0.4 metre offset taken from an offshore benchmark; the gap between that dashed line and the true MLLW surface reaches about one metre at the estuary head, which is the error a constant offset injects. The caption reads Z_MLLW = Z_NAVD88 minus sep(x,y). elevation open shelf estuary head NAVD88 MLLW (true) single benchmark constant (0.4 m) 0.4 m 1.4 m ~1 m error Z_MLLW = Z_NAVD88 − sep(x, y)

Root cause: the separation is a surface, not a number

NAVD88 is fixed to the land by geodetic levelling and the geoid; MLLW is fixed to the water by 19 years of tidal observation over the National Tidal Datum Epoch. The two references drift apart as tidal range changes — amplified up estuaries, in convergent bays, and behind barrier islands. The vertical relationship at a point is

ZMLLW(x,y)=ZNAVD88(x,y)−sep⁡(x,y),Z_\text{MLLW}(x,y) = Z_\text{NAVD88}(x,y) - \operatorname{sep}(x,y),

where sep⁡(x,y)\operatorname{sep}(x,y) is the NAVD88-minus-MLLW separation delivered by NOAA VDatum. Note the sign: where MLLW sits below NAVD88 (the usual coastal case), sep⁡\operatorname{sep} is positive and subtracting it makes the MLLW elevation more negative — a deeper chart depth. Because sep⁡\operatorname{sep} varies, a single constant cc pulled from one benchmark produces a residual that grows with distance:

ε(x,y)=sep⁡(x,y)−c,\varepsilon(x,y) = \operatorname{sep}(x,y) - c,

which is exactly the line-to-line mismatch surveyors see when two overlapping sheets were reduced with different tide stations. Applying one number looks deceptively fine:

import numpy as np

# Elevations already on NAVD88 (metres), spanning shelf to estuary head.
z_navd88 = np.array([-5.10, -4.02, -3.15, -1.88])

# WRONG: one offset from an offshore benchmark applied everywhere.
CONSTANT_OFFSET_M = 0.40
z_mllw_bad = z_navd88 - CONSTANT_OFFSET_M
# Offshore points are ~right; the estuary-head point is off by ~1 m because
# the true separation there is ~1.4 m, not 0.4 m. No error is raised.

The array comes back full of plausible numbers, and nothing flags that the last value is a metre wrong. The fix is to sample the separation per point from the VDatum grid. Getting there first requires a trustworthy horizontal CRS on the points, which is the job of the CRS alignment workflow and, for the code that defines EPSG identity, Geodetic Datums & EPSG Management.

Step-by-step fix with production code

1. Obtain the VDatum separation grid

Pull the NAVD88-to-MLLW separation surface for your region from NOAA VDatum, either as the offline regional grids or via the VDatum REST API. Record the grid’s TSS/geoid model version alongside it — a separation produced under one geoid realisation must not be mixed with elevations tied to another. Treat the grid as a versioned pipeline input, not an ambient download.

2. Load the grid and sample it per point

Open the grid with rasterio, transform the survey points into the grid’s CRS, and bilinearly sample the separation at each location. The sampler raises on empty input and flags points that fall on NoData.

import logging

import numpy as np
import rasterio
from rasterio.enums import Resampling
from pyproj import Transformer

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s | %(levelname)s | %(message)s",
    datefmt="%Y-%m-%dT%H:%M:%S",
)
logger = logging.getLogger(__name__)


def sample_separation(
    grid_path: str,
    lon: np.ndarray,
    lat: np.ndarray,
    points_crs: str = "EPSG:4326",
) -> np.ndarray:
    """Bilinearly sample a VDatum separation grid at lon/lat points.

    Returns an array of separations in metres with np.nan where a point
    falls on the grid's NoData (i.e. outside VDatum coverage). Raises on
    shape mismatch so a truncated coordinate stream cannot pass silently.
    """
    lon = np.asarray(lon, dtype="float64")
    lat = np.asarray(lat, dtype="float64")
    if lon.shape != lat.shape:
        raise ValueError(f"lon/lat shape mismatch: {lon.shape} vs {lat.shape}")
    if lon.size == 0:
        raise ValueError("No points supplied to sample_separation")

    with rasterio.open(grid_path) as src:
        transformer = Transformer.from_crs(points_crs, src.crs, always_xy=True)
        gx, gy = transformer.transform(lon, lat)
        nodata = src.nodata if src.nodata is not None else np.nan

        # rasterio.sample expects (x, y) pairs in the raster CRS.
        samples = np.fromiter(
            (v[0] for v in src.sample(zip(gx, gy), resampling=Resampling.bilinear)),
            dtype="float64",
            count=lon.size,
        )

    sep = np.where(samples == nodata, np.nan, samples)
    n_missing = int(np.isnan(sep).sum())
    if n_missing:
        logger.warning("%d of %d points fell outside VDatum coverage",
                       n_missing, sep.size)
    return sep

3. Apply the datum shift with an explicit sign convention

The separation direction is a declared input, never an assumption. Encode which way the grid points, apply ZMLLW=ZNAVD88−sep⁡Z_\text{MLLW} = Z_\text{NAVD88} - \operatorname{sep}, and refuse to emit points with an unknown separation.

def navd88_to_mllw(
    z_navd88: np.ndarray,
    separation: np.ndarray,
    *,
    sep_is_navd88_minus_mllw: bool = True,
    allow_uncovered: bool = False,
) -> np.ndarray:
    """Convert NAVD88 elevations to MLLW using a per-point separation.

    Parameters
    ----------
    z_navd88 : Elevations on NAVD88, metres, positive up.
    separation : Per-point separation from `sample_separation`, metres.
    sep_is_navd88_minus_mllw : True if the grid stores (NAVD88 - MLLW);
        set False for a grid published as (MLLW - NAVD88).
    allow_uncovered : If False, raise when any point has a NaN separation
        (outside VDatum coverage) instead of emitting NaN elevations.
    """
    z_navd88 = np.asarray(z_navd88, dtype="float64")
    separation = np.asarray(separation, dtype="float64")
    if z_navd88.shape != separation.shape:
        raise ValueError(
            f"shape mismatch: z {z_navd88.shape} vs sep {separation.shape}"
        )

    uncovered = np.isnan(separation)
    if uncovered.any() and not allow_uncovered:
        raise ValueError(
            f"{int(uncovered.sum())} points outside VDatum coverage; "
            "supply a wider grid or pass allow_uncovered=True to flag them."
        )

    signed = separation if sep_is_navd88_minus_mllw else -separation
    z_mllw = z_navd88 - signed
    logger.info("Converted %d points to MLLW (uncovered flagged as NaN=%d)",
                z_mllw.size, int(uncovered.sum()))
    return z_mllw

4. Wire it together for a survey 

lon = np.array([-75.20, -75.14, -75.09, -75.03])
lat = np.array([38.62, 38.70, 38.78, 38.86])
z_navd88 = np.array([-5.10, -4.02, -3.15, -1.88])

sep = sample_separation("data/vdatum/navd88_mllw_MAgrid.tif", lon, lat)
z_mllw = navd88_to_mllw(z_navd88, sep, sep_is_navd88_minus_mllw=True)

For raster-to-raster conversion of an entire DEM rather than a point list, hand the aligned separation surface to the chunked pipeline in Applying MLLW to Coastal Survey Data, which keeps memory bounded for terabyte-scale grids.

Verification and acceptance test

Confirm the conversion at an independent tidal benchmark — a NOAA station where the published NAVD88-to-MLLW relationship is known — and assert the residual is within tolerance. A residual roughly twice the local separation is the fingerprint of an inverted sign.

import pytest


def test_conversion_matches_tide_gauge_benchmark():
    # Lewes, DE (station 8557380): published NAVD88->MLLW separation ~1.40 m,
    # so a benchmark at NAVD88 0.00 m should sit at MLLW ~ -1.40 m.
    bench_lon, bench_lat = np.array([-75.1200]), np.array([38.7827])
    z_navd88 = np.array([0.00])
    expected_mllw = -1.40
    tol = 0.05  # metres

    sep = sample_separation("data/vdatum/navd88_mllw_MAgrid.tif",
                            bench_lon, bench_lat)
    z_mllw = navd88_to_mllw(z_navd88, sep, sep_is_navd88_minus_mllw=True)

    residual = float(z_mllw[0] - expected_mllw)
    assert abs(residual) <= tol, (
        f"benchmark residual {residual:+.3f} m exceeds {tol} m — "
        "check separation grid version or sign convention"
    )


def test_uncovered_points_raise_by_default():
    z = np.array([-3.0])
    sep = np.array([np.nan])          # point outside VDatum coverage
    with pytest.raises(ValueError, match="outside VDatum coverage"):
        navd88_to_mllw(z, sep)

Edge cases and gotchas

  • Outside VDatum coverage. VDatum stops at its hydrodynamic model boundary, so points offshore of the shelf model or up ungauged tidal creeks return NoData. Never let that NaN propagate into published depths — flag those points, clip them, or extend coverage with a documented regional model, but do not silently blank valid soundings.
  • TSS/geoid model version mismatch. The separation depends on the tidal surface (TSS) and geoid realisation VDatum was built with. Mixing a separation grid from one release with elevations tied to a different geoid injects a systematic bias. Pin the VDatum/geoid version as a pipeline input and record it in the output metadata.
  • Sign and depth-vs-elevation conventions. Hydrographic products are frequently depth-positive (down is positive) while topographic elevations are elevation-positive (up is positive). Decide one convention, convert at ingestion, and keep the separation sign consistent with it — the sep_is_navd88_minus_mllw flag exists precisely so this is never guessed.
  • IGLD and other datums. In the Great Lakes the reference is IGLD 85, not MLLW, and estuarine projects sometimes need LWD or a local low-water datum instead. VDatum can target those too, but the separation grid — and therefore the subtraction — is different; do not reuse an MLLW grid for an IGLD conversion. Resolving which vertical datum a deliverable truly needs is part of Geodetic Datums & EPSG Management.

Up: Tidal Datum Transformations in Python