Operational response prioritisation

# --- Reserved parameters (injected by the Notebook Factory; defaults keep the notebook runnable standalone) ---
run_id = ""
inputs_dir = ""
output_dir = "outputs"
area = {}                      # {id, code, name, level, country, country_name, bbox}
area_geojson = ""              # FeatureCollection with the target area
subdivisions_geojson = ""      # FeatureCollection with the aggregation subdivisions
event = {}                     # Montandon STAC item (dict)
event_json = ""
montandon_stac_url = "https://montandon-eoapi-stage.ifrc.org/stac"
hazard = ""

# --- User parameters ---
weighting = "geometric"
w_severity = 0.4
w_vulnerability = 0.3
w_history = 0.2
w_gap = 0.1
top_n = 10
extent_radius_km = 60

# Parameters
run_id = "84481c51-1dc4-4c3f-87bc-28e4b0a5d9b9"
inputs_dir = "/app/data/work/84481c51-1dc4-4c3f-87bc-28e4b0a5d9b9/inputs"
output_dir = "/app/data/work/84481c51-1dc4-4c3f-87bc-28e4b0a5d9b9/outputs"
area = {"id": 16, "code": "NPL.4.2_1", "name": "Karnali", "level": 2, "country": "NPL", "country_name": "Nepal", "bbox": [81.258, 28.7219, 83.6825, 30.447]}
area_geojson = "/app/data/work/84481c51-1dc4-4c3f-87bc-28e4b0a5d9b9/inputs/area.geojson"
subdivisions_geojson = "/app/data/work/84481c51-1dc4-4c3f-87bc-28e4b0a5d9b9/inputs/subdivisions.geojson"
event = {"id": "example-event-npl-earthquake-2026", "bbox": [84.73, 28.23, 84.73, 28.23], "type": "Feature", "links": [{"rel": "self", "href": "https://montandon-eoapi-stage.ifrc.org/stac/collections/reference-events/items/example-event-npl-earthquake-2026", "type": "application/geo+json"}], "assets": {}, "geometry": {"type": "Point", "coordinates": [84.73, 28.23]}, "collection": "reference-events", "properties": {"roles": ["event"], "title": "Earthquake in Nepal (Gorkha region) \u2014 example event", "source": "example", "updated": "2026-09-01T08:00:00Z", "datetime": "2026-09-01T06:11:26Z", "description": "Illustrative M6.8 earthquake near Gorkha used to exercise the event-triggered pipeline.", "end_datetime": "2026-09-01T06:11:26Z", "monty:corr_id": "20260901-NPL-GEO-EAR-GRO-1-GCDB", "start_datetime": "2026-09-01T06:11:26Z", "monty:hazard_codes": ["nat-geo-ear-gro"], "monty:country_codes": ["NPL"], "monty:episode_number": 1}, "stac_version": "1.0.0"}
event_json = "/app/data/work/84481c51-1dc4-4c3f-87bc-28e4b0a5d9b9/inputs/event.json"
montandon_stac_url = "https://montandon-eoapi-stage.ifrc.org/stac"
hazard = ""
top_n = 10
w_gap = 0.1
w_history = 0.2
weighting = "geometric"
w_severity = 0.4
w_vulnerability = 0.3
extent_radius_km = 60.0

import json, math, os, warnings
from pathlib import Path
import numpy as np
import pandas as pd
import geopandas as gpd
import shapely
from shapely.geometry import box

warnings.filterwarnings("ignore", message=".*geographic CRS.*")
Path(output_dir).mkdir(parents=True, exist_ok=True)

def centroids(geoms):
    """Centroids computed in an equal-area projection (avoids the geographic-CRS warning)."""
    return geoms.to_crs(6933).centroid.to_crs(4326)

def load_area_and_subdivisions():
    """Read the boundaries prepared by the app, or build a small demo area when run standalone."""
    if area_geojson and Path(area_geojson).exists():
        area_gdf = gpd.read_file(area_geojson)
        subs = gpd.read_file(subdivisions_geojson) if subdivisions_geojson and Path(subdivisions_geojson).exists() else area_gdf.copy()
        return area_gdf, subs
    demo = gpd.GeoDataFrame({"name": ["Demo area"], "code": ["DEMO"], "level": [1]}, geometry=[box(83.5, 26.5, 88.0, 30.5)], crs=4326)
    cells, names = [], []
    minx, miny, maxx, maxy = demo.total_bounds
    for i in range(3):
        for j in range(3):
            cells.append(box(minx + i * (maxx - minx) / 3, miny + j * (maxy - miny) / 3, minx + (i + 1) * (maxx - minx) / 3, miny + (j + 1) * (maxy - miny) / 3))
            names.append(f"Demo district {i * 3 + j + 1}")
    subs = gpd.GeoDataFrame({"name": names, "code": [f"DEMO.{k}" for k in range(9)], "level": [2] * 9}, geometry=cells, crs=4326)
    return demo, subs

area_gdf, subdivisions = load_area_and_subdivisions()
area_name = area.get("name") or area_gdf.iloc[0].get("name", "Demo area")
print(f"Area: {area_name} · {len(subdivisions)} subdivisions")
subdivisions[["code", "name", "level"]].head()

if not event and event_json and Path(event_json).exists():
    event = json.loads(Path(event_json).read_text())
props = event.get("properties", {}) if event else {}
corr_id = props.get("monty:corr_id", "demo-event")
event_title = props.get("title", "Standalone demo event")
print("Event:", event_title, "·", corr_id)

# ---------------------------------------------------------------------------------------------
# ILLUSTRATIVE DATA. The hazard and population grids below are synthetic, generated so the
# pipeline can be exercised end-to-end without large downloads. Replace `hazard_grid()` and
# `population_grid()` with real loaders (STORM / Aqueduct / GEM / GlobFire rasters, WorldPop / GHSL)
# and the rest of the notebook works unchanged.
# ---------------------------------------------------------------------------------------------
def make_grid(bounds, n=120):
    minx, miny, maxx, maxy = bounds
    step = max((maxx - minx), (maxy - miny)) / n
    xs = np.arange(minx + step / 2, maxx, step)
    ys = np.arange(miny + step / 2, maxy, step)
    X, Y = np.meshgrid(xs, ys)
    return X, Y, step

def hazard_grid(X, Y, seed=0):
    """Synthetic hazard intensity field in the units of the chosen hazard."""
    rng = np.random.default_rng(seed)
    cx, cy = rng.uniform(X.min(), X.max()), rng.uniform(Y.min(), Y.max())
    d = np.hypot(X - cx, Y - cy) / max(X.max() - X.min(), 1e-6)
    ridges = 0.5 + 0.5 * np.sin(X * 9 + rng.uniform(0, 6)) * np.cos(Y * 7 + rng.uniform(0, 6))
    field = np.clip(1.0 - 1.6 * d, 0, 1) * 0.6 + ridges * 0.4
    return np.clip(field, 0, 1)   # 0..1

def population_grid(X, Y, subs, seed=0):
    """Synthetic population: denser near subdivision centroids."""
    rng = np.random.default_rng(seed + 1)
    pop = np.zeros_like(X)
    for geom in subs.geometry:
        c = geom.centroid
        pop += np.exp(-((X - c.x) ** 2 + (Y - c.y) ** 2) / (0.02 + rng.uniform(0, 0.05)))
    pop *= rng.uniform(0.6, 1.4, size=X.shape)
    return pop / pop.sum() * 250_000 * len(subs)   # people per cell

def cells_to_gdf(X, Y, step, **cols):
    geoms = [box(x - step / 2, y - step / 2, x + step / 2, y + step / 2) for x, y in zip(X.ravel(), Y.ravel())]
    return gpd.GeoDataFrame({k: v.ravel() for k, v in cols.items()}, geometry=geoms, crs=4326)

def zonal_sum(cells, zones, value_cols):
    """Sum cell values per zone using cell centroids (fast enough for illustrative grids)."""
    pts = cells.copy()
    pts["geometry"] = centroids(pts.geometry)
    joined = gpd.sjoin(pts, zones[["geometry"]], how="inner", predicate="within")
    agg = joined.groupby("index_right")[value_cols].sum()
    out = zones.copy()
    for c in value_cols:
        out[c] = agg[c].reindex(zones.index).fillna(0).values
    return out

rng = np.random.default_rng(abs(hash(corr_id)) % 10_000)
zones = subdivisions.copy().reset_index(drop=True)
zones["area_km2"] = zones.to_crs(6933).area / 1e6
center = shapely.geometry.shape(event["geometry"]).centroid if event and event.get("geometry") else area_gdf.union_all().centroid
extent = center.buffer(extent_radius_km / 111.0)

# 1 · observed severity: share of the zone inside the extent, weighted by (synthetic) population density
inside = zones.geometry.intersection(extent).to_crs(6933).area / 1e6
density = rng.lognormal(4, 0.6, len(zones))
zones["severity"] = (inside / zones["area_km2"]).clip(0, 1) * (density / density.max())

# 2 · vulnerability & coping capacity (INFORM-style, 0-10 -> 0-1)
zones["vulnerability"] = rng.uniform(2, 8, len(zones)) / 10

# 3 · historical burden: count of past Montandon events in the zone's bbox (falls back to synthetic)
def montandon_history(geom):
    token = os.environ.get("MONTANDON_API_TOKEN")
    if not token:
        return None
    try:
        import httpx
        r = httpx.post(f"{montandon_stac_url.rstrip('/')}/search", headers={"Authorization": f"Bearer {token}"}, timeout=60,
                       json={"collections": ["reference-events"], "bbox": list(geom.bounds), "limit": 1})
        r.raise_for_status()
        return r.json().get("numberMatched") or len(r.json().get("features", []))
    except Exception:
        return None
hist = [montandon_history(g) for g in zones.geometry]
zones["history_source"] = "montandon" if all(h is not None for h in hist) else "synthetic"
zones["history_count"] = [h if h is not None else int(rng.poisson(3)) for h in hist]
zones["history"] = zones["history_count"] / max(zones["history_count"].max(), 1)

# 4 · operational gap & accessibility
zones["gap"] = rng.uniform(0.1, 0.9, len(zones))
zones[["name", "severity", "vulnerability", "history", "gap"]].round(2).head()

w = np.array([w_severity, w_vulnerability, w_history, w_gap]); w = w / w.sum()
comp = zones[["severity", "vulnerability", "history", "gap"]].clip(1e-3, 1).values
if weighting == "geometric":
    zones["priority"] = np.exp((np.log(comp) * w).sum(axis=1))
else:
    zones["priority"] = (comp * w).sum(axis=1)
zones["priority"] = zones["priority"] / zones["priority"].max()
completeness = np.array([1.0, 0.6, 1.0 if zones["history_source"].iloc[0] == "montandon" else 0.5, 0.5])
zones["confidence"] = float((completeness * w).sum())
zones["rank"] = zones["priority"].rank(ascending=False, method="first").astype(int)

def rationale(r):
    drivers = sorted([("observed severity", r.severity), ("vulnerability", r.vulnerability), ("historical burden", r.history), ("operational gap", r.gap)], key=lambda x: -x[1])
    lead = ", ".join(f"{n} ({v:.2f})" for n, v in drivers[:2])
    return f"Ranked #{r['rank']} mainly due to {lead}; {r.history_count} past events on record ({r.history_source})."
zones["rationale"] = zones.apply(rationale, axis=1)
ranked = zones.sort_values("rank")
table = ranked[["rank", "name", "priority", "confidence", "severity", "vulnerability", "history", "gap", "rationale"]].head(top_n)
table.to_csv(Path(output_dir) / "prioritisation.csv", index=False)
ranked.to_file(Path(output_dir) / "prioritisation.geojson", driver="GeoJSON")
table.style.format({"priority": "{:.2f}", "confidence": "{:.2f}", "severity": "{:.2f}", "vulnerability": "{:.2f}", "history": "{:.2f}", "gap": "{:.2f}"}).background_gradient(subset=["priority"], cmap="OrRd")

# Visualisation: lonboard for maps, manywidgets for charts, controls and stats.
# All of these are anywidgets, so they stay interactive in the published (kernel-free) page.
from matplotlib import colormaps
from lonboard import Map, PolygonLayer, ScatterplotLayer
from manywidgets import Chart, Column, Legend, Row, Slider, Stat
from manywidgets.lonboard import LayerToggle

def colors_for(values, cmap="OrRd", vmin=None, vmax=None):
    v = np.asarray(values, dtype=float)
    vmin = np.nanmin(v) if vmin is None else vmin
    vmax = np.nanmax(v) if vmax is None else vmax
    norm = (v - vmin) / (vmax - vmin) if vmax > vmin else np.zeros_like(v)
    rgba = colormaps[cmap](np.nan_to_num(norm))
    return (rgba * 255).astype(np.uint8)

def legend_for(cmap, vmin, vmax, title="", n=5, fmt="{:.0%}"):
    stops = np.linspace(vmin, vmax, n)
    rgb = colors_for(stops, cmap, vmin, vmax)[:, :3]
    return Legend(entries=[[[int(c) for c in rgb[i]], fmt.format(stops[i])] for i in range(n)], title=title)

def bar_chart(names, values, title, y_label, color="#E4572E"):
    chart = Chart(chart_type="bar", title=title, y_label=y_label, height=max(280, 22 * len(names) + 120), legend_enabled=False)
    chart.set_options(scales={"x": {"type": "category", "ticks": {"autoSkip": False, "maxRotation": 60, "minRotation": 30}}})
    chart.add_series(x=[str(n) for n in names], y=[float(v) for v in values], name=y_label, series_type="bar", color=color)
    return chart

def fmt_int(v):
    return f"{float(v):,.0f}"

top = ranked.iloc[0]
Row(
    Stat(label="Areas ranked", value=len(ranked)),
    Stat(label="Top priority", value=top["name"]),
    Stat(label="Top score", value=f"{top['priority']:.2f}"),
    Stat(label="Confidence", value=f"{top['confidence']:.2f}"),
    gap="16px",
)

ranked_layer = PolygonLayer.from_geopandas(ranked[["name", "rank", "priority", "confidence", "rationale", "geometry"]], get_fill_color=colors_for(ranked["priority"], "OrRd", 0, 1), get_line_color=[60, 60, 60], line_width_min_pixels=1, opacity=0.75)
extent_layer = PolygonLayer.from_geopandas(gpd.GeoDataFrame(geometry=[extent], crs=4326), get_fill_color=[46, 111, 176, 40], get_line_color=[46, 111, 176], line_width_min_pixels=2)
Map([ranked_layer, extent_layer])

Row(Column(LayerToggle(layer=ranked_layer, label="Priority index"), LayerToggle(layer=extent_layer, label="Observed extent"), gap="6px"), legend_for("OrRd", 0, 1, title="Composite priority (0–1)", fmt="{:.2f}"), gap="24px", align="start")

bar_chart(table["name"], table["priority"], title=f"Priority index · top {min(top_n, len(ranked))} areas", y_label="composite priority")

from datetime import datetime, timezone
stac_dir = Path(output_dir) / "stac"; (stac_dir / "items").mkdir(parents=True, exist_ok=True)
now = datetime.now(timezone.utc).isoformat()
collection_id = f"response-prioritisation-{corr_id}".lower().replace(" ", "-")
event_link = [{"rel": "derived_from", "href": link["href"], "type": "application/geo+json"} for link in (event.get("links") or []) if link.get("rel") == "self"] if event else []
collection = {
    "type": "Collection", "stac_version": "1.0.0", "id": collection_id,
    "title": f"Response prioritisation · {event_title}",
    "description": f"Ranked priority areas for {event_title}, derived from exposure, observed impact, Montandon history and operational data.",
    "license": "CC-BY-4.0",
    "extent": {"spatial": {"bbox": [list(ranked.total_bounds)]}, "temporal": {"interval": [[now, None]]}},
    "keywords": ["montandon", "response-prioritisation", corr_id],
    "links": event_link,
    "summaries": {"monty:corr_id": [corr_id]},
}
(stac_dir / "collection.json").write_text(json.dumps(collection, indent=1))
for _, r in ranked.head(top_n).iterrows():
    item = {
        "type": "Feature", "stac_version": "1.0.0", "id": f"{collection_id}-{r.get('code', r['rank'])}".replace(".", "-"),
        "collection": collection_id, "geometry": json.loads(shapely.to_geojson(r.geometry)), "bbox": list(r.geometry.bounds),
        "properties": {
            "datetime": now, "title": f"#{r['rank']} {r['name']}", "monty:corr_id": corr_id, "roles": ["impact", "prioritisation"],
            "prioritisation:rank": int(r["rank"]), "prioritisation:score": round(float(r["priority"]), 3), "prioritisation:confidence": round(float(r["confidence"]), 2),
            "prioritisation:components": {"severity": round(float(r.severity), 3), "vulnerability": round(float(r.vulnerability), 3), "history": round(float(r.history), 3), "gap": round(float(r.gap), 3)},
            "prioritisation:rationale": r["rationale"], "admin:code": str(r.get("code", "")), "admin:level": int(r.get("level", 0) or 0),
        },
        "links": event_link, "assets": {},
    }
    (stac_dir / "items" / f"{item['id']}.json").write_text(json.dumps(item, indent=1))
print(f"wrote {collection_id} with {min(top_n, len(ranked))} items to {stac_dir}")