Tidal Datum Transformations in Python #

Automated vertical reference alignment is a non-negotiable prerequisite for coastal modeling, habitat mapping, and regulatory compliance. This workflow sits inside the Marine Spatial Data Fundamentals & Architecture foundation and defines the execution pattern for converting bathymetric and topographic elevations between established vertical reference frames such as NAVD88, MLLW, and MSL. The operational intent is precise: ingest a source elevation raster, align it spatially with a pre-computed tidal offset grid, apply the transformation under strict memory constraints, and serialize a cloud-ready dataset with validated compound CRS metadata. This pattern depends directly on the horizontal-plus-vertical conventions established in CRS alignment for coastal GIS projects — the vertical component of those compound strings is precisely what this transformation rewrites. Vertical datum integrity governs downstream hydrodynamic model boundary conditions and shoreline change detection accuracy; errors introduced here propagate silently through every dependent pipeline stage.

Reference Configuration and Specification #

Before writing a line of transformation code, fix the reference configuration. The table below defines the parameter envelope this workflow targets; values outside these ranges require re-evaluation of chunk strategy and resampling method.

Tidal datum corrections are spatially variable — not scalar offsets. They are derived from harmonic tide models, geoid undulations, and local benchmark surveys, then distributed as rasterized separation surfaces. The transformation this workflow applies is a per-pixel additive shift between the source and target vertical reference planes:

$$Z_{\text{target}}(x, y) = Z_{\text{source}}(x, y) + \Delta_{s \rightarrow t}(x, y)$$

where $\Delta_{s \rightarrow t}(x, y)$ is the signed separation surface sampled at each pixel — positive where the target datum lies above the source datum. Because $\Delta$ varies continuously across the spatial domain, both the source raster and the offset grid must share an identical pixel registration before the addition; the bulk of this workflow is the alignment and validation machinery that guarantees that condition.

The sign convention is the single most error-prone element of the entire transformation. The cross-section below fixes the geometry: each vertical reference plane sits at a different height, and the offset is the signed distance from the source plane to the target plane. Along most of the US Atlantic coast MLLW lies below NAVD88, so the NAVD88→MLLW offset is negative — adding a negative value lowers the elevation onto the MLLW plane.

Selecting the right container format determines pipeline throughput and memory footprint. For multi-dimensional tidal constituent storage, NetCDF carries native temporal axes; for static offset surfaces, GeoTIFF vs NetCDF format selection is the first decision gate before any transformation code runs.

Transformation Flow #

The diagram below shows data flow through the chunked transformation pipeline, from source datum ingestion through compound CRS validation, offset alignment, lazy arithmetic, and COG serialization.

Memory-Constrained Python Implementation #

The following implementation executes a fully chunked, lazy transformation pipeline. It uses rioxarray for spatial alignment, xarray with dask for deferred computation, and pyproj for compound CRS validation. The function signature is typed; every logging call uses the standard library logger rather than print().

import logging
from pathlib import Path

import dask
import rioxarray  # noqa: F401 — registers .rio accessor on xarray
import xarray as xr
from pyproj import CRS

log = logging.getLogger(__name__)


def transform_tidal_datum(
    source_raster: str | Path,
    offset_grid: str | Path,
    output_path: str | Path,
    target_vertical_epsg: str = "5829",
    chunk_size: int = 1024,
) -> None:
    """
    Apply a spatially variable tidal offset to a source elevation raster.

    The transformation is additive:
        Z_target = Z_source + offset
    where *offset* is the signed separation surface between the source datum
    and the target datum (e.g., the NAVD88-to-MLLW separation from VDatum).
    Sign convention: offset is positive where the target datum lies above the
    source datum. Always verify against the offset grid's own metadata.

    Args:
        source_raster: Path to the source elevation GeoTIFF or NetCDF.
        offset_grid: Path to the pre-computed offset surface (source → target).
        output_path: Destination path for the output Cloud-Optimized GeoTIFF.
        target_vertical_epsg: EPSG code for the output vertical datum
            (default "5829" = MLLW).
        chunk_size: Dask tile dimension in pixels; align with storage block size.

    Raises:
        ValueError: If the source raster has no CRS metadata.
        RuntimeError: If post-write CRS validation fails.
    """
    source_raster = Path(source_raster)
    output_path = Path(output_path)

    # ── 1. Load with explicit chunking — keeps full-resolution DEMs off RAM ──
    src = xr.open_dataarray(
        source_raster,
        engine="rasterio",
        chunks={"x": chunk_size, "y": chunk_size},
    )
    offset = xr.open_dataarray(
        offset_grid,
        engine="rasterio",
        chunks={"x": chunk_size, "y": chunk_size},
    )

    # ── 2. Validate horizontal CRS presence ───────────────────────────────────
    if src.rio.crs is None:
        raise ValueError(
            f"Source raster {source_raster} carries no CRS metadata. "
            "Embed the horizontal EPSG before running this pipeline."
        )
    src_crs = CRS.from_wkt(src.rio.crs.to_wkt())

    if offset.rio.crs is not None:
        offset_crs = CRS.from_wkt(offset.rio.crs.to_wkt())
        if src_crs.to_epsg() != offset_crs.to_epsg():
            log.warning(
                "CRS mismatch: source=%s, offset=%s — reprojecting offset grid.",
                src_crs.to_epsg(),
                offset_crs.to_epsg(),
            )
            offset = offset.rio.reproject(src_crs)
    else:
        log.warning("Offset grid carries no CRS; assuming it matches source CRS.")

    # ── 3. Spatial alignment — pixel-perfect resampling before arithmetic ─────
    #   Shape and bounds must match exactly; bilinear preserves smooth gradients
    #   in the offset surface without introducing step artifacts.
    needs_resample = (
        src.rio.shape != offset.rio.shape
        or src.rio.bounds() != offset.rio.bounds()
    )
    if needs_resample:
        log.info("Resampling offset grid to source spatial resolution (bilinear).")
        offset = offset.rio.reproject_match(src, resampling="bilinear")

    # ── 4. Apply vertical shift under Dask lazy evaluation ────────────────────
    #   No compute() call here — the DAG is materialised only at to_raster().
    transformed = (src + offset).rename("elevation_mllw")
    transformed.rio.write_crs(src_crs, inplace=True)

    # ── 5. Serialise as Cloud-Optimized GeoTIFF with compound CRS metadata ────
    transformed.rio.to_raster(
        str(output_path),
        driver="GTiff",
        tiled=True,
        blockxsize=256,
        blockysize=256,
        compress="DEFLATE",
        dtype="float32",
    )

    # ── 6. Post-write CRS smoke-check ─────────────────────────────────────────
    written_crs = xr.open_dataarray(str(output_path), engine="rasterio").rio.crs
    if written_crs is None:
        raise RuntimeError(
            f"Output {output_path} lost its CRS metadata during write. "
            "Inspect GDAL driver options and retry."
        )

    log.info("Transformation complete → %s  (vertical EPSG:%s)", output_path, target_vertical_epsg)

Scheduler configuration for constrained environments #

The pipeline runs safely in containerised workers when the Dask scheduler is configured before the pipeline function is called:

import dask

# Single-threaded: deterministic execution, zero GIL contention, safest for CI.
dask.config.set(scheduler="synchronous")

# Distributed: configure spill-to-disk threshold to match NVMe capacity.
# dask.config.set({"distributed.worker.memory.target": 0.7,
#                  "distributed.worker.memory.spill": 0.85})

For large regional DEMs (>10 GB), chunk_size=512 reduces scheduler overhead while maintaining acceptable cache locality. Never call .values or .compute() before the to_raster() call — doing so forces the entire transformed array into RAM and defeats the lazy evaluation architecture.

Validation Gates and Quality Control #

Every transformation output must pass three mandatory validation gates before it enters a downstream model or archive.

Gate 1 — Benchmark control-point drift #

Compare the transformed surface against a set of independently surveyed tidal benchmark elevations. A mean absolute error above ±0.05 m indicates an offset grid version mismatch or a sign-convention inversion.

import numpy as np
import xarray as xr


def validate_against_benchmarks(
    output_path: str,
    benchmarks: dict[str, tuple[float, float, float]],
    mae_threshold: float = 0.05,
) -> None:
    """
    Sample the transformed raster at benchmark locations and compute MAE.

    Args:
        output_path: Path to the transformed GeoTIFF.
        benchmarks: Mapping of station_id → (easting, northing, known_mllw_elevation).
        mae_threshold: Maximum acceptable mean absolute error in metres.

    Raises:
        AssertionError: If MAE exceeds the threshold, halting pipeline progression.
    """
    da = xr.open_dataarray(output_path, engine="rasterio")
    errors: list[float] = []

    for station_id, (x, y, truth) in benchmarks.items():
        sampled = float(da.sel(x=x, y=y, method="nearest").values)
        err = abs(sampled - truth)
        errors.append(err)
        log.info("  %s  sampled=%.4f m  truth=%.4f m  |err|=%.4f m", station_id, sampled, truth, err)

    mae = float(np.mean(errors))
    log.info("Benchmark MAE: %.4f m  (threshold ±%.3f m)", mae, mae_threshold)
    assert mae <= mae_threshold, (
        f"Benchmark MAE {mae:.4f} m exceeds ±{mae_threshold} m. "
        "Check offset grid version and sign convention."
    )

Gate 2 — NaN propagation audit #

After transformation, NaN coverage must not exceed the source raster’s NoData extent by more than a small tolerance. Unexpected NaN growth signals chunk boundary misalignment or a mismatched offset grid extent.

def check_nan_propagation(source_path: str, output_path: str, max_growth: float = 0.01) -> None:
    """Raise if NaN fraction grows by more than *max_growth* relative to source."""
    src = xr.open_dataarray(source_path, engine="rasterio")
    out = xr.open_dataarray(output_path, engine="rasterio")

    src_nan_frac = float(np.isnan(src).mean().compute())
    out_nan_frac = float(np.isnan(out).mean().compute())
    growth = out_nan_frac - src_nan_frac

    if growth > max_growth:
        raise RuntimeError(
            f"NaN fraction grew from {src_nan_frac:.4f} to {out_nan_frac:.4f} "
            f"(+{growth:.4f}). Check offset grid extent and resampling bounds."
        )
    log.info("NaN propagation check passed (growth=%.4f).", growth)

Gate 3 — Compound CRS metadata completeness #

Use gdalinfo or pyproj to confirm that both the horizontal and vertical components survive the write cycle intact:

import subprocess, json

def assert_crs_complete(output_path: str, expected_horizontal_epsg: int, expected_vertical_epsg: int) -> None:
    result = subprocess.run(
        ["gdalinfo", "-json", output_path],
        capture_output=True, text=True, check=True
    )
    info = json.loads(result.stdout)
    wkt = info.get("coordinateSystem", {}).get("wkt", "")
    assert str(expected_horizontal_epsg) in wkt, f"Horizontal EPSG:{expected_horizontal_epsg} not found in output WKT."
    assert str(expected_vertical_epsg) in wkt, f"Vertical EPSG:{expected_vertical_epsg} not found in output WKT."
    log.info("CRS completeness check passed (EPSG:%d + EPSG:%d).", expected_horizontal_epsg, expected_vertical_epsg)

For inter-datum conversions, consult the NOAA VDatum technical documentation to verify acceptable tolerance thresholds and regional transformation accuracy specifications.

Common Failure Modes and Diagnosis #

Failure 1 — Silent vertical datum stripping during reprojection #

Symptom: Output GeoTIFF passes gdalinfo horizontally but shows no vertical CRS entry. Downstream hydrodynamic solvers treat the surface as MSL regardless of intended datum.

Root cause: rioxarray.reproject() and rasterio.warp.reproject() do not automatically propagate compound CRS strings. The horizontal component is preserved; the vertical component is silently dropped unless explicitly re-written after the reproject call.

Remediation:

# After any reproject call, re-apply the full compound CRS.
da_reprojected = da.rio.reproject("EPSG:26918")
da_reprojected.rio.write_crs("EPSG:26918+5829", inplace=True)

Failure 2 — Offset sign inversion #

Symptom: Transformed elevations are systematically offset by 2× the expected correction in the wrong direction. Benchmark MAE exceeds 1 m.

Root cause: NOAA VDatum offset grids define separation as Z_MLLW = Z_NAVD88 + offset, where offset is typically negative along most of the US Atlantic coast (MLLW lies below NAVD88). Pipelines that invert the sign — computing Z_NAVD88 - offset — introduce a datum inversion error that equals twice the separation magnitude.

Remediation: Always verify offset grid polarity from its accompanying metadata file before production ingestion. Print the offset grid’s 5th and 95th percentile values and cross-reference against NOAA VDatum uncertainty estimates for the region.

Failure 3 — Chunk boundary seam artifacts in the output COG #

Symptom: Visual “grid lines” appear in the transformed DEM at regular intervals matching the Dask chunk size. Seams also appear as step discontinuities in depth contours.

Root cause: The offset grid and source raster were loaded with different chunk sizes, causing reproject_match to produce fractional-pixel misalignment at chunk borders before arithmetic is applied.

Remediation: Load both arrays with the same chunk_size value and verify that offset.rio.shape == src.rio.shape and offset.rio.bounds() == src.rio.bounds() before calling the arithmetic operator. If the shapes still diverge after reproject_match, inspect both datasets for non-standard pixel registrations (area vs. point).

Failure 4 — OOM crash when chunk_size exceeds available memory #

Symptom: Worker pod is OOM-killed during to_raster() with no Python traceback — only a container exit code 137.

Root cause: chunk_size=2048 on a float32 single-band raster consumes approximately 16 MB per chunk. With Dask’s default multi-threaded scheduler spinning up multiple simultaneous tasks, peak RAM can reach 10–20× the per-chunk figure during write.

Remediation: Reduce chunk_size to 256 or 512, or set the Dask scheduler to "synchronous" for single-worker containers. Monitor per-worker RSS via dask.distributed.get_worker().memory_monitor_interval if running on a distributed cluster.

Pipeline Integration and Downstream Handoff #

The tidal datum transformation step is a foundational node in coastal analysis pipelines. Its outputs feed three principal downstream stages:

Hydrodynamic model preprocessors. ADCIRC, SELFE, and Delft3D FM ingestors require MLLW-referenced bathymetric grids as boundary condition inputs. The COG output from this pipeline must be accompanied by a metadata manifest declaring the offset grid version, transformation date, horizontal CRS, and vertical datum EPSG code. Without this provenance record, modellers cannot verify which national tidal datum epoch was applied.

Shoreline change detection. Multi-temporal DEMs compared for shoreline position use MLLW as the consistent reference plane across survey epochs. Any epoch not transformed through this pipeline — or transformed with a different offset grid version — introduces a systematic vertical bias that masquerades as genuine coastal change.

Habitat suitability and inundation indices. Salt marsh, seagrass, and oyster habitat models compute tidal flooding frequency relative to MLLW. Rasters referenced to NAVD88 or MSL produce erroneous inundation extents; vertical misalignment by 0.1 m can shift predicted habitat boundaries by tens of metres in microtidal environments.

When applying this transformation to raw survey point clouds before rasterization, follow the step-by-step workflow in Applying MLLW to Coastal Survey Data, which covers point-level datum shifting, LAS/LAZ file handling, and progressive densification before grid interpolation.

For the broader context of CRS alignment decisions in coastal GIS projects — including compound EPSG string construction and axis-order verification — see the sibling page on CRS architecture, which establishes the pyproj compound string patterns this pipeline relies on.

Long-term archival requires OGC-compliant metadata packaging. Store transformed rasters alongside their source offset grids, transformation logs, and CRS validation reports in immutable object storage. Implement checksum verification (SHA-256) on all outputs and maintain a versioned registry of tidal constituent model releases used for offset grid generation. This ensures full reproducibility when reprocessing historical bathymetry under updated geodetic frameworks or corrected national tidal datum epochs.

---

Related #

• Applying MLLW to Coastal Survey Data — point-cloud and survey-file MLLW alignment before rasterization
• CRS Alignment for Coastal GIS Projects — compound CRS construction, horizontal/vertical separation, pyproj patterns
• Understanding NetCDF vs GeoTIFF for Marine Data — format selection for offset grids and multi-dimensional tidal constituent arrays
• Marine Spatial Data Fundamentals & Architecture — parent section covering the full data architecture