Normalizing EPSG Codes Before a Spatial Merge

Turn mixed, missing, and ambiguous coordinate definitions across branches into one authority code before any geometry is combined — part of CRS Conflict Resolution Across Branches.

Concept & Context

Before two branches of the same layer can be merged safely, they must speak the same coordinate language. In practice they rarely do. One branch carries a clean EPSG:25832 code in a GeoPackage SRS table, another ships a shapefile whose .prj was lost in transit, and a third holds a raw PROJ string pasted from an old processing script. Each of these is a different way of describing — or failing to describe — a reference frame, and a merge that ignores the difference produces the silent corruption described in the parent guide on CRS conflict resolution across branches.

Normalization is the narrow, mechanical step of collapsing all of those forms into a single, unambiguous authority code per layer. It is deliberately conservative: where a definition cannot be resolved with confidence, normalization quarantines the layer rather than guessing. That conservatism is what makes it safe to run automatically. Once every input reports the same authority code and axis order, the actual reprojection and merge become a routine operation, and the downstream check for detecting silent reprojection errors in pull requests has a stable baseline to compare against.

Core Normalization Pipeline

The pipeline runs five steps, each with a strict contract so that an ambiguous input can never leak into the merged output disguised as a confident one.

  1. Collect the declared CRS. Read the reference system from each layer using its native mechanism — the .prj sidecar for shapefiles, the SRS table for a GeoPackage, the embedded metadata for GeoParquet. A missing definition is recorded as None, never silently defaulted.
  2. Resolve to an authority code. Feed WKT, ESRI WKT, and PROJ-string definitions through pyproj and request the authority code with a minimum confidence. A clean match returns something like EPSG:25832.
  3. Match custom strings to the nearest EPSG. When no exact code is returned, score candidate EPSG systems with pyproj’s identification matching and accept the top match only if it clears a confidence floor.
  4. Normalize axis order. Record and standardize the latitude-longitude versus longitude-latitude convention so that two layers agreeing on a code but disagreeing on axis order are reconciled to one order.
  5. Emit a normalized manifest. Produce one row per layer listing the resolved authority code, the axis order, the confidence, and a status. Anything undefined or below the floor is marked quarantine for human assignment.

Working Implementation

"""Normalize the CRS of every layer in a set of branches to one authority code.

Conservative by design: any layer that cannot be resolved with confidence is
marked for quarantine rather than assigned a guessed code.
"""
import geopandas as gpd
import pyproj
from dataclasses import dataclass, asdict
from typing import Optional

CONFIDENCE_FLOOR = 70          # reject nearest-EPSG matches below this score
ALWAYS_XY = True               # fix on longitude-latitude for the normalized output


@dataclass
class NormalizedCRS:
    path: str
    authority: Optional[str]   # e.g. "EPSG:25832", or None when unresolved
    axis_order: Optional[str]  # "xy" (lon,lat) or "yx" (lat,lon)
    confidence: int            # 0–100; 100 for an exact authority match
    status: str                # "ok" | "quarantine"


def resolve_authority(crs: pyproj.CRS) -> tuple[Optional[str], int]:
    """Return (authority_code, confidence). Tries an exact code first, then a
    confidence-scored nearest-EPSG match for raw PROJ / ESRI definitions."""
    exact = crs.to_authority(min_confidence=100)
    if exact is not None:
        return ":".join(exact), 100

    # No exact code: ask PROJ for scored EPSG candidates (best confidence first).
    candidates = crs.list_authority(auth_name="EPSG", min_confidence=CONFIDENCE_FLOOR)
    if candidates:
        best = candidates[0]
        return f"{best.auth_name}:{best.code}", int(best.confidence)
    return None, 0


def axis_convention(crs: pyproj.CRS) -> str:
    """Classify the stored axis order as 'xy' (lon/east first) or 'yx'."""
    first = crs.axis_info[0].abbrev.lower() if crs.axis_info else ""
    return "yx" if first in {"lat", "y", "n"} else "xy"


def normalize_layer(path: str) -> NormalizedCRS:
    gdf = gpd.read_file(path, rows=1)          # header only — CRS without geometry

    # Step 1: a missing definition is never defaulted.
    if gdf.crs is None:
        return NormalizedCRS(path, None, None, 0, "quarantine")

    crs = pyproj.CRS.from_user_input(gdf.crs)

    # Steps 2 & 3: resolve to an authority code, exact or nearest.
    authority, confidence = resolve_authority(crs)
    if authority is None or confidence < CONFIDENCE_FLOOR:
        return NormalizedCRS(path, None, axis_convention(crs), confidence, "quarantine")

    # Step 4: record the axis convention for the resolved system.
    order = axis_convention(crs)
    return NormalizedCRS(path, authority, order, confidence, "ok")


def normalize_branches(paths: list[str]) -> list[dict]:
    """Produce a normalized manifest across all layers of all branches."""
    manifest = [asdict(normalize_layer(p)) for p in paths]
    quarantined = [m for m in manifest if m["status"] == "quarantine"]
    if quarantined:
        names = ", ".join(m["path"] for m in quarantined)
        print(f"[normalize] {len(quarantined)} layer(s) need manual CRS assignment: {names}")
    return manifest


def apply_normalized_crs(path: str, target_authority: str) -> gpd.GeoDataFrame:
    """Reproject a resolved layer to the canonical authority code. GeoPandas
    stores coordinates x/y (lon/lat) internally, so writing the output in the
    canonical code yields the deterministic axis convention set by ALWAYS_XY."""
    gdf = gpd.read_file(path)
    if gdf.crs is None:
        raise ValueError(f"{path}: undefined CRS — resolve via the manifest first")
    return gdf.to_crs(target_authority)

The function deliberately separates resolving a code (normalize_layer) from applying one (apply_normalized_crs). Resolution is read-only and safe to run on every branch in a pull request; application mutates coordinates and runs only after the manifest is clean. A layer that reaches apply_normalized_crs has already been proven to carry a confident authority code.

Validation & Output Verification

Confirm the manifest before letting any layer proceed to reprojection and merge.

No layer proceeds while quarantined. The gate is a single assertion over the manifest:

manifest = normalize_branches(branch_layer_paths)
quarantined = [m for m in manifest if m["status"] == "quarantine"]
assert not quarantined, f"{len(quarantined)} layer(s) have unresolved CRS; assign explicitly"

Every resolved layer shares one authority code. After normalization, the set of distinct codes across a mergeable group must be exactly one:

codes = {m["authority"] for m in manifest if m["status"] == "ok"}
assert len(codes) == 1, f"Layers resolved to multiple CRS: {sorted(codes)}"

Axis order is consistent. Two layers that resolved to the same code but different axis conventions must be reconciled before merge:

orders = {m["axis_order"] for m in manifest if m["status"] == "ok"}
assert len(orders) == 1, f"Inconsistent axis order across layers: {orders}"

Round-trip sanity check. After apply_normalized_crs, the layer’s bounds should sit inside the numeric envelope the target CRS implies — coordinates in the tens where hundreds of thousands are expected mean the reprojection never happened.

Failure Modes

  • Missing .prj treated as a default — a shapefile arrives without its sidecar and a tool assumes EPSG:4326. Root cause: defaulting an undefined CRS from coordinate ranges. Fix: record None and quarantine; assign the CRS explicitly with set_crs from documented provenance.

  • Nearest-EPSG match forced below the floor — a genuinely custom PROJ string is coerced into a standard code that is almost right. Root cause: accepting a low-confidence match. Fix: enforce CONFIDENCE_FLOOR and quarantine anything beneath it for human review.

  • Axis transposition survives normalization — two layers agree on EPSG:4326 but one is stored latitude-first, so features end up mirrored. Root cause: normalizing the code but not the axis order. Fix: classify axis_order per layer and reconcile to one convention with always_xy on output.

  • Same system flagged as two codes — one branch stores an EPSG code and another the equivalent raw WKT, so a naive comparison sees a conflict. Root cause: comparing definition text instead of resolved authority. Fix: resolve both to an authority code before comparing, and only then diff.

Back to CRS Conflict Resolution Across Branches