The finding under test (AIDA): in Logroño, the most socially vulnerable residents at pluvial flood risk have less land that could be greened near their homes than the least vulnerable.
For every map with the current green cover (3 climates × 7 return periods) and both risks (well-being, mobility), among the residents at risk:
resident-weighted mean share of land that could be greened within 100 m, most vulnerable fifth (Q5) against least vulnerable fifth (Q1);
a one-sided permutation test (buildings shuffled between Q1 and Q5, 10,000 times) of “Q5 has less”;
the same at 50 m and 200 m, and with the other land covers, as robustness;
a cross-check against the analysis release, which computed the same measure for the planned scenario with ponds from its own snapshot of the Cadastre.
Decision rule, fixed before running: the finding holds if Q5 has less land that could be greened than Q1 in every current-cover comparison at 100 m, and no comparison shows significantly more (p < 0.05 in the other direction). Many comparisons are made and they are not independent, so an isolated small p-value is not read as evidence.
import json
from pathlib import Path
import numpy as np
import pandas as pd
CLEAN = Path("../data/clean")
RESULTS = Path("../results")
RESULTS.mkdir(parents=True, exist_ok=True)
sources = json.loads(Path("../data/raw/sources.json").read_text())
RESULTS_A = Path(sources["analysis"]["path"]) / "results"
RADII = [50, 100, 200]
RISKS = ["PFRWB", "PFRMA"]
RNG = np.random.default_rng(20261011)
N_PERM = 10_000
room = pd.read_csv(CLEAN / "room_to_green.csv").set_index("ID")
buildings = pd.read_csv(CLEAN / "buildings.csv")
buildings = buildings[buildings["ID"].isin(room.index)]def compare(ids1: pd.Index, ids5: pd.Index, weights: pd.Series, share: pd.Series) -> dict:
"""Resident-weighted mean share, Q5 minus Q1, and one-sided permutation p-values."""
pooled = ids1.append(ids5)
w = weights.reindex(pooled).to_numpy()
ws = w * share.reindex(pooled).to_numpy()
in5 = np.argsort(RNG.random((N_PERM, len(pooled))), axis=1) < len(ids5)
w5, ws5 = in5 @ w, in5 @ ws
null = ws5 / w5 - (ws.sum() - ws5) / (w.sum() - w5)
m1 = float(np.average(share.reindex(ids1), weights=weights.reindex(ids1)))
m5 = float(np.average(share.reindex(ids5), weights=weights.reindex(ids5)))
d = m5 - m1
return {"buildings_q1": len(ids1), "buildings_q5": len(ids5),
"residents_q1": float(weights.reindex(ids1).sum()), "residents_q5": float(weights.reindex(ids5).sum()),
"mean_share_q1": m1, "mean_share_q5": m5, "difference": d,
"p_q5_less": float((1 + (null <= d).sum()) / (1 + N_PERM)),
"p_q5_more": float((1 + (null >= d).sum()) / (1 + N_PERM))}
rows = []
for (lc, cl, rp), g in buildings.groupby(["land_cover", "climate", "rp"]):
g = g.set_index("ID")
for r in RISKS:
at = g[g[f"at_risk_{r}"]]
ids1 = at.index[at["vulnerability_group"] == 1]
ids5 = at.index[at["vulnerability_group"] == 5]
if len(ids1) == 0 or len(ids5) == 0:
continue
for rad in RADII:
rows.append({"land_cover": lc, "climate": cl, "rp": rp, "risk": r, "radius_m": rad,
**compare(ids1, ids5, at["R"], room[f"greenable_share_{rad}m"])})
comp = pd.DataFrame(rows)
comp.to_csv(RESULTS / "room_comparisons.csv", index=False)
main = comp[(comp["land_cover"] == "Baseline") & (comp["radius_m"] == 100)]
print(main[["climate", "rp", "risk", "residents_q1", "residents_q5", "mean_share_q1", "mean_share_q5",
"p_q5_less", "p_q5_more"]].round(3).to_string(index=False))climate rp risk residents_q1 residents_q5 mean_share_q1 mean_share_q5 p_q5_less p_q5_more
CC_Q50 2 PFRWB 6118.803 2461.942 0.189 0.055 0.0 1.0
CC_Q50 2 PFRMA 9319.513 4432.994 0.192 0.045 0.0 1.0
CC_Q50 5 PFRWB 7554.706 2955.211 0.195 0.056 0.0 1.0
CC_Q50 5 PFRMA 11154.467 5087.823 0.186 0.045 0.0 1.0
CC_Q50 10 PFRWB 7848.324 3008.609 0.192 0.055 0.0 1.0
CC_Q50 10 PFRMA 11644.268 5506.585 0.181 0.042 0.0 1.0
CC_Q50 25 PFRWB 8765.185 3733.093 0.190 0.046 0.0 1.0
CC_Q50 25 PFRMA 12670.360 5872.460 0.173 0.039 0.0 1.0
CC_Q50 50 PFRWB 9049.899 3733.093 0.192 0.046 0.0 1.0
CC_Q50 50 PFRMA 12946.087 5887.253 0.174 0.039 0.0 1.0
CC_Q50 100 PFRWB 9301.528 3816.078 0.187 0.045 0.0 1.0
CC_Q50 100 PFRMA 13186.019 5951.379 0.174 0.039 0.0 1.0
CC_Q50 200 PFRWB 9323.835 3858.342 0.187 0.045 0.0 1.0
CC_Q50 200 PFRMA 13256.754 6047.568 0.173 0.038 0.0 1.0
CC_Q90 2 PFRWB 7038.883 2566.710 0.186 0.057 0.0 1.0
CC_Q90 2 PFRMA 10120.010 4543.288 0.195 0.047 0.0 1.0
CC_Q90 5 PFRWB 7639.320 3008.609 0.194 0.055 0.0 1.0
CC_Q90 5 PFRMA 11246.360 5506.585 0.185 0.042 0.0 1.0
CC_Q90 10 PFRWB 8064.492 3486.146 0.195 0.049 0.0 1.0
CC_Q90 10 PFRMA 11872.946 5540.727 0.181 0.041 0.0 1.0
CC_Q90 25 PFRWB 9049.899 3733.093 0.192 0.046 0.0 1.0
CC_Q90 25 PFRMA 12946.087 5887.253 0.174 0.039 0.0 1.0
CC_Q90 50 PFRWB 9301.528 3816.078 0.187 0.045 0.0 1.0
CC_Q90 50 PFRMA 13232.241 5951.379 0.173 0.039 0.0 1.0
CC_Q90 100 PFRWB 9438.684 3901.302 0.186 0.044 0.0 1.0
CC_Q90 100 PFRMA 13256.754 6047.568 0.173 0.038 0.0 1.0
CC_Q90 200 PFRWB 9899.978 3906.562 0.185 0.044 0.0 1.0
CC_Q90 200 PFRMA 13596.719 6050.198 0.172 0.038 0.0 1.0
present 2 PFRWB 5664.806 2375.818 0.173 0.047 0.0 1.0
present 2 PFRMA 8594.264 4341.610 0.189 0.041 0.0 1.0
present 5 PFRWB 7440.555 2746.470 0.196 0.057 0.0 1.0
present 5 PFRMA 10830.012 4789.906 0.190 0.047 0.0 1.0
present 10 PFRWB 7639.320 3008.609 0.194 0.055 0.0 1.0
present 10 PFRMA 11246.360 5487.806 0.185 0.042 0.0 1.0
present 25 PFRWB 8064.492 3436.342 0.195 0.050 0.0 1.0
present 25 PFRMA 11823.936 5506.585 0.181 0.042 0.0 1.0
present 50 PFRWB 8785.032 3733.093 0.190 0.046 0.0 1.0
present 50 PFRMA 12690.207 5872.460 0.173 0.039 0.0 1.0
present 100 PFRWB 9049.899 3733.093 0.192 0.046 0.0 1.0
present 100 PFRMA 12946.087 5887.253 0.174 0.039 0.0 1.0
present 200 PFRWB 9301.528 3816.078 0.187 0.045 0.0 1.0
present 200 PFRMA 13186.019 5945.243 0.174 0.039 0.0 1.0
Robustness: radius and land cover¶
robust = (comp.assign(q5_less=comp["difference"] < 0, sig_less=comp["p_q5_less"] < 0.05,
sig_more=comp["p_q5_more"] < 0.05)
.groupby(["land_cover", "radius_m"])
.agg(comparisons=("q5_less", "size"), q5_less=("q5_less", "sum"),
significant_less=("sig_less", "sum"), significant_more=("sig_more", "sum"),
mean_share_q1=("mean_share_q1", "mean"), mean_share_q5=("mean_share_q5", "mean"))
.reset_index())
robust.to_csv(RESULTS / "room_robustness.csv", index=False)
print(robust.round(3).to_string(index=False)) land_cover radius_m comparisons q5_less significant_less significant_more mean_share_q1 mean_share_q5
Allgreen 50 40 40 40 0 0.127 0.030
Allgreen 100 40 40 40 0 0.187 0.044
Allgreen 200 40 40 40 0 0.256 0.058
Baseline 50 42 42 42 0 0.131 0.030
Baseline 100 42 42 42 0 0.185 0.045
Baseline 200 42 42 42 0 0.248 0.060
NBSplanned 50 40 40 40 0 0.131 0.030
NBSplanned 100 40 40 40 0 0.184 0.045
NBSplanned 200 40 40 40 0 0.248 0.060
NBSplanned_plus 50 42 42 42 0 0.132 0.030
NBSplanned_plus 100 42 42 42 0 0.185 0.045
NBSplanned_plus 200 42 42 42 0 0.248 0.060
Cross-check against the analysis release (planned scenario with ponds, 100 m)¶
The release’s notebook 06 measured the same share from its own copy of the Cadastre,
with at-risk sets from the counts variant; the at-risk sets are the same, so the means
differ only through Cadastre updates between the two downloads.
ref = pd.read_csv(RESULTS_A / "where_vulnerable.csv")
ref = ref.pivot_table(index=["climate", "rp", "risk"], columns="group", values="greenable_share")
ref.columns = [f"release_q{c}" for c in ref.columns]
ours = comp[(comp["land_cover"] == "NBSplanned_plus") & (comp["radius_m"] == 100)].set_index(["climate", "rp", "risk"])
check = ours[["mean_share_q1", "mean_share_q5"]].join(ref, how="inner")
check["abs_diff_q1"] = (check["mean_share_q1"] - check["release_q1"]).abs()
check["abs_diff_q5"] = (check["mean_share_q5"] - check["release_q5"]).abs()
check.reset_index().to_csv(RESULTS / "crosscheck_release.csv", index=False)
print(f"{len(check)} comparisons; largest absolute difference "
f"{check[['abs_diff_q1', 'abs_diff_q5']].max().max():.4f}")40 comparisons; largest absolute difference 0.0000
Context, not part of the finding: does room to green predict who greening takes out of risk?¶
Within each fifth, residents at risk with the current green cover are split into those the all-green scenario takes out of risk (neither risk index above zero) and those who remain; the resident-weighted mean share at 100 m of the two sets is compared. If more room to green went with leaving risk, the taken-out set would have more in most cases.
def at_any(g: pd.DataFrame) -> pd.Series:
return g["at_risk_PFRWB"] | g["at_risk_PFRMA"]
share = room["greenable_share_100m"]
link = []
for (cl, rp), g in buildings[buildings["land_cover"].isin(["Baseline", "Allgreen"])].groupby(["climate", "rp"]):
cur = g[g["land_cover"] == "Baseline"].set_index("ID")
grn = g[g["land_cover"] == "Allgreen"].set_index("ID").reindex(cur.index)
if grn["R"].isna().all():
continue
out = at_any(cur) & ~at_any(grn).fillna(True).astype(bool)
for k in [1, 5]:
m = at_any(cur) & (cur["vulnerability_group"] == k)
o, s_ = cur.index[m & out], cur.index[m & ~out]
if len(o) == 0 or len(s_) == 0:
continue
link.append({"climate": cl, "rp": rp, "group": k,
"share_taken_out": float(np.average(share.reindex(o), weights=cur.loc[o, "R"])),
"share_remaining": float(np.average(share.reindex(s_), weights=cur.loc[s_, "R"]))})
link = pd.DataFrame(link)
link["taken_out_has_more"] = link["share_taken_out"] > link["share_remaining"]
link.to_csv(RESULTS / "room_vs_taken_out.csv", index=False)
link_summary = link.groupby("group").agg(cases=("taken_out_has_more", "size"),
taken_out_has_more=("taken_out_has_more", "sum"),
share_taken_out=("share_taken_out", "mean"),
share_remaining=("share_remaining", "mean"))
print(link_summary.round(3).to_string()) cases taken_out_has_more share_taken_out share_remaining
group
1 20 0 0.144 0.188
5 20 9 0.063 0.040
Does the finding hold?¶
n_main = len(main)
claim = {
"aida": "In Logroño, the most socially vulnerable residents at pluvial flood risk have less land that "
"could be greened near their homes than the least vulnerable.",
"comparisons_current_cover_100m": n_main,
"q5_less_in": int((main["difference"] < 0).sum()),
"significant_q5_less_p_lt_0.05": int((main["p_q5_less"] < 0.05).sum()),
"significant_q5_more_p_lt_0.05": int((main["p_q5_more"] < 0.05).sum()),
"mean_share_q1_range": [float(main["mean_share_q1"].min()), float(main["mean_share_q1"].max())],
"mean_share_q5_range": [float(main["mean_share_q5"].min()), float(main["mean_share_q5"].max())],
"all_comparisons_all_radii_all_covers": int(len(comp)),
"q5_less_in_all": int((comp["difference"] < 0).sum()),
"buildings_matched_in_cadastre": int(len(room)),
"crosscheck_max_abs_diff": float(check[["abs_diff_q1", "abs_diff_q5"]].max().max()),
"context_room_vs_taken_out": {f"q{k}": {"cases": int(r.cases), "taken_out_has_more": int(r.taken_out_has_more)}
for k, r in link_summary.iterrows()},
}
claim["holds"] = (claim["q5_less_in"] == n_main and claim["significant_q5_more_p_lt_0.05"] == 0)
(RESULTS / "claim_test.json").write_text(json.dumps(claim, indent=2, ensure_ascii=False))
claim{'aida': 'In Logroño, the most socially vulnerable residents at pluvial flood risk have less land that could be greened near their homes than the least vulnerable.',
'comparisons_current_cover_100m': 42,
'q5_less_in': 42,
'significant_q5_less_p_lt_0.05': 42,
'significant_q5_more_p_lt_0.05': 0,
'mean_share_q1_range': [0.17218577632828036, 0.19580073516944715],
'mean_share_q5_range': [0.03798285840831739, 0.05745331553164259],
'all_comparisons_all_radii_all_covers': 492,
'q5_less_in_all': 492,
'buildings_matched_in_cadastre': 4540,
'crosscheck_max_abs_diff': 9.020562075079397e-17,
'context_room_vs_taken_out': {'q1': {'cases': 20, 'taken_out_has_more': 0},
'q5': {'cases': 20, 'taken_out_has_more': 9}},
'holds': True}