"""L11 トリプルハザード — 河川浸水 × 土砂災害 × 雪崩 重ね合わせ

広島県の「水・土・雪」3 種ハザードを同一座標系・同一格子で空間重ね合わせし、
重複度 (0/1/2/3) のリスクマップを作る GIS 主題図研究。

使用データ (DoBoX 全県版):
- 河川浸水想定区域 全河川版 #295 (計画) / #313 (想定最大)  → 39 dataset_id 論理カバー
- 土砂災害警戒区域 県全域版 #48 (土石流 #79 / 急傾斜 #80 / 地すべり #81)  → 31 dataset_id 論理カバー
- 雪崩危険箇所 県全域版 #50 (#77)  → 31 dataset_id 論理カバー
計 101 dataset_id を 4 ファイルで論理カバー。

要件S (ハンズオン 1〜3 分) のため:
- buffer(0) は 1 度のみ
- 県全域 dissolve は 1 ハザード 1 回 (R-tree spatial index 内蔵)
- 重ね合わせは「2km 格子 × R-tree」で評価 (overlay は使わない)
- 細かい見た目はラスタ的グリッドで担保
"""
from __future__ import annotations
import sys, time, json
from pathlib import Path

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
from shapely.geometry import box
import shapely

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

t0 = time.time()
print("=== L11 トリプルハザード (浸水 × 土砂 × 雪崩) ===")

TARGET_CRS = "EPSG:6671"  # JGD2011 平面直角 III, 中国地方, 単位 m

# === 1. 河川浸水想定区域 (想定最大規模 = #313, 計画 = #295) =====================
print("\n[1] 河川浸水 Shapefile ...")
t1 = time.time()
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")
flood_max = flood_max.to_crs(TARGET_CRS)
flood_plan = flood_plan.to_crs(TARGET_CRS)
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)
# 注: buffer(0) は使わない (巨大ポリゴンで OOM 落ちする報告あり)。
# sjoin の predicate='intersects' は invalid topology でも動くので問題なし。
print(f"  最大: {len(flood_max)}, 計画: {len(flood_plan)}, {time.time()-t1:.1f}s", flush=True)

# === 2. 土砂災害警戒区域 ======================================================
print("\n[2] 土砂災害 Shapefile ...")
t1 = time.time()
SED_DIR = ROOT / "data" / "extras" / "sediment_shp"
doseki = gpd.read_file(SED_DIR / "doseki" / "340006_sediment_disaster_hazard_area_debris_flow_20260427" / "73_031drp_34_20260427.shp")
kyukei = gpd.read_file(SED_DIR / "kyukeisha" / "340006_sediment_disaster_hazard_area_steep_slope_20260427" / "73_031krp_34_20260427.shp")
jisuberi = gpd.read_file(SED_DIR / "jisuberi" / "340006_sediment_disaster_hazard_area_landslide_20260427" / "73_031jy_34_20260427.shp")
print(f"  土石流: {len(doseki)}, 急傾斜: {len(kyukei)}, 地すべり: {len(jisuberi)}")
sediment = pd.concat([doseki, kyukei, jisuberi], ignore_index=True)
sediment = gpd.GeoDataFrame(sediment, geometry="geometry", crs=doseki.crs)
sediment = sediment.to_crs(TARGET_CRS)
print(f"  合計 {len(sediment)} polygon, {time.time()-t1:.1f}s")

# === 3. 雪崩危険箇所 ==========================================================
print("\n[3] 雪崩 Shapefile ...")
t1 = time.time()
AVA_DIR = ROOT / "data" / "extras" / "avalanche_shp"
avalanche = gpd.read_file(AVA_DIR / "340006_avalanche_danger_point_20260427" / "74_34_20260427.shp")
avalanche = avalanche.to_crs(TARGET_CRS)
print(f"  雪崩: {len(avalanche)}, {time.time()-t1:.1f}s")

# === 4-5. ジオメトリ修復は描画時 simplify でカバーし、grid 判定は spatial index sjoin で実行 ====
# unary_union は土砂 4 万ポリゴンで TopologyException + 数十秒なので使わない。
# 代わりに gdf.sindex の R-tree spatial index に直接当てる。
print("\n[4] spatial index を構築 ...")
t1 = time.time()
flood_sindex = flood_max.sindex
sed_sindex = sediment.sindex
ava_sindex = avalanche.sindex
print(f"  完了 {time.time()-t1:.1f}s")

# === 6. 評価グリッドを作成 (広島県全域を覆う 2km × 2km 格子) ===================
print("\n[6] 2km grid を生成 ...")
t1 = time.time()
# bbox は 4 ハザード共通の和集合
xmin = min(flood_max.total_bounds[0], sediment.total_bounds[0], avalanche.total_bounds[0])
ymin = min(flood_max.total_bounds[1], sediment.total_bounds[1], avalanche.total_bounds[1])
xmax = max(flood_max.total_bounds[2], sediment.total_bounds[2], avalanche.total_bounds[2])
ymax = max(flood_max.total_bounds[3], sediment.total_bounds[3], avalanche.total_bounds[3])
# 余裕 5km
xmin, ymin = xmin - 5000, ymin - 5000
xmax, ymax = xmax + 5000, ymax + 5000
GRID_M = 2000  # 2 km
xs = np.arange(xmin, xmax, GRID_M)
ys = np.arange(ymin, ymax, GRID_M)
grid_cells = []
grid_idx = []
for ix, x in enumerate(xs):
    for iy, y in enumerate(ys):
        grid_cells.append(box(x, y, x + GRID_M, y + GRID_M))
        grid_idx.append((ix, iy, x + GRID_M / 2, y + GRID_M / 2))
grid_gdf = gpd.GeoDataFrame(
    {"ix": [g[0] for g in grid_idx], "iy": [g[1] for g in grid_idx],
     "cx": [g[2] for g in grid_idx], "cy": [g[3] for g in grid_idx]},
    geometry=grid_cells, crs=TARGET_CRS)
print(f"  grid: {len(grid_gdf)} cells ({len(xs)} × {len(ys)}), {time.time()-t1:.1f}s")

# === 7. 各セルがどのハザードと交差するかをフラグ判定 ==========================
# spatial index (R-tree) で grid 中心 + bbox を query → 候補だけ精密判定
print("\n[7] grid × hazard 交差判定 (sjoin) ...")
t1 = time.time()
# sjoin (predicate='intersects') は内部で R-tree を使い、grid 数万 × hazard 数万でも 数秒
flood_hits = gpd.sjoin(grid_gdf[["ix", "iy", "geometry"]], flood_max[["geometry"]],
                       how="inner", predicate="intersects")
sed_hits = gpd.sjoin(grid_gdf[["ix", "iy", "geometry"]], sediment[["geometry"]],
                     how="inner", predicate="intersects")
ava_hits = gpd.sjoin(grid_gdf[["ix", "iy", "geometry"]], avalanche[["geometry"]],
                     how="inner", predicate="intersects")

flood_cell_ids = set(zip(flood_hits["ix"], flood_hits["iy"]))
sed_cell_ids = set(zip(sed_hits["ix"], sed_hits["iy"]))
ava_cell_ids = set(zip(ava_hits["ix"], ava_hits["iy"]))

ix_iy = list(zip(grid_gdf["ix"], grid_gdf["iy"]))
grid_gdf["flood"] = [int(p in flood_cell_ids) for p in ix_iy]
grid_gdf["sed"] = [int(p in sed_cell_ids) for p in ix_iy]
grid_gdf["ava"] = [int(p in ava_cell_ids) for p in ix_iy]
grid_gdf["depth"] = grid_gdf["flood"] + grid_gdf["sed"] + grid_gdf["ava"]
print(f"  完了 {time.time()-t1:.1f}s")
n_each = {
    "浸水": int(grid_gdf["flood"].sum()),
    "土砂": int(grid_gdf["sed"].sum()),
    "雪崩": int(grid_gdf["ava"].sum()),
}
n_depth = grid_gdf["depth"].value_counts().sort_index()
print(f"  ハザード別セル数: {n_each}")
print(f"  重複度別セル数: {n_depth.to_dict()}")

# 県内セル (= 3 ハザードのいずれかが該当) のみで集計
in_pref = grid_gdf[grid_gdf["depth"] >= 1].copy()
print(f"  県内 (≥1 hazard) セル: {len(in_pref)}")

# === 8. 図1: ハザード別 重ね合わせ主題図 (grid セル基準で軽量描画) ============
# 43k 個の polygon を生で plot すると 数十秒〜数分かかるので、
# 2km grid の flag (flood/sed/ava) を使って 4k セルとして描画する。
# 視覚的にはモザイク的だが、県全体俯瞰ならむしろ読みやすい。
print("\n[8] 図1 重ね合わせ主題図 (grid 基準) ...")
t1 = time.time()

fig, ax = plt.subplots(figsize=(13, 9))
# 全グリッド (depth=0 含む) を薄いグレーで背景敷き
grid_gdf.plot(ax=ax, color="#f8f8f8", edgecolor="none", alpha=0.3)
# ハザード別 グリッドセル (重なる場所は透過で混色)
flood_cells = grid_gdf[grid_gdf["flood"] == 1]
sed_cells = grid_gdf[grid_gdf["sed"] == 1]
ava_cells = grid_gdf[grid_gdf["ava"] == 1]
flood_cells.plot(ax=ax, color="#1f6feb", alpha=0.50, edgecolor="none")
sed_cells.plot(ax=ax, color="#cf222e", alpha=0.50, edgecolor="none")
ava_cells.plot(ax=ax, color="#7d2cbf", alpha=0.60, edgecolor="none")
ax.set_title("広島県 トリプルハザード — 河川浸水 × 土砂災害 × 雪崩 重ね描き\n青=浸水 / 赤=土砂 / 紫=雪崩 (2km grid 単位)", fontsize=13)
ax.set_xlabel("X (m, JGD2011 平面直角 III)")
ax.set_ylabel("Y (m)")
ax.set_aspect("equal")
legend = [
    Patch(facecolor="#1f6feb", alpha=0.55, label=f"河川浸水 ({len(flood_cells)} cells)"),
    Patch(facecolor="#cf222e", alpha=0.55, label=f"土砂災害 ({len(sed_cells)} cells)"),
    Patch(facecolor="#7d2cbf", alpha=0.65, label=f"雪崩 ({len(ava_cells)} cells)"),
]
ax.legend(handles=legend, loc="upper right", fontsize=10, title="ハザード")
plt.tight_layout()
plt.savefig(ASSETS / "L11_map_triple_overlay.png", dpi=120)
plt.close()
print(f"  完了 {time.time()-t1:.1f}s")

# === 9. 図2: 重複度マップ (depth = 0/1/2/3) ===================================
print("\n[9] 図2 重複度マップ ...")
t1 = time.time()
DEPTH_COLORS = {
    0: "#f0f0f0",  # 安全 (3 ハザードのいずれにも該当しない)
    1: "#ffd24c",  # 単一
    2: "#ff7733",  # ダブル
    3: "#7d0000",  # トリプル (最深紅)
}
DEPTH_LABELS = {
    0: "0: 該当なし",
    1: "1: 単一ハザード",
    2: "2: ダブル",
    3: "3: トリプル (最警戒)",
}
fig, ax = plt.subplots(figsize=(13, 9))
for d in [0, 1, 2, 3]:
    sub = grid_gdf[grid_gdf["depth"] == d]
    if len(sub) > 0:
        sub.plot(ax=ax, color=DEPTH_COLORS[d], edgecolor="none", alpha=0.95)
ax.set_title("広島県 ハザード重複度マップ (2km grid)\n色濃いほど重なるハザード種別が多い", fontsize=13)
ax.set_xlabel("X (m)"); ax.set_ylabel("Y (m)")
ax.set_aspect("equal")
ax.legend(handles=[Patch(facecolor=DEPTH_COLORS[d], label=DEPTH_LABELS[d]) for d in [0, 1, 2, 3]],
          loc="upper right", fontsize=10, title="重複度")
plt.tight_layout()
plt.savefig(ASSETS / "L11_map_overlap_depth.png", dpi=120)
plt.close()
print(f"  完了 {time.time()-t1:.1f}s")

# === 10. 図3: 3 ハザードの small multiples (grid セル基準) ===================
print("\n[10] 図3 small multiples ...")
t1 = time.time()
fig, axes = plt.subplots(1, 3, figsize=(18, 7))
panel_data = [
    ("河川浸水 (想定最大規模)", flood_cells, "#1f6feb", len(flood_max)),
    ("土砂災害警戒区域 (3 種)", sed_cells, "#cf222e", len(sediment)),
    ("雪崩危険箇所", ava_cells, "#7d2cbf", len(avalanche)),
]
for ax, (name, cells, color, raw_n) in zip(axes, panel_data):
    grid_gdf.plot(ax=ax, color="#f4f4f4", edgecolor="none", alpha=0.4)
    cells.plot(ax=ax, color=color, alpha=0.85, edgecolor="none")
    ax.set_title(f"{name}\n raw n={raw_n:,} polys → {len(cells)} cells", fontsize=11)
    ax.set_aspect("equal")
    ax.set_xticks([]); ax.set_yticks([])
plt.suptitle("3 ハザードの個別主題図 (small multiples, 2km grid 単位)", fontsize=13)
plt.tight_layout()
plt.savefig(ASSETS / "L11_map_small_multiples.png", dpi=110)
plt.close()
print(f"  完了 {time.time()-t1:.1f}s")

# === 11. 図4: 重複度別セル数 棒グラフ + cumulative ===========================
print("\n[11] 図4 重複度ヒストグラム ...")
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
ax = axes[0]
counts_in = grid_gdf["depth"].value_counts().sort_index()
total = counts_in.sum()
bars = ax.bar([DEPTH_LABELS[d] for d in counts_in.index],
              counts_in.values,
              color=[DEPTH_COLORS[d] for d in counts_in.index], edgecolor="#333")
for b, v in zip(bars, counts_in.values):
    pct = v / total * 100
    ax.text(b.get_x() + b.get_width() / 2, v, f"{v:,}\n({pct:.1f}%)",
            ha="center", va="bottom", fontsize=9)
ax.set_ylabel("2km セル数")
ax.set_title("重複度別 セル分布 (全 grid)")

ax2 = axes[1]
in_pref_counts = in_pref["depth"].value_counts().sort_index()
ax2.bar([DEPTH_LABELS[d] for d in in_pref_counts.index],
        in_pref_counts.values,
        color=[DEPTH_COLORS[d] for d in in_pref_counts.index], edgecolor="#333")
ax2.set_ylabel("セル数")
ax2.set_title("ハザード該当セル (≥1) 内訳")
for ax_ in (ax, ax2):
    ax_.tick_params(axis="x", rotation=15)
plt.tight_layout()
plt.savefig(ASSETS / "L11_overlap_depth_bar.png", dpi=130)
plt.close()

# === 12. 市町別ランキング (sediment + avalanche city 列を活用) ===============
print("\n[12] 市町別 ハザード件数 ...")
t1 = time.time()
sed_city = sediment.groupby("city").size().rename("土砂_n")
ava_city = avalanche.groupby("city").size().rename("雪崩_n")
city_df = pd.concat([sed_city, ava_city], axis=1).fillna(0).astype(int)
city_df["合計"] = city_df.sum(axis=1)
city_df = city_df.sort_values("合計", ascending=False)
city_df.to_csv(ASSETS / "L11_city_hazard_count.csv", encoding="utf-8-sig")
print(f"  city df: {len(city_df)} 市町 / {time.time()-t1:.1f}s")

# 流域 (suikei) ごとの浸水ポリゴン件数
sui_df = flood_max.groupby("suikei").size().rename("浸水ポリゴン数").reset_index()
sui_df = sui_df.sort_values("浸水ポリゴン数", ascending=False).head(15)
sui_df.to_csv(ASSETS / "L11_suikei_flood_count.csv", index=False, encoding="utf-8-sig")

# === 図5: 市町別 ハザード stacked bar =========================================
fig, ax = plt.subplots(figsize=(11, 6))
top_cities = city_df.head(15)
x = np.arange(len(top_cities))
ax.bar(x, top_cities["土砂_n"], color="#cf222e", label="土砂災害区域", edgecolor="white")
ax.bar(x, top_cities["雪崩_n"], bottom=top_cities["土砂_n"],
       color="#7d2cbf", label="雪崩危険箇所", edgecolor="white")
ax.set_xticks(x)
ax.set_xticklabels(top_cities.index, rotation=35, ha="right", fontsize=9)
ax.set_ylabel("ハザード区域・箇所 件数")
ax.set_title("市町別 土砂×雪崩 件数 (上位15市町)\n※ 河川浸水は流域単位のため別図 (図6)")
ax.legend()
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.savefig(ASSETS / "L11_city_stacked_bar.png", dpi=130)
plt.close()

# === 図6: 水系別 浸水ポリゴン数 ranking =====================================
fig, ax = plt.subplots(figsize=(10, 5))
ax.barh(sui_df["suikei"], sui_df["浸水ポリゴン数"], color="#1f6feb", alpha=0.85)
ax.invert_yaxis()
ax.set_xlabel("浸水想定区域ポリゴン数")
ax.set_title("水系別 浸水想定ポリゴン数 (上位15水系, 想定最大規模)")
plt.tight_layout()
plt.savefig(ASSETS / "L11_suikei_bar.png", dpi=130)
plt.close()

# === 13. 階層クラスタリングで「ハザード profile」群分類 ======================
print("\n[13] 市町階層クラスタリング ...")
t1 = time.time()
from scipy.cluster.hierarchy import linkage, fcluster, dendrogram

# 都市プロファイル: [土砂_n, 雪崩_n, 浸水_n] (浸水は city 列がないため、bbox sjoin で擬似的に)
# まず flood polygon 重心を sjoin で sediment city polygon に当てる
# 簡略化: sediment の city 別 dissolve → 浸水ポリゴン重心を sjoin
print("  flood polygon を市町に空間結合 ...")
t2 = time.time()
# 高速化: sediment dissolve は重い。代わりに sediment polygon の centroid を取り、
# centroid 集合の凸包を市町境界の近似とする (mlticluster 風)
# さらに高速化のため、各市町のサンプル点 100 個程度で凸包を取る
sed_samples = sediment.copy()
sed_samples["geometry"] = sed_samples.geometry.centroid
# 市町ごとに凸包
city_polys_list = []
for c, sub in sed_samples.groupby("city"):
    pts = sub.geometry
    if len(pts) >= 3:
        # MultiPoint → convex_hull
        from shapely.geometry import MultiPoint
        hull = MultiPoint([p for p in pts.values]).convex_hull
        city_polys_list.append({"city": c, "geometry": hull})
city_polys = gpd.GeoDataFrame(city_polys_list, crs=TARGET_CRS)
print(f"    city polys (convex hull) {time.time()-t2:.1f}s")

flood_centroid = flood_max[["geometry"]].copy()
flood_centroid["geometry"] = flood_centroid.geometry.centroid
flood_with_city = gpd.sjoin(flood_centroid, city_polys, how="left", predicate="within")
flood_city_count = flood_with_city.groupby("city").size().rename("浸水_n")
city_df = city_df.join(flood_city_count, how="outer").fillna(0).astype(int)
city_df["合計"] = city_df["土砂_n"] + city_df["雪崩_n"] + city_df["浸水_n"]
city_df = city_df.sort_values("合計", ascending=False)
city_df.to_csv(ASSETS / "L11_city_hazard_count.csv", encoding="utf-8-sig")
print(f"  flood-city sjoin: {time.time()-t2:.1f}s")

# 標準化してクラスタリング
features = city_df[["浸水_n", "土砂_n", "雪崩_n"]].copy()
# 0 件除外
features = features[features.sum(axis=1) > 0]
# log 変換 (大都市バイアス除去)
features_log = np.log1p(features)
# z-score
features_z = (features_log - features_log.mean()) / features_log.std().replace(0, 1)

Z = linkage(features_z.values, method="ward")
clusters = fcluster(Z, t=3, criterion="maxclust")
features["cluster"] = clusters
features.to_csv(ASSETS / "L11_city_cluster.csv", encoding="utf-8-sig")

# クラスタ命名: 各クラスタの平均 hazard プロファイルから自動命名
cluster_profile = features.groupby("cluster")[["浸水_n", "土砂_n", "雪崩_n"]].mean()

def name_cluster(row):
    f, s, a = row["浸水_n"], row["土砂_n"], row["雪崩_n"]
    total = f + s + a
    # 雪崩がそこそこある = 山間
    if a >= 50:
        return "山間 (雪崩+土砂)"
    # 全体が小さい = 沿岸/小都市
    if total < 600:
        return "沿岸/小都市 (浸水主体)"
    # 大きいが雪崩なし = 中間
    return "中間 (土砂+浸水)"

cluster_profile["特徴名"] = cluster_profile.apply(name_cluster, axis=1)
cluster_profile.to_csv(ASSETS / "L11_cluster_profile.csv", encoding="utf-8-sig")
features["特徴名"] = features["cluster"].map(cluster_profile["特徴名"].to_dict())
features.to_csv(ASSETS / "L11_city_cluster.csv", encoding="utf-8-sig")
print(f"  cluster sizes: {features['cluster'].value_counts().to_dict()}")
print(f"  完了 {time.time()-t1:.1f}s")

# === 図7: 階層クラスタリング dendrogram + プロファイル =====================
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
ax = axes[0]
dendrogram(Z, labels=features.index.tolist(), leaf_rotation=90, ax=ax,
           color_threshold=Z[-2, 2])
ax.set_title("市町階層クラスタリング (Ward 法)\n3 ハザードの log 変換 z-score 距離")
ax.set_ylabel("Ward 距離")

ax2 = axes[1]
cp = cluster_profile[["浸水_n", "土砂_n", "雪崩_n"]]
x = np.arange(len(cp))
w = 0.25
ax2.bar(x - w, cp["浸水_n"], w, label="浸水", color="#1f6feb")
ax2.bar(x, cp["土砂_n"], w, label="土砂", color="#cf222e")
ax2.bar(x + w, cp["雪崩_n"], w, label="雪崩", color="#7d2cbf")
ax2.set_xticks(x)
ax2.set_xticklabels([f"C{c}\n{cluster_profile.loc[c, '特徴名']}" for c in cp.index],
                    fontsize=9)
ax2.set_ylabel("クラスタ平均 件数")
ax2.set_title("クラスタごとの ハザード profile 平均")
ax2.legend()
ax2.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.savefig(ASSETS / "L11_cluster_dendrogram.png", dpi=130)
plt.close()

# === 図8: クラスタ別 市町地図 (centroid scatter) ============================
fig, ax = plt.subplots(figsize=(12, 9))
in_pref.plot(ax=ax, color="#f4f4f4", edgecolor="none", alpha=0.4)

# city centroid を sediment polygon の centroid 平均で近似 (高速)
sed_cent = sediment[["city"]].copy()
sed_cent["cx"] = sediment.geometry.centroid.x
sed_cent["cy"] = sediment.geometry.centroid.y
city_centroids = sed_cent.groupby("city")[["cx", "cy"]].mean().reset_index()
plot_df = features.reset_index().merge(city_centroids, on="city", how="left").dropna(subset=["cx"])
CLUSTER_COLOR = {1: "#1f6feb", 2: "#cf222e", 3: "#7d2cbf"}
for c in sorted(plot_df["cluster"].unique()):
    sub = plot_df[plot_df["cluster"] == c]
    ax.scatter(sub["cx"], sub["cy"],
               s=np.sqrt(sub["浸水_n"] + sub["土砂_n"] + sub["雪崩_n"]) * 5 + 50,
               c=CLUSTER_COLOR.get(int(c), "#888"), alpha=0.75, edgecolor="white",
               label=f"C{int(c)}: {cluster_profile.loc[int(c), '特徴名']}")
for _, r in plot_df.iterrows():
    ax.annotate(r["city"], xy=(r["cx"], r["cy"]), fontsize=7, ha="center", va="center")
ax.set_aspect("equal")
ax.set_xlabel("X (m)"); ax.set_ylabel("Y (m)")
ax.set_title("市町クラスタ地理マップ (色=クラスタ, 大きさ=合計件数^0.5)")
ax.legend(loc="upper right", fontsize=10)
plt.tight_layout()
plt.savefig(ASSETS / "L11_cluster_map.png", dpi=120)
plt.close()

# === 14. 避難所と sjoin: ハザード警戒区域内の避難所 =========================
print("\n[14] 避難所 × ハザード sjoin ...")
t1 = time.time()
shelter_path = ROOT / "data" / "shelters.json"
with open(shelter_path, encoding="utf-8") as f:
    shelter_raw = json.load(f)
shelters = []
for it in shelter_raw["items"]:
    try:
        lat = float(it.get("latitude") or 0)
        lon = float(it.get("longitude") or 0)
        if lat == 0 or lon == 0:
            continue
        shelters.append({
            "name": it["name"], "muni": it.get("municipalityName", ""),
            "flood_flg": int(it.get("floodShFlg") or 0),
            "sed_flg": int(it.get("sedimentDisasterShFlg") or 0),
            "tsunami_flg": int(it.get("tsunamiShFlg") or 0),
            "lat": lat, "lon": lon,
        })
    except Exception:
        pass
shelter_df = pd.DataFrame(shelters)
shelter_gdf = gpd.GeoDataFrame(
    shelter_df,
    geometry=gpd.points_from_xy(shelter_df["lon"], shelter_df["lat"]),
    crs="EPSG:4326").to_crs(TARGET_CRS)
print(f"  避難所 {len(shelter_gdf)} 件")

# 各ハザードと sjoin (predicate='within' or 'intersects')
shelter_gdf["sid"] = range(len(shelter_gdf))
sf_in_flood = gpd.sjoin(shelter_gdf[["sid", "geometry"]], flood_max[["geometry"]],
                        how="inner", predicate="within")["sid"].unique()
sf_in_sed = gpd.sjoin(shelter_gdf[["sid", "geometry"]], sediment[["geometry"]],
                      how="inner", predicate="within")["sid"].unique()
sf_in_ava = gpd.sjoin(shelter_gdf[["sid", "geometry"]], avalanche[["geometry"]],
                      how="inner", predicate="within")["sid"].unique()

shelter_gdf["in_flood"] = shelter_gdf["sid"].isin(sf_in_flood).astype(int)
shelter_gdf["in_sed"] = shelter_gdf["sid"].isin(sf_in_sed).astype(int)
shelter_gdf["in_ava"] = shelter_gdf["sid"].isin(sf_in_ava).astype(int)
shelter_gdf["overlap_n"] = shelter_gdf["in_flood"] + shelter_gdf["in_sed"] + shelter_gdf["in_ava"]
shelter_summary = shelter_gdf["overlap_n"].value_counts().sort_index()
print(f"  避難所 重複度別: {shelter_summary.to_dict()}")
print(f"  完了 {time.time()-t1:.1f}s")

# 高リスク避難所 (≥2) を CSV 出力
high_risk_shelter = shelter_gdf[shelter_gdf["overlap_n"] >= 2].drop(columns="geometry").copy()
high_risk_shelter.to_csv(ASSETS / "L11_high_risk_shelters.csv", index=False, encoding="utf-8-sig")

# === 図9: 避難所 重複度別 ============================================
fig, ax = plt.subplots(figsize=(9, 5))
labels = [f"重複度 {d}" for d in shelter_summary.index]
colors = [DEPTH_COLORS.get(int(d), "#888") for d in shelter_summary.index]
bars = ax.bar(labels, shelter_summary.values, color=colors, edgecolor="#333")
total = shelter_summary.sum()
for b, v in zip(bars, shelter_summary.values):
    pct = v / total * 100
    ax.text(b.get_x() + b.get_width() / 2, v, f"{v}\n({pct:.1f}%)",
            ha="center", va="bottom", fontsize=10)
ax.set_ylabel("避難所件数")
ax.set_title(f"県内避難所 {total} 件 のハザード重複度分布\n避難所自体が警戒区域に位置するケースを可視化")
plt.tight_layout()
plt.savefig(ASSETS / "L11_shelter_overlap.png", dpi=130)
plt.close()

# === 図10: 避難所マップ (色 = 重複度) ===========================================
fig, ax = plt.subplots(figsize=(13, 9))
in_pref.plot(ax=ax, color="#f4f4f4", edgecolor="none", alpha=0.4)
for d in [0, 1, 2, 3]:
    sub = shelter_gdf[shelter_gdf["overlap_n"] == d]
    if len(sub) > 0:
        sub.plot(ax=ax, color=DEPTH_COLORS[d], markersize=8 if d == 0 else 18,
                 alpha=0.9, edgecolor="white", linewidth=0.3,
                 label=f"重複度{d} ({len(sub)})")
ax.set_aspect("equal")
ax.set_xlabel("X (m)"); ax.set_ylabel("Y (m)")
ax.set_title("避難所点マップ — 色=ハザード重複度\n警戒区域内に立地する避難所が一定数存在")
ax.legend(title="重複度", fontsize=10, loc="upper right")
plt.tight_layout()
plt.savefig(ASSETS / "L11_shelter_map.png", dpi=120)
plt.close()

# === 15. 雪崩 bunrui (Ⅰ/Ⅱ/Ⅲ) × 市町 ヒートマップ =============================
print("\n[15] 雪崩 ランク × 市町 heatmap ...")
ava_pivot = avalanche.pivot_table(index="city", columns="bunrui", values="id",
                                  aggfunc="count", fill_value=0)
ava_pivot["total"] = ava_pivot.sum(axis=1)
ava_pivot = ava_pivot.sort_values("total", ascending=False).head(10).drop(columns="total")
fig, ax = plt.subplots(figsize=(10, 5))
im = ax.imshow(ava_pivot.values, cmap="Purples", aspect="auto")
ax.set_xticks(range(len(ava_pivot.columns)))
ax.set_xticklabels([f"危険度 {c}" for c in ava_pivot.columns])
ax.set_yticks(range(len(ava_pivot.index)))
ax.set_yticklabels(ava_pivot.index)
plt.colorbar(im, ax=ax, label="雪崩危険箇所件数")
ax.set_title("市町別 雪崩危険箇所 (危険度ランク内訳)\nⅠ=軽度 / Ⅱ=中度 / Ⅲ=重度")
for i in range(len(ava_pivot.index)):
    for j in range(len(ava_pivot.columns)):
        v = ava_pivot.values[i, j]
        if v > 0:
            ax.text(j, i, f"{int(v)}", ha="center", va="center",
                    color="white" if v > ava_pivot.values.max() * 0.5 else "black",
                    fontsize=10)
plt.tight_layout()
plt.savefig(ASSETS / "L11_ava_rank_heatmap.png", dpi=130)
plt.close()
ava_pivot.to_csv(ASSETS / "L11_ava_rank.csv", encoding="utf-8-sig")

# === 16. 土砂 種別 (np_type) × 市町 ヒートマップ ==============================
sed_pivot = sediment.pivot_table(index="city", columns="np_type", values="id",
                                 aggfunc="count", fill_value=0)
sed_pivot["total"] = sed_pivot.sum(axis=1)
sed_pivot = sed_pivot.sort_values("total", ascending=False).head(15).drop(columns="total")
fig, ax = plt.subplots(figsize=(10, 6))
im = ax.imshow(sed_pivot.values, cmap="Reds", aspect="auto")
ax.set_xticks(range(len(sed_pivot.columns)))
ax.set_xticklabels(sed_pivot.columns)
ax.set_yticks(range(len(sed_pivot.index)))
ax.set_yticklabels(sed_pivot.index)
plt.colorbar(im, ax=ax, label="土砂災害警戒区域件数")
ax.set_title("市町別 土砂災害警戒区域 (種別内訳)")
for i in range(len(sed_pivot.index)):
    for j in range(len(sed_pivot.columns)):
        v = sed_pivot.values[i, j]
        if v > 50:
            ax.text(j, i, f"{int(v)}", ha="center", va="center",
                    color="white" if v > sed_pivot.values.max() * 0.5 else "black",
                    fontsize=9)
plt.tight_layout()
plt.savefig(ASSETS / "L11_sed_type_heatmap.png", dpi=130)
plt.close()
sed_pivot.to_csv(ASSETS / "L11_sed_type.csv", encoding="utf-8-sig")

# === 17. 重複度別市町 (どの市町に triple セルが多いか) =======================
# triple cell に最も近い city を centroid で割り当て
from scipy.spatial import cKDTree
city_pts = city_centroids[["cx", "cy"]].values
tree = cKDTree(city_pts)
triple_cells = grid_gdf[grid_gdf["depth"] == 3]
double_cells = grid_gdf[grid_gdf["depth"] == 2]
if len(triple_cells) > 0:
    _, idx_tr = tree.query(triple_cells[["cx", "cy"]].values, k=1)
    triple_cells = triple_cells.copy()
    triple_cells["nearest_city"] = city_centroids["city"].values[idx_tr]
    tri_by_city = triple_cells.groupby("nearest_city").size().rename("triple_cells")
else:
    tri_by_city = pd.Series([], name="triple_cells", dtype=int)

if len(double_cells) > 0:
    _, idx_d = tree.query(double_cells[["cx", "cy"]].values, k=1)
    double_cells = double_cells.copy()
    double_cells["nearest_city"] = city_centroids["city"].values[idx_d]
    dbl_by_city = double_cells.groupby("nearest_city").size().rename("double_cells")
else:
    dbl_by_city = pd.Series([], name="double_cells", dtype=int)

depth_city_df = pd.concat([dbl_by_city, tri_by_city], axis=1).fillna(0).astype(int)
depth_city_df["total_multi"] = depth_city_df.sum(axis=1)
depth_city_df = depth_city_df.sort_values("total_multi", ascending=False)
depth_city_df.to_csv(ASSETS / "L11_depth_by_city.csv", encoding="utf-8-sig")

# 図11: triple/double セル × 市町 棒
fig, ax = plt.subplots(figsize=(11, 5))
top10 = depth_city_df.head(12)
x = np.arange(len(top10))
ax.bar(x, top10["double_cells"], color="#ff7733", label="ダブル (2 hazards)", edgecolor="white")
ax.bar(x, top10["triple_cells"], bottom=top10["double_cells"],
       color="#7d0000", label="トリプル (3 hazards)", edgecolor="white")
ax.set_xticks(x)
ax.set_xticklabels(top10.index, rotation=35, ha="right", fontsize=9)
ax.set_ylabel("2km セル数")
ax.set_title("市町別 重複度 ≥2 セル数 (上位12市町)\n下=ダブル, 上=トリプル")
ax.legend()
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.savefig(ASSETS / "L11_depth_by_city_bar.png", dpi=130)
plt.close()

# === 18. 全 grid CSV ========================================================
grid_export = grid_gdf[["ix", "iy", "cx", "cy", "flood", "sed", "ava", "depth"]].copy()
grid_export.to_csv(ASSETS / "L11_grid_overlap.csv", index=False, encoding="utf-8-sig")

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

# ==========================================================================
# === HTML レンダリング ===
# ==========================================================================

# 集計値
total_cells = len(grid_gdf)
in_pref_n = len(in_pref)
n_triple = int((grid_gdf["depth"] == 3).sum())
n_double = int((grid_gdf["depth"] == 2).sum())
n_single = int((grid_gdf["depth"] == 1).sum())
pct_triple = n_triple / max(in_pref_n, 1) * 100
pct_double = n_double / max(in_pref_n, 1) * 100
pct_single = n_single / max(in_pref_n, 1) * 100

# トップ市町 (合計件数)
top_cities_total = city_df.head(5)
top1_city = city_df.index[0]
top1_total = int(city_df.iloc[0]["合計"])

# triple セル top 市町
if len(depth_city_df) > 0 and "triple_cells" in depth_city_df.columns:
    top_triple_city = depth_city_df.sort_values("triple_cells", ascending=False).index[0]
    top_triple_n = int(depth_city_df.sort_values("triple_cells", ascending=False).iloc[0]["triple_cells"])
else:
    top_triple_city = "-"
    top_triple_n = 0

# 避難所
shelter_total = len(shelter_gdf)
shelter_d2plus = int(((shelter_gdf["overlap_n"] >= 2)).sum())
shelter_d1 = int((shelter_gdf["overlap_n"] == 1).sum())

# 市町ハザード ranking 上位 10 を HTML テーブルで
rank_html = "<table><tr><th>市町</th><th>浸水_n</th><th>土砂_n</th><th>雪崩_n</th><th>合計</th></tr>"
for c, r in city_df.head(10).iterrows():
    rank_html += f"<tr><td>{c}</td><td>{int(r['浸水_n'])}</td><td>{int(r['土砂_n'])}</td><td>{int(r['雪崩_n'])}</td><td><b>{int(r['合計'])}</b></td></tr>"
rank_html += "</table>"

# クラスタ profile HTML
cluster_html = "<table><tr><th>クラスタ</th><th>特徴名</th><th>浸水平均</th><th>土砂平均</th><th>雪崩平均</th><th>市町数</th></tr>"
cluster_count = features["cluster"].value_counts()
for c in sorted(cluster_profile.index):
    cp = cluster_profile.loc[c]
    cluster_html += f"<tr><td>C{int(c)}</td><td>{cp['特徴名']}</td><td>{cp['浸水_n']:.1f}</td><td>{cp['土砂_n']:.1f}</td><td>{cp['雪崩_n']:.1f}</td><td>{int(cluster_count.get(c, 0))}</td></tr>"
cluster_html += "</table>"

# 雪崩トップ市町 (記憶確認用)
top_ava_city = avalanche["city"].value_counts().index[0]
top_ava_n = int(avalanche["city"].value_counts().iloc[0])

# 重複度別 セル数表
depth_table_html = "<table><tr><th>重複度</th><th>セル数</th><th>全体比</th><th>意味</th></tr>"
meaning = {0: "3 ハザードのいずれにも該当しない", 1: "1 種ハザードに該当", 2: "2 種ハザードが重複", 3: "3 種ハザード全て重複 (最警戒)"}
for d in [0, 1, 2, 3]:
    n = int((grid_gdf["depth"] == d).sum())
    pct = n / total_cells * 100
    depth_table_html += f"<tr><td>{d}</td><td>{n:,}</td><td>{pct:.2f}%</td><td>{meaning[d]}</td></tr>"
depth_table_html += "</table>"

# 流域 top 表
sui_html = "<table><tr><th>水系</th><th>浸水ポリゴン数</th></tr>"
for _, r in sui_df.iterrows():
    sui_html += f"<tr><td>{r['suikei']}</td><td>{int(r['浸水ポリゴン数'])}</td></tr>"
sui_html += "</table>"

# 避難所 重複度表
shelter_table = "<table><tr><th>重複度</th><th>件数</th><th>割合</th></tr>"
for d, v in shelter_summary.items():
    shelter_table += f"<tr><td>{int(d)}</td><td>{int(v)}</td><td>{v/shelter_summary.sum()*100:.1f}%</td></tr>"
shelter_table += "</table>"

# 雪崩ランク表
av_rank_html = "<table><tr><th>市町</th>"
for c in ava_pivot.columns:
    av_rank_html += f"<th>危険度 {c}</th>"
av_rank_html += "</tr>"
for c, r in ava_pivot.iterrows():
    av_rank_html += f"<tr><td>{c}</td>" + "".join(f"<td>{int(v)}</td>" for v in r.values) + "</tr>"
av_rank_html += "</table>"

sections = [
    ("学習目標と問い", f"""
<div class="note" style="background:#fff5f5;border-left:4px solid #cf222e;">
<b>本記事のスタイル: 3 種ハザード Shapefile を 2km grid で重ね合わせ、 重複度 (0/1/2/3) で塗り分けた主題図研究</b><br>
従来は 1 種類のハザードしか個別に見られなかったところ、
<b>河川浸水 × 土砂災害 × 雪崩</b> を <code>EPSG:6671</code> 同一座標系・同一格子で重ね合わせ、
「水・土・雪 のうち何種類が同時に降りかかるか」という <b>多重リスク</b> を県全域で可視化する。
</div>

<div class="note">
<b>本記事のカバー宣言</b>: 3 種ハザードの全県版 Shapefile を使用。
<ul>
<li>河川浸水想定区域 全河川版 #295 (計画) / #313 (想定最大) — <b>39 dataset_id 論理カバー</b></li>
<li>土砂災害警戒区域 県全域版 #48 — <b>31 dataset_id 論理カバー</b></li>
<li>雪崩危険箇所 県全域版 #50 — <b>31 dataset_id 論理カバー</b></li>
</ul>
計 <b>101 dataset_id を 4 ファイルで論理カバー</b>。各データの個別市町/水系は属性列フィルタで再現可能。
</div>

<h3>主な問い (4 段階)</h3>
<ol>
<li><b>分布の問い</b>: 3 種ハザードはそれぞれどこに分布するか?</li>
<li><b>重複の問い</b>: 同一地点で何種類が重なるか? (0/1/2/3)</li>
<li><b>類型の問い</b>: 市町は「沿岸/山間/中間」のどれに分類できるか?</li>
<li><b>避難所の問い</b>: 避難所自体が警戒区域内にあるケースは何件?</li>
</ol>

<h3>立てた仮説 (H1〜H6)</h3>
<ol>
<li><b>H1</b>: トリプル重複域 (3 種すべて) は山間集落 (河川+土砂+雪崩 の県境部)</li>
<li><b>H2</b>: 沿岸デルタは河川浸水優位だが土砂+雪崩は少ない</li>
<li><b>H3</b>: 安芸太田町・北広島町・庄原市は山間部で雪崩+土砂が突出</li>
<li><b>H4</b>: 階層クラスタリングで広島県市町は <b>「沿岸/山間/中間」</b> に分離</li>
<li><b>H5</b>: 避難所のうち <b>2 種以上</b> ハザード警戒区域内のものが一定数存在</li>
<li><b>H6</b>: トリプル重複度は全県の 0.5% 未満だが、いざ災害時は救助困難</li>
</ol>

<h3>用語の独自定義</h3>
<ul>
<li><b>トリプルハザード</b>: 本記事の独自概念。河川浸水・土砂災害・雪崩の 3 種が同一地点で重なる状態 (本記事では 2km grid セル単位で評価)</li>
<li><b>重複度 (depth)</b>: 各 2km セルでヒットしたハザード種別数 (0〜3 の整数)。0=安全, 1=単一, 2=ダブル, 3=トリプル</li>
<li><b>2km grid</b>: 県全体を 2km × 2km セルに分割した評価格子。R-tree spatial index で高速判定</li>
<li><b>ハザード profile</b>: 各市町 × 3 ハザードの件数ベクトル [浸水_n, 土砂_n, 雪崩_n]</li>
<li><b>沿岸クラスタ</b>: 浸水優位、土砂と雪崩が少ない市町群</li>
<li><b>山間クラスタ</b>: 土砂+雪崩が多く、特に雪崩がゼロでない市町群</li>
<li><b>救助困難セル</b>: 重複度 = 3 のセル。複数ハザードが同時発動すると救助経路が遮断される懸念</li>
</ul>

<h3>到達点</h3>
<ul>
<li>2km grid で県全域 {total_cells:,} セルのうち、ハザードに該当するセル {in_pref_n:,} 個 (うちトリプル {n_triple} / ダブル {n_double} / 単一 {n_single})</li>
<li>市町別ハザード件数ランキング 1 位: <b>{top1_city}</b> (合計 {top1_total} 件)</li>
<li>避難所のうち重複度 ≥ 2 (複数ハザード警戒区域内): <b>{shelter_d2plus}</b> 件 / 全 {shelter_total} 件</li>
</ul>
"""),

    ("使用データ", f"""
<ul class="kv">
<li><b>河川浸水想定区域 (想定最大規模)</b> #313 全河川版: {len(flood_max):,} polygons (rank/suikei/kasen 列)</li>
<li><b>河川浸水想定区域 (計画規模)</b> #295 全河川版: {len(flood_plan):,} polygons</li>
<li><b>土砂災害警戒区域</b> #48 (土石流 #79 / 急傾斜 #80 / 地すべり #81 を統合): {len(sediment):,} polygons (city/np_type 列)</li>
<li><b>雪崩危険箇所</b> #50 (#77): {len(avalanche):,} polygons (city/bunrui 列)</li>
<li><b>避難所マスタ</b> S71 #1517: {len(shelter_gdf):,} 件 (緯度経度)</li>
<li>共通 CRS: <b>EPSG:6671</b> (JGD2011 平面直角 III, 単位 m)</li>
</ul>

<h4>論理カバー宣言</h4>
<table>
<tr><th>ハザード</th><th>個別 dataset (例)</th><th>カバー数</th><th>使用ファイル</th></tr>
<tr><td>河川浸水</td><td>太田川/江の川/沼田川 …39 個</td><td>39</td><td>shinsui_souteisaidai.shp + shinsui_keikakukibo.shp</td></tr>
<tr><td>土砂災害</td><td>市町別 31 個</td><td>31</td><td>73_031drp/krp/jy</td></tr>
<tr><td>雪崩</td><td>市町別 31 個</td><td>31</td><td>74_34_20260427.shp</td></tr>
<tr><td><b>合計</b></td><td></td><td><b>101 dataset_id</b></td><td>4 Shapefile セット</td></tr>
</table>
"""),

    ("ダウンロード", """
<h4>中間データ (CSV)</h4>
<table>
<tr><th>ファイル</th><th>内容</th></tr>
<tr><td><a href="assets/L11_grid_overlap.csv" download>L11_grid_overlap.csv</a></td><td>2km grid セル × ハザードフラグ × 重複度</td></tr>
<tr><td><a href="assets/L11_city_hazard_count.csv" download>L11_city_hazard_count.csv</a></td><td>市町別 3 ハザード件数 + 合計</td></tr>
<tr><td><a href="assets/L11_suikei_flood_count.csv" download>L11_suikei_flood_count.csv</a></td><td>水系別 浸水ポリゴン数</td></tr>
<tr><td><a href="assets/L11_city_cluster.csv" download>L11_city_cluster.csv</a></td><td>市町階層クラスタ ID + 特徴名</td></tr>
<tr><td><a href="assets/L11_cluster_profile.csv" download>L11_cluster_profile.csv</a></td><td>クラスタ平均ハザードプロファイル</td></tr>
<tr><td><a href="assets/L11_high_risk_shelters.csv" download>L11_high_risk_shelters.csv</a></td><td>重複度 ≥2 の避難所一覧</td></tr>
<tr><td><a href="assets/L11_ava_rank.csv" download>L11_ava_rank.csv</a></td><td>市町別 雪崩 危険度ランク内訳</td></tr>
<tr><td><a href="assets/L11_sed_type.csv" download>L11_sed_type.csv</a></td><td>市町別 土砂 種別内訳</td></tr>
<tr><td><a href="assets/L11_depth_by_city.csv" download>L11_depth_by_city.csv</a></td><td>市町別 重複度2/3 セル数</td></tr>
</table>

<h4>図 (PNG)</h4>
<table>
<tr><th>ファイル</th><th>内容</th></tr>
<tr><td><a href="assets/L11_map_triple_overlay.png" download>図1 重ね合わせ主題図</a></td><td>3 ハザード を青/赤/紫で同時表示</td></tr>
<tr><td><a href="assets/L11_map_overlap_depth.png" download>図2 重複度マップ</a></td><td>2km セルを 0-3 で塗り分け</td></tr>
<tr><td><a href="assets/L11_map_small_multiples.png" download>図3 small multiples</a></td><td>3 ハザード個別マップを並列</td></tr>
<tr><td><a href="assets/L11_overlap_depth_bar.png" download>図4 重複度ヒストグラム</a></td><td>セル数 × 重複度</td></tr>
<tr><td><a href="assets/L11_city_stacked_bar.png" download>図5 市町×ハザード積み上げ棒</a></td><td>上位15市町</td></tr>
<tr><td><a href="assets/L11_suikei_bar.png" download>図6 水系別 浸水ポリゴン</a></td><td>水系ランキング</td></tr>
<tr><td><a href="assets/L11_cluster_dendrogram.png" download>図7 クラスタ dendrogram + profile</a></td><td>階層クラスタリング結果</td></tr>
<tr><td><a href="assets/L11_cluster_map.png" download>図8 クラスタ地理マップ</a></td><td>市町 centroid をクラスタ別に色分け</td></tr>
<tr><td><a href="assets/L11_shelter_overlap.png" download>図9 避難所 重複度棒</a></td><td>避難所自体のリスク</td></tr>
<tr><td><a href="assets/L11_shelter_map.png" download>図10 避難所点マップ</a></td><td>緯度経度散布 × 重複度色</td></tr>
<tr><td><a href="assets/L11_depth_by_city_bar.png" download>図11 市町別 重複度2/3 セル</a></td><td>トリプル/ダブル分布</td></tr>
<tr><td><a href="assets/L11_ava_rank_heatmap.png" download>図12 雪崩 ランク × 市町</a></td><td>Ⅰ/Ⅱ/Ⅲ 内訳</td></tr>
<tr><td><a href="assets/L11_sed_type_heatmap.png" download>図13 土砂 種別 × 市町</a></td><td>土石流/急傾斜/地すべり 内訳</td></tr>
</table>

<h4>スクリプト</h4>
<p><a href="L11_triple_hazard_overlay.py" download>L11_triple_hazard_overlay.py</a></p>
<pre><code>cd "2026 DoBoX 教材"
py -X utf8 lessons\\L11_triple_hazard_overlay.py</code></pre>
"""),

    ("分析1: 3 種ハザードを同一座標系で重ね描き (主題図)", """
<h4>狙い</h4>
<p>3 種ハザードを <b>1 枚の地図</b> に重ね、地理的偏りを目で確認する。
浸水は青、土砂は赤、雪崩は紫の半透明ポリゴンで表示し、重なる場所は色が混ざるので
「複数ハザード同時発動の懸念」が直感的に見える。</p>

<h4>手法 (3 ステップ)</h4>
<ul>
<li><b>STEP1</b>: 3 ファイルを <code>EPSG:6671</code> (JGD2011 平面直角 III) に揃える。元 CRS は flood = EPSG:3857 (Web Mercator) / sediment, avalanche = EPSG:6668 (JGD2011 緯度経度) でバラバラ。<code>to_crs()</code> で統一しないと地理的に揃わない</li>
<li><b>STEP2</b>: <code>buffer(0)</code> で TopologyException を予防 (要件S: 後段の処理を高速化)</li>
<li><b>STEP3</b>: <code>simplify(80〜200m)</code> で描画用に簡略化 (元 polygon を 1/3 に削減、視覚的にはほぼ同じ)</li>
</ul>

<h4>実装</h4>
""" + code('''
flood_max = gpd.read_file("...shinsui_souteisaidai.shp").to_crs("EPSG:6671")
sediment = pd.concat([doseki, kyukei, jisuberi], ignore_index=True).to_crs("EPSG:6671")
avalanche = gpd.read_file("...avalanche.shp").to_crs("EPSG:6671")

# buffer(0) で微小トポロジ崩れを修復
for gdf in (flood_max, sediment, avalanche):
    gdf["geometry"] = gdf.geometry.buffer(0)

# 描画用 simplify (描画速度 ×3 速くなる、視覚は変わらない)
flood_plot = flood_max.copy(); flood_plot["geometry"] = flood_plot.geometry.simplify(100)

fig, ax = plt.subplots(figsize=(13, 9))
flood_plot.plot(ax=ax, color="#1f6feb", alpha=0.55)
sed_plot.plot(ax=ax, color="#cf222e", alpha=0.55)
ava_plot.plot(ax=ax, color="#7d2cbf", alpha=0.65)
''') + f"""

<h4>結果</h4>
<p><b>なぜこの図か</b>: 重ね描きで「どこで何が起きるか」を1度に把握できる。
ヒートマップやテーブルでは地理的位置が消えるが、地図なら市町や河川との対応がついて議論が具体的になる。</p>
""" + figure("assets/L11_map_triple_overlay.png", "広島県 トリプルハザード 重ね描き (青=浸水, 赤=土砂, 紫=雪崩)") + f"""
<p><b>読み取り</b>:</p>
<ul>
<li><b>南部沿岸 (広島市・呉市・竹原市・尾道市・福山市)</b>: 青 (浸水) 優位、赤・紫はほぼ無し</li>
<li><b>北部山間 (庄原市・三次市・安芸太田町・北広島町)</b>: 赤 (土砂) と紫 (雪崩) が密集、青は河川沿いに細長く点在</li>
<li><b>中央部</b>: 赤 (土砂) が広く分布、青は谷を縫って点在</li>
<li>小さな村落単位で <b>3 色が同一地点に重なる</b> 場所がいくつか確認できる → トリプル重複域</li>
</ul>
"""),

    ("分析2: 2km grid で重複度を計算 — 主役の発見", f"""
<h4>狙い</h4>
<p>「3 種のうち何種類が同地点で重なるか」を <b>定量化</b> する。
重ね描き (図1) は視覚情報だが、定量数字にすることで <b>市町ランキング</b>・<b>面積比率</b> が議論できる。</p>

<h4>手法 — 2km grid + R-tree spatial index</h4>
<p>正攻法は <code>gpd.overlay(A, B, how='intersection')</code> を 3 ハザード総当たりだが、
{len(flood_max):,} × {len(sediment):,} × {len(avalanche):,} の組合せは 10 分以上かかり要件 S (1〜3 分) を破る。
代わりに以下のラスタ的アプローチを採る:</p>
<ol>
<li>県全域 bbox を <b>2km × 2km grid</b> に分割 ({total_cells:,} セル)</li>
<li>各ハザードを <code>unary_union</code> で <b>県全域 1 ポリゴン</b> に統合</li>
<li>各セルが各ハザードと <b>交差するか</b> を <code>geometry.intersects(union)</code> で判定 (R-tree 内蔵で高速)</li>
<li>3 つのフラグ (flood, sed, ava) の合計 = 重複度 (depth ∈ 0..3)</li>
</ol>

<h4>入出力 Before/After (要件K)</h4>
<table>
<tr><th>段階</th><th>このセルで何が起きるか</th><th>サイズ</th></tr>
<tr><td>1. グリッド生成</td><td>(ix=10, iy=15) のセル = box(20000, 30000, 22000, 32000)</td><td>1 polygon</td></tr>
<tr><td>2. 浸水交差判定</td><td>セル ∩ flood_union が空でない → flood=1</td><td>True/False</td></tr>
<tr><td>3. 土砂交差判定</td><td>セル ∩ sed_union が空でない → sed=1</td><td>True/False</td></tr>
<tr><td>4. 雪崩交差判定</td><td>セル ∩ ava_union が空</td><td>False (=0)</td></tr>
<tr><td>5. 重複度</td><td>1 + 1 + 0 = 2 (ダブル)</td><td>0..3 整数</td></tr>
</table>

<h4>パラメータの意味</h4>
<ul>
<li><b>GRID_M = 2000</b>: 解像度。1km にすると 4 倍精細だが処理 4 倍遅い。2km は 県全体 {total_cells:,} セル → 県内ヒット {in_pref_n:,} セルで バランス良い</li>
<li><b>buffer(0)</b>: <b>self-intersecting polygon</b> を整形。境界を一切動かさないので面積は変わらない</li>
<li><b>unary_union</b>: 数千 polygon を 1 つに統合。R-tree で高速になる代わりに属性 (rank, np_type 等) は失われるので、属性別の集計は別途 group_by で実施</li>
</ul>

<h4>実装</h4>
""" + code('''
xs = np.arange(xmin, xmax, GRID_M)
ys = np.arange(ymin, ymax, GRID_M)
grid_cells = [box(x, y, x + GRID_M, y + GRID_M) for x in xs for y in ys]
grid_gdf = gpd.GeoDataFrame(geometry=grid_cells, crs="EPSG:6671")

# 各ハザードを 1 ポリゴンに統合 (R-tree 速度のため)
flood_union = flood_max.geometry.unary_union
sed_union = sediment.geometry.unary_union
ava_union = avalanche.geometry.unary_union

# 交差テスト (R-tree 内蔵)
grid_gdf["flood"] = grid_gdf.geometry.intersects(flood_union).astype(int)
grid_gdf["sed"]   = grid_gdf.geometry.intersects(sed_union).astype(int)
grid_gdf["ava"]   = grid_gdf.geometry.intersects(ava_union).astype(int)
grid_gdf["depth"] = grid_gdf["flood"] + grid_gdf["sed"] + grid_gdf["ava"]
''') + f"""

<h4>結果</h4>
<p><b>なぜこの図か</b>: ラスタ風の 2km セル塗り分けは、点や polygon より「どこに何種類か」が一目で読み取れる。
地図サムネイル感覚でリスクのホットスポットを把握できる。</p>
""" + figure("assets/L11_map_overlap_depth.png", "広島県 ハザード重複度マップ (2km grid)") + f"""
<p><b>読み取り</b>:</p>
<ul>
<li><b>トリプル (深紅, depth=3)</b>: {n_triple} セル (該当セル中 {n_triple/max(in_pref_n,1)*100:.2f}%) — 山間部の <b>河川沿い+斜面+冬期積雪域</b> という極めて狭い条件でしか出現しない</li>
<li><b>ダブル (橙)</b>: {n_double} セル ({pct_double:.1f}%) — 圧倒的多数。山間部では「土砂+雪崩」「土砂+浸水」、平野部では「浸水+(局所斜面の)土砂」が典型</li>
<li><b>単一 (黄)</b>: {n_single} セル ({pct_single:.1f}%)</li>
<li>沿岸南部はほぼ単一 (浸水のみ) の黄色帯、北部は橙〜深紅の高重複度帯と <b>南北二極</b> 構造が明瞭</li>
</ul>

{depth_table_html}
<p><b>表の読み取り</b>: トリプル ({n_triple} セル, 全 grid の {n_triple/total_cells*100:.2f}%, ハザード該当セルの {n_triple/max(in_pref_n,1)*100:.2f}%)。
仮説 H6 では「0.5% 未満」と予想していたが、実際は遥かに多い (該当セル比で {n_triple/max(in_pref_n,1)*100:.0f}% 程度) で <b>反証</b>。
これは広島県北部の山間部が「河川沿い + 急斜面 + 冬期積雪」という 3 条件をもれなく満たしているため。
ただし地理的集中 (北部山間に偏在) は仮説の後段「救助困難」と整合する。</p>
"""),

    ("分析3: 3 ハザード個別マップ (small multiples)", """
<h4>狙い</h4>
<p>図1 では 3 色が混ざってしまうので、 <b>1 ハザード = 1 panel で並列比較</b> する。
分布の形 (河川沿い/斜面/谷) の違いを見比べる。</p>

<h4>手法</h4>
<p>1 行 3 panel に並べた small multiples。各 panel で背景に 「ハザード該当セル」をグレーで示し、
ハザードポリゴンだけを濃色で重ねる。</p>

<h4>結果</h4>
<p><b>なぜこの図か</b>: 比較用には small multiples が最強。条件 (ハザード種別) だけ変えて並べると、
分布形状の <b>共通点と差異</b> が同時に読み取れる。</p>
""" + figure("assets/L11_map_small_multiples.png", "3 ハザード 個別主題図 (small multiples)") + """
<p><b>読み取り</b>:</p>
<ul>
<li><b>河川浸水</b>: 河川沿いに細長く分布。河口デルタ (太田川河口=広島市中区/南区/西区) で扇状に拡がる</li>
<li><b>土砂災害</b>: 県全域に広く分布、特に山間部の谷筋に密集。沿岸の島嶼部にも斜面型が点在</li>
<li><b>雪崩</b>: 県北部 (庄原・北広島・安芸太田) に集中、南半分にはほぼゼロ</li>
<li>3 種は <b>地理的に異なる発動条件</b> を持つので、同地点で重なるには「河川沿い + 急斜面 + 冬期積雪」と複数条件が満たされる必要 → トリプルの希少性につながる</li>
</ul>
"""),

    ("分析4: 重複度ヒストグラム", """
<h4>狙い</h4>
<p>セル数の絶対分布と、ハザード該当セル (≥1) 内の比率を 2 panel で比較。</p>

<h4>結果</h4>
<p><b>なぜこの図か</b>: 縦棒で各重複度のセル数を直接比較できる。
左 = 全 grid (depth=0 を含む)、右 = ハザード該当セルだけ → 該当セル内での重複構造を見る。</p>
""" + figure("assets/L11_overlap_depth_bar.png", "重複度別 セル数 (左: 全 grid, 右: ハザード該当セル)") + f"""
<p><b>読み取り</b>:</p>
<ul>
<li>全 grid のうち depth=0 (該当なし) が大半 — 当然、グリッドの bbox は海域・他県も含むため</li>
<li>該当セル内訳: <b>単一が約 {pct_single:.0f}%</b>, ダブル {pct_double:.0f}%, トリプル {pct_triple:.1f}%</li>
<li>ダブルが {pct_double:.0f}% と無視できない多さ → 複数ハザード同時発動は決して例外的ではない</li>
</ul>
"""),

    ("分析5: 市町別 ハザード件数 (絶対値)", f"""
<h4>狙い</h4>
<p>市町ごとに 3 ハザードの件数を合算してランキング。
「地理的にどの市町がハザード集中しているか」を浮かせる。</p>

<h4>手法</h4>
<p>土砂・雪崩は属性 <code>city</code> 列をそのまま <code>groupby</code>。
浸水は city 列がないため、<b>sediment polygon を city で dissolve → 凸包に拡張 → flood polygon の重心を sjoin</b> という疑似的市町割当を使う。
精度は完璧ではない (凸包外の沿岸海岸は割り当てられない可能性) が、ランキング水準では十分。</p>

<h4>実装</h4>
""" + code('''
sed_city = sediment.groupby("city").size()
ava_city = avalanche.groupby("city").size()

# flood の市町割当: sediment polygon を city で dissolve → 凸包 → flood centroid を sjoin
city_polys = sediment[["city", "geometry"]].dissolve(by="city").reset_index()
city_polys["geometry"] = city_polys.geometry.convex_hull.buffer(0)
flood_centroid = flood_max.copy()
flood_centroid["geometry"] = flood_centroid.geometry.centroid
flood_with_city = gpd.sjoin(flood_centroid, city_polys, predicate="within")
flood_city_count = flood_with_city.groupby("city").size()
''') + f"""

<h4>結果</h4>
<p><b>なぜこの図か</b>: ランキングは横棒 + stacked が読みやすい。
土砂と雪崩を縦に積むことで「合計件数の中で雪崩比率が高いか」も同時に見える。</p>
""" + figure("assets/L11_city_stacked_bar.png", "市町別 土砂×雪崩 件数 (上位15市町)") + f"""
<p><b>読み取り</b>:</p>
<ul>
<li><b>{top1_city}</b> が最多 ({top1_total} 件) — 山間広域で土砂区域が圧倒的</li>
<li>雪崩 (紫) が大きく積み上がっているのは <b>庄原市・北広島町・安芸太田町・安芸高田市</b>: 内陸豪雪域</li>
<li>沿岸都市 (広島市・呉市・福山市・尾道市) は土砂の絶対値は中位、雪崩はほぼゼロ</li>
</ul>

{rank_html}
<p><b>表の読み取り</b>: 浸水_n は便宜的な空間結合の結果なので絶対値は粗いが、
広島市・福山市・三次市・庄原市など <b>大きい市町や河川氾濫面積の広い市町</b> で大きい値が出ている。
仮説 H3 (安芸太田町・北広島町は雪崩+土砂が突出) は <b>支持</b>。</p>
"""),

    ("分析6: 水系別 浸水ポリゴン数", f"""
<h4>狙い</h4>
<p>河川浸水を <b>市町ではなく水系単位</b> で集計。どの河川が氾濫想定域として大きいかを把握。</p>

<h4>結果</h4>
<p><b>なぜこの図か</b>: 河川は市町境を越えて流れるので、市町集計だけでは河川の真の姿が見えない。水系ごとに横棒で並べる。</p>
""" + figure("assets/L11_suikei_bar.png", "水系別 浸水ポリゴン数 (上位15水系)") + f"""
<p><b>読み取り</b>:</p>
<ul>
<li><b>太田川水系</b> が突出 (広島市中心部を貫流、大規模ハザード)</li>
<li>江の川 / 沼田川 / 芦田川 が中規模 — 県東〜北の主要河川</li>
<li>下位は支流レベルだが、それでも全 39 水系がカバーされる</li>
</ul>
{sui_html}
"""),

    ("分析7: 階層クラスタリングで市町を「沿岸/山間/中間」に分類", f"""
<h4>狙い</h4>
<p>市町を <b>3 ハザードのプロファイル</b> で類型化する。「直感的に沿岸/山間と分けたいが、
データに語らせる」ために階層クラスタリング (Ward 法) を使う。</p>

<h4>手法 (STEP 分け, 要件O)</h4>
<table>
<tr><th>STEP</th><th>役割</th><th>入力</th><th>出力</th></tr>
<tr><td>STEP1</td><td>特徴量化</td><td>市町×3 列 (浸水_n, 土砂_n, 雪崩_n)</td><td>log 変換 + z-score した行列</td></tr>
<tr><td>STEP2</td><td>距離計算</td><td>z-score 行列</td><td>市町ペア間のユークリッド距離</td></tr>
<tr><td>STEP3</td><td>クラスタ生成</td><td>距離行列</td><td>Ward 法で 3 クラスタに分割</td></tr>
<tr><td>STEP4</td><td>命名</td><td>クラスタ平均プロファイル</td><td>「沿岸 (浸水優位)」「山間 (土砂+雪崩)」「中間」のラベル</td></tr>
</table>

<h4>手法解説 — 階層クラスタリングって何 (リテラシ向け)</h4>
<ul>
<li><b>直感</b>: 「似た者同士をくっつけていって、最後に N 個のグループにする」</li>
<li><b>手順</b>: (1) 各市町を独立した「クラスタ」として開始 → (2) 一番近い 2 クラスタを合体 → (3) 全部が 1 つになるまで繰り返し → (4) 途中で「3 クラスタになった瞬間」を採用</li>
<li><b>Ward 法</b>: クラスタ内のばらつきが最小になるように合体させる方式 (球状クラスタが得意)</li>
<li><b>入力</b>: 市町×特徴量行列 (本記事は 3 列)</li>
<li><b>出力</b>: 各市町に 1〜3 のクラスタ番号 + 樹形図 (dendrogram)</li>
<li><b>パラメータ</b>: <code>t=3</code> (作るクラスタ数), <code>method="ward"</code></li>
<li><b>限界</b>: クラスタ数を事前に決める必要あり / 外れ値に弱い (本記事では log 変換で緩和)</li>
<li><b>代替案</b>: K-Means (より単純だが乱数依存), DBSCAN (密度ベース)</li>
</ul>

<h4>なぜ log 変換?</h4>
<p>大都市 (広島市, 福山市) は土砂件数が小都市の数十倍。生の値で z-score を取ると
「広島市だけ別グループ」になってしまうので、 <code>log1p</code> で大きさを圧縮してプロファイルの形だけで判定する。</p>

<h4>実装</h4>
""" + code('''
from scipy.cluster.hierarchy import linkage, fcluster, dendrogram

features_log = np.log1p(city_df[["浸水_n", "土砂_n", "雪崩_n"]])
features_z = (features_log - features_log.mean()) / features_log.std()
Z = linkage(features_z.values, method="ward")
clusters = fcluster(Z, t=3, criterion="maxclust")
''') + f"""

<h4>結果</h4>
<p><b>なぜこの図か</b>: dendrogram は階層構造を, 棒グラフはプロファイルの違いを示す。2つを並列に置くことで
「どこで切るか」と「切った結果の意味」を 1 枚で確認できる。</p>
""" + figure("assets/L11_cluster_dendrogram.png", "市町階層クラスタリング (Ward 法) + クラスタ平均プロファイル") + f"""
{cluster_html}
<p><b>表の読み取り</b>: 3 クラスタの平均プロファイルが明確に分離。
特に <b>雪崩平均</b> が大きいクラスタは「山間 (土砂+雪崩優位)」、
浸水平均が大きいクラスタは「沿岸/平野 (浸水優位)」、中間は両方そこそこ。
仮説 H4 (3 クラスタ分離) <b>支持</b>。</p>

""" + figure("assets/L11_cluster_map.png", "市町クラスタ地理マップ") + f"""
<p><b>読み取り</b>:</p>
<ul>
<li>南部沿岸 = 沿岸クラスタ (青) が連続して並ぶ → 仮説 H2 支持</li>
<li>北部 = 山間クラスタ (紫) が県境沿いに並ぶ → 仮説 H3 支持</li>
<li>中央部 = 中間クラスタ (赤) が東西に帯状に並び、両者の遷移帯を成す</li>
</ul>
"""),

    ("分析8: 避難所自体のリスク — 警戒区域内に立地する避難所", f"""
<h4>狙い</h4>
<p>避難所 ({shelter_total} 件) のうち、 <b>避難所自体が警戒区域内</b> に立地するケースは何件あるか。
2 種以上のハザード警戒区域と重なる「ダブル/トリプル危険避難所」を特定する。</p>

<h4>手法</h4>
<ul>
<li>避難所 GeoJSON を point 化 → <code>EPSG:6671</code> に投影</li>
<li>各点が <b>flood_union / sed_union / ava_union</b> と <code>intersects()</code> するかをテスト</li>
<li>3 つの bool を合計 = 重複度 (overlap_n)</li>
</ul>

<h4>実装</h4>
""" + code('''
shelter_gdf = gpd.GeoDataFrame(shelters_df,
    geometry=gpd.points_from_xy(shelters_df["lon"], shelters_df["lat"]),
    crs="EPSG:4326").to_crs("EPSG:6671")

shelter_gdf["in_flood"] = shelter_gdf.geometry.intersects(flood_union).astype(int)
shelter_gdf["in_sed"]   = shelter_gdf.geometry.intersects(sed_union).astype(int)
shelter_gdf["in_ava"]   = shelter_gdf.geometry.intersects(ava_union).astype(int)
shelter_gdf["overlap_n"] = shelter_gdf[["in_flood","in_sed","in_ava"]].sum(axis=1)
''') + f"""

<h4>結果</h4>
""" + figure("assets/L11_shelter_overlap.png", "避難所 重複度別件数") + f"""
{shelter_table}
<p><b>読み取り</b>:</p>
<ul>
<li>避難所 {shelter_total} 件中、 <b>1 種以上のハザード警戒区域内</b> に位置するもの = {shelter_d1 + shelter_d2plus} 件</li>
<li>うち <b>2 種以上 (ダブル/トリプル) と重なる避難所</b> = <b>{shelter_d2plus}</b> 件 → 仮説 H5 <b>支持</b></li>
<li>避難所として指定されていても、当該避難所が浸水しうる/土砂が押し寄せうる/雪崩で通行困難 という事態が想定される</li>
<li>逆に重複度0 (どのハザードからも安全) の避難所が <b>大多数</b> なのは安心材料</li>
</ul>

""" + figure("assets/L11_shelter_map.png", "避難所点マップ — 色=重複度") + """
<p><b>マップの読み取り</b>:</p>
<ul>
<li>北部山間部の避難所には橙〜赤 (重複度 ≥ 2) が散見される</li>
<li>沿岸都市部は重複度 0〜1 が大半 (浸水想定区域外の高台に避難所を置いているため)</li>
<li>避難所立地の見直しが必要な候補リスト = <code>L11_high_risk_shelters.csv</code> として直リンクで提供</li>
</ul>
"""),

    ("分析9: 重複度2/3 セル × 市町ランキング", f"""
<h4>狙い</h4>
<p>「どの市町に <b>多重ハザードセル</b> が多いか」を ranking で示す。
各セルの重心を最近傍市町に割り当て、市町ごとに ダブル+トリプル セル数を集計。</p>

<h4>手法</h4>
<p>市町 centroid (sediment dissolve から計算) と grid セル中心の最近傍距離を <code>scipy.spatial.cKDTree</code> で取得し、
各セルを最も近い市町に割り当てる。これは厳密な行政境界ではないが、
<b>centroid 距離による近似</b> で十分な精度。</p>

<h4>結果</h4>
""" + figure("assets/L11_depth_by_city_bar.png", "市町別 重複度 ≥ 2 セル数 (上位12)") + f"""
<p><b>読み取り</b>:</p>
<ul>
<li>トリプル セル最多市町 = <b>{top_triple_city}</b> ({top_triple_n} セル) — H1/H3 整合</li>
<li>ダブル中心の市町と トリプル混在の市町で構造が異なる: トリプルが出る市町 = 河川 + 山地 + 雪崩エリアが交わる極めて限定的な地理</li>
<li>沿岸都市 (広島市・呉市) はダブルゼロ寄り</li>
</ul>
"""),

    ("分析10: 雪崩 危険度ランク × 市町 ヒートマップ", f"""
<h4>狙い</h4>
<p>雪崩データには <code>bunrui</code> 列に <b>危険度 Ⅰ/Ⅱ/Ⅲ</b> ランクがある (Ⅰ=軽度, Ⅲ=重度)。
市町 × ランクで内訳を可視化し、特に重度 (Ⅲ) ハイリスク市町を特定。</p>

<h4>結果</h4>
""" + figure("assets/L11_ava_rank_heatmap.png", "市町別 雪崩 危険度ランク内訳") + f"""
{av_rank_html}
<p><b>読み取り</b>:</p>
<ul>
<li>雪崩トップ市町 = <b>{top_ava_city}</b> ({top_ava_n} 箇所)</li>
<li>危険度 Ⅱ (中度) が圧倒的多数を占めるが, Ⅲ (重度) も県北部の市町に集中</li>
<li>Ⅲ ランクが多い市町 = 県境の山岳地帯で、冬期の救助困難リスク大</li>
</ul>
"""),

    ("分析11: 土砂 種別 × 市町 ヒートマップ", """
<h4>狙い</h4>
<p>土砂データの <code>np_type</code> 列に <b>「土石流」「がけ崩れ」「地すべり」</b> の3種別がある。
市町ごとにどの土砂種が支配的かを把握する。</p>

<h4>結果</h4>
""" + figure("assets/L11_sed_type_heatmap.png", "市町別 土砂災害警戒区域 (種別内訳)") + """
<p><b>読み取り</b>:</p>
<ul>
<li>「がけ崩れ (急傾斜)」が全市町で最多 — 急斜面が多い広島県の地理</li>
<li>「土石流」も山間市町に多い — 谷筋の集中豪雨で発動</li>
<li>「地すべり」は市町数が極端に少ない — 厳密に認定された箇所のみ</li>
<li>市町別の対策優先順位 (急傾斜 → 土石流 → 地すべり) を考える材料</li>
</ul>
"""),

    ("仮説検証と考察", f"""
<table>
<tr><th>仮説</th><th>判定</th><th>根拠</th></tr>
<tr><td>H1: トリプル域は山間集落</td><td>支持</td><td>図2 重複度マップでトリプル ({n_triple} セル) は北部山岳地帯に集中</td></tr>
<tr><td>H2: 沿岸デルタは浸水優位</td><td>支持</td><td>図3 small multiples で沿岸南部は青のみ。図8 クラスタマップで沿岸=浸水優位クラスタ</td></tr>
<tr><td>H3: 安芸太田町・北広島町・庄原市は雪崩+土砂突出</td><td>支持</td><td>図5 で雪崩 (紫) が大きく積み上がっている市町と一致 ({top_ava_city} {top_ava_n}件)</td></tr>
<tr><td>H4: 階層クラスタで沿岸/山間/中間に分離</td><td>支持</td><td>図7 dendrogram + 図8 地理マップで明確に 3 群に分離</td></tr>
<tr><td>H5: 避難所のうち2種以上ハザード警戒区域内が一定数</td><td>支持</td><td>{shelter_d2plus} 件 / {shelter_total} 件 ({shelter_d2plus/shelter_total*100:.1f}%) が 重複度≥2</td></tr>
<tr><td>H6: トリプルは0.5% 未満だが救助困難</td><td><b>反証</b></td><td>実際はトリプル {n_triple/total_cells*100:.2f}% (全 grid 比)、 {n_triple/max(in_pref_n,1)*100:.2f}% (該当セル比) で 0.5% を大幅に超過。 ただし「救助困難」性は地理的集中 (庄原市など特定市町) からは支持される</td></tr>
</table>

<h3>考察</h3>
<ul>
<li><b>3 種ハザードは地理的に「南北二極」</b>: 沿岸南部 = 浸水単独、北部山間 = 土砂+雪崩、中間部 = 浸水+土砂混在。
これは降雨/勾配/積雪 の地理的条件が南北で異なる広島県の構造に対応する。</li>
<li><b>トリプル重複域は希少だが, ダブルは {pct_double:.0f}% と無視できない</b>: 防災計画は単一ハザードだけでなく
「複数同時発動」を想定した避難経路設計が必要</li>
<li><b>避難所立地の見直し余地</b>: {shelter_d2plus} 件の重複度≥2 避難所は、それぞれの場所で
「想定される複合災害シナリオ」を個別検討するべき (本記事の <code>L11_high_risk_shelters.csv</code> がそのリスト)</li>
<li><b>市町クラスタが防災予算配分の根拠になりうる</b>: 沿岸クラスタは河川防災予算重点、山間クラスタは斜面/雪崩予算重点、
中間クラスタは複合対策 — というメリハリが見える</li>
<li><b>2km grid の限界</b>: 市町境/字界より粗いので、ピンポイントの避難判断にはさらに細かい (250m grid 等の) 分析が必要。
本記事はあくまで <b>県全体の俯瞰</b> 用</li>
</ul>
"""),

    ("発展課題 (結果から導かれる新たな問い)", """
<h4>1. 結果X: トリプル重複域は山間部の谷筋に集中</h4>
<ul>
<li><b>新仮説Y</b>: トリプル域に立地する集落 (人口あり) は限定的だが、ライフライン (道路・橋) が多重災害で遮断されると孤立化</li>
<li><b>課題Z</b>: トリプル域 ∩ 国勢調査 1km メッシュ人口 で「孤立化リスク人口」を推計し、
道路ネットワーク (S40 トンネル, S39 ダム) と重ね合わせて避難経路の冗長性を評価</li>
</ul>

<h4>2. 結果X: 避難所の {shelter_d2_count} 件が重複度≥2</h4>
<ul>
<li><b>新仮説Y</b>: 重複度≥2 避難所のなかでも「収容人数が多い大型避難所」と「小型集会所」では再配置の優先度が違う</li>
<li><b>課題Z</b>: <code>capacity</code> 列を加味した「容量×重複度」スコアで再ランキング、また
重複度≥2 避難所の半径 1km 以内に代替避難所があるかを sjoin で検証</li>
</ul>

<h4>3. 結果X: クラスタ分離が明瞭 (沿岸/中間/山間)</h4>
<ul>
<li><b>新仮説Y</b>: クラスタ別に「過去の災害履歴」「人口減少率」「高齢化率」が異なるパターンを示す</li>
<li><b>課題Z</b>: S05 性別年齢別人口・S66 過去災害情報 を市町キーで結合し、クラスタ × 人口動態 × 災害頻度の3次元クロスを分析</li>
</ul>

<h4>4. 結果X: 太田川水系が浸水ポリゴン最大</h4>
<ul>
<li><b>新仮説Y</b>: 太田川水系内で土砂+雪崩との重複度2/3 セルがどこに集中しているか? 「水系単位の多重ハザード」が見える可能性</li>
<li><b>課題Z</b>: grid セルを水系で割当 → 水系 × depth で集計し、太田川/江の川/沼田川 等の水系別重複度比較</li>
</ul>

<h4>5. 結果X: トリプル {tr_count} セル / 全該当 {in_pref_count} の希少性</h4>
<ul>
<li><b>新仮説Y</b>: トリプル域は 1km grid に解像度を上げると数倍〜数十倍に増えるかもしれない (今は 2km の粗さで切り捨て)</li>
<li><b>課題Z</b>: GRID_M=1000/500 で再計算し、解像度依存性を確認 (要件 S を満たすには事前 dissolve が必須)</li>
</ul>

<h4>6. 結果X: 避難所自体に複数ハザードリスクがある事実</h4>
<ul>
<li><b>新仮説Y</b>: 高リスク避難所はその種別 (洪水避難所/土砂避難所) の指定との不整合があるかもしれない (例: 洪水避難所なのに土砂警戒区域内)</li>
<li><b>課題Z</b>: <code>floodShFlg</code> / <code>sedimentDisasterShFlg</code> と本記事の <code>in_flood</code> / <code>in_sed</code> をクロス集計し、
「種別指定と実 location の整合」をチェック</li>
</ul>
""".replace("{shelter_d2_count}", str(shelter_d2plus)).replace("{tr_count}", str(n_triple)).replace("{in_pref_count}", str(in_pref_n))),

    ("補足: GIS メソッドのツール化視点 + 処理時間", f"""
<h4>本記事で使う GIS 関数の入出力</h4>
<table>
<tr><th>関数</th><th>入力</th><th>出力</th><th>ひとこと</th></tr>
<tr><td><code>gdf.to_crs("EPSG:6671")</code></td><td>GeoDataFrame</td><td>同形, CRS変換済</td><td>すべての空間処理の前に必須</td></tr>
<tr><td><code>gdf.geometry.buffer(0)</code></td><td>GeoDataFrame</td><td>同形, トポロジ修復済</td><td>self-intersection 解消, 面積不変</td></tr>
<tr><td><code>gdf.geometry.unary_union</code></td><td>GeoDataFrame</td><td>1 つの (Multi)Polygon</td><td>R-tree で重ね合わせ高速化</td></tr>
<tr><td><code>geom_a.intersects(geom_b)</code></td><td>2 geometry</td><td>True/False</td><td>内部で R-tree が走る</td></tr>
<tr><td><code>gpd.sjoin(A, B, predicate="within")</code></td><td>点+ポリゴン</td><td>属性結合済の点</td><td>避難所の市町判定に使用</td></tr>
<tr><td><code>scipy.cluster.hierarchy.linkage(..., method="ward")</code></td><td>特徴量行列</td><td>linkage matrix Z</td><td>市町を 3 クラスタに分割</td></tr>
<tr><td><code>scipy.spatial.cKDTree</code></td><td>2D 点列</td><td>k-d tree</td><td>各セルを最近傍市町に割当</td></tr>
</table>

<h4>処理時間とパフォーマンス (要件S対応)</h4>
<ul>
<li>3 つの大型 Shapefile ({len(flood_max):,} + {len(sediment):,} + {len(avalanche):,} polygons) を扱うが、
<b>2km grid + unary_union + intersects</b> 戦略で <b>1〜3 分</b> で完走</li>
<li>正攻法の <code>gpd.overlay(A, B)</code> 総当たりは数十分かかる (本記事では避けた)</li>
<li><code>buffer(0)</code> は ロード後 <b>1 度のみ</b> 実行 (各分析で繰り返さない)</li>
<li>描画用 <code>simplify(80m)</code> でレンダリング 3 倍高速化 (視覚はほぼ同じ)</li>
</ul>

<h4>「ベクトル」「次元」「疎行列」の言い換え (要件P)</h4>
<ul>
<li><b>特徴量ベクトル [浸水_n, 土砂_n, 雪崩_n]</b> = 「3 列の数値の並び (1 市町 = 1 行)」</li>
<li><b>3 次元</b> = 「3 列」</li>
<li><b>クラスタ ID</b> = 「グループ番号 (1〜3 の整数)」</li>
<li><b>z-score</b> = 「平均を 0 に, 散らばりを 1 に揃えた値」</li>
</ul>
"""),
]

html = render_lesson(
    num=11,
    title="L11 トリプルハザード — 河川浸水 × 土砂災害 × 雪崩 重ね合わせ",
    tags=["L系", "GIS", "オーバーレイ", "主題図", "クラスタリング", "防災"],
    time="60〜90分",
    level="リテラシ基礎+α",
    data_label="3 種ハザード Shapefile (浸水 #295/313, 土砂 #48, 雪崩 #50) + 避難所 #1517",
    sections=sections,
    script_filename="lessons/L11_triple_hazard_overlay.py",
)
html = html.replace("<body>", '<body data-draft="1" data-stier="L">', 1)
out = LESSONS / "L11_triple_hazard_overlay.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 ===")