Coastal & Marine Spatial Analysis Pipelines

Applying MLLW to Coastal Survey Data

Coastal survey deliverables almost never arrive on the datum the charting or coastal-engineering team actually needs. Multibeam soundings come out of acquisition software on a local benchmark, airborne lidar arrives on NAVD88, and merged products drift onto Mean Sea Level — none of which line up with the Mean Lower Low Water (MLLW) reference that nautical charts, dredging volumes, and shoreline-change baselines are computed against. The task here is narrow and operational: take a raster already referenced to a known source datum and re-reference it to MLLW by adding a precomputed separation surface, doing it under strict memory limits so terabyte-scale bathymetry does not blow the worker’s RAM. This page is one step within the broader tidal datum transformation workflow; read that parent first for the full configuration envelope, library pins, and the harmonic-to-surface theory behind the offset grids consumed below.

Why the Naive Approach Corrupts the Output

MLLW is a tidal datum: it is the arithmetic mean of the lower of the two daily low waters over the 19-year National Tidal Datum Epoch. Because it is hydrodynamically defined rather than gravity-defined, the separation between MLLW and a geometric datum like NAVD88 is not a constant — it can swing by more than a metre between the head of an estuary and the open shelf only kilometres away. That single fact rules out the most common mistake: subtracting a scalar “tide correction” pulled from a single benchmark.

Operational pipelines never recompute MLLW from harmonic constituents at runtime. They consume a precomputed signed separation surface — a NOAA VDatum grid, a regional tidal model output, or an agency MLLW-to-NAVD88 raster — and apply it per pixel. The transformation is additive:

Z_MLLW = Z_raw + Offset

where Offset is the signed separation from the source datum to MLLW. The second silent corruption vector lives in that sign. A grid published as “MLLW relative to NAVD88” and a grid published as “NAVD88 relative to MLLW” differ by a global sign flip; apply the wrong one and every elevation is wrong by roughly twice the local tidal offset, with no error raised. The naive call that produces a plausible-but-wrong raster looks like this:

import rioxarray as rxr

survey = rxr.open_rasterio("raw_survey_navd88.tif")
offset = rxr.open_rasterio("navd88_to_mllw.tif")

# WRONG on three counts: no chunking (OOM on large extents),
# no CRS/grid alignment (broadcast against mismatched coords -> all-NaN),
# and no check of the offset's sign convention.
mllw = survey + offset

If the two grids do not share an exact coordinate index, xarray aligns on the intersection and the result is mostly NaN; if they do happen to overlap, the unverified sign quietly doubles the error. Both failures pass a casual min/max glance. The fix below makes alignment explicit, keeps memory bounded, and forces the sign convention to be a declared input.

Chunked MLLW Datum-Shift Data Flow Two lazy, Dask-chunked inputs — a NAVD88 survey raster and a signed VDatum separation surface — converge on a reproject_match alignment stage. The aligned grids feed an additive vertical shift Z_MLLW = Z_raw + sign × Offset, computed under lazy evaluation. The output is tagged with the MLLW vertical CRS (EPSG:5829) and serialized as a Cloud-Optimized GeoTIFF or NetCDF4. A verification stage samples tidal benchmark control points: if the mean residual is within ±0.05 m the run is accepted, otherwise the offset sign is flipped and the shift re-runs. Survey raster NAVD88 · chunked, lazy Dask array Separation surface VDatum · signed source→MLLW offset reproject_match offset → survey grid: CRS · res · extent Additive shift Z_MLLW = Z_raw + sign × Offset Tag + serialize MLLW CRS (EPSG:5829) → COG / NetCDF4 Verify benchmark MAE ≤ 0.05 m? yes Accept output MLLW-referenced grid no → flip offset_sign, re-run

Step-by-Step Fix With Production Code

1. Confirm the offset grid’s sign and source datum

Before any code, read the offset grid’s documentation and confirm two things: the source datum it converts from, and the direction the values point. Encode that as an explicit offset_sign so the pipeline cannot run on an assumption. For VDatum-derived grids, the published convention is documented in the NOAA VDatum documentation; authoritative vertical CRS codes (MLLW is EPSG:5829) come from the EPSG Geodetic Parameter Dataset.

2. Open both grids as chunked, lazy arrays

Loading a full coastal extent eagerly triggers an immediate OOM in any containerised worker. Open both rasters with explicit Dask chunks aligned to the storage block size so only the active window ever resides in RAM. This mirrors the memory-bounded ingestion pattern used across the tidal datum transformation workflow.

3. Validate the horizontal CRS, then reproject_match the offset

The offset surface must be coerced onto the survey grid’s exact CRS, resolution, extent, and origin before the arithmetic — otherwise the broadcast produces edge artifacts and NaN propagation. If the survey arrives on the wrong horizontal CRS entirely, reproject it first using the CRS alignment workflow rather than letting this pipeline guess.

4. Apply the additive shift and tag the vertical CRS

The shift itself is one lazy expression; the work is in attaching correct metadata so downstream GIS engines, hydrodynamic solvers, and charting systems read the output as MLLW rather than re-interpreting it as the source datum.

The complete, runnable pipeline:

import logging
import pathlib

import dask
import rioxarray  # noqa: F401 — registers the .rio accessor
import xarray as xr
from rasterio.crs import CRS
from rasterio.enums import Resampling

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


def apply_mllw_transform(
    survey_path: pathlib.Path,
    offset_path: pathlib.Path,
    output_path: pathlib.Path,
    *,
    chunk_size: int = 512,
    target_horizontal_crs: str = "EPSG:26918",
    target_vertical_crs: str = "EPSG:5829",  # MLLW
    offset_sign: int = 1,  # +1 if grid stores (MLLW - source); -1 otherwise
    resampling_method: Resampling = Resampling.bilinear,
) -> None:
    """Re-reference a coastal survey raster to MLLW with a precomputed offset surface.

    The offset surface must encode the signed separation between the survey's source
    datum and MLLW. Set `offset_sign` from the grid's documentation; it is never inferred.

    Parameters
    ----------
    survey_path : Raw bathymetry/elevation grid (GeoTIFF/NetCDF) on the source datum.
    offset_path : Precomputed source-to-MLLW separation surface.
    output_path : Destination for the MLLW-referenced output (.tif or .nc).
    chunk_size  : Square Dask chunk dimension, in pixels; tune to storage block size.
    target_horizontal_crs : Expected horizontal CRS of the survey grid.
    target_vertical_crs   : Target vertical EPSG (5829 = MLLW).
    offset_sign : +1 if the grid stores (MLLW - source), -1 for (source - MLLW).
    resampling_method : Resampling strategy for offset-grid alignment.
    """
    if offset_sign not in (1, -1):
        raise ValueError(f"offset_sign must be +1 or -1, got {offset_sign!r}")

    logger.info("Initializing MLLW transform: %s", survey_path.name)

    # Synchronous scheduler gives predictable peak memory on a single node.
    dask.config.set(scheduler="synchronous")

    survey_da = xr.open_dataarray(
        survey_path, engine="rasterio",
        chunks={"x": chunk_size, "y": chunk_size},
    )
    survey_crs = CRS.from_wkt(survey_da.rio.crs.to_wkt())

    if not survey_crs.equals(CRS.from_string(target_horizontal_crs)):
        raise ValueError(
            f"Survey CRS EPSG:{survey_crs.to_epsg()} does not match target "
            f"{target_horizontal_crs}. Reproject before applying the datum shift."
        )
    logger.info("Survey grid: shape=%s, CRS=EPSG:%s", survey_da.shape, survey_crs.to_epsg())

    offset_da = xr.open_dataarray(
        offset_path, engine="rasterio",
        chunks={"x": chunk_size, "y": chunk_size},
    )
    offset_crs = CRS.from_wkt(offset_da.rio.crs.to_wkt())

    # Coerce the offset onto the survey grid: same CRS, resolution, extent, origin.
    if not offset_crs.equals(survey_crs):
        logger.warning("Offset CRS mismatch; reprojecting offset to survey CRS.")
        offset_da = offset_da.rio.reproject(survey_crs, resampling=resampling_method)

    if (survey_da.rio.shape != offset_da.rio.shape
            or survey_da.rio.bounds() != offset_da.rio.bounds()):
        logger.info("Resampling offset surface to match survey grid.")
        offset_da = offset_da.rio.reproject_match(survey_da, resampling=resampling_method)

    # Apply the signed additive shift lazily: Z_MLLW = Z_raw + sign * Offset.
    logger.info("Applying chunked vertical datum shift (offset_sign=%+d).", offset_sign)
    mllw_da = (survey_da + offset_sign * offset_da).rename("elevation_mllw")

    mllw_da.attrs.update({
        "units": "meters",
        "long_name": "Elevation relative to Mean Lower Low Water",
        "vertical_datum": "MLLW",
        "vertical_datum_epsg": target_vertical_crs.split(":")[1],
        "source_datum_epsg": str(survey_crs.to_epsg()),
        "offset_grid": offset_path.name,
        "offset_sign": offset_sign,
    })
    mllw_da.rio.write_crs(survey_crs, inplace=True)

    logger.info("Serializing output to %s", output_path)
    if str(output_path).endswith(".nc"):
        mllw_da.to_netcdf(output_path, mode="w", format="NETCDF4", engine="netcdf4")
    else:
        mllw_da.rio.to_raster(
            output_path,
            driver="GTiff",
            tiled=True,
            blockxsize=256,
            blockysize=256,
            compress="DEFLATE",
            dtype="float32",
        )
    logger.info("MLLW transform complete.")


if __name__ == "__main__":
    apply_mllw_transform(
        survey_path=pathlib.Path("data/raw_survey_navd88.tif"),
        offset_path=pathlib.Path("data/offsets/navd88_to_mllw.tif"),
        output_path=pathlib.Path("data/processed/survey_mllw.tif"),
        offset_sign=1,
    )

The choice between Cloud-Optimized GeoTIFF and NetCDF4 output is itself a downstream-handoff decision; see NetCDF vs GeoTIFF for marine data for which to emit for charting versus model-forcing consumers.

Verification and Acceptance Test

Plausible output is not correct output. Confirm the shift against independent tidal benchmark control points — stations where the published MLLW-to-source separation is known — and assert the residual stays within the workflow’s acceptance threshold (±0.05 m against benchmarks):

import numpy as np
import rioxarray as rxr

mllw = rxr.open_rasterio("data/processed/survey_mllw.tif").squeeze()

# control_points: list of (x, y, expected_mllw_elevation) at benchmark stations
control_points = [(412350.0, 4516120.0, -3.214), (418900.0, 4519880.0, -1.087)]

residuals = []
for x, y, expected in control_points:
    sampled = float(mllw.sel(x=x, y=y, method="nearest"))
    residuals.append(sampled - expected)

mae = float(np.mean(np.abs(residuals)))
logger.info("Benchmark residuals (m): %s | MAE=%.4f", [round(r, 4) for r in residuals], mae)
assert mae <= 0.05, f"MLLW shift failed acceptance: MAE {mae:.4f} m exceeds 0.05 m"

# Guard against the silent alignment failure: the output must not be mostly NaN.
valid_fraction = float((~np.isnan(mllw)).mean())
assert valid_fraction > 0.95, f"Output {1 - valid_fraction:.1%} NaN — check offset alignment"

A large MAE that is roughly twice the local separation is the fingerprint of an inverted offset_sign; flip it and re-run. A high NaN fraction means the offset grid never aligned to the survey grid.

Edge Cases and Gotchas

  • Offset coverage gaps blank valid bathymetry. VDatum and regional tidal models stop at their model boundary, so the separation surface carries NaN offshore or up narrow channels. Adding NaN blanks otherwise-valid soundings. Fill small gaps with a nearest-neighbour extrapolation of the offset before the arithmetic, and flag pixels outside model coverage rather than silently emitting them.
  • NoData sentinels masquerade as elevations. If the source raster uses a sentinel like -9999 instead of a declared NoData mask, reproject_match bilinear resampling smears it into neighbouring real values. Call survey_da.rio.write_nodata(..., inplace=True) and convert to NaN before resampling so the sentinel never enters the arithmetic.
  • Vertical CRS is dropped on round-trip. GeoTIFF only carries the horizontal CRS in its standard tags; the MLLW vertical reference lives in the array attributes and any compound-CRS WKT you embed. Tools that re-open the file and ignore attributes will treat it as the source datum again. For regulatory exchange where the vertical datum must travel with the data, prefer NetCDF4 with a grid_mapping variable, and never assume EPSG:5829 survives a careless gdal_translate.

Up: Tidal Datum Transformations in Python