Best Practices for CRS Transformations in Movement Data

import geopandas as gpd
import pandas as pd
from pyproj import Transformer, CRS
from typing import Optional
import numpy as np


def transform_movement_trajectories(
    gdf: gpd.GeoDataFrame,
    target_epsg: int,
    z_col: Optional[str] = None,
    t_col: Optional[str] = None,
) -> gpd.GeoDataFrame:
    """
    Transform a GeoDataFrame of trajectory points to a metric projected CRS.

    Parameters
    ----------
    gdf : gpd.GeoDataFrame
        Input points. Must have .crs defined and Point geometry.
    target_epsg : int
        EPSG code for the target projected (metric) CRS, e.g. 32632 for UTM zone 32N.
    z_col : str, optional
        Column name holding elevation values. If supplied, Z is transformed
        as a third coordinate dimension (vertical datum shift applied if applicable).
    t_col : str, optional
        Column name holding timestamps. If supplied, rows are re-sorted by this
        column after transformation to guarantee temporal ordering.

    Returns
    -------
    gpd.GeoDataFrame
        Copy of input with updated geometry and CRS. Original index is reset.

    Raises
    ------
    ValueError
        If gdf has no CRS, or if target_epsg resolves to a geographic (non-metric) CRS.
    """
    if gdf is None or len(gdf) == 0:
        # Return empty GeoDataFrame with target CRS already set
        empty = gpd.GeoDataFrame(columns=gdf.columns if gdf is not None else [])
        return empty

    if gdf.crs is None:
        raise ValueError(
            "Input GeoDataFrame has no CRS defined. "
            "Call gdf.set_crs(epsg=<SOURCE_EPSG>) before transformation."
        )

    target_crs = CRS.from_epsg(target_epsg)
    if target_crs.is_geographic:
        raise ValueError(
            f"EPSG:{target_epsg} is a geographic CRS (degrees). "
            "Target must be a projected, metric CRS (e.g. UTM, Lambert, State Plane). "
            "Distance and speed calculations must use metric units — never raw WGS84."
        )

    # always_xy=True: force (lon, lat) → (easting, northing) regardless of CRS
    # axis-order definition. Without this flag, WKT2-defined CRS such as EPSG:4326
    # expect (lat, lon) input and will silently transpose your coordinates.
    transformer = Transformer.from_crs(
        gdf.crs,
        target_crs,
        always_xy=True,
    )

    x_src = gdf.geometry.x.values
    y_src = gdf.geometry.y.values

    # --- 3D transform (elevation-aware) ---
    if z_col is not None and z_col in gdf.columns:
        z_src = gdf[z_col].values.astype(float)
        # Replace NaN elevations with 0.0 so the transform call does not fail;
        # NaN positions are flagged in validation rather than dropped here.
        z_src_clean = np.where(np.isnan(z_src), 0.0, z_src)
        try:
            x_new, y_new, z_new = transformer.transform(x_src, y_src, z_src_clean)
            gdf = gdf.copy()
            # Restore NaN mask on elevation output
            gdf[z_col] = np.where(np.isnan(z_src), np.nan, z_new)
        except Exception:
            # Vertical datum transform failed (mismatched vertical CRS pair).
            # Fall back to 2D transform; pass Z through unchanged with a warning.
            import warnings
            warnings.warn(
                f"3D transform failed for EPSG:{target_epsg}. "
                "Falling back to 2D transform; elevation values are NOT reprojected.",
                stacklevel=2,
            )
            x_new, y_new = transformer.transform(x_src, y_src)
            gdf = gdf.copy()
    else:
        x_new, y_new = transformer.transform(x_src, y_src)
        gdf = gdf.copy()

    # Rebuild Point geometry from transformed coordinates
    gdf["geometry"] = gpd.points_from_xy(x_new, y_new)
    gdf = gdf.set_crs(target_crs, allow_override=True)

    # Re-sort by timestamp if supplied: vectorized transforms do not guarantee
    # row order when input is not already sorted by time.
    if t_col is not None and t_col in gdf.columns:
        gdf = gdf.sort_values(t_col).reset_index(drop=True)

    return gdf

from pyproj import CRS


def utm_epsg_for_longitude(lon_deg: float) -> int:
    """Return the WGS84 UTM zone EPSG code for a given longitude."""
    zone = int((lon_deg + 180) / 6) + 1
    # Northern hemisphere: 32600 + zone; southern: 32700 + zone
    return 32600 + zone  # adjust to 32700 + zone for latitudes < 0


# Usage: compute centroid of each trajectory batch, then transform that batch only
centroid_lon = gdf.geometry.x.mean()
target_epsg = utm_epsg_for_longitude(centroid_lon)
transformed = transform_movement_trajectories(gdf, target_epsg=target_epsg, t_col="timestamp")

import numpy as np


def validate_transformed_trajectories(
    gdf: gpd.GeoDataFrame,
    expected_crs_epsg: int,
    max_speed_ms: float = 83.0,  # ~300 km/h; adjust per transport mode
    t_col: str = "timestamp",
) -> dict:
    """
    Post-transformation sanity checks. Returns a dict of check results.
    Raises AssertionError if the CRS check fails (hard requirement).
    """
    results = {}

    # 1. CRS matches target
    assert gdf.crs is not None and gdf.crs.to_epsg() == expected_crs_epsg, (
        f"CRS mismatch: expected EPSG:{expected_crs_epsg}, "
        f"got {gdf.crs.to_epsg() if gdf.crs else 'None'}"
    )
    results["crs_ok"] = True

    # 2. No NaN coordinates (transformation failed silently for some rows)
    nan_x = gdf.geometry.x.isna().sum()
    nan_y = gdf.geometry.y.isna().sum()
    results["nan_coords"] = int(nan_x + nan_y)
    if results["nan_coords"] > 0:
        import warnings
        warnings.warn(f"{results['nan_coords']} NaN coordinate(s) after transformation.")

    # 3. Coordinates are finite (Inf values indicate out-of-domain projection)
    inf_count = (~np.isfinite(gdf.geometry.x.values)).sum()
    inf_count += (~np.isfinite(gdf.geometry.y.values)).sum()
    results["inf_coords"] = int(inf_count)

    # 4. Velocity bounds check (requires at least 2 points and a timestamp column)
    if t_col in gdf.columns and len(gdf) >= 2:
        gdf_sorted = gdf.sort_values(t_col)
        dx = np.diff(gdf_sorted.geometry.x.values)
        dy = np.diff(gdf_sorted.geometry.y.values)
        dt = np.diff(gdf_sorted[t_col].astype("int64")) / 1e9  # ns → seconds
        dt_safe = np.where(dt == 0, np.nan, dt)
        speeds = np.sqrt(dx**2 + dy**2) / dt_safe
        overspeed = int(np.nansum(speeds > max_speed_ms))
        results["overspeed_segments"] = overspeed
    else:
        results["overspeed_segments"] = None

    return results


# Example usage
checks = validate_transformed_trajectories(transformed_gdf, expected_crs_epsg=32632)
print(checks)
# {'crs_ok': True, 'nan_coords': 0, 'inf_coords': 0, 'overspeed_segments': 0}