#!/usr/bin/env python3
"""Probe ground elevation across the demo block with the USGS 3DEP point service.

Usage:
    Geo/scripts/elevation_probe.py Geo/raw/osm-api-map-2026-11-06.osm

Samples a grid along the boulevard (centreline and both sidewalk centrelines),
the hero footprint corners, and the crossing midpoints, then reports the relief
across the demo and whether a separate elevation source is needed.

The decision rule follows `AGENTS.md`: acquire separately licensed elevation data
unless the DEM shows less than 1 m of height variation across the demo bounds.

Writes the raw responses to `Geo/raw/epqs-<date>.json` so the numbers can be
re-checked without re-querying. Service: USGS Elevation Point Query Service,
3DEP 1 m seamless DEM (CONUS). Standard library only.
"""

import argparse
import concurrent.futures
import datetime
import json
import os
import time
import urllib.parse
import urllib.request

import osm_common as oc

EPQS = "https://epqs.nationalmap.gov/v1/json"
UA = "gig-line-geo/0.1 (local development, elevation check)"


def query(lat: float, lon: float, retries: int = 5) -> dict:
    params = urllib.parse.urlencode(
        {"x": f"{lon:.7f}", "y": f"{lat:.7f}", "wkid": "4326", "units": "Meters", "includeDate": "false"}
    )
    url = f"{EPQS}?{params}"
    for attempt in range(retries):
        req = urllib.request.Request(url, headers={"User-Agent": UA})
        try:
            with urllib.request.urlopen(req, timeout=40) as resp:
                payload = json.loads(resp.read().decode("utf-8"))
            value = payload.get("value")
            if value in (None, "", "No Data"):
                raise ValueError(f"no elevation returned: {payload}")
            return {
                "lat": round(lat, 7),
                "lon": round(lon, 7),
                "elevation_m": float(value),
                "raster_id": payload.get("rasterId"),
                "resolution_m": payload.get("resolution"),
            }
        except Exception as exc:  # noqa: BLE001
            if attempt == retries - 1:
                # a single dead point must not kill the whole probe
                print(f"    ! no elevation at {lat:.6f},{lon:.6f}: {exc}", flush=True)
                return {"lat": round(lat, 7), "lon": round(lon, 7), "elevation_m": None,
                        "raster_id": None, "resolution_m": None}
            time.sleep(3 * (attempt + 1))
    return {}


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("osm")
    ap.add_argument("--blvd", default="Hollywood Boulevard")
    ap.add_argument("--west", default="North Orange Drive")
    ap.add_argument("--east", default="North Highland Avenue")
    ap.add_argument("--hero", type=int, default=424024147, help="hero building OSM way id")
    ap.add_argument("--workers", type=int, default=4, help="parallel requests; the service is slow")
    ap.add_argument("--pause", type=float, default=0.0)
    args = ap.parse_args()

    nodes, ways = oc.load(args.osm)
    line, lat0, lon0, lengths, block_len = oc.block_extent(nodes, ways, args.blvd, args.west, args.east)

    samples = []
    grid = [block_len * i / 4 for i in range(5)]
    for s in grid:
        for label, offset in (("south sidewalk", -12.6), ("centreline", 0.0), ("north sidewalk", 13.6)):
            lat, lon = oc.point_at(line, lengths, s, offset, lat0, lon0)
            samples.append({"label": f"s={s:5.1f} m {label}", "lat": lat, "lon": lon})

    hero_refs = next((refs for wid, _, refs in ways if str(wid) == str(args.hero)), None)
    if hero_refs:
        hero_pts = oc.way_points(nodes, hero_refs)
        lats = [p[0] for p in hero_pts]
        lons = [p[1] for p in hero_pts]
        corners = [
            ("near NW", min(lats), min(lons)),
            ("near NE", min(lats), max(lons)),
            ("far SW", max(lats), min(lons)),
            ("far SE", max(lats), max(lons)),
        ]
        for label, lat, lon in corners:
            samples.append({"label": f"hero {label}", "lat": lat, "lon": lon})

    for wid, tags, refs in ways:
        if tags.get("footway") != "crossing":
            continue
        pts = oc.way_points(nodes, refs)
        if not pts:
            continue
        mid = (sum(p[0] for p in pts) / len(pts), sum(p[1] for p in pts) / len(pts))
        dist, s, _ = oc.project(mid, line, lat0, lon0, lengths)
        # only crossings of THIS boulevard, not side-street crossings inside the window
        if 0 <= s <= block_len and dist <= 4.0:
            samples.append({"label": f"crossing way {wid} at s={s:.0f} m", "lat": mid[0], "lon": mid[1]})

    def fetch(sample):
        r = query(sample["lat"], sample["lon"])
        r["label"] = sample["label"]
        return r

    started = time.time()
    with concurrent.futures.ThreadPoolExecutor(max_workers=max(1, args.workers)) as pool:
        futures = {pool.submit(fetch, s): s for s in samples}
        results = []
        for done, future in enumerate(concurrent.futures.as_completed(futures), start=1):
            r = future.result()
            results.append(r)
            value = "no data" if r["elevation_m"] is None else f"{r['elevation_m']:7.2f} m"
            print(f"  [{done:02d}/{len(samples)}] {r['label']:28} {value:>10}", flush=True)
            if args.pause:
                time.sleep(args.pause)
    order = {s["label"]: i for i, s in enumerate(samples)}
    results.sort(key=lambda r: order.get(r["label"], 999))
    print(f"query time {time.time() - started:.0f} s for {len(samples)} samples")

    walkable = [r for r in results if r["elevation_m"] is not None and not r["label"].startswith("hero")]
    building = [r for r in results if r["elevation_m"] is not None and r["label"].startswith("hero")]
    values = [r["elevation_m"] for r in walkable]
    lo, hi = min(values), max(values)
    relief = hi - lo

    print("\nsamples by group (walkable surfaces only below; building interiors are diagnostics):")
    for key in ("crossing", "hero"):
        group = [r for r in results if r["elevation_m"] is not None and key in r["label"]]
        if group:
            g = [r["elevation_m"] for r in group]
            print(f"  {key:12} n={len(group):2}  min {min(g):.2f}  max {max(g):.2f}  spread {max(g)-min(g):.2f} m")

    print("\ncross-slope check (south sidewalk minus north sidewalk):")
    for s in [block_len * i / 4 for i in range(5)]:
        south = next(r for r in results if r["label"] == f"s={s:5.1f} m south sidewalk")
        north = next(r for r in results if r["label"] == f"s={s:5.1f} m north sidewalk")
        drop = south["elevation_m"] - north["elevation_m"]
        print(f"  s={s:5.1f} m: south {south['elevation_m']:.2f} m, north {north['elevation_m']:.2f} m, "
              f"difference {drop:+.2f} m over ~26 m")

    if building:
        print("\nwarning: points inside the hero footprint are not ground reference; 3DEP can return")
        print("         building or void values there. They are excluded from the relief figure.")
    print(f"\nelevation range on walkable surfaces: {lo:.2f} to {hi:.2f} m, relief {relief:.2f} m")
    first = next(r for r in walkable)
    print(f"source: USGS 3DEP point service, raster {first['raster_id']}, resolution {first['resolution_m']} m")
    if relief < 1.0:
        print("decision: relief under 1 m, a separate elevation source is NOT required (AGENTS.md rule)")
    else:
        print("decision: relief at or above 1 m, a separately licensed elevation source IS required")
    print("datum: 3DEP CONUS heights are orthometric NAVD88; treat absolute values as such, not ellipsoidal")

    today = datetime.date.today().isoformat()
    out_path = os.path.join("Geo", "raw", f"epqs-{today}.json")
    os.makedirs(os.path.dirname(out_path), exist_ok=True)
    payload = {
        "service": EPQS,
        "processed_utc": datetime.datetime.now(datetime.timezone.utc).isoformat(timespec="seconds"),
        "source_dem": "USGS 3DEP seamless 1 m (CONUS)",
        "vertical_datum": "NAVD88 orthometric (3DEP CONUS)",
        "query": {
            "blvd": args.blvd,
            "west": args.west,
            "east": args.east,
            "hero_way": args.hero,
            "block_len_m": round(block_len, 2),
        },
        "summary": {"min_m": lo, "max_m": hi, "relief_m": round(relief, 3), "samples": len(results)},
        "samples": results,
    }
    with open(out_path, "w", encoding="utf-8") as handle:
        json.dump(payload, handle, indent=2)
    print(f"wrote {out_path}")


if __name__ == "__main__":
    main()