1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
"""How good are the models? Everything is scored on data the model did not see.
1. Risk model on the full grid of each held-out AOI (real flood share, not the balanced training
sample): ROC AUC, average precision, and precision / recall when flagging the top 5 % / 10 % of
the area. Compared with the EFAS RP100 map and with the CEMS-river variant.
2. Building damage: do depth (EFAS + DEM check) and terrain risk separate damaged from undamaged
buildings in the CEMS grading maps (builtUpP damage_gra)?
3. Collects the earlier checks (roads, DEM depth check, in-between model) into one table.
Writes outputs/eval/summary.csv and per-part CSVs.
"""
import argparse
import geopandas as gpd, numpy as np, pandas as pd, shapely
from sklearn.metrics import average_precision_score, roc_auc_score
import risk_model
from main import DATA
from paths import CACHE, OUTPUTS
from risk_model import FEATURES, Area, new_model
from road_model import sample_raster
OUT = OUTPUTS / "eval"
DAMAGED = {"Damaged", "Destroyed", "Possibly damaged"}
def top_share(score, truth, share):
"""Precision and recall when the highest-scoring `share` of cells is flagged."""
k = max(1, int(len(score) * share))
flag = np.zeros(len(score), bool)
flag[np.argsort(-score)[:k]] = True
return (truth & flag).sum() / k, (truth & flag).sum() / max(truth.sum(), 1)
def grid_scores(name, score, truth):
row = dict(area=name, flooded_pct=100 * truth.mean(), auc=roc_auc_score(truth, score),
avg_precision=average_precision_score(truth, score), baseline_precision=truth.mean())
for s in (0.05, 0.10):
row[f"precision_top{int(s * 100)}"], row[f"recall_top{int(s * 100)}"] = top_share(score, truth, s)
return row
def risk_eval(rivers, features=FEATURES):
df, grids = pd.read_pickle(CACHE / f"risk_features_{rivers}.pkl")
rows = []
for name, g in grids.items():
truth = np.asarray(g["flooded"], bool)
if truth.sum() < 10 or (~truth).sum() < 10:
continue
model = new_model().fit(df.loc[df.area != name, features], df.loc[df.area != name, "flooded"])
rows.append(dict(model=f"risk ({rivers} rivers)", **grid_scores(name, model.predict_proba(g["feats"][features])[:, 1], truth)))
rows.append(dict(model="EFAS RP100 depth", **grid_scores(name, g["feats"]["rp100"].to_numpy(), truth)))
return pd.DataFrame(rows)
def damage_eval():
"""CEMS builtUpP only lists damaged buildings, so two tests:
(a) hit: graded buildings vs GHS-OBAT buildings within 1 km of a graded one but > 20 m from any;
(b) severity: damaged / destroyed vs possibly damaged, among graded buildings."""
from access_model import buildings_in
from scipy.spatial import cKDTree
rows = []
for folder in sorted(DATA.glob("ems_*/*_GRA_PRODUCT_*_vector")):
path = next(iter(folder.glob("*_builtUpP_*.shp")), None)
if path is None:
continue
g = gpd.read_file(path).to_crs(3035)
g = g[g.damage_gra.isin(DAMAGED)]
area = Area(folder)
gxy = shapely.get_coordinates(g.geometry.representative_point().values)
_, oxy = buildings_in(area.aoi)
d = cKDTree(gxy).query(oxy)[0]
nxy = oxy[(d > 20) & (d <= 1000)]
tif = risk_model.OUT / f"risk_{area.name}.tif"
def scores(xy):
out = {"EFAS RP100 depth (DEM-checked)": area.depth("rp100", xy),
"EFAS RP100 depth (raw)": area.depth("rp100", xy, dem_check=False),
"inside observed 2021 extent": shapely.contains_xy(area.flood, xy[:, 0], xy[:, 1]).astype(float)}
if tif.exists():
out["terrain risk (held-out map)"] = np.nan_to_num(sample_raster(tif, xy))
return out
tests = {"hit": (np.vstack([gxy, nxy]), np.r_[np.ones(len(gxy)), np.zeros(len(nxy))].astype(bool)),
"severity": (gxy, g.damage_gra.isin({"Damaged", "Destroyed"}).to_numpy())}
for test, (xy, y) in tests.items():
if y.sum() < 5 or (~y).sum() < 5:
continue
for k, sc in scores(xy).items():
rows.append(dict(test=test, area=area.name, score=k, n=len(y), positive_pct=100 * y.mean(), auc=roc_auc_score(y, sc)))
print(f"{area.name}: {len(gxy)} graded, {len(nxy)} ungraded nearby", flush=True)
del area
return pd.DataFrame(rows)
if __name__ == "__main__":
p = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
p.parse_args()
OUT.mkdir(parents=True, exist_ok=True)
pd.set_option("display.width", 200)
risk = pd.concat([risk_eval(r) for r in ("osm", "cems") if (CACHE / f"risk_features_{r}.pkl").exists()])
risk = risk.drop_duplicates(["model", "area"])
risk.round(4).to_csv(OUT / "risk_grid.csv", index=False)
med = risk.groupby("model")[["auc", "avg_precision", "baseline_precision", "precision_top5", "recall_top5",
"precision_top10", "recall_top10"]].median()
print("\n1. Flood extent on the full held-out AOI grid (medians over AOIs)\n" + med.round(3).to_string())
dmg = damage_eval()
dmg.round(3).to_csv(OUT / "building_damage.csv", index=False)
print("\n2. Building damage, ROC AUC (hit: graded vs ungraded nearby; severity: damaged/destroyed vs possibly)\n"
+ dmg.pivot_table(index=["test", "area"], columns="score", values="auc").round(3).to_string())
summary = [dict(model=m, check="flood extent, held-out AOI grid", metric="median ROC AUC / avg precision",
value=f"{r.auc:.3f} / {r.avg_precision:.3f} (base rate {r.baseline_precision:.3f})") for m, r in med.iterrows()]
summary += [dict(model=sc, check=f"building damage, {t} (CEMS grades)", metric="median ROC AUC", value=f"{v:.3f}")
for (t, sc), v in dmg.groupby(["test", "score"]).auc.median().items()]
for f, check in [("roads/damage_validation.csv", "road damage (CEMS grades)"), ("eval/depth_check.csv", "DEM depth check")]:
if (OUTPUTS / f).exists():
summary.append(dict(model=f.split("/")[0], check=check, metric="see file", value=str(OUTPUTS / f)))
summary.append(dict(model="interp_model", check="predict AOI15 18/07 from 16 + 20/07", metric="IoU",
value="0.14 (shape morphing 0.007, nearest image 0.02)"))
pd.DataFrame(summary).to_csv(OUT / "summary.csv", index=False)
print(f"\nWrote {OUT}/summary.csv, risk_grid.csv, building_damage.csv")
|
