"""Reproducible scenario analysis for SHWRP-2026-016. The program screens the June 2026 national sidewalk inventory for Keelung and links district-level March 2026 population data. It does not identify actual planting sites, existing trees, utilities, land ownership, wind exposure, soil volume, or measured cooling effects. All tree counts are planning scenarios. """ from __future__ import annotations import csv import io import json import math from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd ROOT = Path(__file__).resolve().parent REPO = ROOT.parent.parent TMP = REPO / "tmp" / "research016" OUT = ROOT / "outputs" FIG = OUT / "figures" OUT.mkdir(parents=True, exist_ok=True) FIG.mkdir(parents=True, exist_ok=True) REPORT = "SHWRP-2026-016" DATA_DATE_SIDEWALK = "2026-06" DATA_DATE_POPULATION = "2026-03" MIN_CLEAR = 1.5 PLANTING_STRIP = 1.0 MODELED_SPACING = 8.0 def find_keelung_geojson() -> Path: candidates = list((TMP / "sidewalk").glob("*WGS84.geojson")) # The downloaded ZIP may expose mojibake file names on some Windows tools; # identify the file from its official county attribute instead of its name. for path in candidates: if path.stat().st_size > 5_000_000: continue data = json.loads(path.read_text(encoding="utf-8")) features = data.get("features", []) if features and features[0].get("properties", {}).get("COUNTY_NA") == "基隆市": return path raise FileNotFoundError("Keelung WGS84 sidewalk GeoJSON was not found") def load_sidewalks() -> tuple[pd.DataFrame, list[dict]]: path = find_keelung_geojson() geo = json.loads(path.read_text(encoding="utf-8")) rows = [] for idx, feature in enumerate(geo["features"], start=1): p = feature["properties"] total = pd.to_numeric(p.get("SW_WTH"), errors="coerce") clear = pd.to_numeric(p.get("SWW_WTH"), errors="coerce") length = pd.to_numeric(p.get("SW_LENG"), errors="coerce") facility = total - clear if pd.notna(total) and pd.notna(clear) else np.nan valid = pd.notna(length) and length > 0 and pd.notna(total) and total > 0 and pd.notna(clear) and clear >= 0 and clear <= total length_anomaly = bool(valid and length > 5000) if length_anomaly: scenario = "非植樹優先／資料待查" elif valid and clear >= MIN_CLEAR and facility >= PLANTING_STRIP: scenario = "直接容納候選" elif valid and total >= MIN_CLEAR + PLANTING_STRIP: scenario = "需道路空間重分配" else: scenario = "非植樹優先／資料待查" rows.append({ "segment_id": f"KLSW-{idx:04d}", "district": p.get("VILL_NAME") or "未標示", "road_name": p.get("NAME") or "", "start": p.get("PSTART") or "", "end": p.get("PEND") or "", "length_m": float(length) if pd.notna(length) else np.nan, "length_quality_flag": "待查核異常值：單筆長度超過5公里" if length_anomaly else "通過研究規則", "total_width_m": float(total) if pd.notna(total) else np.nan, "clear_width_m": float(clear) if pd.notna(clear) else np.nan, "facility_strip_proxy_m": float(facility) if pd.notna(facility) else np.nan, "ramps_recorded": p.get("SW_RAMP"), "screening_class": scenario, "modeled_tree_positions_8m": math.floor(float(length) / MODELED_SPACING) if scenario != "非植樹優先／資料待查" and pd.notna(length) else 0, "source_date": DATA_DATE_SIDEWALK, }) return pd.DataFrame(rows), geo["features"] def load_population() -> pd.DataFrame: path = TMP / "keelung_population_202603.csv" rows = list(csv.reader(io.StringIO(path.read_text(encoding="cp950")))) header = rows[2] output = [] # The official sheet lists total/male/female for city and each district. for row in rows[3::3]: if len(row) < len(header): continue area = row[0].strip() if not area or area == "基隆市": continue total = int(row[2]) age65 = sum(int(value or 0) for value in row[68:104]) age0_14 = sum(int(value or 0) for value in row[3:18]) output.append({ "district": area.replace("　", "").strip(), "population": total, "age_65_plus": age65, "share_65_plus": age65 / total, "age_0_14": age0_14, "share_0_14": age0_14 / total, "source_date": DATA_DATE_POPULATION, }) return pd.DataFrame(output) def sensitivity(sidewalks: pd.DataFrame) -> pd.DataFrame: rows = [] for clear in [1.2, 1.5, 1.8]: for strip in [0.8, 1.0, 1.2]: valid = ( (sidewalks["clear_width_m"] >= clear) & (sidewalks["facility_strip_proxy_m"] >= strip) & (sidewalks["length_m"] > 0) & (sidewalks["length_quality_flag"] == "通過研究規則") ) candidate_length = sidewalks.loc[valid, "length_m"].sum() rows.append({ "minimum_clear_width_m": clear, "modeled_planting_strip_m": strip, "segments": int(valid.sum()), "candidate_length_km": candidate_length / 1000, "tree_positions_6m": int(np.floor(sidewalks.loc[valid, "length_m"] / 6).sum()), "tree_positions_8m": int(np.floor(sidewalks.loc[valid, "length_m"] / 8).sum()), "tree_positions_10m": int(np.floor(sidewalks.loc[valid, "length_m"] / 10).sum()), }) return pd.DataFrame(rows) def district_summary(sidewalks: pd.DataFrame, population: pd.DataFrame) -> pd.DataFrame: group = sidewalks.groupby(["district", "screening_class"], dropna=False).agg( segments=("segment_id", "count"), length_m=("length_m", "sum"), tree_positions_8m=("modeled_tree_positions_8m", "sum"), ).reset_index() pivot = group.pivot(index="district", columns="screening_class", values=["segments", "length_m", "tree_positions_8m"]).fillna(0) pivot.columns = [f"{a}_{b}" for a, b in pivot.columns] pivot = pivot.reset_index() total = sidewalks.groupby("district").agg(recorded_segments=("segment_id", "count"), recorded_length_m=("length_m", "sum")).reset_index() result = total.merge(pivot, on="district", how="left").merge(population, on="district", how="left") direct = result.get("length_m_直接容納候選", pd.Series(0, index=result.index)) redesign = result.get("length_m_需道路空間重分配", pd.Series(0, index=result.index)) result["screened_candidate_length_m"] = direct + redesign result["candidate_length_per_1000_people_m"] = result["screened_candidate_length_m"] / result["population"] * 1000 # Coarse equity priority proxy: elderly residents per km of screened candidate. result["elderly_pressure_proxy"] = result["age_65_plus"] / (result["screened_candidate_length_m"] / 1000).replace(0, np.nan) return result.sort_values("district").reset_index(drop=True) def scenario_summary(sidewalks: pd.DataFrame) -> pd.DataFrame: records = [] for label, predicate in [ ("S0 官方人行道紀錄（含待查核值）", pd.Series(True, index=sidewalks.index)), ("S1 直接容納候選", sidewalks["screening_class"] == "直接容納候選"), ("S2 加入道路重分配候選", sidewalks["screening_class"] != "非植樹優先／資料待查"), ]: subset = sidewalks[predicate] records.append({ "scenario": label, "segments": len(subset), "length_km": subset["length_m"].sum() / 1000, "modeled_tree_positions_6m": int(np.floor(subset["length_m"] / 6).sum()) if not label.startswith("S0") else np.nan, "modeled_tree_positions_8m": int(np.floor(subset["length_m"] / 8).sum()) if not label.startswith("S0") else np.nan, "modeled_tree_positions_10m": int(np.floor(subset["length_m"] / 10).sum()) if not label.startswith("S0") else np.nan, }) return pd.DataFrame(records) def setup_font() -> None: font_file = REPO / "public" / "fonts" / "NotoSansTC-Regular.ttf" if font_file.exists(): from matplotlib import font_manager font_manager.fontManager.addfont(str(font_file)) plt.rcParams["font.family"] = font_manager.FontProperties(fname=str(font_file)).get_name() else: plt.rcParams["font.family"] = ["Microsoft JhengHei", "DejaVu Sans"] plt.rcParams["axes.unicode_minus"] = False plt.rcParams["figure.facecolor"] = "#fbf7f0" def make_figures(sidewalks: pd.DataFrame, features: list[dict], districts: pd.DataFrame, sens: pd.DataFrame, scenarios: pd.DataFrame) -> None: setup_font() brown, green, gold, gray, red = "#6f3519", "#315c4c", "#c28b42", "#9a938c", "#a34a3a" colors = {"直接容納候選": green, "需道路空間重分配": gold, "非植樹優先／資料待查": gray} fig, ax = plt.subplots(figsize=(10.8, 8.2)) for row, feature in zip(sidewalks.itertuples(), features): geom = feature.get("geometry") or {} polygons = geom.get("coordinates", []) for poly in polygons: if not poly or not poly[0]: continue xy = np.asarray(poly[0]) ax.plot(xy[:, 0], xy[:, 1], color=colors[row.screening_class], linewidth=.7, alpha=.62) for label, color in colors.items(): ax.plot([], [], color=color, linewidth=4, label=label) ax.set_title("圖1　基隆市人行道紀錄的幾何初篩", fontsize=17, color=brown, weight="bold") ax.set_xlabel("經度"); ax.set_ylabel("緯度"); ax.legend(frameon=False, loc="best") ax.text(.01, -.11, "僅依2026年6月人行道寬度欄位分類；不是樹位清冊，也未查核地下管線、產權、坡度與風場。", transform=ax.transAxes, fontsize=9, color="#625b55") fig.tight_layout(); fig.savefig(FIG / "fig1_geometric_screening_map.png", dpi=220, bbox_inches="tight"); plt.close(fig) fig, ax = plt.subplots(figsize=(11, 6.4)) vals = sidewalks[["total_width_m", "clear_width_m", "facility_strip_proxy_m"]].rename(columns={"total_width_m":"人行道總寬", "clear_width_m":"淨寬", "facility_strip_proxy_m":"設施帶代理寬度"}) ax.hist([vals[c].clip(upper=8).dropna() for c in vals], bins=np.arange(0, 8.25, .25), label=list(vals.columns), color=[gold, green, brown], alpha=.68) ax.axvline(MIN_CLEAR, color=red, linestyle="--", label="研究基準：1.5 m淨寬") ax.set_title("圖2　人行道寬度欄位分布（8公尺以上截尾顯示）", fontsize=17, color=brown, weight="bold") ax.set_xlabel("公尺"); ax.set_ylabel("線段筆數"); ax.legend(frameon=False); ax.grid(axis="y", alpha=.2) fig.tight_layout(); fig.savefig(FIG / "fig2_width_distribution.png", dpi=220); plt.close(fig) cats = ["直接容納候選", "需道路空間重分配", "非植樹優先／資料待查"] data = sidewalks.groupby(["district", "screening_class"])["length_m"].sum().unstack(fill_value=0) / 1000 data = data.reindex(columns=cats, fill_value=0).sort_index() fig, ax = plt.subplots(figsize=(11.5, 6.8)) data.plot(kind="bar", stacked=True, color=[colors[c] for c in cats], ax=ax) ax.set_title("圖3　各行政區人行道紀錄長度與情境分類", fontsize=17, color=brown, weight="bold") ax.set_xlabel(""); ax.set_ylabel("公里"); ax.tick_params(axis="x", rotation=0); ax.legend(frameon=False); ax.grid(axis="y", alpha=.2) fig.tight_layout(); fig.savefig(FIG / "fig3_district_scenarios.png", dpi=220); plt.close(fig) matrix = sens.pivot(index="minimum_clear_width_m", columns="modeled_planting_strip_m", values="candidate_length_km") fig, ax = plt.subplots(figsize=(9.5, 6.6)) image = ax.imshow(matrix.values, cmap="YlGn", aspect="auto") for i in range(matrix.shape[0]): for j in range(matrix.shape[1]): ax.text(j, i, f"{matrix.iloc[i,j]:.1f} km", ha="center", va="center", fontsize=12, weight="bold") ax.set_xticks(range(matrix.shape[1]), [f"{x:.1f}" for x in matrix.columns]); ax.set_yticks(range(matrix.shape[0]), [f"{x:.1f}" for x in matrix.index]) ax.set_xlabel("假設設施／植栽帶寬度（m）"); ax.set_ylabel("最低淨寬（m）") ax.set_title("圖4　直接容納候選長度的門檻敏感度", fontsize=17, color=brown, weight="bold") fig.colorbar(image, ax=ax, label="候選長度（km）"); fig.tight_layout(); fig.savefig(FIG / "fig4_threshold_sensitivity.png", dpi=220); plt.close(fig) d = districts.sort_values("elderly_pressure_proxy", ascending=True) fig, ax = plt.subplots(figsize=(11, 6.8)) bars = ax.barh(d["district"], d["candidate_length_per_1000_people_m"], color=green, alpha=.82) ax.set_xlabel("每千人可初篩候選長度（m）"); ax.set_title("圖5　行政區人口、老年人口與候選長度的粗粒度公平檢視", fontsize=16, color=brown, weight="bold") ax2 = ax.twiny(); ax2.plot(d["share_65_plus"]*100, d["district"], "o-", color=gold, linewidth=2); ax2.set_xlabel("65歲以上人口占比（%）", color=gold) ax.grid(axis="x", alpha=.2); ax.text(.01, -.13, "僅為行政區尺度，不足以辨識里、街廓或個人的熱暴露與需求。", transform=ax.transAxes, fontsize=9, color="#625b55") fig.tight_layout(); fig.savefig(FIG / "fig5_equity_screen.png", dpi=220, bbox_inches="tight"); plt.close(fig) fig, ax = plt.subplots(figsize=(10.8, 6.5)) plot = scenarios.iloc[1:].copy() x = np.arange(len(plot)); width=.24 for offset, col, label, color in [(-width,"modeled_tree_positions_6m","6 m間距",gold),(0,"modeled_tree_positions_8m","8 m間距",green),(width,"modeled_tree_positions_10m","10 m間距",brown)]: ax.bar(x+offset, plot[col], width, label=label, color=color) ax.set_xticks(x, [s.replace("S1 ","").replace("S2 ","") for s in plot["scenario"]]); ax.set_ylabel("情境樹位數（個）") ax.set_title("圖6　不同間距下的規劃樹位情境", fontsize=17, color=brown, weight="bold"); ax.legend(frameon=False); ax.grid(axis="y", alpha=.2) ax.text(.01, -.14, "樹位數＝候選線段長度除以假設間距後向下取整；不是現有或核准植樹數。", transform=ax.transAxes, fontsize=9, color="#625b55") fig.tight_layout(); fig.savefig(FIG / "fig6_tree_position_scenarios.png", dpi=220, bbox_inches="tight"); plt.close(fig) def main() -> None: sidewalks, features = load_sidewalks() population = load_population() sens = sensitivity(sidewalks) districts = district_summary(sidewalks, population) scenarios = scenario_summary(sidewalks) sidewalks.to_csv(OUT / "keelung_sidewalk_screening.csv", index=False, encoding="utf-8-sig") population.to_csv(OUT / "district_demographics.csv", index=False, encoding="utf-8-sig") districts.to_csv(OUT / "district_summary.csv", index=False, encoding="utf-8-sig") sens.to_csv(OUT / "threshold_sensitivity.csv", index=False, encoding="utf-8-sig") scenarios.to_csv(OUT / "scenario_summary.csv", index=False, encoding="utf-8-sig") make_figures(sidewalks, features, districts, sens, scenarios) direct = sidewalks[sidewalks["screening_class"] == "直接容納候選"] combined = sidewalks[sidewalks["screening_class"] != "非植樹優先／資料待查"] summary = { "report_number": REPORT, "analysis_type": "scenario-based geometric screening", "sidewalk_data_date": DATA_DATE_SIDEWALK, "population_data_date": DATA_DATE_POPULATION, "recorded_segments": int(len(sidewalks)), "recorded_length_km_including_flagged_values": round(float(sidewalks["length_m"].sum()/1000), 3), "quality_screened_length_km": round(float(sidewalks.loc[sidewalks["length_quality_flag"] == "通過研究規則", "length_m"].sum()/1000), 3), "flagged_length_records": int((sidewalks["length_quality_flag"] != "通過研究規則").sum()), "direct_fit_segments": int(len(direct)), "direct_fit_length_km": round(float(direct["length_m"].sum()/1000), 3), "combined_candidate_segments": int(len(combined)), "combined_candidate_length_km": round(float(combined["length_m"].sum()/1000), 3), "modeled_tree_positions_8m_direct": int(np.floor(direct["length_m"]/8).sum()), "modeled_tree_positions_8m_combined": int(np.floor(combined["length_m"]/8).sum()), "districts": int(districts["district"].nunique()), "measured_cooling_effect": False, "field_audit_completed": False, "important_limit": "width-only screening; existing trees, utilities, ownership, soil, slope, wind and sightlines are unobserved", } (OUT / "analysis_summary.json").write_text(json.dumps(summary, ensure_ascii=False, indent=2), encoding="utf-8") print(json.dumps(summary, ensure_ascii=False, indent=2)) if __name__ == "__main__": main()