Handling timezone shifts in cross-border mobility data

import numpy as np
import pandas as pd
import geopandas as gpd
from zoneinfo import ZoneInfo
from pathlib import Path
import warnings

warnings.filterwarnings("ignore", category=UserWarning)


def haversine_km(
    lat1: np.ndarray,
    lon1: np.ndarray,
    lat2: np.ndarray,
    lon2: np.ndarray,
) -> np.ndarray:
    """
    Vectorized haversine distance in kilometres.

    NOTE: Uses raw WGS84 degrees — suitable only for velocity anomaly
    detection. For metric-precision distance analytics, project to a
    local UTM zone (e.g. EPSG:32632) before computing distances.
    """
    R = 6371.0
    lat1r, lon1r, lat2r, lon2r = map(np.radians, [lat1, lon1, lat2, lon2])
    dlat = lat2r - lat1r
    dlon = lon2r - lon1r
    a = np.sin(dlat / 2) ** 2 + np.cos(lat1r) * np.cos(lat2r) * np.sin(dlon / 2) ** 2
    return 2.0 * R * np.arctan2(np.sqrt(a), np.sqrt(1.0 - a))


def normalize_cross_border_timestamps(
    raw_df: pd.DataFrame,
    tz_boundaries_path: str | Path,
    velocity_ceiling_kmh: float = 1200.0,
) -> pd.DataFrame:
    """
    Ingest cross-border mobility records, resolve missing IANA timezones
    spatially, normalize to UTC, and flag physically impossible pings.

    Parameters
    ----------
    raw_df : DataFrame with columns:
        device_id   (str)   — unique device or asset identifier
        local_time  (str)   — ISO-8601 or locale timestamp string
        lat         (float) — WGS84 latitude
        lon         (float) — WGS84 longitude
        timezone    (str)   — IANA name (e.g. 'Europe/Paris'); may be NaN
    tz_boundaries_path : path to a GeoPackage or Shapefile whose 'tzid'
        column contains IANA timezone names (e.g. timezone-boundary-builder)
    velocity_ceiling_kmh : consecutive pings implying speed above this
        value are flagged as anomalies (default 1 200 km/h covers all
        surface and subsonic air transport)

    Returns
    -------
    DataFrame with additional columns:
        utc_timestamp  — UTC-normalized pd.Timestamp
        tz_source      — 'provided' | 'spatial' | 'missing'
        is_anomaly     — bool, True when implied velocity exceeds ceiling
    """
    if raw_df.empty:
        return raw_df.copy()

    required = {"device_id", "local_time", "lat", "lon", "timezone"}
    missing_cols = required - set(raw_df.columns)
    if missing_cols:
        raise ValueError(f"Input DataFrame missing columns: {missing_cols}")

    df = raw_df.copy()
    df["local_time"] = pd.to_datetime(df["local_time"], format="ISO8601", utc=False)
    df["tz_source"] = np.where(df["timezone"].notna(), "provided", "missing")

    # ------------------------------------------------------------------ #
    # Step 1: Spatially resolve missing IANA timezones                    #
    # ------------------------------------------------------------------ #
    missing_mask = df["timezone"].isna()
    if missing_mask.any():
        gdf_missing = gpd.GeoDataFrame(
            df[missing_mask].copy(),
            geometry=gpd.points_from_xy(
                df.loc[missing_mask, "lon"],
                df.loc[missing_mask, "lat"],
            ),
            crs="EPSG:4326",
        )
        tz_world = gpd.read_file(tz_boundaries_path)
        joined = gpd.sjoin(
            gdf_missing,
            tz_world[["tzid", "geometry"]],
            how="left",
            predicate="within",
        )
        resolved = joined["tzid"].values
        df.loc[missing_mask, "timezone"] = resolved

        # Fallback: nearest-neighbour for points in polygon gaps (maritime,
        # rural borders). Only triggered when sjoin still returns NaN.
        still_missing = missing_mask & df["timezone"].isna()
        if still_missing.any():
            tz_centroids = tz_world.copy()
            tz_centroids["geometry"] = tz_world.geometry.centroid
            for idx in df.index[still_missing]:
                pt = gpd.GeoSeries(
                    [gpd.points_from_xy([df.at[idx, "lon"]], [df.at[idx, "lat"]])[0]],
                    crs="EPSG:4326",
                ).to_crs("EPSG:3857")
                centroids_proj = tz_centroids.to_crs("EPSG:3857")
                nearest_idx = centroids_proj.geometry.distance(pt.iloc[0]).idxmin()
                df.at[idx, "timezone"] = tz_world.at[nearest_idx, "tzid"]

        df.loc[missing_mask & df["timezone"].notna(), "tz_source"] = "spatial"

    # ------------------------------------------------------------------ #
    # Step 2: Localize to IANA zone → convert to UTC                      #
    # DST rules:                                                           #
    #   ambiguous='NaT'      — fall-back duplicates flagged, not guessed  #
    #   nonexistent='shift_forward' — spring-forward gaps advanced        #
    # ------------------------------------------------------------------ #
    def _to_utc(row: pd.Series) -> pd.Timestamp:
        if pd.isna(row["timezone"]):
            return pd.NaT
        ts = row["local_time"]
        try:
            tz = ZoneInfo(row["timezone"])
        except Exception:
            return pd.NaT
        if ts.tzinfo is None:
            ts = ts.tz_localize(
                tz,
                ambiguous="NaT",
                nonexistent="shift_forward",
            )
        if ts is pd.NaT or ts is None:
            return pd.NaT
        return ts.tz_convert(ZoneInfo("UTC"))

    df["utc_timestamp"] = df.apply(_to_utc, axis=1)

    # Drop rows where DST ambiguity made localization impossible; these
    # should be quarantined and reviewed with the originating data owner.
    pre_drop = len(df)
    df = df.dropna(subset=["utc_timestamp"]).copy()
    if len(df) < pre_drop:
        import logging
        logging.getLogger(__name__).warning(
            "Dropped %d rows with unresolvable DST ambiguity.",
            pre_drop - len(df),
        )

    # ------------------------------------------------------------------ #
    # Step 3: Kinematic anomaly detection                                  #
    # ------------------------------------------------------------------ #
    df = df.sort_values(["device_id", "utc_timestamp"]).reset_index(drop=True)

    grp = df.groupby("device_id")
    prev_lat = grp["lat"].shift(1)
    prev_lon = grp["lon"].shift(1)
    prev_ts  = grp["utc_timestamp"].shift(1)

    dist_km  = haversine_km(
        prev_lat.values, prev_lon.values,
        df["lat"].values, df["lon"].values,
    )
    elapsed_h = (df["utc_timestamp"] - prev_ts).dt.total_seconds() / 3600.0
    velocity  = np.where(elapsed_h > 0.0, dist_km / elapsed_h, np.nan)

    df["is_anomaly"] = velocity > velocity_ceiling_kmh
    return df

import logging

def validate_normalized_mobility(df: pd.DataFrame) -> None:
    """Post-run sanity checks for UTC-normalized mobility data."""
    assert "utc_timestamp" in df.columns, "utc_timestamp column missing"
    assert df["utc_timestamp"].dt.tz is not None, "utc_timestamp must be tz-aware"

    # All timestamps should be UTC
    assert str(df["utc_timestamp"].dt.tz) == "UTC", "Non-UTC timezone found"

    # Timestamps must be monotonically increasing per device after sort
    is_monotone = (
        df.groupby("device_id")["utc_timestamp"]
        .apply(lambda s: s.is_monotonic_increasing)
        .all()
    )
    assert is_monotone, "Non-monotonic timestamps detected per device"

    anomaly_rate = df["is_anomaly"].mean()
    if anomaly_rate > 0.05:
        logging.warning(
            "Anomaly rate %.1f%% exceeds 5%% — check source timezone metadata.",
            anomaly_rate * 100,
        )

    # tz_source distribution provides a quick data-quality signal
    print(df["tz_source"].value_counts().to_string())
    print(f"Rows: {len(df):,} | Anomalies: {df['is_anomaly'].sum():,}")