Enforcing CRS Contracts in a Generation Pipeline

Coordinate reference system (CRS) metadata that is assumed rather than asserted lets a generation pipeline reproject, mislabel, or strip a layer’s CRS silently, so the geometries are numerically wrong long before any test that looks at their values.

Part of Scoping Rules & Data Contracts: that page defines the contract every synthetic artifact owes its consumers — schema, nullability, units, CRS; this page hardens the single most silently violated clause, the CRS, by asserting it at each boundary where a layer changes hands instead of trusting whatever metadata happens to ride along.

Root Cause: An Assumed CRS Is Not an Enforced CRS

A generation pipeline is a chain of components — a sampler, a reprojection step, a tessellator, a writer — and every boundary between them is a place where CRS metadata can be lost or quietly changed. Three mechanisms cause the drift, and none of them raises an error on its own:

  1. Missing CRS treated as a default. A GeoDataFrame built from raw coordinate arrays has crs = None. The next component that needs a CRS often assumes EPSG:4326 rather than refusing to proceed, so metre-valued projected coordinates get labelled as degrees and every downstream distance is wrong by orders of magnitude.
  2. Implicit reprojection on join or overlay. When two layers with different CRSs meet in an overlay or spatial join, some code paths reproject one to match the other automatically. The operation succeeds, the result looks plausible, and the output CRS is whichever layer happened to be on the left — a decision no one made deliberately.
  3. Axis-order confusion. EPSG:4326 is formally latitude-longitude, but most tooling stores longitude-latitude. A boundary that swaps the two produces coordinates that are still inside a valid range, so no bounds check fires, yet every point is transposed across the diagonal of the globe.

The common thread is that CRS is carried as advisory metadata the pipeline trusts, when it should be a contract the pipeline asserts. The fix is not a better default; defaults are the disease. The fix is to declare the expected CRS at every boundary, assert the incoming layer matches it, and refuse to guess. A cheap, powerful cross-check is the bounding box: geographic coordinates live in [-180, 180] × [-90, 90], while a projected CRS such as EPSG:6933 produces values in the millions of metres. If a layer declares EPSG:4326 but its coordinates exceed 180, the declaration is a lie and the pipeline should stop. This is the CRS-specific instance of the same fail-loud discipline that governs automated fallbacks for missing spatial attributes: a missing value must trigger an explicit, logged decision, never a silent guess.

CRS contract gates asserting the declared reference system at every generation-pipeline boundary Top row, left to right: a generator stage emitting points declared as EPSG 4326; a contract gate; a reproject stage converting to the metric CRS EPSG 6933; a contract gate; a tessellate stage producing polygons; a contract gate; and a write stage producing a GeoPackage. Each diamond-shaped gate asserts that the incoming layer's declared CRS matches the expected CRS and that its bounding box falls in the valid range for that CRS, and passes the layer along a primary-coloured path when both hold. A separate lower path in the accent colour illustrates the failure being prevented: a layer arriving with a null CRS is implicitly assumed to be EPSG 4326, crosses a boundary that has no gate, and reaches a corrupted output box marked with a warning symbol, showing that the absence of a gate is what lets silent drift through. Generator declares EPSG:4326 Reproject → EPSG:6933 (m) Tessellate polygons Write GeoPackage ✓ ✓ ✓ CRS contract gate: assert declared CRS == expected · bbox in range Layer with crs = None implicitly assumed 4326 ⚠ Corrupted output degrees labelled as metres no gate on this boundary → silent drift
A contract gate on every boundary asserts the declared CRS and range before a layer moves on; the lower path shows what an ungated boundary lets through — a null-CRS layer assumed to be geographic and written out as nonsense.

Minimal Reproducer: Watch a Layer Drift Silently

Reproduce the failure in isolation. A layer built from projected metre coordinates but never assigned a CRS gets assumed to be geographic by the next operation, and nothing complains.

python
# requirements.txt — pin majors; let patch releases float
geopandas==0.14.*        # CRS-aware layers; GDAL 3.x backend
shapely==2.*             # GEOS 3.11+ geometry engine
pyproj==3.*              # CRS objects and transforms

python
import geopandas as gpd
from shapely.geometry import Point

# Points are actually in metres (EPSG:6933), but the CRS was never set.
pts = gpd.GeoDataFrame(geometry=[Point(-8_600_000, 5_400_000),
                                 Point(-8_601_000, 5_401_000)])
print(pts.crs)                       # None  <- the root cause, and it is not an error

# A downstream step "helpfully" assumes geographic coordinates:
pts = pts.set_crs("EPSG:4326", allow_override=True)   # WRONG, but silent
print(pts.total_bounds)              # [-8601000  5400000 -8600000  5401000]
# Longitude of -8.6 million "degrees" is absurd, yet no exception was raised.

The bounding box is the tell: values in the millions cannot be degrees. The reproducer confirms both the missing-CRS root cause and the cheap check that catches it.

Fix: Declare, Assert, and Fail Loud at Every Boundary

Centralize the contract in one assertion and call it on the way into and out of every component. The assertion refuses to guess: a None CRS is a hard failure, a mismatched CRS is a hard failure, and coordinates outside the declared CRS’s envelope are a hard failure. This is the discipline the CI/CD integration for spatial data stage relies on to keep drift out of a release.

Step 1 — A single contract assertion. One function encodes the whole clause: exact CRS identity plus a bounding-box sanity check keyed to whether the CRS is geographic or projected.

python
import geopandas as gpd
from pyproj import CRS

class CRSContractError(AssertionError):
    """Raised when a layer violates its declared CRS contract."""

def assert_crs(gdf: gpd.GeoDataFrame, expected: str, boundary: str) -> gpd.GeoDataFrame:
    if gdf.crs is None:
        raise CRSContractError(f"[{boundary}] layer has no CRS; refusing to assume one")
    if not CRS.from_user_input(gdf.crs).equals(CRS.from_user_input(expected)):
        raise CRSContractError(f"[{boundary}] CRS {gdf.crs} != contract {expected}")
    minx, miny, maxx, maxy = gdf.total_bounds
    if CRS.from_user_input(expected).is_geographic:          # degrees
        ok = (-180 <= minx and maxx <= 180 and -90 <= miny and maxy <= 90)
    else:                                                    # projected: reject degree-range values
        ok = max(abs(minx), abs(maxx), abs(miny), abs(maxy)) > 180
    if not ok:
        raise CRSContractError(f"[{boundary}] bbox {gdf.total_bounds} implausible for {expected}")
    return gdf

Step 2 — Reproject explicitly, never implicitly. A boundary that changes CRS must assert the input CRS, call to_crs deliberately, and assert the output. No operation is allowed to reproject as a side effect.

python
def reproject(gdf: gpd.GeoDataFrame, src: str, dst: str) -> gpd.GeoDataFrame:
    assert_crs(gdf, src, boundary="reproject:in")     # input must already match the contract
    out = gdf.to_crs(dst)                              # the ONLY place a CRS may change
    return assert_crs(out, dst, boundary="reproject:out")

Step 3 — Wrap every component boundary. Guard the generator’s output, the tessellation input and output, and the writer. Where two layers meet, assert both share the CRS before the overlay so no path can reproject one silently. The tessellation contract is the same one enforced by the polygon tessellation algorithms stage, which assumes its input is already in a metric CRS.

python
def tessellate_stage(points: gpd.GeoDataFrame, work_crs: str = "EPSG:6933"):
    points = assert_crs(points, work_crs, boundary="tessellate:in")   # must be metres
    cells = build_voronoi(points)                     # your tessellator; inherits input CRS
    cells = cells.set_crs(work_crs, allow_override=False)  # never override a real CRS
    return assert_crs(cells, work_crs, boundary="tessellate:out")

def write_stage(gdf: gpd.GeoDataFrame, path: str, out_crs: str = "EPSG:4326"):
    gdf = reproject(gdf, src="EPSG:6933", dst=out_crs)   # explicit final reprojection
    gdf.to_file(path, driver="GPKG")                     # GeoPackage persists the CRS

Verification Step: Gate Implicit Reprojection in CI

Assert the contract as pytest cases so a regression that drops or swaps a CRS fails the build. The key positive test is that an unlabelled or mislabelled layer is rejected, not silently corrected.

python
import geopandas as gpd, pytest
from shapely.geometry import Point

def test_crs_contract_rejects_missing_and_mismatched():
    # 1. A layer with no CRS is refused, not defaulted.
    bare = gpd.GeoDataFrame(geometry=[Point(0, 0)])
    with pytest.raises(CRSContractError, match="no CRS"):
        assert_crs(bare, "EPSG:4326", boundary="test")

    # 2. Metre coordinates mislabelled as degrees fail the bbox check.
    mislabeled = gpd.GeoDataFrame(geometry=[Point(-8_600_000, 5_400_000)],
                                  crs="EPSG:4326")
    with pytest.raises(CRSContractError, match="implausible"):
        assert_crs(mislabeled, "EPSG:4326", boundary="test")

    # 3. A correct geographic layer passes and round-trips through reprojection.
    good = gpd.GeoDataFrame(geometry=[Point(-8.6, 43.2)], crs="EPSG:4326")
    metres = reproject(good, "EPSG:4326", "EPSG:6933")
    assert metres.crs.to_epsg() == 6933
    assert (reproject(metres, "EPSG:6933", "EPSG:4326")
            .geometry.iloc[0].x == pytest.approx(-8.6, abs=1e-6))

The rejection tests are the load-bearing pair: it is not enough that valid data flows through — the gate must actively stop the null-CRS and mislabelled cases the pipeline would otherwise wave past.

Edge Cases & Gotchas

Axis order at EPSG:4326. The authority definition of EPSG:4326 is latitude-longitude, but GeoJSON, shapely, and most in-memory tooling use longitude-latitude. A layer read through a driver that honours authority axis order arrives transposed, and because a swapped point often still lands inside [-90, 90] for both axes, the bounding-box check can pass. Assert the expected axis order explicitly with pyproj and, for regions with |lat| < 90 < |lon| impossible under a swap, add a domain-specific bounds check tighter than the global envelope.

Antimeridian-spanning bounding box. A layer that straddles ±180° longitude reports a total_bounds of roughly [-180, ..., 180, ...] — a box that appears to wrap the whole planet even though the data occupies a narrow strip across the Pacific. A naive extent check reads this as valid but a downstream clip against it fails. Detect the wrap (min near −180 and max near +180 with an empty middle) and either split the layer at the seam or shift into a region-local projected CRS before the extent assertion.

UTM zone boundaries. Pipelines that pick a UTM zone from a layer’s centroid silently choose different zones for layers whose centroids fall on opposite sides of a 6°-wide boundary, so two datasets that should share a CRS end up in EPSG:32630 and EPSG:32631. The contract must pin one explicit projected CRS for the whole run — a single equal-area CRS such as EPSG:6933, or one named UTM zone — rather than deriving it per-layer from geometry.

Frequently Asked Questions

Why not just default a missing CRS to EPSG:4326 and move on?
Because a missing CRS is a signal that something upstream lost metadata, and defaulting hides that failure while producing wrong numbers when the data was actually projected. A degree of longitude is not a metre, so mislabelling metre coordinates as degrees corrupts every distance, area, and join downstream. The contract refuses to guess and forces the upstream producer to declare its CRS explicitly.
Which CRS should the internal working stages use?
Use a metric CRS for any stage that computes distances, areas, or tessellations, because those operations are meaningless in degrees. A global equal-area CRS such as EPSG:6933 works for worldwide data; a single named UTM zone is tighter for a regional dataset. Reproject to EPSG:4326 only at the final write boundary if consumers expect geographic coordinates, and assert the CRS on both sides of that reprojection.
The bounding-box check passes but coordinates still look wrong. What else can drift?
Axis order. EPSG:4326 is formally latitude-longitude while most tooling uses longitude-latitude, and a transposed point frequently still falls inside the valid global range, so the extent check passes. Assert the expected axis order explicitly through pyproj and add a domain-tight bounds check for your region, where a swap would push a coordinate outside the expected local envelope even though it stays inside the global one.