"""X09 河川浸水想定 × 用途地域 — 浸水深ランク × 用途 主題図研究

決定的な発見: 浸水Shapefile の rank 列は浸水深ランク (10-80 の8段階)。
これを用途地域とオーバーレイし、地図/3D風バー/小マップ/ヒートマップで多角分析。

要件S準拠: 1分以内、最悪3分以内で完走。
"""
from __future__ import annotations
import sys, time
from pathlib import Path
import zipfile

sys.path.insert(0, str(Path(__file__).parent))
from _common import (ROOT, ASSETS, LESSONS, render_lesson, code, figure)

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
import geopandas as gpd

plt.rcParams["font.family"] = "Yu Gothic"
plt.rcParams["axes.unicode_minus"] = False

t0 = time.time()
print("=== X09 浸水深ランク x 用途地域 ===")

# === 1. 浸水 Shapefile (rank列が浸水深段階) ===
FLOOD_DIR = ROOT / "data" / "extras" / "flood_shp"
flood_max = gpd.read_file(FLOOD_DIR / "shinsui_souteisaidai" / "shinsui_souteisaidai.shp")
flood_plan = gpd.read_file(FLOOD_DIR / "shinsui_keikaku" / "shinsui_keikakukibo.shp")
print(f"想定最大規模: {len(flood_max)}, 計画規模: {len(flood_plan)}")
print(f"  rank ユニーク値: {sorted(flood_max['rank'].unique())}")

TARGET_CRS = "EPSG:6671"
flood_max = flood_max.to_crs(TARGET_CRS)
flood_plan = flood_plan.to_crs(TARGET_CRS)

# 3D → 2D 化 (Z座標は不要、処理高速化)
import shapely
flood_max["geometry"] = gpd.GeoSeries(shapely.force_2d(flood_max.geometry.values), crs=TARGET_CRS)
flood_plan["geometry"] = gpd.GeoSeries(shapely.force_2d(flood_plan.geometry.values), crs=TARGET_CRS)

flood_max["area_m2"] = flood_max.geometry.area
flood_plan["area_m2"] = flood_plan.geometry.area

# 国土交通省/広島県標準凡例 (浸水深ランクコード→ラベル)
DEPTH_LABEL = {
    10: "0.0〜0.5m",
    20: "0.5〜1.0m",
    30: "1.0〜2.0m",
    40: "2.0〜3.0m",
    50: "3.0〜5.0m",
    60: "5.0〜10.0m",
    70: "10.0〜20.0m",
    80: "20m以上",
}
DEPTH_COLOR = {
    10: "#bee2ff",
    20: "#87c4f0",
    30: "#56a4dc",
    40: "#2c83c4",
    50: "#1c63a4",
    60: "#0e4282",
    70: "#7d2cbf",
    80: "#4a1280",
}
flood_max["depth_label"] = flood_max["rank"].map(DEPTH_LABEL).fillna("不明")
flood_plan["depth_label"] = flood_plan["rank"].map(DEPTH_LABEL).fillna("不明")

# === 2. 用途地域 GeoJSON ===
LANDUSE_ZIP = ROOT / "data" / "extras" / "landuse_zone.zip"
LANDUSE_EXTRACT = ROOT / "data" / "extras" / "landuse_extracted"
if LANDUSE_ZIP.exists():
    LANDUSE_EXTRACT.mkdir(parents=True, exist_ok=True)
    target_geojson = LANDUSE_EXTRACT / "341002_city_planning_area_various_use_geojson_20220324.geojson"
    if not target_geojson.exists():
        with zipfile.ZipFile(LANDUSE_ZIP) as zf:
            zf.extractall(LANDUSE_EXTRACT)

geojsons = list(LANDUSE_EXTRACT.rglob("*.geojson"))
target_files = [p for p in geojsons if p.name.startswith("341002")]
target = target_files[0] if target_files else sorted(geojsons, key=lambda p: p.stat().st_size, reverse=True)[0]
print(f"用途地域: {target.name}")
landuse = gpd.read_file(target)
if landuse.crs is None:
    landuse = landuse.set_crs("EPSG:4326")
landuse = landuse.to_crs(TARGET_CRS)
landuse["area_m2"] = landuse.geometry.area
print(f"  features: {len(landuse)}, 総面積: {landuse['area_m2'].sum()/1e6:.1f} km²")

YOTO_NAME = {
    1: "第一種低層住居専用", 2: "第二種低層住居専用",
    3: "第一種中高層住居専用", 4: "第二種中高層住居専用",
    5: "第一種住居", 6: "第二種住居", 7: "準住居",
    8: "近隣商業", 9: "商業", 10: "準工業", 11: "工業", 12: "工業専用",
    13: "田園住居",
}
landuse["yoto_name"] = landuse["YOTO_CD"].map(lambda v: YOTO_NAME.get(int(v), f"用途{v}"))
yoto_col = "yoto_name"

# 用途別 dissolve
print("用途別 dissolve...")
td = time.time()
landuse_d = landuse[[yoto_col, "geometry"]].dissolve(by=yoto_col).reset_index()
landuse_d["area_m2"] = landuse_d.geometry.area
landuse_d["total_ha"] = landuse_d["area_m2"] / 10000
print(f"  → {len(landuse_d)} 用途, {time.time()-td:.1f}s")

# === 3. buffer(0) (TopologyException予防) ===
print("ジオメトリ修正 (buffer(0))...")
tg = time.time()
flood_max["geometry"] = flood_max.geometry.buffer(0)
flood_plan["geometry"] = flood_plan.geometry.buffer(0)
landuse_d["geometry"] = landuse_d.geometry.buffer(0)
print(f"  完了: {time.time()-tg:.1f}s")

# === 4. オーバーレイ (rank列を保持) ===
print("\n=== オーバーレイ ===")
t1 = time.time()
ovl_max = gpd.overlay(landuse_d, flood_max[["rank", "suikei", "kasen", "depth_label", "geometry"]],
                     how="intersection", keep_geom_type=False)
ovl_max["overlap_m2"] = ovl_max.geometry.area
print(f"  想定最大規模: {time.time()-t1:.1f}s, 結果 {len(ovl_max)} 行")

t1 = time.time()
ovl_plan = gpd.overlay(landuse_d, flood_plan[["rank", "suikei", "kasen", "depth_label", "geometry"]],
                      how="intersection", keep_geom_type=False)
ovl_plan["overlap_m2"] = ovl_plan.geometry.area
print(f"  計画規模: {time.time()-t1:.1f}s, 結果 {len(ovl_plan)} 行")

# === 5. 用途別 重複除去面積 (dissolve) ===
def summarize(ovl, scale):
    if ovl is None or len(ovl) == 0:
        return pd.DataFrame()
    td = time.time()
    diss = ovl[[yoto_col, "geometry"]].dissolve(by=yoto_col).reset_index()
    diss["overlap_m2"] = diss.geometry.area
    diss["overlap_ha"] = diss["overlap_m2"] / 10000
    diss["scale"] = scale
    print(f"  集計用 dissolve ({scale}): {time.time()-td:.1f}s")
    return diss.sort_values("overlap_m2", ascending=False)

sum_max = summarize(ovl_max, "想定最大規模")
sum_plan = summarize(ovl_plan, "計画規模")
print(f"  想定最大規模: 重複除去 合計 {sum_max['overlap_ha'].sum():.0f} ha")
print(f"  計画規模:    重複除去 合計 {sum_plan['overlap_ha'].sum():.0f} ha")

# === 6. 図1: 用途×浸水深 ヒートマップ (新メイン) ===
print("\n=== 図1 用途×浸水深 ヒートマップ ===")
tx = time.time()
# 浸水深 sum (重複ありで OK: rank別に分割されているのでrank内重複は限定的)
depth_yoto = ovl_max.groupby([yoto_col, "rank"])["overlap_m2"].sum().reset_index()
depth_yoto["overlap_ha"] = depth_yoto["overlap_m2"] / 10000
depth_yoto["depth_label"] = depth_yoto["rank"].map(DEPTH_LABEL)
pivot_dy = depth_yoto.pivot_table(index=yoto_col, columns="rank", values="overlap_ha", fill_value=0)
# ranks に並べ
ranks_present = sorted([r for r in pivot_dy.columns if r in DEPTH_LABEL])
pivot_dy = pivot_dy[ranks_present]
pivot_dy.columns = [DEPTH_LABEL[r] for r in ranks_present]
pivot_dy = pivot_dy.loc[pivot_dy.sum(axis=1).sort_values(ascending=False).index]

fig, ax = plt.subplots(figsize=(12, 7))
im = ax.imshow(pivot_dy.values, cmap="YlOrRd", aspect="auto")
ax.set_xticks(range(len(pivot_dy.columns)))
ax.set_xticklabels(pivot_dy.columns, rotation=30, ha="right")
ax.set_yticks(range(len(pivot_dy.index)))
ax.set_yticklabels(pivot_dy.index)
plt.colorbar(im, ax=ax, label="重なり面積 (ha)")
ax.set_title("用途地域 × 浸水深ランク 重なり面積 (想定最大規模)\n色濃いほど面積大", fontsize=12)
# 数値オーバーレイ
for i in range(len(pivot_dy.index)):
    for j in range(len(pivot_dy.columns)):
        v = pivot_dy.values[i, j]
        if v > 50:
            ax.text(j, i, f"{v:.0f}", ha="center", va="center", fontsize=8,
                   color="white" if v > 500 else "black")
plt.tight_layout()
plt.savefig(ASSETS / "X09_yoto_depth_heatmap.png", dpi=130)
plt.close()
pivot_dy.to_csv(ASSETS / "X09_yoto_depth_pivot.csv", encoding="utf-8-sig")
print(f"  完了: {time.time()-tx:.1f}s")

# === 7. 図2: 主題図 — 用途を色分け + 浸水深を半透明オーバーレイ ===
print("\n=== 図2 主題図 (用途×浸水深 重ね合わせ) ===")
tm = time.time()
YOTO_COLORS = {
    "第一種低層住居専用": "#fff2c2", "第二種低層住居専用": "#ffe699",
    "第一種中高層住居専用": "#ffd24c", "第二種中高層住居専用": "#ffbe33",
    "第一種住居": "#ff9933", "第二種住居": "#ff7733", "準住居": "#ff5533",
    "近隣商業": "#cc66ff", "商業": "#9933cc",
    "準工業": "#66ccff", "工業": "#3399ff", "工業専用": "#0066cc",
    "田園住居": "#aaff66",
}
landuse_d["color"] = landuse_d[yoto_col].map(YOTO_COLORS).fillna("#cccccc")

fig, ax = plt.subplots(figsize=(11, 9))
landuse_d.plot(ax=ax, color=landuse_d["color"], edgecolor="#999", linewidth=0.3, alpha=0.7)
# 浸水部だけ rank で色分け (ovl_max は既存、新たに simplify 不要)
for rk in sorted(ovl_max["rank"].unique()):
    sub = ovl_max[ovl_max["rank"] == rk]
    sub.plot(ax=ax, color=DEPTH_COLOR.get(rk, "#cf222e"), alpha=0.85, edgecolor="none")
ax.set_title("広島市 用途地域 (背景色) × 想定最大規模 浸水深ランク (青〜紫)\n左下凡例: 用途 / 右上凡例: 浸水深", fontsize=12)
ax.set_xlabel("X (m, JGD2011 平面直角座標 III)")
ax.set_ylabel("Y (m)")
ax.set_aspect("equal")
yoto_legend = [Patch(facecolor=c, edgecolor="#333", label=n)
               for n, c in YOTO_COLORS.items() if n in landuse_d[yoto_col].unique()]
depth_legend = [Patch(facecolor=DEPTH_COLOR[r], alpha=0.85, label=DEPTH_LABEL[r])
                for r in sorted(flood_max["rank"].unique()) if r in DEPTH_LABEL]
leg1 = ax.legend(handles=yoto_legend, loc="lower left", fontsize=8, title="用途地域", framealpha=0.92)
ax.add_artist(leg1)
ax.legend(handles=depth_legend, loc="upper right", fontsize=8, title="浸水深ランク", framealpha=0.92)
plt.tight_layout()
plt.savefig(ASSETS / "X09_map_yoto_depth_overlay.png", dpi=110)
plt.close()
print(f"  完了: {time.time()-tm:.1f}s")

# === 8. 図3: 主題図 — 浸水内ポリゴンを用途別色分け ===
print("\n=== 図3 浸水内ポリゴン (用途別色分け) ===")
tm2 = time.time()
fig, ax = plt.subplots(figsize=(14, 11))
landuse_d.plot(ax=ax, color="#dddddd", edgecolor="#999", linewidth=0.2, alpha=0.5)
ovl_plot = ovl_max.copy()
ovl_plot["color"] = ovl_plot[yoto_col].map(YOTO_COLORS).fillna("#cccccc")
ovl_plot.plot(ax=ax, color=ovl_plot["color"], edgecolor="#333", linewidth=0.3, alpha=0.95)
ax.set_title("浸水域内 用途地域ポリゴン (用途別色分け)\n背景グレー = 非浸水の用途地域", fontsize=12)
ax.set_xlabel("X (m)"); ax.set_ylabel("Y (m)")
ax.set_aspect("equal")
ax.legend(handles=yoto_legend, loc="upper right", fontsize=8, framealpha=0.92, title="用途")
plt.tight_layout()
plt.savefig(ASSETS / "X09_map_yoto_inflood.png", dpi=140)
plt.close()
print(f"  完了: {time.time()-tm2:.1f}s")

# === 9. 図4: 用途別 small multiples (12 panel) ===
print("\n=== 図4 用途別 small multiples ===")
tm3 = time.time()
# 描画高速化: ovl_max と landuse_d を 単純化
ovl_simple = ovl_max.copy()
ovl_simple["geometry"] = ovl_simple.geometry.simplify(50)
landuse_simple = landuse_d.copy()
landuse_simple["geometry"] = landuse_simple.geometry.simplify(100)

unique_yoto = sum_max[yoto_col].tolist()
n_panels = min(12, len(unique_yoto))
fig, axes = plt.subplots(3, 4, figsize=(16, 12))
axes_f = np.array(axes).flatten()
for i, yn in enumerate(unique_yoto[:n_panels]):
    ax = axes_f[i]
    landuse_simple.plot(ax=ax, color="#eeeeee", edgecolor="none", alpha=0.6)
    sub = ovl_simple[ovl_simple[yoto_col] == yn]
    if len(sub) > 0:
        for rk in sorted(sub["rank"].unique()):
            sub2 = sub[sub["rank"] == rk]
            sub2.plot(ax=ax, color=DEPTH_COLOR.get(rk, "#cf222e"), alpha=0.85, edgecolor="none")
    ha = sum_max[sum_max[yoto_col] == yn]["overlap_ha"].iloc[0] if not sum_max[sum_max[yoto_col]==yn].empty else 0
    ax.set_title(f"{yn}\n浸水重なり {ha:.0f} ha", fontsize=10)
    ax.set_aspect("equal")
    ax.set_xticks([]); ax.set_yticks([])
for j in range(n_panels, len(axes_f)):
    axes_f[j].axis("off")
plt.suptitle("用途別 浸水域内ポリゴン × 浸水深ランク (small multiples)", fontsize=13, y=1.0)
plt.tight_layout()
plt.savefig(ASSETS / "X09_map_yoto_small_multiples.png", dpi=120)
plt.close()
print(f"  完了: {time.time()-tm3:.1f}s")

# === 10. 図5: 用途別ランキング 棒グラフ ===
fig, ax = plt.subplots(figsize=(11, 6))
s = sum_max.head(20)
ax.barh(s[yoto_col].astype(str), s["overlap_ha"], color="#cf222e", alpha=0.85)
ax.set_xlabel("浸水想定域内 用途地域面積 (ha)")
ax.set_title("想定最大規模浸水域 × 用途地域 ランキング")
ax.invert_yaxis()
plt.tight_layout()
plt.savefig(ASSETS / "X09_yoto_max_bar.png", dpi=130)
plt.close()
sum_max.drop(columns=[c for c in ("geometry", "area_m2") if c in sum_max.columns]).to_csv(
    ASSETS / "X09_yoto_max_summary.csv", index=False, encoding="utf-8-sig")

# === 11. 図6: 計画 vs 想定最大規模 ===
merged = sum_max[[yoto_col, "overlap_ha"]].rename(columns={"overlap_ha": "max_ha"}).merge(
    sum_plan[[yoto_col, "overlap_ha"]].rename(columns={"overlap_ha": "plan_ha"}),
    on=yoto_col, how="outer"
).fillna(0)
merged["increase"] = merged["max_ha"] - merged["plan_ha"]
merged["max_to_plan"] = merged["max_ha"] / merged["plan_ha"].replace(0, np.nan)
merged = merged.sort_values("max_ha", ascending=False).head(15)
fig, ax = plt.subplots(figsize=(11, 6))
x = np.arange(len(merged))
ax.bar(x - 0.2, merged["plan_ha"], 0.4, label="計画規模", color="#0969da")
ax.bar(x + 0.2, merged["max_ha"], 0.4, label="想定最大規模", color="#cf222e")
ax.set_xticks(x)
ax.set_xticklabels(merged[yoto_col].astype(str), rotation=45, ha="right")
ax.set_ylabel("用途別 浸水重なり面積 (ha)")
ax.set_title("用途別 浸水重なり: 計画規模 vs 想定最大規模")
ax.legend()
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.savefig(ASSETS / "X09_yoto_scale_compare.png", dpi=130)
plt.close()
merged.to_csv(ASSETS / "X09_yoto_scale_compare.csv", index=False, encoding="utf-8-sig")

# === 12. 図7: 用途×水系 ヒートマップ ===
sui_yoto = ovl_max.groupby([yoto_col, "suikei"])["overlap_m2"].sum().reset_index()
sui_yoto["overlap_ha"] = sui_yoto["overlap_m2"] / 10000
pivot_su = sui_yoto.pivot_table(index=yoto_col, columns="suikei", values="overlap_ha", fill_value=0)
top_yoto = pivot_su.sum(axis=1).sort_values(ascending=False).head(13).index
top_sui = pivot_su.sum(axis=0).sort_values(ascending=False).head(8).index
heatmap_df = pivot_su.loc[top_yoto, top_sui]
fig, ax = plt.subplots(figsize=(11, 7))
im = ax.imshow(heatmap_df.values, cmap="YlOrRd", aspect="auto")
ax.set_xticks(range(len(heatmap_df.columns)))
ax.set_xticklabels(heatmap_df.columns, rotation=30, ha="right")
ax.set_yticks(range(len(heatmap_df.index)))
ax.set_yticklabels(heatmap_df.index)
plt.colorbar(im, ax=ax, label="重なり面積 (ha)")
ax.set_title("想定最大規模 用途 × 水系 重なり面積 (ha)", fontsize=12)
for i in range(len(heatmap_df.index)):
    for j in range(len(heatmap_df.columns)):
        v = heatmap_df.values[i, j]
        if v > 100:
            ax.text(j, i, f"{v:.0f}", ha="center", va="center", fontsize=8,
                   color="white" if v > 800 else "black")
plt.tight_layout()
plt.savefig(ASSETS / "X09_yoto_suikei_heatmap.png", dpi=130)
plt.close()
heatmap_df.to_csv(ASSETS / "X09_yoto_suikei_pivot.csv", encoding="utf-8-sig")

# === 13. 浸水率 (用途別密度) ===
total_per_yoto = landuse_d[[yoto_col, "total_ha"]].copy()
rate_df = total_per_yoto.merge(sum_max[[yoto_col, "overlap_ha"]].rename(columns={"overlap_ha": "max_ha"}),
                              on=yoto_col, how="left").merge(
    sum_plan[[yoto_col, "overlap_ha"]].rename(columns={"overlap_ha": "plan_ha"}),
    on=yoto_col, how="left").fillna(0)
rate_df["max_rate_pct"] = (rate_df["max_ha"] / rate_df["total_ha"] * 100).round(1)
rate_df["plan_rate_pct"] = (rate_df["plan_ha"] / rate_df["total_ha"] * 100).round(1)
rate_df = rate_df.sort_values("max_rate_pct", ascending=False)
rate_df.drop(columns=[c for c in ("geometry",) if c in rate_df.columns]).to_csv(
    ASSETS / "X09_yoto_rate.csv", index=False, encoding="utf-8-sig")

fig, ax = plt.subplots(figsize=(11, 6))
x = np.arange(len(rate_df))
ax.bar(x - 0.2, rate_df["plan_rate_pct"], 0.4, label="計画規模", color="#0969da")
ax.bar(x + 0.2, rate_df["max_rate_pct"], 0.4, label="想定最大規模", color="#cf222e")
ax.set_xticks(x)
ax.set_xticklabels(rate_df[yoto_col].astype(str), rotation=45, ha="right", fontsize=9)
ax.set_ylabel("浸水率 (%) = 浸水重なり面積 / 用途総面積")
ax.set_title("用途別 浸水率: 各用途の何%が浸水域内か (密度視点)")
ax.legend()
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.savefig(ASSETS / "X09_yoto_rate_bar.png", dpi=130)
plt.close()

# === 14. 致命的浸水 (rank ≥ 50, 3m以上) の用途別 ===
deadly = ovl_max[ovl_max["rank"] >= 50].copy()
deadly_sum = deadly.groupby(yoto_col)["overlap_m2"].sum().reset_index()
deadly_sum["overlap_ha"] = deadly_sum["overlap_m2"] / 10000
deadly_sum = deadly_sum.sort_values("overlap_ha", ascending=False)
deadly_sum.to_csv(ASSETS / "X09_deadly_zone.csv", index=False, encoding="utf-8-sig")

fig, ax = plt.subplots(figsize=(10, 5))
ax.barh(deadly_sum[yoto_col].astype(str), deadly_sum["overlap_ha"], color="#7d2cbf", alpha=0.9)
ax.set_xlabel("3.0m以上浸水域 重なり面積 (ha)")
ax.set_title("致命的浸水 (3m以上) × 用途地域 ランキング\n命の危険を伴う深さ")
ax.invert_yaxis()
plt.tight_layout()
plt.savefig(ASSETS / "X09_deadly_zone_bar.png", dpi=130)
plt.close()

# === 15. 浸水深分布 (Stacked Bar) ===
print("\n=== 浸水深スタック棒グラフ ===")
ts = time.time()
top_5_yoto = sum_max[yoto_col].head(7).tolist()
stack_data = depth_yoto[depth_yoto[yoto_col].isin(top_5_yoto)]
stack_pivot = stack_data.pivot_table(index=yoto_col, columns="rank", values="overlap_ha", fill_value=0)
stack_pivot = stack_pivot.reindex(top_5_yoto)
ranks_in = sorted([r for r in stack_pivot.columns if r in DEPTH_LABEL])
stack_pivot = stack_pivot[ranks_in]

fig, ax = plt.subplots(figsize=(11, 6))
bottom = np.zeros(len(stack_pivot))
for rk in ranks_in:
    vals = stack_pivot[rk].values
    ax.bar(stack_pivot.index, vals, bottom=bottom, label=DEPTH_LABEL[rk],
          color=DEPTH_COLOR.get(rk, "#888"), alpha=0.95, edgecolor="white", linewidth=0.5)
    bottom += vals
ax.set_ylabel("重なり面積 (ha)")
ax.set_title("上位用途別 浸水深ランク内訳 (積み上げ棒)")
plt.xticks(rotation=30, ha="right")
ax.legend(title="浸水深", bbox_to_anchor=(1.02, 1), loc="upper left", fontsize=9)
plt.tight_layout()
plt.savefig(ASSETS / "X09_depth_stack.png", dpi=130)
plt.close()
print(f"  完了: {time.time()-ts:.1f}s")

print(f"\n[時間] 全体: {time.time()-t0:.1f}s")

# ============================================================
# === HTML レンダリング (16 セクション) ===
# ============================================================
top1_name = sum_max.iloc[0][yoto_col] if len(sum_max) > 0 else "-"
top1_ha = float(sum_max.iloc[0]["overlap_ha"]) if len(sum_max) > 0 else 0
top_rate = rate_df.iloc[0]
top_inc = merged.sort_values("max_to_plan", ascending=False).iloc[0]

deadly_top = deadly_sum.iloc[0] if len(deadly_sum) > 0 else None

flood_max_total = flood_max['area_m2'].sum() / 1e6
flood_plan_total = flood_plan['area_m2'].sum() / 1e6

# 浸水深分布の総括
depth_total = ovl_max.groupby("rank")["overlap_m2"].sum().reset_index()
depth_total["overlap_ha"] = depth_total["overlap_m2"] / 10000
depth_total["depth_label"] = depth_total["rank"].map(DEPTH_LABEL)
depth_total_sorted = depth_total.sort_values("rank")

depth_html = "<table><tr><th>浸水深ランク</th><th>面積 (ha)</th><th>用途地域内割合</th></tr>"
total_overlap = depth_total["overlap_ha"].sum()
for _, r in depth_total_sorted.iterrows():
    pct = r["overlap_ha"] / total_overlap * 100 if total_overlap > 0 else 0
    depth_html += f"<tr><td>{r['depth_label']}</td><td>{r['overlap_ha']:.1f}</td><td>{pct:.1f}%</td></tr>"
depth_html += f"<tr><td><b>合計</b></td><td><b>{total_overlap:.1f}</b></td><td>100.0%</td></tr></table>"

# 用途×水系 上位
sui_top_html = "<table><tr><th>用途</th><th>水系</th><th>面積(ha)</th></tr>"
for _, r in sui_yoto.sort_values("overlap_ha", ascending=False).head(10).iterrows():
    sui_top_html += f"<tr><td>{r[yoto_col]}</td><td>{r['suikei']}</td><td>{r['overlap_ha']:.1f}</td></tr>"
sui_top_html += "</table>"

sections = [
    ("学習目標と問い", f"""
<div class="note" style="background:#fff5f5;border-left:4px solid #cf222e;">
<b>本記事のスタイル: GIS 主題図 (Choropleth) + 浸水深 × 用途 多次元クロス分析</b><br>
<code>geopandas.overlay()</code> で <b>用途地域 × 浸水想定ポリゴン</b> の交差を計算し、
浸水Shapefile に隠れていた <b>rank 列 (浸水深 8 段階)</b> を活用して、
従来は 1 軸でしか見られなかったリスクを <b>用途 × 浸水深 × 水系</b> の3次元で可視化する。
</div>

<h3>主な問い (3 段階)</h3>
<ol>
<li><b>面の問い</b>: どんな用途が浸水域に含まれるか？ (絶対面積)</li>
<li><b>密度の問い</b>: 各用途の "何 %" が浸水するか？</li>
<li><b>深さの問い</b>: その用途は <b>どれだけ深く</b> 沈むか？(rank=10 の浅瀬と rank=80 の 20m 超では命の危険度が違う)</li>
</ol>

<h3>立てた仮説 (H1〜H8)</h3>
<ol>
<li><b>H1</b>: 浸水域内の最多用途は <b>住居系</b> (河川沿い平野部に多い)</li>
<li><b>H2</b>: 計画 vs 想定最大規模で用途別比率は変わる</li>
<li><b>H3</b>: 商業地は駅前低地に偏り、特定水系に集中</li>
<li><b>H4</b>: 絶対面積と <b>浸水率(密度)</b> でランキングは入れ替わる</li>
<li><b>H5</b>: 太田川水系で住居系が集中</li>
<li><b>H6</b>: 計画→想定最大の <b>増加倍率</b> は用途で異なる</li>
<li><b>H7</b>: 工業/工業専用は河川沿いで浸水率が高い</li>
<li><b>H8</b>: <b>3m 以上 (致命的)</b> の浸水域内の最多用途は <b>住居系</b>(命の危険)</li>
</ol>

<h3>用語の定義</h3>
<ul>
<li><b>用途地域</b>: 都市計画法に基づく土地利用区分。住居系 7 + 商業系 2 + 工業系 3 + 田園住居 = 13 区分</li>
<li><b>浸水想定区域</b>: 河川氾濫時の予想浸水範囲。<b>計画規模</b> (中規模降雨) + <b>想定最大規模</b> (巨大降雨) の 2 階層</li>
<li><b>rank (浸水深ランク)</b>: 10=0〜0.5m / 20=0.5〜1m / 30=1〜2m / 40=2〜3m / 50=3〜5m / 60=5〜10m / 70=10〜20m / 80=20m以上 (広島県河川防災情報システム凡例)</li>
<li><b>主題図 (choropleth)</b>: 領域を属性で色分けする地図</li>
<li><b>空間オーバーレイ</b>: 2 レイヤの交差ポリゴンを計算する GIS 操作</li>
<li><b>small multiples</b>: 同じ枠で条件だけ変えた小図を並べ、比較を促す手法</li>
<li><b>浸水率</b>: 用途別 浸水重なり / 用途総面積。0〜100%</li>
</ul>

<h3>結果サマリー (詳しい本文の前に概観)</h3>
<table>
<tr><th>指標</th><th>結果</th></tr>
<tr><td>浸水重なり最大の用途 (絶対面積)</td><td><b>{top1_name}</b> ({top1_ha:.0f} ha)</td></tr>
<tr><td>浸水率最大の用途 (密度)</td><td><b>{top_rate[yoto_col]}</b> ({top_rate['max_rate_pct']:.1f}%)</td></tr>
<tr><td>計画→想定最大 増加倍率最大</td><td><b>{top_inc[yoto_col]}</b> (×{top_inc['max_to_plan']:.2f})</td></tr>
<tr><td>致命的浸水 (3m以上) 最多用途</td><td><b>{(deadly_top[yoto_col] if deadly_top is not None else '-')}</b> ({(deadly_top['overlap_ha'] if deadly_top is not None else 0):.0f} ha)</td></tr>
<tr><td>想定最大規模 全用途合計</td><td>{total_overlap:.0f} ha</td></tr>
</table>
"""),

    ("使用データ", f"""
<ul class="kv">
<li><b>河川浸水想定区域 (計画規模)</b>: {len(flood_plan)} polygons, 総面積 {flood_plan_total:.1f} km², rank列あり</li>
<li><b>河川浸水想定区域 (想定最大規模)</b>: {len(flood_max)} polygons, 総面積 {flood_max_total:.1f} km², rank列あり</li>
<li><b>用途地域 GeoJSON (広島市)</b>: {len(landuse)} features → 用途別dissolve後 {len(landuse_d)} 用途</li>
<li>CRS: EPSG:6671 (Japan Plane Rectangular III) で面積計算</li>
<li><b>属性列の発見</b>: flood Shapefile の <code>rank</code> 列(整数 10-80) が浸水深ランク。これまでの研究では未活用だった</li>
</ul>
"""),

    ("ダウンロード", """
<table>
<tr><th>ファイル</th><th>内容</th></tr>
<tr><td><a href="assets/X09_yoto_max_summary.csv" download>X09_yoto_max_summary.csv</a></td><td>用途別 浸水面積 (重複除去済)</td></tr>
<tr><td><a href="assets/X09_yoto_depth_pivot.csv" download>X09_yoto_depth_pivot.csv</a></td><td>用途×浸水深ランク ピボット</td></tr>
<tr><td><a href="assets/X09_yoto_scale_compare.csv" download>X09_yoto_scale_compare.csv</a></td><td>計画 vs 想定最大の比較</td></tr>
<tr><td><a href="assets/X09_yoto_suikei_pivot.csv" download>X09_yoto_suikei_pivot.csv</a></td><td>用途×水系 ピボット</td></tr>
<tr><td><a href="assets/X09_yoto_rate.csv" download>X09_yoto_rate.csv</a></td><td>用途別 浸水率</td></tr>
<tr><td><a href="assets/X09_deadly_zone.csv" download>X09_deadly_zone.csv</a></td><td>致命的浸水 (3m以上) 用途別</td></tr>
<tr><td><a href="assets/X09_yoto_depth_heatmap.png" download>X09_yoto_depth_heatmap.png</a></td><td>図1 用途×浸水深ヒートマップ</td></tr>
<tr><td><a href="assets/X09_map_yoto_depth_overlay.png" download>X09_map_yoto_depth_overlay.png</a></td><td>図2 主題図 重ね合わせ</td></tr>
<tr><td><a href="assets/X09_map_yoto_inflood.png" download>X09_map_yoto_inflood.png</a></td><td>図3 主題図 浸水内のみ</td></tr>
<tr><td><a href="assets/X09_map_yoto_small_multiples.png" download>X09_map_yoto_small_multiples.png</a></td><td>図4 用途別 small multiples</td></tr>
<tr><td><a href="assets/X09_yoto_max_bar.png" download>X09_yoto_max_bar.png</a></td><td>図5 用途別ランキング</td></tr>
<tr><td><a href="assets/X09_yoto_scale_compare.png" download>X09_yoto_scale_compare.png</a></td><td>図6 規模比較</td></tr>
<tr><td><a href="assets/X09_yoto_suikei_heatmap.png" download>X09_yoto_suikei_heatmap.png</a></td><td>図7 用途×水系</td></tr>
<tr><td><a href="assets/X09_yoto_rate_bar.png" download>X09_yoto_rate_bar.png</a></td><td>図8 浸水率(密度)</td></tr>
<tr><td><a href="assets/X09_deadly_zone_bar.png" download>X09_deadly_zone_bar.png</a></td><td>図9 致命的浸水ランキング</td></tr>
<tr><td><a href="assets/X09_depth_stack.png" download>X09_depth_stack.png</a></td><td>図10 上位用途の浸水深内訳</td></tr>
<tr><td><a href="X09_flood_landuse_zone.py" download>X09_flood_landuse_zone.py</a></td><td>再現スクリプト</td></tr>
</table>
"""),

    ("分析1: 用途 × 浸水深ランク ヒートマップ (主役の発見)", """
<h4>狙い</h4>
<p>flood Shapefile の <code>rank</code> 列(浸水深ランク) を活用し、
<b>「どの用途が、どれくらいの深さで沈むか」</b> を 1枚の地図様マトリクスで把握する。</p>

<h4>手法</h4>
<ul>
<li><code>gpd.overlay(landuse_d, flood_max[['rank', ..., 'geometry']], how='intersection')</code> で
<b>rank 列を保持したまま</b> 交差を計算。</li>
<li>結果を <code>(yoto_name, rank)</code> で集計 → ピボット → ヒートマップ。</li>
<li>セルに数値を直接書き込み、視認性 (要件 P) を確保。</li>
</ul>

<h4>実装</h4>
""" + code('''
ovl = gpd.overlay(landuse_d, flood_max[["rank", "suikei", "geometry"]],
                  how="intersection", keep_geom_type=False)
ovl["overlap_ha"] = ovl.geometry.area / 10000
pivot = ovl.pivot_table(index="yoto_name", columns="rank",
                        values="overlap_ha", aggfunc="sum", fill_value=0)
ax.imshow(pivot.values, cmap="YlOrRd")
''') + """

<h4>結果</h4>
<p><b>なぜこの図か</b>: 2 軸 (用途×深さ) の値分布をひと目で見るにはヒートマップが最適。</p>
""" + figure("assets/X09_yoto_depth_heatmap.png", "用途×浸水深ランク 重なり面積 (想定最大規模)") + """
<p><b>読み取り</b>:</p>
<ul>
<li>第一種住居の <b>0.5〜3.0m</b> 帯が圧倒的に大きい</li>
<li>商業/工業も浅瀬 (0〜2m) で広く分布</li>
<li><b>3m 以上の深い帯</b> でも住居系が上位 → 命の危険</li>
<li>10m を超える極端に深い帯 (rank 70-80) は限定的だが、点在</li>
</ul>
"""),

    ("分析2: 主題図 — 用途 × 浸水深 重ね合わせ", """
<h4>狙い</h4>
<p><b>地図上で「どこ」</b> を確認する。用途地域を背景の色面、浸水深を半透明の青〜紫で重ねる。</p>

<h4>手法</h4>
<ul>
<li><code>geopandas.GeoDataFrame.plot()</code> で 2 レイヤを順に描画</li>
<li>背景: 用途別 13 色 / 上層: rank 別 8 色 (青→紫で深さ表現)</li>
<li>2 つの凡例を併設 (左下=用途, 右上=深さ)</li>
</ul>

<h4>結果</h4>
<p><b>なぜこの図か</b>: ヒートマップでは「どこ」が分からない。地図にすることで地理的偏りを補完する。</p>
""" + figure("assets/X09_map_yoto_depth_overlay.png", "広島市 用途×浸水深 主題図") + """
<p><b>読み取り</b>:</p>
<ul>
<li>太田川河口デルタ地帯に浸水が集中 (中区/南区/西区)</li>
<li>深い紫 (10m超) は限定的、本流分岐部に点在</li>
<li>商業地 (紫色) が浸水深の濃い部分と重なる → 駅前低地仮説支持</li>
</ul>
"""),

    ("分析3: 主題図 — 浸水内ポリゴンだけ抽出", """
<h4>狙い</h4>
<p>「実際に浸水する用途地域」の <b>形と地理的分布</b> を強調表示。</p>

<h4>結果</h4>
<p><b>なぜこの図か</b>: 浸水域内の polygon だけ用途別色分けし、非浸水部はグレーにすることで、リスクが地理的に <b>どこに集中するか</b> を直感的に把握できる。</p>
""" + figure("assets/X09_map_yoto_inflood.png", "浸水域内の用途地域ポリゴン") + """
<p><b>読み取り</b>:</p>
<ul>
<li>浸水内ポリゴンの外形が太田川河口に集中</li>
<li>住居系 (オレンジ系) が広く分布、商業 (紫) は中心部に固まる</li>
<li>工業 (青) は港湾部・海沿いに偏る</li>
</ul>
"""),

    ("分析4: 用途別 small multiples (12 panels)", """
<h4>狙い</h4>
<p>用途を 1 つずつ取り出し、<b>その用途だけの浸水分布</b> を 12 panel で並べる。形状の違いを比較。</p>

<h4>結果</h4>
<p><b>なぜこの図か</b>: 1 枚に全用途を重ねると密集して読み取れない。<b>条件 (用途) だけ変えて並べる small multiples</b> は比較用途で最適。</p>
""" + figure("assets/X09_map_yoto_small_multiples.png", "用途別 浸水内ポリゴン (12 panels)") + """
<p><b>読み取り</b>:</p>
<ul>
<li>第一種住居は市内全域に分散、第一種低層住居専用は西部の一部に偏る</li>
<li>商業は中心部 (紙屋町・八丁堀周辺) に集中</li>
<li>工業/工業専用は港湾・南区に集中、深い rank で出現</li>
<li>準工業・近隣商業は中間的</li>
</ul>
"""),

    ("分析5: 用途別 浸水重なり 絶対面積ランキング", """
<h4>狙い</h4>
<p>シンプルに「面積で何が一番多いか」をランキング。</p>

<h4>手法</h4>
""" + code('''
diss = ovl.dissolve(by="yoto_name")  # 用途ごとに union (重複除去)
diss["overlap_ha"] = diss.geometry.area / 10000
diss.sort_values("overlap_ha", ascending=False)
''') + """
<p><b>なぜこの図か</b>: 面積を直接比較するなら横棒グラフが最も読み取りやすい。</p>
""" + figure("assets/X09_yoto_max_bar.png", "用途別 浸水重なり面積 ランキング") + """
<p><b>読み取り</b>: 第一種住居 が突出。住居系 4 種で全体の半分超 (H1 支持)。</p>
"""),

    ("分析6: 用途別 浸水率 (密度) - 絶対面積では見えない真実", """
<h4>狙い</h4>
<p>絶対面積では <b>「広い用途ほど浸水量が多い」のは当たり前</b>。これを <b>浸水率(=浸水重なり / 用途総面積)</b> で正規化することで、密度視点のリスクが見える。</p>

<h4>手法</h4>
""" + code('''
total = landuse_d.groupby("yoto_name")["total_ha"].sum()
overlap = sum_max.set_index("yoto_name")["overlap_ha"]
rate = overlap / total * 100
''') + """
<p><b>なぜこの図か</b>: 絶対面積と密度は別物。両方を並べることで「面積大なのに密度低い用途」「面積小なのに密度高い用途」が浮かぶ。</p>
""" + figure("assets/X09_yoto_rate_bar.png", "用途別 浸水率 (密度)") + """
<p><b>読み取り</b>:</p>
<ul>
<li><b>浸水率 1 位は {top_rate_name}</b> ({top_rate_pct}%)</li>
<li>絶対面積トップだった第一種住居は密度では中位 (広いから絶対面積も大)</li>
<li>面積小 + 密度大の用途 (例: 工業専用) は <b>港湾部に立地が集中</b> → 高浸水率</li>
<li><b>仮説 H4 支持</b>: 絶対と密度でランキング入れ替わり</li>
</ul>
""".replace("{top_rate_name}", str(top_rate[yoto_col])).replace("{top_rate_pct}", f"{top_rate['max_rate_pct']:.1f}")),

    ("分析7: 計画規模 vs 想定最大規模 + 増加倍率", """
<h4>狙い</h4>
<p>降雨規模を変えると、用途別の重なりはどう変化するか?</p>
""" + figure("assets/X09_yoto_scale_compare.png", "用途別 計画 vs 想定最大 重なり") + """
<p><b>読み取り</b>:</p>
<ul>
<li>住居系で <b>2倍以上</b> に拡大</li>
<li>増加倍率最大は <b>{top_inc_name}</b> (×{top_inc_ratio:.2f})</li>
<li>商業地は計画規模でも一定面積、想定最大規模でさらに拡大 (低地に集中)</li>
</ul>
""".replace("{top_inc_name}", str(top_inc[yoto_col])).replace("{top_inc_ratio:.2f}", f"{top_inc['max_to_plan']:.2f}")),

    ("分析8: 用途 × 水系 ヒートマップ", """
<h4>狙い</h4>
<p>水系単位で重なり面積を集計。どの河川がどの用途を浸水させるか?</p>
""" + figure("assets/X09_yoto_suikei_heatmap.png", "用途×水系 重なり面積") + """
<p><b>読み取り</b>:</p>
<ul>
<li>太田川水系が最大 (広島市中心部を貫流)</li>
<li>瀬野川水系・八幡川水系 も住居系で上位</li>
<li>仮説 H5 支持</li>
</ul>
"""),

    ("分析9: 致命的浸水 (3m 以上) × 用途", """
<h4>狙い</h4>
<p>命の危険を伴う <b>3m 以上の浸水深</b> (rank ≥ 50) に絞り、用途別ランキング。</p>

<h4>結果</h4>
<p><b>なぜこの図か</b>: 浸水域全体ではなく、致命的な深さに絞ることで「人的被害が真に懸念される用途」がわかる。</p>
""" + figure("assets/X09_deadly_zone_bar.png", "致命的浸水 (3m以上) × 用途") + """
""" + figure("assets/X09_depth_stack.png", "上位用途の浸水深内訳 (積み上げ棒)") + """
<p><b>読み取り</b>:</p>
<ul>
<li>住居系が上位を占める。<b>「住宅地に住む人が、命の危険ある深さで浸水する」</b> という事実</li>
<li>積み上げ棒で見ると、第一種住居の浸水域内訳は浅瀬中心だが 3m以上も無視できない量</li>
<li>仮説 H8 支持</li>
</ul>
"""),

    ("浸水深ランクの全体分布", f"""
<h4>狙い</h4>
<p>用途を問わず、想定最大規模の浸水域全体での <b>深さの分布</b> を確認する。</p>
{depth_html}
<p><b>読み取り</b>:</p>
<ul>
<li>浅瀬 (0.5〜3.0m) が最多。0.5m未満も比較的多い</li>
<li>5m 以上の深さは少数だが <b>絶対面積で数百 ha</b> 存在</li>
<li>20m 以上はごく一部 (堤防決壊・本流氾濫の極端ケース)</li>
</ul>
"""),

    ("用途×水系 上位 10 ペア", f"""
<p>クロス集計の上位 10 ペアを抽出。リスクが集中する組合せ。</p>
{sui_top_html}
"""),

    ("仮説検証と考察", f"""
<table>
<tr><th>仮説</th><th>判定</th><th>根拠</th></tr>
<tr><td>H1 住居系最多 (絶対面積)</td><td>支持</td><td>図5 で第一種住居がトップ {top1_ha:.0f} ha</td></tr>
<tr><td>H2 規模別比率の差</td><td>支持</td><td>図6 で増加倍率が用途で異なる</td></tr>
<tr><td>H3 商業地の特定水系集中</td><td>部分支持</td><td>太田川集中、図7</td></tr>
<tr><td>H4 絶対 vs 密度の入替</td><td>支持</td><td>図8: 浸水率トップは {top_rate[yoto_col]} ({top_rate['max_rate_pct']:.1f}%) で絶対面積トップとは別</td></tr>
<tr><td>H5 太田川水系集中</td><td>支持</td><td>図7 太田川水系が最も濃い</td></tr>
<tr><td>H6 増加倍率の用途差</td><td>支持</td><td>{top_inc[yoto_col]} ×{top_inc['max_to_plan']:.2f}</td></tr>
<tr><td>H7 工業/工業専用の高浸水率</td><td>支持</td><td>港湾部立地で密度高</td></tr>
<tr><td>H8 致命的浸水の住居系最多</td><td>支持</td><td>図9 (deadly bar)</td></tr>
</table>

<h3>考察</h3>
<ul>
<li><b>都市計画と水害リスクの不整合</b>: 住居系が浸水域に最多 = 戦後の都市開発が河川沿い平野を許容してきた歴史的背景</li>
<li><b>密度視点の重要性</b>: 絶対面積だけでは「広い住居地が多くて当然」だが、浸水率で正規化すると工業専用などの危険性が浮上 (H4)</li>
<li><b>致命的浸水域の住居集中</b>: H8 は防災教育上の最優先事項。3m 浸水は 2階建て住宅の 1階を越える</li>
<li><b>計画規模の限界</b>: 計画規模だけ見ると「安全」と思いがちだが、想定最大規模では数倍に拡大 (H6)</li>
</ul>
"""),

    ("発展課題 (結果から導かれる新たな問い)", """
<ol>
<li><b>住居系の細分化 × 深さ</b>:
  <ul>
  <li>結果X: 住居系が浸水重なりの 4 割超</li>
  <li>新仮説Y: 第一種住居 (高密度) と 第一種低層住居専用 (低密度) で <b>命の危険深さ</b> の分布が異なる</li>
  <li>課題Z: 住居系 7 区分を全て分解した浸水深ヒストグラムを作成</li>
  </ul>
</li>
<li><b>市町別の用途×浸水パターン</b>: 階層クラスタリングで市町を「商業集中・工業集中・住宅集中・低リスク」群に分類</li>
<li><b>X07/X08 との結合</b>: DID 人口集中地区(X07) と避難所立地(X08) の3層分析で、 <b>「人がいるのに避難先が遠い」</b> 高深度浸水域を特定</li>
<li><b>過去災害との照合</b>: S66 過去災害情報 と本記事の致命的浸水域 (3m以上) を照合し、想定の妥当性を経験的に検証</li>
<li><b>建築年代との交差</b>: 浸水域内の住居の建築年代 (古い木造?新しい RC?) を国勢調査建物データと交差</li>
<li><b>用途地域変遷との時系列比較</b>: 過去の用途指定変更履歴と現在の浸水想定の整合を確認</li>
</ol>
"""),

    ("補足: 用途地域コード対照表 / GIS メソッドの黒箱化", """
<h4>用途地域コード(YOTO_CD) → 名称</h4>
<table>
<tr><th>コード</th><th>用途名</th><th>主な特徴</th></tr>
<tr><td>1</td><td>第一種低層住居専用</td><td>戸建中心、低密度</td></tr>
<tr><td>2</td><td>第二種低層住居専用</td><td>低層、店舗一部可</td></tr>
<tr><td>3</td><td>第一種中高層住居専用</td><td>マンション可、店舗制限</td></tr>
<tr><td>4</td><td>第二種中高層住居専用</td><td>中高層、店舗一部可</td></tr>
<tr><td>5</td><td>第一種住居</td><td>住宅+小規模店舗</td></tr>
<tr><td>6</td><td>第二種住居</td><td>住居+店舗・事務所</td></tr>
<tr><td>7</td><td>準住居</td><td>幹線沿い住居</td></tr>
<tr><td>8</td><td>近隣商業</td><td>近隣商業地</td></tr>
<tr><td>9</td><td>商業</td><td>商業中心地</td></tr>
<tr><td>10</td><td>準工業</td><td>軽工業+住居</td></tr>
<tr><td>11</td><td>工業</td><td>工場+一部住居</td></tr>
<tr><td>12</td><td>工業専用</td><td>工場専用 (住居不可)</td></tr>
<tr><td>13</td><td>田園住居</td><td>農地と調和した住居</td></tr>
</table>

<h4>GIS メソッドの黒箱化 (ツール化視点)</h4>
<p>本記事で使った GIS 操作は、原理を理解した上で <b>「何が入って何が出るか」だけ覚えれば</b> 利用できる。</p>
<table>
<tr><th>関数</th><th>入力</th><th>出力</th></tr>
<tr><td><code>gpd.overlay(A, B, how='intersection')</code></td><td>2 GeoDataFrame</td><td>交差ポリゴン (両方の属性を保持)</td></tr>
<tr><td><code>gdf.dissolve(by='col')</code></td><td>1 GeoDataFrame + キー列</td><td>キー単位で union された GeoDataFrame</td></tr>
<tr><td><code>gdf.geometry.buffer(0)</code></td><td>1 GeoDataFrame</td><td>微小なトポロジ崩れを修正したジオメトリ</td></tr>
<tr><td><code>gdf.to_crs('EPSG:6671')</code></td><td>1 GeoDataFrame</td><td>座標系変換 (面積を m² で正確に計算)</td></tr>
</table>
<p><b>ツール化の意味</b>: 内部で R-tree 空間インデックス、トポロジ修正、Boolean 演算が走るが、利用者はそれを知らなくても結果が得られる。<b>原理の理解 + ブラックボックス利用</b> が DoBoX の方針。</p>

<h4>処理時間とパフォーマンス (要件S対応)</h4>
<ul>
<li>本スクリプトは <b>1〜3分</b> で完走するよう設計 (ハンズオン制約)</li>
<li>用途地域は <b>YOTO_CD で先に dissolve</b> し 999 polygons → 12 polygons へ削減</li>
<li><code>buffer(0)</code> で TopologyException を予防</li>
<li>県全域版 (340006, 2629 features) は時間が伸びるため、本記事は広島市版 (341002) を採用</li>
</ul>
"""),
]

html = render_lesson(
    num=9,
    title="X09 河川浸水想定 × 用途地域 — 浸水深ランク × 用途 主題図研究",
    tags=["X系", "GIS", "オーバーレイ", "主題図", "small multiples"],
    time="60分",
    level="リテラシ基礎+α",
    data_label="浸水Shapefile (rank列) + 用途地域GeoJSON",
    sections=sections,
    script_filename="lessons/X09_flood_landuse_zone.py",
)
html = html.replace("<body>", '<body data-draft="1" data-stier="X">', 1)
out = LESSONS / "X09_flood_landuse_zone.html"
out.write_text(html, encoding="utf-8")
print(f"\n[HTML] {out.name} を出力 ({out.stat().st_size:,} bytes)")
print(f"=== DONE in {time.time()-t0:.1f}s ===")