Reprojecting Mixed-CRS Datasets with PyProj

You have inherited a folder of layers from several sources — a county parcel export in a state plane system, a utilities layer in Universal Transverse Mercator (UTM), and a basemap someone saved as WGS 84 degrees — and you need every one of them aligned to a single coordinate reference system (CRS) before any validation rule runs. This page shows how to inventory the declared CRS of each layer, detect where the declared CRS contradicts the actual coordinates, select an accurate datum-shift pipeline with pyproj, and batch-reproject everything to one canonical CRS. It sits within CRS Precision Standards, and the normalization it produces is the input every geometry and topology check depends on.

Prerequisites

• Python 3.10+ — required for the | union type hints used below.
• pyproj 3.6+ (pip install pyproj>=3.6) — bundles PROJ 9.x, which exposes TransformerGroup, per-operation accuracy metadata, and only_best transformation selection. Confirm with import pyproj; print(pyproj.proj_version_str).
• GeoPandas 0.14+ and Shapely 2.0+ — for reading layers and applying to_crs.
• PROJ data grids installed. Datum shifts such as NADCON5 (North American Datum conversion) and NTv2 (National Transformation version 2) depend on grid files. Run pyproj.datadir.get_data_dir() and confirm .tif/.gsb grids are present, or set PROJ_NETWORK=ON to fetch them from the PROJ content delivery network on demand. Without grids, transformations silently degrade to a lower-accuracy fallback.
• A chosen canonical CRS. Pick one projected CRS whose zone spans the whole dataset extent and whose linear unit matches your tolerance matrix. The examples target EPSG:6543 (NAD83(2011) / North Carolina, US survey feet); substitute your own.

Step-by-Step Procedure

Step 1 — Inventory the declared CRS of every layer

Read the declared CRS from each source and record its EPSG code, unit, and datum. This manifest is the ground truth for both mismatch detection and pipeline selection.

from pathlib import Path
import geopandas as gpd
from pyproj import CRS

def inventory_crs(paths: list[str | Path]) -> list[dict]:
    """Record the declared CRS of every layer for downstream reprojection."""
    manifest = []
    for path in paths:
        gdf = gpd.read_file(path, rows=1)  # read one row: metadata only, cheap
        crs = gdf.crs
        entry = {"path": str(path), "declared": None, "epsg": None,
                 "unit": None, "datum": None, "is_geographic": None}
        if crs is not None:
            pyproj_crs = CRS.from_user_input(crs)
            entry.update(
                declared=pyproj_crs.name,
                epsg=pyproj_crs.to_epsg(),
                unit=pyproj_crs.axis_info[0].unit_name,
                datum=pyproj_crs.datum.name if pyproj_crs.datum else None,
                is_geographic=pyproj_crs.is_geographic,
            )
        manifest.append(entry)
    return manifest

Verification: Print the manifest and confirm every layer has a non-null epsg. Any row where declared is None is an undefined-CRS layer — it cannot be reprojected until a CRS is assigned by hand, so route it to manual review rather than guessing.

Step 2 — Detect declared-versus-actual CRS mismatches

A file can declare one CRS while its coordinates were captured in another. The cheapest reliable detector compares coordinate magnitude against the valid axis range of the declared CRS: geographic degrees must fall within ±180/±90, while projected values are typically far larger.

from pyproj import CRS

def detect_crs_mismatch(gdf: gpd.GeoDataFrame) -> dict:
    """Flag layers whose coordinate magnitude contradicts their declared CRS."""
    crs = CRS.from_user_input(gdf.crs)
    minx, miny, maxx, maxy = gdf.total_bounds
    max_abs = max(abs(minx), abs(miny), abs(maxx), abs(maxy))

    suspect = False
    reason = "consistent"
    if crs.is_geographic and max_abs > 360:
        suspect, reason = True, "declared geographic but coordinates are projected-scale"
    if crs.is_projected and max_abs <= 360:
        suspect, reason = True, "declared projected but coordinates look like degrees"

    # Cross-check against the CRS area-of-use bounds when available
    aou = crs.area_of_use
    if aou and crs.is_geographic:
        if not (aou.west - 1 <= minx <= aou.east + 1):
            suspect, reason = True, "longitude outside declared CRS area of use"

    return {"declared_epsg": crs.to_epsg(), "max_abs_coord": float(max_abs),
            "suspect": suspect, "reason": reason}

Verification: Run against a layer you deliberately mislabel — for example a UTM layer with its CRS overwritten to EPSG:4326. The function must return suspect=True. A clean layer returns suspect=False. Never reproject a suspect layer with to_crs, because to_crs trusts the declared CRS; instead correct the declaration with set_crs(allow_override=True) first.

Step 3 — Select an explicit transformation pipeline

Between two datums there are usually several candidate operations with different accuracies. TransformerGroup enumerates them so you can choose deliberately instead of accepting whatever PROJ defaults to.

from pyproj.transformer import TransformerGroup

def choose_transformer(source_epsg: int, target_epsg: int):
    """Return the highest-accuracy available transformation and its metadata."""
    group = TransformerGroup(source_epsg, target_epsg, always_xy=True)
    if not group.transformers:
        raise RuntimeError(f"No transformation from {source_epsg} to {target_epsg}")

    if group.unavailable_operations:
        # A more accurate grid-based operation exists but its grid is missing
        print(f"WARNING: {len(group.unavailable_operations)} higher-accuracy "
              f"operation(s) unavailable — install PROJ grids or set PROJ_NETWORK=ON")

    best = group.transformers[0]  # transformers are sorted best-first
    accuracy = group.operations[0].accuracy  # metres, or None if unknown
    print(f"Selected: {group.operations[0].name} (accuracy ~{accuracy} m)")
    return best, accuracy

Verification: For a NAD27-to-NAD83 pair (EPSG:4267 to EPSG:4269), a correctly provisioned environment reports a NADCON/NTv2 grid operation with sub-metre accuracy. If unavailable_operations is non-empty, the selected transform is the ballpark fallback and every reprojected coordinate carries a one-to-two-metre datum error.

Step 4 — Batch-reproject every layer to the canonical CRS

With the manifest and a target CRS, reproject each layer and record the operation actually used per layer so the run is auditable.

import geopandas as gpd

CANONICAL_EPSG = 6543  # NAD83(2011) / North Carolina (ftUS)

def reproject_all(manifest: list[dict], target_epsg: int = CANONICAL_EPSG
                  ) -> tuple[list[gpd.GeoDataFrame], list[dict]]:
    """Reproject every inventoried layer to one canonical CRS with audit metadata."""
    layers, audit = [], []
    for entry in manifest:
        if entry["epsg"] is None:
            audit.append({"path": entry["path"], "status": "skipped_undefined_crs"})
            continue
        gdf = gpd.read_file(entry["path"])
        if entry["epsg"] == target_epsg:
            out = gdf
            op_name, accuracy = "identity", 0.0
        else:
            _, accuracy = choose_transformer(entry["epsg"], target_epsg)
            out = gdf.to_crs(epsg=target_epsg)
            op_name = f"{entry['epsg']}->{target_epsg}"
        layers.append(out)
        audit.append({"path": entry["path"], "status": "reprojected",
                      "operation": op_name, "accuracy_m": accuracy,
                      "target_epsg": target_epsg})
    return layers, audit

Verification: Assert every output layer shares the target CRS: assert all(l.crs.to_epsg() == CANONICAL_EPSG for l in layers). The audit list should carry an accuracy_m value for each reprojected layer — store it alongside the dataset hash so downstream consumers can trace positional accuracy.

Step 5 — Verify alignment with control points

Coordinate counts and CRS labels can match while the data is still misaligned. Confirm that a known control point lands where it should, within the accuracy budget from Step 3.

from pyproj import Transformer

def verify_control_point(source_epsg: int, target_epsg: int,
                         src_xy: tuple[float, float],
                         expected_xy: tuple[float, float],
                         budget_m: float = 0.5) -> bool:
    """Transform a known control point and confirm it lands within the budget."""
    t = Transformer.from_crs(source_epsg, target_epsg, always_xy=True)
    got_x, got_y = t.transform(*src_xy)
    dx, dy = got_x - expected_xy[0], got_y - expected_xy[1]
    offset = (dx ** 2 + dy ** 2) ** 0.5
    ok = offset <= budget_m
    print(f"Control point offset: {offset:.3f} (budget {budget_m}) -> "
          f"{'PASS' if ok else 'FAIL'}")
    return ok

Verification: Use a published geodetic control monument whose coordinates are known in both systems. A failing offset larger than the budget means the wrong pipeline was chosen or a datum grid is missing — go back to Step 3 before trusting the batch.

Interpreting Results

The audit records from Step 4 tell you not just that reprojection ran, but how trustworthy each layer now is. Map the accuracy field to an action:

A layer that reprojects without error can still be wrong: to_crs faithfully applies whatever transformation PROJ selected, and a missing grid produces a clean-looking but silently displaced result. The accuracy column is the difference between a defensible reprojection and one that will surface as a boundary dispute later. Once every layer shares the canonical CRS and its accuracy is recorded, the dataset is ready for Shapely geometry checks and the precision audits described in the parent standard.

Gotchas & Edge Cases

always_xy matters more than any other flag. PROJ honours each CRS’s declared axis order, so geographic systems default to latitude-longitude, not longitude-latitude. Construct every Transformer and TransformerGroup with always_xy=True to force the x-then-y convention GeoPandas and Shapely expect. Omitting it swaps easting and northing and produces coordinates that look plausible but are transposed.

A missing grid degrades accuracy without raising. When a required NADCON or NTv2 grid is absent, PROJ does not error — it selects a lower-accuracy null transformation and continues. Always inspect TransformerGroup.unavailable_operations; an empty declared-CRS pipeline list is not the same as a correct one.

Reprojecting geographic degrees as the canonical CRS breaks metric tolerances. If you normalize everything to EPSG:4326, area and distance thresholds expressed in metres cannot be applied cleanly, because a degree is not a fixed ground distance. Choose a projected canonical CRS whose zone covers the extent, and confirm decimal-precision expectations against Setting Decimal Precision for Survey Boundaries.

State plane feet and metres are easy to confuse. Many US datasets exist in both US survey feet and international feet variants of the same zone. EPSG:2264 (feet) and EPSG:32119 (metres) cover the same ground with different units; mixing them scales coordinates by roughly 3.28 and passes magnitude checks. Compare the unit_name in the manifest, not just the numeric range.

Anti-meridian and pole-crossing extents defeat UTM. A dataset spanning more than one UTM zone or crossing 180° longitude cannot be reprojected to a single UTM zone without distortion at the edges. Use a wider projected system (a national grid or an equal-area projection) as the canonical CRS for continental or global extents.

When to Escalate

The single-process pyproj and GeoPandas approach here is right for tens of layers and datasets up to a few hundred thousand features. Escalate when:

• Reprojection must run inside the database. For millions of rows already in PostGIS, ST_Transform inside a spatial-indexed table avoids the serialization round-trip. Provision the same PROJ grids on the database host so accuracy matches the Python path, then apply the precision controls from CRS Precision Standards after transformation.
• The layer count outgrows a single process. Hundreds of layers arriving continuously belong in a batched, parallel workflow rather than a serial loop. Partition by source directory and reproject partitions independently.
• Datum epochs and plate motion matter. For high-accuracy geodesy across time (for example NAD83(2011) at a specific epoch versus a current ITRF realization), a static transformation is insufficient; you need time-dependent, epoch-aware pipelines and should involve a geodesist.
• Mismatch detection needs to feed a rule engine. When declared-versus-actual checks must run as versioned, auditable rules rather than an ad hoc script, wire them into a validation harness such as the one in Implementing Shapely Geometry Checks in Python.

Frequently Asked Questions

Related:

• CRS Precision Standards — the tolerance and rounding controls that run once every layer shares one CRS
• Setting Decimal Precision for Survey Boundaries — how many decimals to retain after reprojection for legally defensible parcel vertices
• Implementing Shapely Geometry Checks in Python — the validity pass that consumes a CRS-normalized dataset

Back to CRS Precision Standards