CRS Alignment for Coastal GIS Projects

Operational intent: This page defines the end-to-end CRS alignment workflow that standardises heterogeneous marine and coastal datasets into a single project-specific spatial reference before automated spatial analysis. It sits within the Marine Spatial Data Fundamentals & Architecture framework, where CRS alignment is the mandatory pipeline gate that precedes every downstream operation — from bathymetric differencing to hydrodynamic model forcing. Coastal pipelines routinely ingest bathymetric rasters, vessel telemetry, habitat polygons, and tidal gauge records. When those inputs carry divergent datums, ellipsoids, or projection definitions, spatial operations fail silently: intersection geometries drift, raster overlays misregister, and spatial indexes return false negatives.

---

Reference Configuration and Specification Table #

Marine spatial operations require explicit compound coordinate reference systems that decouple horizontal positioning from vertical referencing. Relying on standalone EPSG codes without explicit datum context introduces vertical offsets that invalidate coastal elevation models, flood inundation simulations, and bathymetric differencing operations.

The table below captures the canonical configuration for a US East Coast coastal project. Substitute UTM zone, geoid model, and vertical datum as your regional footprint dictates.

A compound CRS is the mechanism that makes this work: it pairs an independent horizontal reference (the projected grid that fixes easting/northing) with an independent vertical reference (the datum that fixes the meaning of zero). The diagram below shows why a bare horizontal EPSG code is insufficient for coastal work — the same ellipsoidal height resolves to wildly different “depths” depending on which vertical datum is layered beneath it.

Metadata extraction strategies differ significantly between formats. The NetCDF vs GeoTIFF format guide details how to interrogate each container’s CRS metadata before feeding it into the aligner below.

Memory-Constrained Python Implementation #

Production pipelines cannot afford to load multi-gigabyte bathymetric mosaics or decade-long telemetry archives into RAM. The implementation below enforces strict memory boundaries through windowed raster processing and chunked vector transformation. It uses pyproj for deterministic datum shifts and rasterio for tile-based I/O. Each engineering decision is annotated inline.

"""
crs_aligner.py — production CRS alignment for coastal marine pipelines.

Dependencies:
    pip install pyproj>=3.6 rasterio>=1.3 geopandas>=0.14 pyogrio shapely>=2.0
"""

from __future__ import annotations

import json
import logging
import time
from dataclasses import dataclass, field
from pathlib import Path
from typing import Optional

import numpy as np
import rasterio
from rasterio.crs import CRS
from rasterio.warp import Resampling, calculate_default_transform, reproject
import geopandas as gpd
from pyproj import Transformer
from pyproj import CRS as ProjCRS
from shapely.validation import make_valid

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

DRIFT_HARD_LIMIT_M: float = 0.05   # reject threshold
DRIFT_WARN_LIMIT_M: float = 0.025  # warning threshold (pyproj accuracy estimate)


@dataclass
class AlignmentManifest:
    """Provenance record written alongside every aligned output."""
    source_crs: str
    target_crs: str
    transformation_timestamp: str
    pyproj_version: str
    rasterio_version: str
    max_drift_m: float
    validation_passed: bool
    notes: list[str] = field(default_factory=list)

    def write(self, path: Path) -> None:
        path.write_text(json.dumps(self.__dict__, indent=2))
        logger.info("Manifest written: %s", path)


class CRSAligner:
    """
    Aligns heterogeneous coastal datasets to a single compound CRS.

    Parameters
    ----------
    target_h_crs:
        EPSG string or PROJ string for the target horizontal CRS.
        Example: 'EPSG:32618' for WGS 84 / UTM zone 18N.
    chunk_size:
        Tile height/width in pixels for windowed raster reprojection.
        Default 4096 keeps peak RSS under ~512 MB per band at float32.
    resampling:
        rasterio Resampling kernel. Use Resampling.bilinear for raw
        bathymetry; Resampling.lanczos for slope/curvature derivatives.
    """

    def __init__(
        self,
        target_h_crs: str,
        chunk_size: int = 4096,
        resampling: Resampling = Resampling.bilinear,
    ) -> None:
        self.target_h: CRS = CRS.from_string(target_h_crs)
        self.chunk: int = chunk_size
        self.resampling = resampling
        # Source CRS is detected dynamically; this transformer is rebuilt
        # per dataset inside the transform methods using Transformer.from_crs.
        self._target_proj: ProjCRS = ProjCRS.from_user_input(target_h_crs)

    # ------------------------------------------------------------------
    # Vector alignment
    # ------------------------------------------------------------------

    def transform_vector(
        self,
        src_path: Path,
        dst_path: Path,
        control_points: Optional[np.ndarray] = None,
    ) -> AlignmentManifest:
        """
        Reproject a vector dataset and validate coordinate drift.

        Parameters
        ----------
        src_path:   Input GeoPackage / Shapefile path.
        dst_path:   Output GeoPackage path.
        control_points:
            Optional (N, 4) array of [src_x, src_y, ref_x, ref_y] in
            source and target CRS units respectively. Omit for datasets
            without known control.
        """
        if src_path.suffix not in {".gpkg", ".shp", ".geojson", ".fgb"}:
            raise ValueError(f"Unsupported vector format: {src_path.suffix}")

        gdf = gpd.read_file(src_path, engine="pyogrio")
        if gdf.crs is None:
            raise ValueError(f"Source vector lacks CRS definition: {src_path}")

        source_crs_str = gdf.crs.to_wkt()
        aligned = gdf.to_crs(self.target_h)

        # Repair topology after datum shift — high-curvature coastal
        # boundaries occasionally produce self-intersections post-reproject.
        aligned["geometry"] = aligned["geometry"].apply(make_valid)

        invalid_count = (~aligned.geometry.is_valid).sum()
        if invalid_count > 0:
            raise RuntimeError(
                f"{invalid_count} geometries remain invalid after make_valid(). "
                "Inspect source topology before proceeding."
            )

        aligned.to_file(dst_path, driver="GPKG", engine="pyogrio")
        logger.info("Vector aligned: %s → %s", src_path, dst_path)

        max_drift = self._validate_control_points(gdf.crs, control_points)

        return self._build_manifest(
            source_crs_str,
            max_drift,
            notes=[f"Repaired {invalid_count} invalid geometries via make_valid()"],
        )

    # ------------------------------------------------------------------
    # Raster alignment (windowed, memory-bounded)
    # ------------------------------------------------------------------

    def transform_raster_chunked(
        self,
        src_path: Path,
        dst_path: Path,
        control_points: Optional[np.ndarray] = None,
    ) -> AlignmentManifest:
        """
        Windowed raster reprojection. Memory ceiling ≈ chunk_size² × bands × 4 bytes.

        rasterio.warp.reproject handles the affine transform and resampling
        kernel natively. We compute output dimensions with
        calculate_default_transform so the output grid is pixel-aligned to
        the target CRS without rounding bias.
        """
        with rasterio.open(src_path) as src:
            if src.crs is None:
                raise ValueError(f"Source raster lacks CRS definition: {src_path}")

            source_crs_str = src.crs.to_wkt()
            src_crs = src.crs   # retain for post-write drift validation

            dst_transform, dst_width, dst_height = calculate_default_transform(
                src.crs, self.target_h, src.width, src.height, *src.bounds
            )

            profile = src.profile.copy()
            profile.update(
                crs=self.target_h,
                transform=dst_transform,
                width=dst_width,
                height=dst_height,
                dtype=rasterio.float32,
                compress="lzw",
                predictor=2,   # horizontal differencing — effective for float32
                tiled=True,
                blockxsize=self.chunk,
                blockysize=self.chunk,
            )

            with rasterio.open(dst_path, "w", **profile) as dst:
                for band_idx in range(1, src.count + 1):
                    reproject(
                        source=rasterio.band(src, band_idx),
                        destination=rasterio.band(dst, band_idx),
                        src_transform=src.transform,
                        src_crs=src.crs,
                        dst_transform=dst_transform,
                        dst_crs=self.target_h,
                        resampling=self.resampling,
                    )
                    logger.info(
                        "Raster band %d/%d reprojected: %s", band_idx, src.count, src_path
                    )

        # Reuse the source CRS captured above — do not reopen the file
        # (a second rasterio.open would leak an unclosed dataset handle).
        max_drift = self._validate_control_points(src_crs, control_points)
        return self._build_manifest(source_crs_str, max_drift)

    # ------------------------------------------------------------------
    # Validation helpers
    # ------------------------------------------------------------------

    def _validate_control_points(
        self,
        src_crs: CRS,
        control_points: Optional[np.ndarray],
    ) -> float:
        """
        Compute Euclidean drift (metres) between transformed control points
        and published reference coordinates.

        Returns max drift across all provided control points, or 0.0 if none.
        """
        if control_points is None or len(control_points) == 0:
            return 0.0

        transformer = Transformer.from_crs(src_crs, self.target_h, always_xy=True)
        src_x, src_y = control_points[:, 0], control_points[:, 1]
        ref_x, ref_y = control_points[:, 2], control_points[:, 3]

        tx, ty = transformer.transform(src_x, src_y)
        drift = np.sqrt((tx - ref_x) ** 2 + (ty - ref_y) ** 2)
        rms = float(np.sqrt(np.mean(drift ** 2)))

        logger.info(
            "Control-point drift — max: %.4f m, rms: %.4f m, n=%d",
            drift.max(), rms, len(drift),
        )

        if drift.max() > DRIFT_HARD_LIMIT_M:
            raise RuntimeError(
                f"Coordinate drift {drift.max():.4f} m exceeds hard limit "
                f"{DRIFT_HARD_LIMIT_M} m. Check datum grid installation and "
                "PROJ_DATA path. Pipeline halted."
            )
        return float(drift.max())

    def _build_manifest(
        self,
        source_crs: str,
        max_drift: float,
        notes: Optional[list[str]] = None,
    ) -> AlignmentManifest:
        import pyproj
        return AlignmentManifest(
            source_crs=source_crs,
            target_crs=self.target_h.to_wkt(),
            transformation_timestamp=time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()),
            pyproj_version=pyproj.__version__,
            rasterio_version=rasterio.__version__,
            max_drift_m=max_drift,
            validation_passed=max_drift <= DRIFT_HARD_LIMIT_M,
            notes=notes or [],
        )

Validation Gates and Quality Control #

Transformation without validation introduces silent errors. The pipeline enforces a three-stage hard gate before any aligned dataset proceeds to downstream analysis.

Stage 1: Control-Point Drift Verification #

Extract known survey markers, tide gauge benchmarks, or GNSS control stations from the National Geodetic Survey (NGS) or NOAA CO-OPS databases. Apply the pipeline transformer to their source coordinates and compute the planar offset against published target coordinates.

For each control point $i$, the drift is the Euclidean distance between the transformed coordinate $(\hat{x}_i, \hat{y}_i)$ and the published reference $(x_i^{\text{ref}}, y_i^{\text{ref}})$ in target-CRS metres:

$$d_i = \sqrt{(\hat{x}_i - x_i^{\text{ref}})^2 + (\hat{y}_i - y_i^{\text{ref}})^2}$$

The pipeline gates on the maximum drift rather than the mean, so that a single outlier benchmark cannot be masked by an otherwise tight fit:

$$d_{\max} = \max_{1 \le i \le N} d_i \le 0.05\ \text{m}$$

It also reports the root-mean-square drift as a population-level summary of transformation quality:

$$d_{\text{rms}} = \sqrt{\frac{1}{N} \sum_{i=1}^{N} d_i^{\,2}}$$

The _validate_control_points method above implements the $d_{\max}$ gate; load your NGS control as a NumPy array:

import numpy as np

# Columns: [src_lon, src_lat, ref_east_m, ref_north_m]
# Source in WGS84 geographic; reference (published NGS) in target projected CRS.
ngs_control = np.array([
    [-75.1234, 38.5678, 412345.12, 4271890.34],
    [-75.2345, 38.6789, 402100.88, 4283210.67],
])

aligner = CRSAligner("EPSG:32618")
manifest = aligner.transform_raster_chunked(
    Path("bathy_raw.tif"),
    Path("bathy_aligned.tif"),
    control_points=ngs_control,
)
manifest.write(Path("bathy_aligned_manifest.json"))

Stage 2: Vertical Offset Audit #

Compare transformed raster elevations against NOAA VDatum reference surfaces at tide gauge benchmark locations. NOAA VDatum uncertainty guidelines define acceptable tolerances for coastal zones; deviations greater than 0.10 m trigger a pipeline halt. The pyproj library exposes estimated transformation accuracy via Transformer.get_last_used_operation():

from pyproj import Transformer

t = Transformer.from_crs("EPSG:4326", "EPSG:32618", always_xy=True)
t.transform(-75.1234, 38.5678)
op = t.get_last_used_operation()
logger.info("Transformer accuracy estimate: %s m", op.accuracy)

if op.accuracy and op.accuracy > DRIFT_WARN_LIMIT_M:
    logger.warning(
        "Transformer accuracy %.4f m exceeds warning threshold. "
        "Verify proj-data datum grids are installed.",
        op.accuracy,
    )

Missing datum grids are the most common cause of unexpectedly large accuracy estimates. Install them with:

pip install proj-data
# Or for PROJ native:
projsync --source-id us_noaa

Stage 3: Topology Integrity Check #

Run shapely.is_valid on transformed vector geometries before writing to the output GeoPackage. Datum shifts occasionally produce self-intersecting polygons or collapsed linestrings in high-curvature coastal boundaries — especially for polygons representing bay shorelines or inlet channels. The implementation above calls make_valid() then re-checks and raises if any geometry remains invalid.

For interactive debugging of CRS mismatches during development, fixing CRS mismatches in QGIS and GeoPandas covers diagnostic patterns that complement this programmatic approach.

Common Failure Modes and Diagnosis #

1. Axis-Order Inversion (Lat/Lon vs Lon/Lat) #

Symptom: Transformed coordinates appear in the ocean or on a different continent; drift exceeds tens of thousands of metres.

Root cause: EPSG definitions for geographic CRSs follow the ISO 19111 standard and place latitude before longitude. Older pyproj versions and GDAL bindings sometimes honour this; newer ones do not unless always_xy=True is set. A single missing flag silently transposes every coordinate.

Fix: Always construct transformers with always_xy=True. Verify with a known point:

t = Transformer.from_crs("EPSG:4326", "EPSG:32618", always_xy=True)
east, north = t.transform(-75.0, 38.5)   # lon, lat order
assert 300_000 < east < 800_000, f"Axis order suspect: east={east}"

2. Missing Datum Grid (Geoid or Tidal Shift File) #

Symptom: pyproj raises ProjError: Cannot find proj.db or returns NaN coordinates; Transformer.get_last_used_operation().accuracy is -1 (unknown).

Root cause: The vertical datum shift from ellipsoidal height to NAVD88 or MLLW requires grid files (.tif or legacy .gtx) in the PROJ data directory. If those files are absent, PROJ either silently falls back to an approximate transformation or raises an error depending on the error_if_best_transformation_not_available flag.

Fix:

from pyproj import Transformer

# Raise immediately if the best available transformation cannot be used
t = Transformer.from_crs(
    "EPSG:4326+5703",
    "EPSG:32618+5703",
    always_xy=True,
    error_if_best_transformation_not_available=True,
)

Then install the missing grids: projsync --source-id us_noaa.

3. CRS-Less Legacy Datasets #

Symptom: src.crs is None (rasterio) or gdf.crs is None (geopandas); the pipeline raises ValueError immediately.

Root cause: Older bathymetric surveys (pre-2000) and some agency data transfers omit the .prj sidecar or GEOGCS tag in GeoTIFF metadata.

Fix: Assign the CRS from known provenance metadata before running the aligner. Do not guess:

import rasterio
from rasterio.crs import CRS

with rasterio.open("bathy_no_crs.tif") as src:
    profile = src.profile
    data = src.read()

profile["crs"] = CRS.from_epsg(4326)   # confirmed from data provider metadata
with rasterio.open("bathy_crs_assigned.tif", "w", **profile) as dst:
    dst.write(data)

Document the provenance of the CRS assignment in the output manifest’s notes field.

4. Silent Integer Overflow in Band Dtype Mismatch #

Symptom: Reprojected bathymetry shows large blocks of NoData or extreme depth spikes; drift validation passes because control points are correct.

Root cause: Source rasters stored as int16 carrying depth values in centimetres (e.g., −3200 cm) are read by rasterio as signed 16-bit integers. The profile.update(dtype=rasterio.float32) conversion in the aligner handles this correctly — but if the profile update is skipped and Resampling.bilinear is applied to integer data, bilinear interpolation can produce float values that are then truncated back to int16, corrupting border pixels.

Fix: Always cast to float32 before reprojection. The aligner above does this unconditionally in profile.update.

Pipeline Integration and Downstream Handoff #

Aligned datasets are staged in a standardised directory hierarchy with the following structure:

pipeline/
├── aligned/
│   ├── bathy_aligned.tif
│   ├── bathy_aligned_manifest.json
│   ├── habitat_aligned.gpkg
│   └── habitat_aligned_manifest.json
└── raw/
    ├── bathy_raw.tif
    └── habitat_raw.gpkg

Each output file is accompanied by an AlignmentManifest JSON record (schema above) containing source CRS, target CRS, transformation timestamp, library versions, maximum observed coordinate drift, and validation status. Downstream pipeline stages must read the manifest and abort if validation_passed is false.

CF-convention NetCDF outputs produced by the tidal datum transformation workflow must include crs_wkt and geoid_model global attributes. GeoTIFFs must carry a GDAL_METADATA domain block documenting the exact PROJ string and grid files used. This ensures reproducibility when hydrodynamic models are retrained or regulatory baselines shift.

Vessel telemetry streams — such as those parsed from raw AIS feeds by the AIS NMEA sentence parser — require special handling. Transform positional fixes in batch windows rather than per-message to maintain transformer state consistency and reduce computational overhead. A single Transformer instance is not thread-safe; create one per worker when parallelising across vessels.

After CRS alignment, datasets enter the spatial analysis queue. Surface models proceed to bathymetric DEM interpolation; habitat polygons proceed to overlay analysis; vessel trajectories proceed to behaviour segmentation. The aligned CRS guarantees that all these operations share a coherent spatial reference, eliminating the silent misregistration errors that otherwise corrupt their outputs.

---

Related #

• Fixing CRS Mismatch in QGIS and GeoPandas — step-by-step debug guide for interactive and programmatic CRS diagnostics
• Tidal Datum Transformations in Python — applying MLLW, NAVD88, and LAT vertical references to coastal survey data
• Parsing AIS NMEA Sentences with Python — sanitising positional fixes from vessel telemetry before rasterisation
• Understanding NetCDF vs GeoTIFF for Marine Data — format-specific CRS metadata extraction and storage patterns

Parent: Marine Spatial Data Fundamentals & Architecture