Population & infrastructure risk exposure

# --- 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 ---
hazard = "flood"
intensity_threshold = 0.5
return_period = 100
climate_scenario = "current"
aggregation_level = 2
population_source = "worldpop_100m"
age_group = "all"
infrastructure_layers = ["schools", "hospitals"]

# Parameters
run_id = "cac7f4a1-35ad-4e4f-87c7-018d2ee74051"
inputs_dir = "/app/data/work/cac7f4a1-35ad-4e4f-87c7-018d2ee74051/inputs"
output_dir = "/app/data/work/cac7f4a1-35ad-4e4f-87c7-018d2ee74051/outputs"
area = {"id": 20, "code": "NPL.5.3_1", "name": "Lumbini", "level": 2, "country": "NPL", "country_name": "Nepal", "bbox": [82.6981, 27.33, 84.4328, 28.2706]}
area_geojson = "/app/data/work/cac7f4a1-35ad-4e4f-87c7-018d2ee74051/inputs/area.geojson"
subdivisions_geojson = "/app/data/work/cac7f4a1-35ad-4e4f-87c7-018d2ee74051/inputs/subdivisions.geojson"
event = {}
event_json = ""
montandon_stac_url = "https://montandon-eoapi-stage.ifrc.org/stac"
hazard = "flood"
age_group = "all"
return_period = 100
climate_scenario = "current"
aggregation_level = 2
population_source = "worldpop_100m"
intensity_threshold = 0.5
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()

# ---------------------------------------------------------------------------------------------
# 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((area.get("code", "DEMO"), hazard, return_period, climate_scenario))) % 10_000
X, Y, step = make_grid(area_gdf.total_bounds)
intensity = hazard_grid(X, Y, seed)
if climate_scenario != "current":
    intensity = np.clip(intensity * (1.15 if climate_scenario == "rcp45" else 1.3), 0, 1)
intensity *= 1 + math.log10(return_period / 100) * 0.2 if return_period else 1
mask = (intensity >= intensity_threshold).astype(int)
print(f"{hazard} · threshold {intensity_threshold} · {return_period}-yr · {climate_scenario}: {mask.mean():.1%} of grid cells exposed")

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
exposed = population * mask
cells = cells_to_gdf(X, Y, step, intensity=intensity, mask=mask, population=population, exposed=exposed)
cells = cells[cells.intersects(area_gdf.union_all())]
print(f"Population ({population_source}, {age_group}): {population.sum():,.0f} · exposed: {exposed.sum():,.0f}")

stats = zonal_sum(cells, subdivisions, ["population", "exposed"])
stats["exposed_share"] = (stats["exposed"] / stats["population"].replace(0, np.nan)).fillna(0)

# Facilities: synthetic points per subdivision (replace with OSM / HOT exports)
rng = np.random.default_rng(seed)
fac_rows = []
for layer in infrastructure_layers:
    per_zone = {"schools": 25, "hospitals": 4, "airports": 1, "roads": 10}.get(layer, 5)
    for idx, geom in subdivisions.geometry.items():
        minx, miny, maxx, maxy = geom.bounds
        n = 0
        while n < per_zone:
            p = shapely.Point(rng.uniform(minx, maxx), rng.uniform(miny, maxy))
            if geom.contains(p):
                fac_rows.append({"layer": layer, "zone": idx, "geometry": p}); n += 1
facilities = gpd.GeoDataFrame(fac_rows, crs=4326) if fac_rows else gpd.GeoDataFrame(columns=["layer", "zone", "geometry"], crs=4326)
if len(facilities):
    hit = gpd.sjoin(facilities, cells[cells["mask"] == 1][["geometry"]], how="left", predicate="within")
    facilities["exposed"] = hit["index_right"].notna().groupby(level=0).any().reindex(facilities.index).fillna(False).astype(bool)
    for layer in infrastructure_layers:
        f = facilities[facilities.layer == layer]
        stats[f"{layer}_total"] = f.groupby("zone").size().reindex(stats.index).fillna(0).astype(int)
        stats[f"{layer}_exposed"] = f[f.exposed].groupby("zone").size().reindex(stats.index).fillna(0).astype(int)

cols = ["name", "population", "exposed", "exposed_share"] + [c for c in stats.columns if c.endswith("_exposed") or c.endswith("_total")]
table = stats[cols].sort_values("exposed", ascending=False)
table["population"] = table["population"].round(0); table["exposed"] = table["exposed"].round(0)
table.to_csv(Path(output_dir) / "exposure_by_subdivision.csv", index=False)
stats.to_file(Path(output_dir) / "exposure_by_subdivision.geojson", driver="GeoJSON")
table.style.format({"population": "{:,.0f}", "exposed": "{:,.0f}", "exposed_share": "{:.1%}"}).background_gradient(subset=["exposed_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}"

n_fac = int(sum(stats[c].sum() for c in stats.columns if c.endswith("_exposed"))) if infrastructure_layers else 0
Row(
    Stat(label=f"Population ({age_group})", value=fmt_int(stats["population"].sum())),
    Stat(label="Exposed population", value=fmt_int(stats["exposed"].sum())),
    Stat(label="Share exposed", value=f"{stats['exposed'].sum() / max(stats['population'].sum(), 1):.1%}"),
    Stat(label="Facilities exposed", value=n_fac),
    gap="16px",
)

vmax = max(float(stats["exposed_share"].max()), 0.01)
zones_layer = PolygonLayer.from_geopandas(stats[["name", "population", "exposed", "exposed_share", "geometry"]], get_fill_color=colors_for(stats["exposed_share"], "OrRd", 0, vmax), get_line_color=[60, 60, 60], line_width_min_pixels=1, opacity=0.75)
layers = [zones_layer]
fac_layer = None
if len(facilities) and facilities.exposed.any():
    fx = facilities[facilities.exposed]
    fac_layer = ScatterplotLayer.from_geopandas(fx[["layer", "geometry"]], get_fill_color=[31, 36, 48], radius_min_pixels=3)
    layers.append(fac_layer)
Map(layers)

toggles = [LayerToggle(layer=zones_layer, label="Exposed share")] + ([LayerToggle(layer=fac_layer, label="Exposed facilities")] if fac_layer is not None else [])
Row(Column(*toggles, gap="6px"), legend_for("OrRd", 0, vmax, title="Exposed population share"), gap="24px", align="start")

bar_chart(table["name"], table["exposed"], title=f"{area_name} · {hazard} ≥ {intensity_threshold} · {return_period}-yr · {climate_scenario}", y_label=f"Exposed population ({age_group})")