CRS Normalization Across Mixed Datasets

pip install "geopandas>=0.14.0" "rasterio>=1.3.0" "pyproj>=3.6.0" "shapely>=2.0.0" "numpy>=1.24.0"

import pyproj
from pyproj import CRS

print(f"PROJ Data Dir: {pyproj.datadir.get_data_dir()}")
print(f"Network Enabled: {pyproj.network.is_network_enabled()}")

if not pyproj.network.is_network_enabled():
    pyproj.network.set_network_enabled(True)
    print("PROJ network enabled — datum grids will be downloaded on demand")

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


def extract_crs_metadata(source_path: str | Path) -> dict:
    """
    Extract CRS metadata from vector or raster source without loading full data.
    Returns a dict with 'crs', 'epsg', 'is_geographic', 'units', 'source'.
    """
    path = Path(source_path)
    result: dict = {"source": str(path), "crs": None, "epsg": None,
                    "is_geographic": None, "units": None}

    try:
        # Raster path: GeoTIFF, IMG, VRT, etc.
        if path.suffix.lower() in {".tif", ".tiff", ".img", ".vrt", ".nc"}:
            with rasterio.open(path) as src:
                if src.crs is None:
                    return result  # CRS absent — caller must decide
                crs = CRS.from_wkt(src.crs.to_wkt())
        else:
            # Vector path: read one row to avoid loading full dataset
            gdf = gpd.read_file(path, rows=1)
            crs = gdf.crs

        if crs is None:
            return result

        result["crs"] = crs
        result["epsg"] = crs.to_epsg()                 # None if no EPSG code
        result["is_geographic"] = crs.is_geographic
        result["units"] = crs.axis_info[0].unit_name if crs.axis_info else "unknown"

    except Exception as exc:
        result["error"] = str(exc)

    return result


def audit_dataset_crs(paths: list[str | Path]) -> list[dict]:
    """Inventory CRS across a list of sources and flag problems."""
    inventory = [extract_crs_metadata(p) for p in paths]
    for item in inventory:
        item["needs_normalization"] = (
            item["crs"] is None
            or item["epsg"] is None
            or item["is_geographic"]  # geographic CRS unsuitable for metric ops
        )
    return inventory

def validate_target_crs(epsg_code: int) -> CRS:
    """
    Load and validate a projected CRS suitable for metric spatial operations.
    Raises ValueError with a diagnostic message on failure.
    """
    crs = CRS.from_epsg(epsg_code)

    if not crs.is_projected:
        raise ValueError(
            f"EPSG:{epsg_code} is a geographic CRS ({crs.name}). "
            "Metric operations require a projected CRS."
        )

    units = crs.axis_info[0].unit_name
    if units != "metre":
        raise ValueError(
            f"EPSG:{epsg_code} uses '{units}' as linear unit. "
            "Pipeline requires metric (metre) units."
        )

    return crs


# UTM Zone 33N — suitable for Central Europe
target_crs = validate_target_crs(32633)

import shapely
import numpy as np


def normalize_vector(
    gdf: gpd.GeoDataFrame,
    target_crs: CRS,
    *,
    repair_invalid: bool = True,
) -> tuple[gpd.GeoDataFrame, dict]:
    """
    Reproject a GeoDataFrame to target_crs and optionally repair invalid geometries.

    Returns the reprojected GeoDataFrame and a diagnostics dict.
    """
    if gdf.crs is None:
        raise ValueError(
            "Source CRS is undefined. Assign the correct CRS before reprojecting."
        )

    # Shapely 2.0 vectorized validity check on the array of geometries
    pre_invalid = int(np.sum(~shapely.is_valid(gdf.geometry.values)))

    normalized = gdf.to_crs(target_crs)

    post_invalid_mask = ~shapely.is_valid(normalized.geometry.values)
    post_invalid = int(np.sum(post_invalid_mask))

    if repair_invalid and post_invalid > 0:
        # make_valid operates on the full geometry array — no .apply() loop
        repaired_geoms = shapely.make_valid(normalized.geometry.values)
        normalized = normalized.copy()
        normalized.geometry = gpd.array.GeometryArray(repaired_geoms)
        post_invalid_after_repair = int(
            np.sum(~shapely.is_valid(normalized.geometry.values))
        )
    else:
        post_invalid_after_repair = post_invalid

    diagnostics = {
        "source_crs": gdf.crs.to_epsg(),
        "target_crs": target_crs.to_epsg(),
        "pre_transform_invalid": pre_invalid,
        "post_transform_invalid": post_invalid,
        "post_repair_invalid": post_invalid_after_repair,
        "features": len(normalized),
    }
    return normalized, diagnostics

import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from pathlib import Path


def normalize_raster(
    src_path: str | Path,
    dst_path: str | Path,
    target_crs: CRS,
    resampling: Resampling = Resampling.bilinear,
) -> dict:
    """
    Reproject a raster to target_crs using the specified resampling algorithm.
    Returns a diagnostics dict with source/destination metadata.
    """
    with rasterio.open(src_path) as src:
        if src.crs is None:
            raise ValueError(f"Raster {src_path} has no CRS — cannot reproject.")

        transform, width, height = calculate_default_transform(
            src.crs, target_crs, src.width, src.height, *src.bounds
        )

        dst_meta = src.meta.copy()
        dst_meta.update({
            "crs": target_crs.to_wkt(),
            "transform": transform,
            "width": width,
            "height": height,
            "nodata": src.nodata if src.nodata is not None else -9999,
        })

        with rasterio.open(dst_path, "w", **dst_meta) 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=transform,
                    dst_crs=target_crs,
                    resampling=resampling,
                )

        return {
            "source_crs": CRS.from_wkt(src.crs.to_wkt()).to_epsg(),
            "target_crs": target_crs.to_epsg(),
            "source_shape": (src.height, src.width),
            "destination_shape": (height, width),
            "bands": src.count,
            "resampling": resampling.name,
        }

import geopandas as gpd
from pyproj import CRS
from shapely.geometry import box


def validate_normalized_vector(
    original: gpd.GeoDataFrame,
    normalized: gpd.GeoDataFrame,
    expected_bounds: tuple[float, float, float, float],
    *,
    area_tolerance_pct: float = 2.0,
) -> dict:
    """
    Run quality gates on a reprojected GeoDataFrame.

    Checks:
    - Total bounds fall within expected regional extent
    - Polygon area is preserved within tolerance (catches datum-shift errors)
    - No invalid geometries remain
    """
    import numpy as np

    issues: list[str] = []

    # Bounds gate: transformed data must fall within expected region
    minx, miny, maxx, maxy = normalized.total_bounds
    exp_minx, exp_miny, exp_maxx, exp_maxy = expected_bounds
    if not (exp_minx <= minx and maxx <= exp_maxx and
            exp_miny <= miny and maxy <= exp_maxy):
        issues.append(
            f"Bounds {(minx, miny, maxx, maxy)} outside expected {expected_bounds}"
        )

    # Area gate: compare polygon areas before and after (projected units only)
    if normalized.crs and normalized.crs.is_projected:
        orig_proj = original.to_crs(normalized.crs)
        orig_area = orig_proj.geometry.area.sum()
        norm_area = normalized.geometry.area.sum()
        if orig_area > 0:
            pct_diff = abs(norm_area - orig_area) / orig_area * 100
            if pct_diff > area_tolerance_pct:
                issues.append(
                    f"Area changed by {pct_diff:.2f}% (threshold {area_tolerance_pct}%)"
                )

    # Validity gate
    invalid_count = int(np.sum(~shapely.is_valid(normalized.geometry.values)))
    if invalid_count:
        issues.append(f"{invalid_count} geometries remain invalid after normalization")

    return {"passed": len(issues) == 0, "issues": issues, "feature_count": len(normalized)}

def infer_crs_from_coordinates(gdf: gpd.GeoDataFrame) -> CRS | None:
    """
    Heuristically infer CRS from coordinate ranges.
    Returns a candidate CRS or None if inference is ambiguous.
    Caller must confirm before assigning.
    """
    bounds = gdf.total_bounds  # (minx, miny, maxx, maxy)
    minx, miny, maxx, maxy = bounds

    # Geographic CRS: coordinates in [-180,180] x [-90,90]
    if -180 <= minx and maxx <= 180 and -90 <= miny and maxy <= 90:
        return CRS.from_epsg(4326)

    # Web Mercator: coordinates in metres up to ~20 million
    if abs(minx) < 20_100_000 and abs(miny) < 20_100_000:
        return CRS.from_epsg(3857)

    # Large positive coordinates suggest a national grid or UTM
    # Cannot reliably distinguish — return None for manual resolution
    return None

from pyproj import Transformer

def check_datum_shift_accuracy(
    source_epsg: int,
    target_epsg: int,
    sample_lon: float,
    sample_lat: float,
) -> dict:
    """
    Transform a sample point and report the accuracy flag from pyproj.
    PROJ reports accuracy=inf when grid-based shift is unavailable.
    """
    transformer = Transformer.from_crs(
        source_epsg, target_epsg, always_xy=True, authority="EPSG"
    )
    x, y = transformer.transform(sample_lon, sample_lat)

    # pyproj 3.4+ exposes best_available flag on the transformer
    has_grids = getattr(transformer, "has_ballpark_transformation", False) is False

    return {
        "source": source_epsg,
        "target": target_epsg,
        "sample_in": (sample_lon, sample_lat),
        "sample_out": (x, y),
        "grid_based_shift_available": has_grids,
        "warning": None if has_grids else (
            "Datum shift uses approximate Helmert parameters. "
            "Enable PROJ network access to download precise grid files."
        ),
    }

def reproject_in_chunks(
    src_path: str | Path,
    target_crs: CRS,
    chunk_size: int = 50_000,
) -> gpd.GeoDataFrame:
    """
    Reproject a large vector file in row chunks to control peak memory usage.
    Returns a single reprojected GeoDataFrame via concat.
    """
    import pyogrio

    total = pyogrio.read_info(src_path)["features"]
    chunks: list[gpd.GeoDataFrame] = []

    for offset in range(0, total, chunk_size):
        chunk = gpd.read_file(
            src_path,
            rows=slice(offset, min(offset + chunk_size, total)),
        )
        chunks.append(chunk.to_crs(target_crs))

    result = gpd.GeoDataFrame(
        gpd.pd.concat(chunks, ignore_index=True), crs=target_crs
    )
    return result

import shapely
import numpy as np

# Batch validity check — operates on C-level geometry array
geom_array = gdf.geometry.values  # numpy array of geometry objects
valid_flags: np.ndarray = shapely.is_valid(geom_array)

# Batch repair — no Python loop, no .apply()
repaired: np.ndarray = shapely.make_valid(geom_array)

# Coordinate transformation without GeoPandas overhead
# Useful for coordinate-only transforms (e.g. rounding, axis swaps)
def swap_xy(geom_arr: np.ndarray) -> np.ndarray:
    """Swap X/Y coordinate order on an array of geometries."""
    return shapely.transform(geom_arr, lambda coords: coords[:, ::-1])

from typing import Protocol
import logging

logger = logging.getLogger(__name__)


class CRSNormalizationHook:
    """Drop-in ETL hook that enforces a unified CRS before downstream processing."""

    def __init__(self, target_epsg: int):
        self.target_crs = validate_target_crs(target_epsg)

    def __call__(self, gdf: gpd.GeoDataFrame, source_id: str) -> gpd.GeoDataFrame:
        if gdf.crs is None:
            inferred = infer_crs_from_coordinates(gdf)
            if inferred is None:
                raise ValueError(
                    f"[{source_id}] CRS is None and cannot be inferred. "
                    "Assign the correct EPSG code before ingestion."
                )
            logger.warning("[%s] CRS inferred as EPSG:%s — verify before production use.",
                           source_id, inferred.to_epsg())
            gdf = gdf.set_crs(inferred)

        if not gdf.crs.equals(self.target_crs):
            normalized, diag = normalize_vector(gdf, self.target_crs)
            logger.info("[%s] Reprojected EPSG:%s → EPSG:%s | %s",
                        source_id, diag["source_crs"], diag["target_crs"], diag)
            return normalized

        return gdf

from pathlib import Path
import json

def normalize_with_dead_letter(
    gdf: gpd.GeoDataFrame,
    target_crs: CRS,
    source_id: str,
    dead_letter_dir: Path,
) -> gpd.GeoDataFrame | None:
    """
    Attempt CRS normalization; write failed datasets to dead-letter storage.
    Returns None on failure so the caller can skip rather than crash.
    """
    try:
        normalized, diag = normalize_vector(gdf, target_crs)
        if diag["post_repair_invalid"] > 0:
            raise ValueError(
                f"{diag['post_repair_invalid']} geometries remain invalid after repair."
            )
        return normalized
    except Exception as exc:
        dead_letter_dir.mkdir(parents=True, exist_ok=True)
        # Persist the failed dataset for manual review
        gdf.to_file(dead_letter_dir / f"{source_id}.gpkg", driver="GPKG")
        (dead_letter_dir / f"{source_id}.error.json").write_text(
            json.dumps({"source_id": source_id, "error": str(exc)}, indent=2)
        )
        logger.error("[%s] CRS normalization failed — routed to dead-letter: %s", source_id, exc)
        return None

# In your CI test suite (pytest)
# tests/test_crs_normalization.py
import geopandas as gpd
from pyproj import CRS
from pathlib import Path

def test_all_fixtures_have_declared_crs():
    fixture_dir = Path("tests/fixtures")
    for path in fixture_dir.glob("**/*.gpkg"):
        gdf = gpd.read_file(path, rows=1)
        assert gdf.crs is not None, f"{path.name} has no CRS declaration"
        assert gdf.crs.to_epsg() is not None, f"{path.name} CRS has no EPSG code"