app/registers/nycha.py

"""nycha_development_exposure — flood-exposure briefing per NYCHA development.

Same pattern as the MTA-entrance specialist, but NYCHA developments are
*polygons* not points, so the metrics shift to overlap fractions:

  - % of footprint inside the 2012 Sandy Inundation Zone (empirical)
  - % of footprint inside DEP Extreme-2080 / Moderate-2050 scenarios
    (modeled, broken out by depth class)
  - Representative-point elevation, HAND, TWI (proxy)
  - Footprint area (km²)
  - Distance from query point to development boundary

Joins:
  - data/nycha.geojson  (NYC Open Data, 218 NYCHA developments)
  - data/sandy_inundation.geojson
  - DEP Stormwater Flood Map polygons (3 scenarios)
  - data/nyc_dem_30m.tif, data/hand.tif

Per queried (lat, lon), returns developments whose centroid is within
the radius (default 2000 m — NYCHA developments are sparser than
subway entrances, so the radius is wider).

Honest scope:
  - This is exposure, not damage forecast. We say "85% of this
    development's footprint is inside the 2012 Sandy zone" — not
    "this development will flood next storm".
  - All overlap fractions are computed in EPSG:2263 (NYC State Plane,
    feet) for accurate area arithmetic in the city.
"""

from __future__ import annotations

import json
import logging
import math
import sys
from dataclasses import dataclass
from functools import lru_cache
from pathlib import Path

_ROOT = Path(__file__).resolve().parents[2]
if str(_ROOT) not in sys.path:
    sys.path.insert(0, str(_ROOT))

log = logging.getLogger("riprap.nycha")

DATA = _ROOT / "data"
NYCHA = DATA / "nycha.geojson"

DEFAULT_RADIUS_M = 2000
DEFAULT_MAX_PER_QUERY = 5


@dataclass
class DevelopmentFinding:
    development: str
    tds_num: str
    borough: str
    centroid_lat: float
    centroid_lon: float
    distance_m: float
    rep_elevation_m: float | None
    rep_hand_m: float | None
    inside_sandy_2012: bool
    dep_extreme_2080_class: int          # 0=outside, 1/2/3 = depth class
    dep_extreme_2080_label: str
    dep_moderate_2050_class: int
    dep_moderate_2050_label: str
    dep_moderate_current_class: int
    dep_moderate_current_label: str


@lru_cache(maxsize=1)
def _load_nycha():
    import geopandas as gpd
    gdf = gpd.read_file(NYCHA).to_crs("EPSG:2263")  # feet, accurate areas
    gdf["centroid_2263"] = gdf.geometry.centroid
    return gdf.reset_index(drop=True)


@lru_cache(maxsize=1)
def _load_sandy_2263():
    """Load the Sandy zone in EPSG:2263 once. Already used by
    app.flood_layers.sandy_inundation but we want the geometry directly
    for overlap-fraction math."""
    import geopandas as gpd
    g = gpd.read_file(DATA / "sandy_inundation.geojson").to_crs("EPSG:2263")
    # Some NYC OEM Sandy polygons have hole-orientation issues that
    # blow up unary_union. buffer(0) fixes self-intersections without
    # changing the footprint at sub-foot precision.
    g["geometry"] = g.geometry.buffer(0)
    return g.geometry.union_all()


@lru_cache(maxsize=4)
def _load_dep_2263(scenario: str):
    """DEP scenario polygons in EPSG:2263, with depth_class column."""
    import geopandas as gpd
    p = DATA / "dep" / f"{scenario}.geojson"
    if not p.exists():
        # Fallback to whatever the existing dep_stormwater module loaded.
        from app.flood_layers import dep_stormwater
        gdf = dep_stormwater.load(scenario)
        return gdf.to_crs("EPSG:2263") if gdf.crs is not None else gdf
    return gpd.read_file(p).to_crs("EPSG:2263")


def _haversine_m(lat1, lon1, lat2, lon2) -> float:
    R = 6371000.0
    p1, p2 = math.radians(lat1), math.radians(lat2)
    dp = math.radians(lat2 - lat1); dl = math.radians(lon2 - lon1)
    a = math.sin(dp / 2) ** 2 + math.cos(p1) * math.cos(p2) * math.sin(dl / 2) ** 2
    return 2 * R * math.asin(math.sqrt(a))


def _sample_raster(raster_path: Path, lat: float, lon: float) -> float | None:
    if not raster_path.exists():
        return None
    try:
        import rasterio
        with rasterio.open(raster_path) as src:
            v = next(src.sample([(lon, lat)]))[0]
            v = float(v)
            if math.isnan(v) or v == src.nodata:
                return None
            return v
    except Exception:
        log.exception("raster sample failed for %s", raster_path)
        return None


def _developments_near(lat: float, lon: float, radius_m: float):
    """Return developments whose centroid is within `radius_m` of
    (lat, lon). Uses haversine on centroids re-projected back to
    EPSG:4326 — the bbox prefilter gets us close, then exact distance."""
    import geopandas as gpd
    gdf = _load_nycha()
    # Re-project centroids to 4326 for haversine
    cents_4326 = gpd.GeoSeries(gdf["centroid_2263"], crs="EPSG:2263").to_crs("EPSG:4326")
    deg = radius_m / 90_000
    cent_lat = cents_4326.y
    cent_lon = cents_4326.x
    mask = ((cent_lat >= lat - deg) & (cent_lat <= lat + deg)
            & (cent_lon >= lon - deg) & (cent_lon <= lon + deg))
    sub = gdf[mask].copy()
    if sub.empty:
        return sub, []
    sub["clat"] = cent_lat[mask].values
    sub["clon"] = cent_lon[mask].values
    sub["distance_m"] = sub.apply(
        lambda r: _haversine_m(lat, lon, r["clat"], r["clon"]),
        axis=1,
    )
    sub = sub[sub["distance_m"] <= radius_m].sort_values("distance_m")
    return sub, sub.index.tolist()


def _overlap_pct(geom_2263, mask_geom_2263) -> float:
    """% of geom_2263's area that intersects mask_geom_2263."""
    if mask_geom_2263 is None or mask_geom_2263.is_empty:
        return 0.0
    inter = geom_2263.intersection(mask_geom_2263)
    if inter.is_empty:
        return 0.0
    return round(100.0 * inter.area / max(geom_2263.area, 1e-9), 2)


def _dep_overlap(geom_2263, scenario: str) -> tuple[float, float]:
    """Return (pct_any_depth, pct_deep_contiguous) of a polygon's area
    inside the DEP scenario."""
    try:
        gdf = _load_dep_2263(scenario)
    except Exception:
        log.exception("DEP load failed for %s", scenario)
        return 0.0, 0.0
    if gdf is None or gdf.empty:
        return 0.0, 0.0
    # Bbox-prefilter the DEP polygons to those near our development.
    minx, miny, maxx, maxy = geom_2263.bounds
    cand = gdf.cx[minx:maxx, miny:maxy]
    if cand.empty:
        return 0.0, 0.0
    # DEP NYC stormwater FGDB uses `Flooding_Category` (int16):
    # 1=nuisance, 2=shallow, 3=deep contiguous (>4 ft).
    cat_col = "Flooding_Category" if "Flooding_Category" in cand.columns else None
    any_geom = cand.geometry.buffer(0).union_all()
    if cat_col:
        deep = cand[cand[cat_col] == 3]
        deep_geom = deep.geometry.buffer(0).union_all() if not deep.empty else None
    else:
        deep_geom = None
    pct_any = _overlap_pct(geom_2263, any_geom)
    pct_deep = _overlap_pct(geom_2263, deep_geom) if deep_geom is not None else 0.0
    return pct_any, pct_deep


_DEPTH_LABEL = {
    0: "outside",
    1: "Nuisance (>4 in to 1 ft)",
    2: "Deep & Contiguous (1-4 ft)",
    3: "Deep Contiguous (>4 ft)",
}


def summary_for_point(lat: float, lon: float,
                       radius_m: float = DEFAULT_RADIUS_M,
                       max_developments: int = DEFAULT_MAX_PER_QUERY) -> dict:
    """Return the N nearest tier-1-3 NYCHA developments to (lat, lon)
    within radius_m, with their pre-computed exposure flags from the
    register catalog at data/registers/nycha.json.

    The catalog is the source of truth for which developments are
    flood-exposed (the bake script ran the polygon-overlap math once,
    citywide). Per-query work is haversine + dict lookup — sub-ms even
    on the HF Space CPU. Developments outside the tier-1-3 catalog
    (truly unexposed inland sites) are intentionally not surfaced;
    "no NYCHA developments at risk within 1 mi" is the honest answer
    for low-exposure queries.
    """
    from app.registers._loader import nearest_n
    hits = nearest_n("nycha", lat, lon, radius_m, max_developments)
    if not hits:
        return {"available": False,
                "n_developments": 0,
                "radius_m": radius_m,
                "developments": []}

    findings: list[DevelopmentFinding] = []
    for distance_m, row in hits:
        snap = row.get("snap") or {}
        dep = snap.get("dep") or {}
        microtopo = snap.get("microtopo") or {}

        def _dep_class(scen: str) -> int:
            d = dep.get(scen) or {}
            return int(d.get("depth_class") or 0)

        c2080 = _dep_class("dep_extreme_2080")
        c2050 = _dep_class("dep_moderate_2050")
        ccur  = _dep_class("dep_moderate_current")

        elev = microtopo.get("point_elev_m")
        hand = microtopo.get("aoi_hand_m") or microtopo.get("hand_m")

        findings.append(DevelopmentFinding(
            development=str(row.get("name", "")),
            tds_num=str(row.get("tds_num", "")),
            borough=str(row.get("borough", "")),
            centroid_lat=round(float(row["lat"]), 5),
            centroid_lon=round(float(row["lon"]), 5),
            distance_m=round(distance_m, 1),
            rep_elevation_m=round(float(elev), 2) if elev is not None else None,
            rep_hand_m=round(float(hand), 2) if hand is not None else None,
            inside_sandy_2012=bool(snap.get("sandy")),
            dep_extreme_2080_class=c2080,
            dep_extreme_2080_label=_DEPTH_LABEL.get(c2080, "outside"),
            dep_moderate_2050_class=c2050,
            dep_moderate_2050_label=_DEPTH_LABEL.get(c2050, "outside"),
            dep_moderate_current_class=ccur,
            dep_moderate_current_label=_DEPTH_LABEL.get(ccur, "outside"),
        ))

    n_in_sandy = sum(1 for f in findings if f.inside_sandy_2012)
    n_in_2080 = sum(1 for f in findings if f.dep_extreme_2080_class > 0)
    return {
        "available": True,
        "n_developments": len(findings),
        "radius_m": radius_m,
        "n_inside_sandy_2012": n_in_sandy,
        "n_in_dep_extreme_2080": n_in_2080,
        "developments": [vars(f) for f in findings],
        "citation": ("Pre-computed from NYC Open Data NYCHA Developments "
                     "(phvi-damg) joined to Sandy 2012 Inundation Zone "
                     "(5xsi-dfpx) + NYC DEP Stormwater Flood Maps + "
                     "USGS 3DEP DEM. See data/registers/nycha.json."),
    }


def main() -> int:
    import argparse
    ap = argparse.ArgumentParser()
    ap.add_argument("--lat", type=float, required=True)
    ap.add_argument("--lon", type=float, required=True)
    ap.add_argument("--radius", type=float, default=DEFAULT_RADIUS_M)
    ap.add_argument("--max", type=int, default=DEFAULT_MAX_PER_QUERY)
    args = ap.parse_args()
    s = summary_for_point(args.lat, args.lon, args.radius, args.max)
    print(json.dumps(s, indent=2, default=str))
    return 0


if __name__ == "__main__":
    sys.exit(main())