6. Hurricane Helene, 48 hours later

import sys, pathlib, json
ROOT = pathlib.Path("..").resolve() if pathlib.Path("..").exists() else pathlib.Path("../..").resolve()
sys.path.insert(0, str(ROOT))

import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import Polygon, Point, shape

import lonboard
from lonboard import Map, ScatterplotLayer, PolygonLayer, BitmapLayer

from widgets.hurricane_dashboard.widget import HurricaneDashboard

DATA = ROOT / "data" / "helene"
print("data dir:", DATA)

def parse_hurdat2(path):
    rows = []
    with open(path) as f:
        cur_id = cur_name = None
        for line in f:
            parts = [p.strip() for p in line.strip(", \n").split(",")]
            if parts and parts[0].startswith("AL") and len(parts) == 3:
                cur_id, cur_name = parts[0], parts[1]
                continue
            if not cur_id or len(parts) < 8:
                continue
            try:
                rows.append({
                    "storm_id": cur_id,
                    "name": cur_name,
                    "date": parts[0],
                    "time": parts[1],
                    "lat": float(parts[4][:-1]) * (1 if parts[4][-1] == "N" else -1),
                    "lon": float(parts[5][:-1]) * (-1 if parts[5][-1] == "W" else 1),
                    "wind_kt": int(parts[6]),
                    "pressure_mb": int(parts[7]) if parts[7] != "-999" else np.nan,
                })
            except (ValueError, IndexError):
                continue
    return pd.DataFrame(rows)

storms = parse_hurdat2(DATA / "hurdat2.txt")
helene = storms[storms.storm_id == "AL092024"].copy()
helene["wind_mph"] = helene["wind_kt"] * 1.151
print(f"{len(helene)} track points; peak wind {helene.wind_mph.max():.0f} mph")
helene.head()

fig, ax = plt.subplots(figsize=(7, 4.2))
sc = ax.scatter(helene.lon, helene.lat, c=helene.wind_mph, cmap="magma_r",
                s=22, edgecolor="white", linewidth=0.4)
ax.plot(helene.lon, helene.lat, color="0.4", linewidth=1, alpha=0.6, zorder=0)
# Asheville reference point
ax.plot(-82.55, 35.59, marker="*", color="#67c1ff", markersize=14, zorder=5)
ax.annotate("Asheville", xy=(-82.55, 35.59), xytext=(-79.5, 36.5),
            color="#67c1ff", fontsize=10, fontweight="bold")
plt.colorbar(sc, ax=ax, label="max wind (mph)")
ax.set_xlabel("longitude"); ax.set_ylabel("latitude")
ax.set_title("Hurricane Helene best track (HURDAT2 / NOAA NHC)")
ax.grid(alpha=0.25)
plt.tight_layout()
plt.show()

def load_osm_ways(path, kind):
    """Return a GeoDataFrame with one row per way (polygon) or node (point)."""
    with open(path) as f:
        data = json.load(f)
    rows = []
    for el in data["elements"]:
        if el["type"] == "way" and "geometry" in el:
            ring = [(g["lon"], g["lat"]) for g in el["geometry"]]
            if len(ring) >= 3:
                try:
                    poly = Polygon(ring)
                    if poly.is_valid and poly.area > 0:
                        rows.append({"id": el["id"], "geometry": poly,
                                     "tags": el.get("tags", {})})
                except Exception:
                    pass
        elif el["type"] == "node":
            rows.append({"id": el["id"], "geometry": Point(el["lon"], el["lat"]),
                         "tags": el.get("tags", {})})
    gdf = gpd.GeoDataFrame(rows, crs="EPSG:4326")
    return gdf

bldgs = load_osm_ways(DATA / "osm_buildings.json", "building")
hosps = load_osm_ways(DATA / "osm_hospitals.json", "hospital")
print(f"{len(bldgs):,} buildings, {len(hosps):,} hospitals")

# Centroids for fast map rendering
bldgs["centroid"] = bldgs.geometry.centroid

flood = gpd.read_file(DATA / "flood_extent.geojson")
print(f"{len(flood)} flood polygons; depths {sorted(flood.depth_m.unique())} m")

# sjoin building centroids onto flood polygons
cent_gdf = gpd.GeoDataFrame(bldgs[["id"]].copy(), geometry=bldgs["centroid"], crs="EPSG:4326")
joined = gpd.sjoin(cent_gdf, flood[["depth_m", "geometry"]], how="left", predicate="within")
# Where a building falls in multiple overlapping polygons, take the deepest
depth_per_building = joined.groupby("id")["depth_m"].max().fillna(0.0)
bldgs["flood_depth_m"] = bldgs["id"].map(depth_per_building).fillna(0.0)

print("flooded buildings:", int((bldgs.flood_depth_m > 0).sum()))
print("max depth assigned:", float(bldgs.flood_depth_m.max()))

def bucket(d):
    if d <= 0:    return "dry"
    if d <= 1:    return "shallow"
    if d <= 3:    return "deep"
    return "drowned"

bldgs["bucket"] = bldgs["flood_depth_m"].map(bucket)
counts = bldgs["bucket"].value_counts().reindex(["dry", "shallow", "deep", "drowned"], fill_value=0)
print(counts)

# Hospitals affected = those whose centroid is in the flood polygon
hosp_pts = gpd.GeoDataFrame(hosps[["id", "tags"]].copy(),
                            geometry=hosps.geometry.centroid,
                            crs="EPSG:4326")
hosp_joined = gpd.sjoin(hosp_pts, flood[["geometry"]], how="left", predicate="within")
hospitals_affected = int(hosp_joined["index_right"].notna().sum())

# Rough population estimate: 2.4 people per residential building (urban avg).
PEOPLE_PER_BLDG = 2.4

# Threshold semantics: 0 = show everything, 1 = ≥shallow, 2 = ≥deep, 3 = drowned only
threshold_layer_map = [
    ["b_dry", "b_shallow", "b_deep", "b_drowned"],
    ["b_shallow", "b_deep", "b_drowned"],
    ["b_deep", "b_drowned"],
    ["b_drowned"],
]

stats_tables = []
for buckets in threshold_layer_map:
    name_to_count = {f"b_{b}": int(counts.get(b, 0)) for b in counts.index}
    visible_buildings = sum(name_to_count.get(k, 0) for k in buckets)
    # Hospitals only count as "affected" if any flooded layer is on
    hosp_count = hospitals_affected if any(k != "b_dry" for k in buckets) else 0
    # People estimate: same buildings × ratio, but only flooded ones contribute
    flooded_visible = sum(name_to_count.get(k, 0) for k in buckets if k != "b_dry")
    stats_tables.append({
        "buildings": visible_buildings,
        "hospitals": hosp_count,
        "people_est": int(flooded_visible * PEOPLE_PER_BLDG),
    })
stats_tables

import base64

def png_to_data_url(path):
    return "data:image/png;base64," + base64.b64encode(path.read_bytes()).decode()

def safe_layer(points, color, radius=14):
    """ScatterplotLayer that tolerates empty buckets by inserting a phantom
    point off the visible map. Lonboard's geoarrow conversion can't handle a
    truly empty geometry list."""
    if len(points) == 0:
        points = gpd.GeoSeries([Point(0, 0)], crs="EPSG:4326")
    geom = gpd.GeoDataFrame(geometry=gpd.GeoSeries(points.values, crs="EPSG:4326"))
    return ScatterplotLayer.from_geopandas(
        geom, get_fill_color=color, get_radius=radius,
        radius_min_pixels=1.5, opacity=0.85,
    )

# Layers (each is a separate lonboard layer with its own model_id)
bitmap_after = BitmapLayer(
    image=png_to_data_url(DATA / "modis_after.png"),
    bounds=[-82.70, 35.45, -82.40, 35.72],  # west, south, east, north
    opacity=0.55,
)

flood_layer = PolygonLayer.from_geopandas(
    flood,
    get_fill_color=[80, 160, 235, 110],
    get_line_color=[80, 160, 235, 220],
    line_width_min_pixels=1,
)

def bucket_centroids(name):
    return bldgs.loc[bldgs.bucket == name, "centroid"]

building_layers = {
    "b_dry":     safe_layer(bucket_centroids("dry"),     [110, 110, 110, 160],  8),
    "b_shallow": safe_layer(bucket_centroids("shallow"), [255, 220,  80, 240], 16),
    "b_deep":    safe_layer(bucket_centroids("deep"),    [255, 140,  60, 240], 20),
    "b_drowned": safe_layer(bucket_centroids("drowned"), [220,  50,  50, 250], 24),
}

hosp_layer = ScatterplotLayer.from_geopandas(
    gpd.GeoDataFrame(geometry=gpd.GeoSeries(hosp_pts.geometry.values, crs="EPSG:4326")),
    get_fill_color=[230, 80, 200, 250],
    get_radius=80,
    radius_min_pixels=6,
    stroked=True,
    get_line_color=[255, 255, 255, 255],
    line_width_min_pixels=2,
)

m = Map(
    layers=[bitmap_after, flood_layer,
            building_layers["b_dry"], building_layers["b_shallow"],
            building_layers["b_deep"], building_layers["b_drowned"],
            hosp_layer],
    view_state={"longitude": -82.555, "latitude": 35.585, "zoom": 12.5, "pitch": 0, "bearing": 0},
    show_tooltip=False,
)

layer_uuids = {
    "modis_after": bitmap_after.model_id,
    "flood": flood_layer.model_id,
    "b_dry": building_layers["b_dry"].model_id,
    "b_shallow": building_layers["b_shallow"].model_id,
    "b_deep": building_layers["b_deep"].model_id,
    "b_drowned": building_layers["b_drowned"].model_id,
    "hospitals": hosp_layer.model_id,
}

dashboard = HurricaneDashboard(
    layer_uuids=layer_uuids,
    initial_visible={
        "modis_after": True,
        "flood": True,
        "hospitals": True,
        # building layers are controlled by the threshold, default = all visible
        "b_dry": True, "b_shallow": True, "b_deep": True, "b_drowned": True,
    },
    toggle_layers=[
        {"key": "modis_after", "label": "Satellite (MODIS post-storm)"},
        {"key": "flood",       "label": "Flood extent"},
        {"key": "hospitals",   "label": "Hospitals"},
    ],
    threshold_layer_map=threshold_layer_map,
    stats_tables=stats_tables,
    threshold=0,
)
print(f"map: {len(m.layers)} layers · dashboard: {len(layer_uuids)} layer bindings")

from IPython.display import display
display(dashboard)
display(m)