scripts/analysis/check_cavity_proximity.py

"""
Real, independent karst cross-check for every real gauge station, using
the actual raw BDCavités data -- not the precomputed distance_to_
nearest_cavity_km column, which turned out to never have been saved
anywhere in the current reach-graph pipeline (confirmed directly: the
column is absent from both real enriched node files, despite the raw
cavité data itself existing on disk).

Built specifically because the current karst/losing-reach flag
(verified_continuous=False on 3 real edges of La Risle's bétoire
stretch) is a hardcoded list based on real station naming evidence
("[amont bétoire]"/"[aval bétoire]"), not a systematic scan -- it can
only ever confirm what was already suspected from names, never surface
a DIFFERENT real station that also sits near a real cavity but never
got a "bétoire"-style name. This script is that systematic scan.

Usage:
    python -m scripts.analysis.check_cavity_proximity --data-root datasets
"""
import argparse
import json
from pathlib import Path
from typing import List, Optional

import numpy as np
import pandas as pd


def _haversine_km_vec(lat1, lon1, lat2, lon2):
    """Real haversine distance, vectorized -- same formula used
    throughout this project's other real distance computations."""
    R = 6371.0
    lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])
    dlat, dlon = lat2 - lat1, lon2 - lon1
    a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
    return R * 2 * np.arcsin(np.sqrt(a))


def load_real_cavity_coordinates(geojson_path: Path, allowed_types: Optional[List[str]] = None) -> np.ndarray:
    """
    Extracts real [lat, lon] pairs directly from the raw BDCavités
    GeoJSON -- point features only, skipping anything malformed rather
    than crashing the whole load on a single bad real record (large
    real government datasets do have occasional genuine data-quality
    issues; skip-and-count is more useful here than a hard failure).

    allowed_types, when given, filters on the real "type_cavite" field.
    Confirmed as a real, necessary filter, not an optional refinement:
    an unfiltered first run showed 24 of 27 real stations within 1km of
    *some* cavity, with counts in the thousands -- because BDCavités
    catalogs every cavity type, and the real breakdown (confirmed
    directly against this project's actual file) is dominated by
    quarries (carrière, 21,813), civil/military structures (ouvrage
    civil/militaire, ~9,900), and cellars (cave, 3,990) -- none
    hydrologically meaningful. Only "naturelle" (8,313) corresponds to
    genuine natural cavities, the category that could actually indicate
    real karst/water-loss connectivity. "indéterminé" (18,349) is
    genuinely ambiguous -- could hide real unclassified natural
    features or unclassified quarries -- so it's deliberately not
    included by default; check separately with it included rather than
    assume either way.
    """
    with open(geojson_path) as f:
        data = json.load(f)

    coords = []
    n_skipped = 0
    n_filtered = 0
    for feature in data.get("features", []):
        geom = feature.get("geometry")
        if geom is None or geom.get("type") != "Point":
            n_skipped += 1
            continue
        if allowed_types is not None:
            real_type = feature.get("properties", {}).get("type_cavite")
            if real_type not in allowed_types:
                n_filtered += 1
                continue
        lon, lat = geom["coordinates"][0], geom["coordinates"][1]
        coords.append((lat, lon))

    filter_msg = f", {n_filtered} filtered out by type" if allowed_types is not None else ""
    print(f"Loaded {len(coords)} real cavity point(s) from {geojson_path.name} "
          f"({n_skipped} non-point/malformed feature(s) skipped{filter_msg})")
    return np.array(coords)


def check_cavity_proximity(stations_df: pd.DataFrame, cavity_coords: np.ndarray, radius_km: float = 20.0) -> pd.DataFrame:
    """
    Real nearest-cavity distance and real cavity count within radius_km,
    per real station -- same two metrics documented for the (never
    actually saved) distance_to_nearest_cavity_km / n_cavities_within_
    20km features, computed fresh here directly against the real raw
    data.
    """
    cav_lat, cav_lon = cavity_coords[:, 0], cavity_coords[:, 1]
    results = []
    for _, row in stations_df.iterrows():
        dists = _haversine_km_vec(row["latitude"], row["longitude"], cav_lat, cav_lon)
        nearest_km = float(dists.min())
        n_within_radius = int((dists <= radius_km).sum())
        results.append({
            "station_code": row["station_code"],
            "distance_to_nearest_cavity_km": nearest_km,
            "n_cavities_within_20km": n_within_radius,
        })
    return pd.DataFrame(results).sort_values("distance_to_nearest_cavity_km")


def main() -> None:
    parser = argparse.ArgumentParser(description="Real cavity-proximity cross-check for every real gauge station")
    parser.add_argument("--data-root", type=Path, default=Path("datasets"))
    parser.add_argument("--stations-path", type=Path, default=None)
    parser.add_argument("--close-threshold-km", type=float, default=1.0,
                         help="Real distance below which a station is flagged as plausibly karst-connected.")
    args = parser.parse_args()

    stations_path = args.stations_path or (args.data_root / "station_elevations.csv")
    stations_df = pd.read_csv(stations_path)

    geojson_path = args.data_root / "bdcavites" / "cavite_localisee.geojson"
    if not geojson_path.exists():
        print(f"No real cavity file found at {geojson_path} -- run "
              f"scripts.data_acquisition.download_bdcavites first.")
        return

    # Two real, confirmed-necessary filtered runs, not one -- an
    # unfiltered check showed 24/27 real stations within 1km of *some*
    # cavity (counts in the thousands), because BDCavités catalogs
    # every cavity type and quarries/civil structures/cellars dominate
    # the real breakdown. "naturelle" is the conservative, defensible
    # filter; "naturelle + indéterminé" is a broader one that also
    # captures any real natural feature that was never properly typed
    # -- reported side by side rather than picking one and hiding the
    # other's disagreement.
    for label, allowed_types in [
        ("naturelle only (conservative)", ["naturelle"]),
        ("naturelle + indéterminé (broader)", ["naturelle", "indéterminé"]),
    ]:
        print()
        print("=" * 70)
        print(f"Filter: {label}")
        print("=" * 70)
        cavity_coords = load_real_cavity_coordinates(geojson_path, allowed_types=allowed_types)
        if len(cavity_coords) == 0:
            print("Zero real cavities matched this filter -- nothing to check against.")
            continue

        result_df = check_cavity_proximity(stations_df, cavity_coords)
        print(result_df.to_string(index=False))

        close = result_df[result_df["distance_to_nearest_cavity_km"] < args.close_threshold_km]
        print()
        if len(close) > 0:
            print(f"Station(s) within {args.close_threshold_km}km of a real, filtered cavity "
                  f"(cross-check against the known bétoire stations -- H605641101, H605022010, "
                  f"H605641401, H605641201 -- anything else here is a real, new candidate):")
            print(close.to_string(index=False))
        else:
            print(f"No real station within {args.close_threshold_km}km of a real, filtered cavity.")

        output_path = args.data_root / f"cavity_proximity_check_{allowed_types[0]}.csv"
        result_df.to_csv(output_path, index=False)
        print(f"Saved to {output_path}")
    print(f"\nSaved to {output_path}")


if __name__ == "__main__":
    main()