Post-disaster impact estimation

# --- 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 ---
extent_radius_km = 80
population_source = "worldpop_100m"
age_group = "all"
infrastructure_layers = ["schools", "hospitals"]
timesteps = 3

# Parameters
run_id = "023eb058-6f9c-485a-aa8f-4b32ba3ffaee"
inputs_dir = "/app/data/work/023eb058-6f9c-485a-aa8f-4b32ba3ffaee/inputs"
output_dir = "/app/data/work/023eb058-6f9c-485a-aa8f-4b32ba3ffaee/outputs"
area = {"id": 4, "code": "NPL.3_1", "name": "Far-Western", "level": 1, "country": "NPL", "country_name": "Nepal", "bbox": [80.0601, 28.3944, 81.8097, 30.2468]}
area_geojson = "/app/data/work/023eb058-6f9c-485a-aa8f-4b32ba3ffaee/inputs/area.geojson"
subdivisions_geojson = "/app/data/work/023eb058-6f9c-485a-aa8f-4b32ba3ffaee/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/023eb058-6f9c-485a-aa8f-4b32ba3ffaee/inputs/event.json"
montandon_stac_url = "https://montandon-eoapi-stage.ifrc.org/stac"
hazard = ""
age_group = "all"
timesteps = 3
extent_radius_km = 80.0
population_source = "worldpop_100m"
infrastructure_layers = ["schools", "hospitals"]

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")
print("Event:", props.get("title", "(standalone demo event)"), "·", corr_id)

extent = None
token = os.environ.get("MONTANDON_API_TOKEN")
if event and token:
    try:
        import httpx
        r = httpx.post(f"{montandon_stac_url.rstrip('/')}/search", headers={"Authorization": f"Bearer {token}"}, timeout=60,
                       json={"limit": 50, "filter-lang": "cql2-json", "filter": {"op": "=", "args": [{"property": "monty:corr_id"}, corr_id]}})
        r.raise_for_status()
        related = [f for f in r.json().get("features", []) if f.get("geometry") and f["geometry"]["type"] in ("Polygon", "MultiPolygon")]
        if related:
            extent = gpd.GeoDataFrame.from_features(related, crs=4326)
            print(f"Using {len(extent)} polygon product(s) linked to {corr_id}")
    except Exception as exc:
        print("Montandon lookup failed, using stand-in extent:", exc)
if extent is None:
    if event and event.get("geometry"):
        center = shapely.geometry.shape(event["geometry"]).centroid
    else:
        center = area_gdf.union_all().centroid
    extent = gpd.GeoDataFrame({"product": ["stand-in buffer"]}, geometry=[center.buffer(extent_radius_km / 111.0)], crs=4326)
    print(f"Stand-in extent: {extent_radius_km} km around {center.x:.2f}, {center.y:.2f}")
impact_zone = extent.union_all()

# ---------------------------------------------------------------------------------------------
# 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

seed = abs(hash(corr_id)) % 10_000
X, Y, step = make_grid(area_gdf.total_bounds, n=220)
share = {"all": 1.0, "under_18": 0.33, "under_5": 0.10, "over_65": 0.06}[age_group]
population = population_grid(X, Y, subdivisions, seed) * share
cells = cells_to_gdf(X, Y, step, population=population)
cells = cells[cells.intersects(area_gdf.union_all())]

# successive products: the extent grows with each time step (illustrative)
frames = []
for t in range(1, timesteps + 1):
    zone_t = impact_zone.buffer((t - timesteps) * 0.15 * impact_zone.length / (2 * math.pi)) if timesteps > 1 else impact_zone
    inside = centroids(cells.geometry).within(zone_t)
    c = cells.copy(); c["affected"] = c["population"].where(inside, 0)
    s = zonal_sum(c, subdivisions, ["population", "affected"])
    s["timestep"] = t
    frames.append(s)
latest = frames[-1]
latest["affected_share"] = (latest["affected"] / latest["population"].replace(0, np.nan)).fillna(0)
evolution = pd.concat(frames).groupby("timestep")["affected"].sum().round(0)
print("Affected population per product time step:"); print(evolution.to_string())

table = latest[["name", "population", "affected", "affected_share"]].sort_values("affected", ascending=False)
table.to_csv(Path(output_dir) / "impact_by_zone.csv", index=False)
latest.drop(columns=["timestep"]).to_file(Path(output_dir) / "impact_by_zone.geojson", driver="GeoJSON")
table.style.format({"population": "{:,.0f}", "affected": "{:,.0f}", "affected_share": "{:.1%}"}).background_gradient(subset=["affected_share"], 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}"

Row(
    Stat(label=f"Population in area ({age_group})", value=fmt_int(latest["population"].sum())),
    Stat(label="Affected (latest product)", value=fmt_int(latest["affected"].sum())),
    Stat(label="Share affected", value=f"{latest['affected'].sum() / max(latest['population'].sum(), 1):.1%}"),
    Stat(label="Products compared", value=timesteps),
    gap="16px",
)

vmax = max(float(latest["affected_share"].max()), 0.01)
zones_layer = PolygonLayer.from_geopandas(latest[["name", "affected", "affected_share", "geometry"]], get_fill_color=colors_for(latest["affected_share"], "OrRd", 0, vmax), get_line_color=[60, 60, 60], line_width_min_pixels=1, opacity=0.7)
extent_layer = PolygonLayer.from_geopandas(extent[["geometry"]], get_fill_color=[46, 111, 176, 50], get_line_color=[46, 111, 176], line_width_min_pixels=2)
Map([zones_layer, extent_layer])

Row(Column(LayerToggle(layer=zones_layer, label="Affected share"), LayerToggle(layer=extent_layer, label="Observed extent"), gap="6px"), legend_for("OrRd", 0, vmax, title="Affected population share", fmt="{:.1%}"), gap="24px", align="start")

chart = Chart(title=f"{props.get('title', 'Event')} · affected population per product", x_label="product time step", y_label="affected population", height=320)
chart.add_series(x=[int(i) for i in evolution.index], y=[float(v) for v in evolution.values], name="affected population", color="#E4572E")
chart