Converting mixed EPSG codes to a unified CRS in Python

import pyproj, geopandas, shapely
print(pyproj.__version__, geopandas.__version__, shapely.__version__)
# Confirm PROJ data directory is populated
print(pyproj.datadir.get_data_dir())

import geopandas as gpd
import pyproj
from pyproj.exceptions import CRSError
import shapely
import logging
import numpy as np
from pathlib import Path
from typing import Iterable, Union
import pandas as pd

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

# Datums known to require grid-based corrections when converting to WGS 84 / ETRS89
_DATUM_SHIFT_DATUMS = frozenset(["NAD27", "OSGB 1936", "Tokyo", "ED50", "Pulkovo 1942"])


def _detect_datum_shift(source_crs: pyproj.CRS, target_crs: pyproj.CRS) -> bool:
    """Return True if the source→target transition crosses a datum that needs grid corrections."""
    try:
        src_datum = source_crs.datum.name if source_crs.datum else ""
        tgt_datum = target_crs.datum.name if target_crs.datum else ""
        return src_datum in _DATUM_SHIFT_DATUMS and src_datum != tgt_datum
    except Exception:
        return False


def unify_mixed_epsg(
    datasets: Iterable[Union[gpd.GeoDataFrame, Path, str]],
    target_epsg: int = 4326,
    strict: bool = True,
    enable_network: bool = True,
) -> gpd.GeoDataFrame:
    """
    Validate and reproject datasets carrying mixed EPSG codes to a single target CRS.

    Parameters
    ----------
    datasets:
        Iterable of file paths (str/Path) or already-loaded GeoDataFrame objects.
    target_epsg:
        EPSG code for the unified output CRS. Default 4326 (WGS 84 geographic).
    strict:
        If True, raise CRSError on undefined or unresolvable CRS. If False, log and
        skip the offending dataset so the pipeline continues.
    enable_network:
        Allow PROJ to download datum-shift grid files on demand. Recommended when
        inputs may carry NAD27, OSGB36, or other legacy datums.

    Returns
    -------
    gpd.GeoDataFrame
        Concatenated frame in target_epsg with a ``_source_epsg`` audit column.
    """
    if enable_network:
        pyproj.network.set_network_enabled(True)
        logger.info("PROJ network access enabled for datum-grid downloads.")

    target_crs = pyproj.CRS.from_epsg(target_epsg)
    logger.info("Target CRS: %s (EPSG:%s)", target_crs.name, target_epsg)

    converted_frames: list[gpd.GeoDataFrame] = []

    for idx, ds in enumerate(datasets):
        # Load from disk or copy an existing frame to avoid mutating caller's data
        gdf: gpd.GeoDataFrame = (
            gpd.read_file(ds) if isinstance(ds, (Path, str)) else ds.copy()
        )

        if gdf.empty:
            logger.warning("Dataset %s is empty — skipping.", idx)
            continue

        # ── 1. Check CRS is defined ──────────────────────────────────────────
        if gdf.crs is None:
            msg = f"Dataset {idx} has undefined CRS."
            if strict:
                raise CRSError(msg + " Cannot proceed in strict mode.")
            logger.warning("%s Assigning target CRS as fallback.", msg)
            gdf = gdf.set_crs(target_crs)
            gdf["_source_epsg"] = None
            converted_frames.append(gdf)
            continue

        # ── 2. Parse and validate against the EPSG registry ─────────────────
        try:
            source_crs = pyproj.CRS.from_user_input(gdf.crs)
            source_epsg = source_crs.to_epsg()
            if source_epsg is None:
                logger.warning(
                    "Dataset %s has no EPSG code (custom/compound CRS). WKT: %s",
                    idx,
                    source_crs.to_wkt()[:120],
                )
            else:
                logger.info("Dataset %s validated: EPSG:%s", idx, source_epsg)
        except CRSError as exc:
            if strict:
                raise CRSError(f"Dataset {idx} has invalid CRS: {exc}") from exc
            logger.error("Dataset %s CRS parse failed: %s — skipping.", idx, exc)
            continue

        # ── 3. Detect datum shifts requiring grid corrections ────────────────
        if _detect_datum_shift(source_crs, target_crs):
            logger.warning(
                "Dataset %s crosses a datum shift (%s → %s). "
                "Ensure PROJ grid files are available for sub-metre accuracy.",
                idx,
                source_crs.datum.name,
                target_crs.datum.name,
            )

        # Skip transform if already in target CRS
        if source_crs.equals(target_crs):
            logger.info("Dataset %s already in target CRS — no transform needed.", idx)
            gdf["_source_epsg"] = source_epsg
            converted_frames.append(gdf)
            continue

        # ── 4. Reproject ─────────────────────────────────────────────────────
        try:
            reprojected = gdf.to_crs(target_crs)
        except Exception as exc:
            if strict:
                raise
            logger.error("Reprojection failed for dataset %s: %s — skipping.", idx, exc)
            continue

        # ── 5. Validate output geometry (Shapely 2 vectorised) ───────────────
        validity_arr = shapely.is_valid(reprojected.geometry.values)
        invalid_count = int(np.sum(~validity_arr))
        if invalid_count:
            logger.warning(
                "Dataset %s: %s geometries invalid after reprojection. "
                "Run geometry repair before downstream joins.",
                idx,
                invalid_count,
            )

        reprojected["_source_epsg"] = source_epsg
        converted_frames.append(reprojected)

    if not converted_frames:
        logger.warning("No valid datasets to concatenate. Returning empty GeoDataFrame.")
        return gpd.GeoDataFrame(geometry=gpd.GeoSeries([], crs=target_crs))

    unified = gpd.GeoDataFrame(
        pd.concat(converted_frames, ignore_index=True),
        crs=target_crs,
    )
    logger.info(
        "Unified GeoDataFrame: %s rows, CRS=%s", len(unified), target_crs.name
    )
    return unified