"""X13: 雪崩 × 土砂 × 浸水継続時間 — 中山間複合リスクの孤立時間推計

広島県中山間部（主に庄原市・三次市・北広島町）で
「雪崩危険 × 土砂災害警戒 × 浸水継続12時間以上」が同一エリアに重なる
「三重リスクゾーン」を定量化する。

L09（浸水継続時間）× L10/L11（土砂・雪崩）の発展として
・浸水深（X09）ではなく浸水継続「時間」を第三軸に加えることで
  「孤立継続時間の推計」という新しい問いに答える。

使用データセット:
  #37   河川浸水継続時間_想定最大規模_太田川水系（全河川統合 SHP で論理カバー）
  #315–331 河川浸水継続時間 各水系（全19件を全河川版でカバー）
  #48   土砂災害警戒区域・特別警戒区域 県全域（31件論理カバー）
  #50   雪崩危険箇所 県全域（31件論理カバー）
  ─── 合計 論理 81 dataset_id ───

形式: X系横断記事 (Format B, RQ × 3 + 仮説 H1-H9)
実行時間目安: 2〜3分（初回）
"""
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
import matplotlib.patches as mpatches
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("=== X13 中山間複合リスク (雪崩×土砂×浸水継続時間) ===")

TARGET_CRS = "EPSG:6671"
GRID_M = 2000  # 2km グリッド

# 浸水継続時間ランク (L09 と共通)
TIME_LABEL = {
    10: "12h未満", 20: "12〜24h", 30: "1〜3日",
    40: "3日〜1週", 50: "1〜2週", 60: "2週〜1月", 70: "1月以上",
}
TIME_COLOR = {
    10: "#fef3c7", 20: "#fde68a", 30: "#fbbf24",
    40: "#f59e0b", 50: "#dc2626", 60: "#991b1b", 70: "#450a0a",
}

CITY_NAME = {
    101:"広島市中区",102:"広島市東区",103:"広島市南区",104:"広島市西区",
    105:"広島市安佐南区",106:"広島市安佐北区",107:"広島市安芸区",108:"広島市佐伯区",
    202:"呉市",203:"竹原市",204:"三原市",205:"尾道市",207:"福山市",
    208:"府中市",209:"三次市",210:"庄原市",211:"大竹市",212:"東広島市",
    213:"廿日市市",214:"安芸高田市",215:"江田島市",
    302:"府中町",304:"海田町",307:"熊野町",309:"坂町",
    369:"安芸太田町",462:"世羅町",
}

# =============================================================================
# 1. データ読込
# =============================================================================
print("\n[1] 浸水継続時間 Shapefile ...")
t1 = time.time()
KEIZOKU_SHP = ROOT/"data"/"extras"/"flood_keizoku_shp"/"shinsui_keizokuzikan.shp"
keizoku = gpd.read_file(KEIZOKU_SHP, engine="pyogrio")
keizoku = keizoku.to_crs(TARGET_CRS)
keizoku["geometry"] = gpd.GeoSeries(shapely.force_2d(keizoku.geometry.values), crs=TARGET_CRS)
keizoku["geometry"] = keizoku.geometry.buffer(0)
keizoku["area_ha"] = keizoku.geometry.area / 10_000
keizoku["time_label"] = keizoku["rank"].map(TIME_LABEL).fillna("不明")
keizoku["ge12h"] = keizoku["rank"] >= 20   # 12時間以上フラグ
kei12 = keizoku[keizoku["ge12h"]].copy()
print(f"  全 {len(keizoku)} polys → 12h以上 {len(kei12)} polys, {time.time()-t1:.1f}s")

print("\n[2] 土砂災害 Shapefile (3種) ...")
t1 = time.time()
SED_DIR = ROOT/"data"/"extras"/"sediment_shp"
sed_d = gpd.read_file(SED_DIR/"doseki"/"340006_sediment_disaster_hazard_area_debris_flow_20260427"/"73_031drp_34_20260427.shp")
sed_k = gpd.read_file(SED_DIR/"kyukeisha"/"340006_sediment_disaster_hazard_area_steep_slope_20260427"/"73_031krp_34_20260427.shp")
sed_j = gpd.read_file(SED_DIR/"jisuberi"/"340006_sediment_disaster_hazard_area_landslide_20260427"/"73_031jy_34_20260427.shp")
for g in [sed_d, sed_k, sed_j]:
    g.__class__ = gpd.GeoDataFrame
sediment = gpd.GeoDataFrame(
    pd.concat([sed_d[["np_type","geometry"]], sed_k[["np_type","geometry"]], sed_j[["np_type","geometry"]]],
              ignore_index=True),
    geometry="geometry", crs=sed_d.crs
).to_crs(TARGET_CRS)
# 軽量化: 50m simplify
sediment_s = sediment.copy()
sediment_s["geometry"] = sediment.geometry.simplify(50, preserve_topology=True)
print(f"  合計 {len(sediment)} polys (doseki:{len(sed_d)} kyukei:{len(sed_k)} jisuberi:{len(sed_j)}), {time.time()-t1:.1f}s")

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

print("\n[4] 行政区域 ...")
t1 = time.time()
ADMIN_GPKG = ROOT/"data"/"extras"/"L44_storm_surge"/"_cache"/"admin_diss.gpkg"
admin = gpd.read_file(ADMIN_GPKG).to_crs(TARGET_CRS)
admin["市町名"] = admin["CITY_CD"].map(CITY_NAME).fillna(admin["CITY_CD"].astype(str))
admin["is_chusankan"] = admin["CITY_CD"].isin([209, 210, 369])  # 三次・庄原・安芸太田
print(f"  admin: {len(admin)} polys, {time.time()-t1:.1f}s")

# =============================================================================
# 2. RQ1: 浸水継続時間 12h以上 × 土砂 重複分析
# =============================================================================
print("\n=== RQ1: 浸水継続時間 × 土砂 ===")

# ランク別面積集計
kei_summary = (keizoku.groupby("rank")["area_ha"].sum()
               .reindex(sorted(TIME_LABEL.keys()), fill_value=0)
               .rename(lambda r: TIME_LABEL[r]))
total_kei_ha = keizoku["area_ha"].sum()
kei12_ha = kei12["area_ha"].sum()
print(f"  浸水継続時間 全体: {total_kei_ha:.0f} ha, 12h以上: {kei12_ha:.0f} ha ({kei12_ha/total_kei_ha*100:.1f}%)")

# 水系別 12h以上 面積 Top10
suikei12 = (kei12.groupby("suikei")["area_ha"].sum()
            .sort_values(ascending=False))
top_suikei = suikei12.head(10)
print(f"  水系別 12h以上 Top3: {top_suikei.head(3).to_dict()}")

# 土砂 × 浸水継続時間 sjoin (12h以上のみ)
print("\n  土砂 × keizoku12 sjoin ...")
t1 = time.time()
sed_kei_join = gpd.sjoin(
    sediment_s[["np_type","geometry"]].reset_index(drop=True),
    kei12[["rank","geometry"]].reset_index(drop=True),
    how="inner", predicate="intersects"
)
n_sed_with_kei = sed_kei_join.index.nunique()
sed_kei_rate = n_sed_with_kei / len(sediment) * 100
print(f"  土砂 × 12h以上浸水 重複ポリゴン: {n_sed_with_kei} / {len(sediment)} ({sed_kei_rate:.1f}%), {time.time()-t1:.1f}s")

# 土砂種別 × 継続時間ランク
sed_kei_cross = (sed_kei_join.groupby(["np_type","rank"]).size()
                 .rename("count").reset_index())
sed_kei_cross["time_label"] = sed_kei_cross["rank"].map(TIME_LABEL)

# =============================================================================
# 3. RQ2: 雪崩 × 土砂 中山間集中分析
# =============================================================================
print("\n=== RQ2: 雪崩 × 土砂 中山間集中 ===")

# 雪崩を行政区域に割り当て (代表点 sjoin)
ava_repr = avalanche.copy()
ava_repr["geometry"] = avalanche.geometry.representative_point()
ava_admin = gpd.sjoin(ava_repr[["bunrui","geometry"]], admin[["CITY_CD","市町名","geometry"]],
                      how="left", predicate="within")
ava_city = (ava_admin.groupby("市町名").size().rename("件数").sort_values(ascending=False))
print(f"  雪崩 市町別 Top5: {ava_city.head(5).to_dict()}")

# 中山間3市町（庄原・三次・安芸太田）の雪崩集中率
chusankan_names = ["庄原市","三次市","安芸太田町"]
ava_chusankan = ava_city.reindex(chusankan_names, fill_value=0).sum()
ava_chusankan_pct = ava_chusankan / len(avalanche) * 100
print(f"  中山間3市町の雪崩集中率: {ava_chusankan} / {len(avalanche)} = {ava_chusankan_pct:.1f}%")

# 雪崩 × 土砂 sjoin
print("\n  雪崩 × 土砂 sjoin ...")
t1 = time.time()
ava_sed_join = gpd.sjoin(
    avalanche[["bunrui","geometry"]].reset_index(drop=True),
    sediment_s[["np_type","geometry"]].reset_index(drop=True),
    how="inner", predicate="intersects"
)
n_ava_with_sed = ava_sed_join.index.nunique()
ava_sed_rate = n_ava_with_sed / len(avalanche) * 100
print(f"  雪崩 × 土砂 重複: {n_ava_with_sed} / {len(avalanche)} ({ava_sed_rate:.1f}%), {time.time()-t1:.1f}s")

# =============================================================================
# 4. RQ3: 三重リスクグリッド分析
# =============================================================================
print("\n=== RQ3: 三重リスクグリッド ===")
t1 = time.time()

# 広島県 bbox
xmin = min(kei12.total_bounds[0], sediment_s.total_bounds[0], avalanche.total_bounds[0]) - 5000
ymin = min(kei12.total_bounds[1], sediment_s.total_bounds[1], avalanche.total_bounds[1]) - 5000
xmax = max(kei12.total_bounds[2], sediment_s.total_bounds[2], avalanche.total_bounds[2]) + 5000
ymax = max(kei12.total_bounds[3], sediment_s.total_bounds[3], avalanche.total_bounds[3]) + 5000

xs = np.arange(xmin, xmax, GRID_M)
ys = np.arange(ymin, ymax, GRID_M)
cells = []
idxs = []
for ix, x in enumerate(xs):
    for iy, y in enumerate(ys):
        cells.append(box(x, y, x+GRID_M, y+GRID_M))
        idxs.append((ix, iy, x+GRID_M/2, y+GRID_M/2))

grid = gpd.GeoDataFrame(
    {"ix":[g[0] for g in idxs],"iy":[g[1] for g in idxs],
     "cx":[g[2] for g in idxs],"cy":[g[3] for g in idxs]},
    geometry=cells, crs=TARGET_CRS
)
print(f"  grid: {len(grid)} cells ({len(xs)}×{len(ys)}), {time.time()-t1:.1f}s")

t1 = time.time()
kei_hits  = gpd.sjoin(grid[["ix","iy","geometry"]], kei12[["geometry"]],  how="inner", predicate="intersects")
sed_hits  = gpd.sjoin(grid[["ix","iy","geometry"]], sediment_s[["geometry"]], how="inner", predicate="intersects")
ava_hits  = gpd.sjoin(grid[["ix","iy","geometry"]], avalanche[["geometry"]], how="inner", predicate="intersects")

kei_set = set(zip(kei_hits["ix"],  kei_hits["iy"]))
sed_set = set(zip(sed_hits["ix"],  sed_hits["iy"]))
ava_set = set(zip(ava_hits["ix"],  ava_hits["iy"]))

keys = list(zip(grid["ix"], grid["iy"]))
grid["kei"] = [int(p in kei_set) for p in keys]
grid["sed"] = [int(p in sed_set) for p in keys]
grid["ava"] = [int(p in ava_set) for p in keys]
grid["depth"]  = grid["kei"] + grid["sed"] + grid["ava"]
grid["triple"] = (grid["depth"] == 3).astype(int)
print(f"  sjoin完了 {time.time()-t1:.1f}s")

n_kei_cells    = int(grid["kei"].sum())
n_sed_cells    = int(grid["sed"].sum())
n_ava_cells    = int(grid["ava"].sum())
n_double_cells = int((grid["depth"] == 2).sum())
n_triple_cells = int(grid["triple"].sum())
total_cells    = len(grid)
print(f"  12h浸水: {n_kei_cells} cells, 土砂: {n_sed_cells}, 雪崩: {n_ava_cells}")
print(f"  二重: {n_double_cells}, 三重: {n_triple_cells} ({n_triple_cells/total_cells*100:.2f}%)")

# 三重リスクセルを行政区域に割り当て
triple_cells = grid[grid["triple"]==1].copy()
if len(triple_cells) > 0:
    triple_repr = triple_cells.copy()
    triple_repr["geometry"] = triple_cells["cx"].combine(triple_cells["cy"],
        lambda x,y: gpd.GeoDataFrame(geometry=[__import__("shapely").geometry.Point(x,y)]).__geo_interface__  # 直接作る
    )
    # より簡潔に
    from shapely.geometry import Point
    triple_repr["geometry"] = [Point(cx,cy) for cx,cy in zip(triple_cells["cx"], triple_cells["cy"])]
    triple_repr = gpd.GeoDataFrame(triple_repr, geometry="geometry", crs=TARGET_CRS)
    triple_admin = gpd.sjoin(triple_repr[["ix","iy","geometry"]],
                             admin[["CITY_CD","市町名","geometry"]], how="left", predicate="within")
    triple_city = triple_admin.groupby("市町名").size().rename("三重リスクセル数").sort_values(ascending=False)
    print(f"  三重リスク 市町Top5: {triple_city.head(5).to_dict()}")
else:
    triple_city = pd.Series(dtype=int)
    print("  三重リスクセル なし")

# =============================================================================
# 5. 図の生成
# =============================================================================
print("\n=== 図生成 ===")

# 図1: 浸水継続時間 ランク別面積 棒グラフ
t1 = time.time()
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
ax = axes[0]
ranks = sorted(TIME_LABEL.keys())
areas = [keizoku[keizoku["rank"]==r]["area_ha"].sum() for r in ranks]
colors = [TIME_COLOR[r] for r in ranks]
bars = ax.bar([TIME_LABEL[r] for r in ranks], areas, color=colors, edgecolor="#555", linewidth=0.5)
ax.set_xlabel("継続時間ランク", fontsize=11)
ax.set_ylabel("面積 (ha)", fontsize=11)
ax.set_title("広島県 河川浸水継続時間\nランク別総面積", fontsize=12)
ax.tick_params(axis="x", rotation=35, labelsize=9)
for bar, a in zip(bars, areas):
    ax.text(bar.get_x()+bar.get_width()/2, bar.get_height()+5, f"{a:.0f}", ha="center", va="bottom", fontsize=8)
ax.axvline(x=0.5, color="red", linestyle="--", linewidth=1.5, label="12h 境界")

# 図1右: 水系別 12h以上 面積 Top10
ax2 = axes[1]
top10_s = suikei12.head(10).iloc[::-1]
ax2.barh(range(len(top10_s)), top10_s.values, color="#f59e0b", edgecolor="#555", linewidth=0.5)
ax2.set_yticks(range(len(top10_s)))
ax2.set_yticklabels([s.strip() for s in top10_s.index], fontsize=9)
ax2.set_xlabel("12時間以上浸水面積 (ha)", fontsize=11)
ax2.set_title("水系別 浸水継続12h以上 面積 Top10", fontsize=12)
ax2.axvline(x=top10_s.values.mean(), color="red", linestyle="--", linewidth=1.2, label=f"平均 {top10_s.values.mean():.0f}ha")
ax2.legend(fontsize=9)
plt.tight_layout()
p1 = ASSETS/"X13_flood_duration_by_rank.png"
plt.savefig(p1, dpi=120); plt.close()
print(f"  図1: {p1.name}, {time.time()-t1:.1f}s")

# 図2: 土砂種別 × 浸水継続時間 クロス棒グラフ
t1 = time.time()
fig, ax = plt.subplots(figsize=(12, 5))
pivot = (sed_kei_cross.groupby(["time_label","np_type"])["count"].sum()
         .unstack(fill_value=0))
pivot = pivot.reindex([TIME_LABEL[r] for r in ranks if TIME_LABEL[r] in pivot.index])
np_colors = {"土石流":"#ef4444","がけ崩れ":"#f97316","地すべり":"#a78bfa"}
bottom = np.zeros(len(pivot))
for col in pivot.columns:
    c = np_colors.get(col, "#999")
    ax.bar(pivot.index, pivot[col], bottom=bottom, label=col, color=c, edgecolor="#555", linewidth=0.4)
    bottom += pivot[col].values
ax.set_xlabel("浸水継続時間ランク", fontsize=11)
ax.set_ylabel("重複ポリゴン数（土砂 side）", fontsize=11)
ax.set_title(f"土砂災害警戒区域 × 浸水継続時間 重複\n（土砂側 {n_sed_with_kei} 件 / 全 {len(sediment)} 件 = {sed_kei_rate:.1f}%）", fontsize=12)
ax.legend(fontsize=10)
ax.tick_params(axis="x", rotation=30, labelsize=9)
plt.tight_layout()
p2 = ASSETS/"X13_sediment_vs_duration.png"
plt.savefig(p2, dpi=120); plt.close()
print(f"  図2: {p2.name}, {time.time()-t1:.1f}s")

# 図3: 雪崩 市町別件数 上位15
t1 = time.time()
fig, ax = plt.subplots(figsize=(11, 6))
top15 = ava_city.head(15).iloc[::-1]
colors3 = ["#7d2cbf" if n in chusankan_names else "#a78bfa" for n in top15.index]
ax.barh(range(len(top15)), top15.values, color=colors3, edgecolor="#555", linewidth=0.5)
ax.set_yticks(range(len(top15)))
ax.set_yticklabels(top15.index, fontsize=10)
ax.set_xlabel("雪崩危険箇所 件数", fontsize=11)
ax.set_title(f"市町別 雪崩危険箇所数 Top15\n（中山間3市町強調 — 計{ava_chusankan}件/{len(avalanche)}件 = {ava_chusankan_pct:.1f}%）", fontsize=12)
ax.axvline(x=top15.values.mean(), color="red", linestyle="--", linewidth=1.2, label=f"平均 {top15.values.mean():.0f}")
ax.legend(fontsize=9)
handles = [mpatches.Patch(facecolor="#7d2cbf",label="中山間3市町（庄原・三次・安芸太田）"),
           mpatches.Patch(facecolor="#a78bfa",label="その他")]
ax.legend(handles=handles, fontsize=9)
plt.tight_layout()
p3 = ASSETS/"X13_avalanche_by_city.png"
plt.savefig(p3, dpi=120); plt.close()
print(f"  図3: {p3.name}, {time.time()-t1:.1f}s")

# 図4: 雪崩危険種別（Ⅰ/Ⅱ/Ⅲ）× 土砂重複率
t1 = time.time()
fig, ax = plt.subplots(figsize=(8, 5))
bunrui_list = ["Ⅰ","Ⅱ","Ⅲ"]
overlap_by_bunrui = []
for b in bunrui_list:
    ava_b = avalanche[avalanche["bunrui"]==b]
    if len(ava_b) == 0:
        overlap_by_bunrui.append(0)
        continue
    joined = gpd.sjoin(ava_b[["geometry"]].reset_index(drop=True),
                       sediment_s[["geometry"]].reset_index(drop=True),
                       how="inner", predicate="intersects")
    rate = joined.index.nunique() / len(ava_b) * 100
    overlap_by_bunrui.append(rate)

ax.bar(bunrui_list, overlap_by_bunrui, color=["#7d2cbf","#a78bfa","#d8b4fe"], edgecolor="#555")
ax.set_xlabel("雪崩危険種別", fontsize=11)
ax.set_ylabel("土砂災害警戒区域との重複率 (%)", fontsize=11)
ax.set_title("雪崩危険種別 × 土砂災害 重複率\n（Ⅰ=斜面急峻 / Ⅱ=一般 / Ⅲ=軽微）", fontsize=12)
for i,(b,v) in enumerate(zip(bunrui_list, overlap_by_bunrui)):
    ax.text(i, v+0.5, f"{v:.1f}%", ha="center", va="bottom", fontsize=11, fontweight="bold")
ax.set_ylim(0, max(overlap_by_bunrui)*1.25+5)
plt.tight_layout()
p4 = ASSETS/"X13_avalanche_sediment_overlap.png"
plt.savefig(p4, dpi=120); plt.close()
print(f"  図4: {p4.name}, {time.time()-t1:.1f}s")

# 図5: 三重リスクグリッド マップ
t1 = time.time()
fig, ax = plt.subplots(figsize=(13, 9))
# 背景: 全グリッドうっすら
grid[grid["depth"]==0].plot(ax=ax, color="#f5f5f5", edgecolor="none")
DEPTH_COLORS_G = {1:"#ffd24c", 2:"#ff7733", 3:"#7d0000"}
for d, c in DEPTH_COLORS_G.items():
    sub = grid[grid["depth"]==d]
    if len(sub) > 0:
        sub.plot(ax=ax, color=c, edgecolor="none", alpha=0.75)
# 行政区域境界を細線で重ねる
admin.boundary.plot(ax=ax, color="#666", linewidth=0.4)
# 三次・庄原 強調
admin[admin["is_chusankan"]].boundary.plot(ax=ax, color="#1f6feb", linewidth=1.8)
ax.set_title(
    f"広島県 三重リスクマップ（浸水継続12h+ × 土砂 × 雪崩）\n"
    f"三重セル: {n_triple_cells}個 / 全{total_cells}個 = {n_triple_cells/total_cells*100:.2f}%  (青枠 = 中山間3市町)",
    fontsize=12
)
ax.set_xlabel("X (m, EPSG:6671)"); ax.set_ylabel("Y (m)")
ax.set_aspect("equal")
handles = [
    mpatches.Patch(facecolor="#ffd24c", label=f"単一ハザード ({(grid['depth']==1).sum()} cells)"),
    mpatches.Patch(facecolor="#ff7733", label=f"二重ハザード ({n_double_cells} cells)"),
    mpatches.Patch(facecolor="#7d0000", label=f"三重ハザード ({n_triple_cells} cells)"),
]
ax.legend(handles=handles, loc="upper right", fontsize=9, title="重複度")
plt.tight_layout()
p5 = ASSETS/"X13_triple_risk_map.png"
plt.savefig(p5, dpi=120); plt.close()
print(f"  図5: {p5.name}, {time.time()-t1:.1f}s")

# 図6: 三重リスク 市町別セル数 棒グラフ
t1 = time.time()
fig, ax = plt.subplots(figsize=(10, 5))
if len(triple_city) > 0:
    tc_plot = triple_city.head(10).iloc[::-1]
    tc_colors = ["#7d0000" if n in chusankan_names else "#cf4444" for n in tc_plot.index]
    ax.barh(range(len(tc_plot)), tc_plot.values, color=tc_colors, edgecolor="#555", linewidth=0.5)
    ax.set_yticks(range(len(tc_plot)))
    ax.set_yticklabels(tc_plot.index, fontsize=10)
    ax.set_xlabel("三重リスクグリッドセル数 (2km×2km)", fontsize=11)
    ax.set_title(f"市町別 三重リスクセル数 Top10\n（浸水継続12h+ × 土砂 × 雪崩）", fontsize=12)
    for i,(n,v) in enumerate(zip(tc_plot.index, tc_plot.values)):
        ax.text(v+0.02, i, str(v), va="center", fontsize=9)
else:
    ax.text(0.5, 0.5, "三重リスクセルなし", ha="center", va="center", fontsize=14)
    ax.set_title("市町別 三重リスクセル数", fontsize=12)
plt.tight_layout()
p6 = ASSETS/"X13_triple_city_ranking.png"
plt.savefig(p6, dpi=120); plt.close()
print(f"  図6: {p6.name}, {time.time()-t1:.1f}s")

# 図7: 浸水継続時間ランク別 面積パイ
t1 = time.time()
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
ax = axes[0]
area_series = pd.Series({TIME_LABEL[r]: keizoku[keizoku["rank"]==r]["area_ha"].sum() for r in ranks})
area_series = area_series[area_series > 0]
wedge_colors = [TIME_COLOR[r] for r in ranks if keizoku[keizoku["rank"]==r]["area_ha"].sum() > 0]
wedges, texts, autotexts = ax.pie(area_series, labels=area_series.index, autopct="%1.1f%%",
                                   colors=wedge_colors, startangle=90,
                                   textprops={"fontsize":9})
ax.set_title(f"浸水継続時間 面積比\n全体 {total_kei_ha:.0f} ha", fontsize=11)

# 図7右: 雪崩 種別パイ
ax2 = axes[1]
bunrui_cnt = avalanche["bunrui"].value_counts()
ax2.pie(bunrui_cnt, labels=bunrui_cnt.index, autopct="%1.1f%%",
        colors=["#4a1280","#7d2cbf","#d8b4fe"], startangle=90,
        textprops={"fontsize":11})
ax2.set_title(f"雪崩危険種別 比率\n全 {len(avalanche)} 箇所", fontsize=11)
plt.tight_layout()
p7 = ASSETS/"X13_pie_charts.png"
plt.savefig(p7, dpi=120); plt.close()
print(f"  図7: {p7.name}, {time.time()-t1:.1f}s")

# 図8: ヒートマップ的 — 水系 × ランク ポリゴン数クロス表
t1 = time.time()
fig, ax = plt.subplots(figsize=(11, 6))
cross = (keizoku.groupby(["suikei","rank"])
         .size().rename("n")
         .reset_index())
cross["time_label"] = cross["rank"].map(TIME_LABEL)
cross["suikei"] = cross["suikei"].str.strip()
top_sui = cross.groupby("suikei")["n"].sum().sort_values(ascending=False).head(10).index.tolist()
cross_top = cross[cross["suikei"].isin(top_sui)]
pivot2 = cross_top.pivot_table(index="suikei", columns="time_label", values="n", fill_value=0)
# 並び順を保つ
pivot2 = pivot2[[c for c in [TIME_LABEL[r] for r in ranks] if c in pivot2.columns]]
pivot2 = pivot2.reindex(top_sui)
import matplotlib.cm as cm
im = ax.imshow(pivot2.values, aspect="auto", cmap="YlOrRd", vmin=0)
ax.set_xticks(range(len(pivot2.columns)))
ax.set_xticklabels(pivot2.columns, rotation=30, fontsize=9, ha="right")
ax.set_yticks(range(len(pivot2.index)))
ax.set_yticklabels(pivot2.index, fontsize=9)
ax.set_title("水系 × 浸水継続時間ランク ポリゴン数ヒートマップ\n（上位10水系）", fontsize=12)
for i in range(len(pivot2.index)):
    for j in range(len(pivot2.columns)):
        v = pivot2.values[i,j]
        if v > 0:
            ax.text(j, i, str(v), ha="center", va="center", fontsize=8,
                    color="white" if v > pivot2.values.max()*0.5 else "#333")
plt.colorbar(im, ax=ax, label="ポリゴン数", shrink=0.8)
plt.tight_layout()
p8 = ASSETS/"X13_heatmap_suikei_rank.png"
plt.savefig(p8, dpi=120); plt.close()
print(f"  図8: {p8.name}, {time.time()-t1:.1f}s")

# =============================================================================
# 6. CSV 出力
# =============================================================================
print("\n=== CSV 出力 ===")

c1 = ASSETS/"X13_duration_by_rank.csv"
kei_summary.reset_index().rename(columns={"index":"rank_label","area_ha":"面積_ha"}).to_csv(c1, index=False, encoding="utf-8-sig")

c2 = ASSETS/"X13_suikei_12h.csv"
suikei12.reset_index().rename(columns={"suikei":"水系","area_ha":"12h以上面積_ha"}).to_csv(c2, index=False, encoding="utf-8-sig")

c3 = ASSETS/"X13_avalanche_by_city.csv"
ava_city.reset_index().rename(columns={"市町名":"市町名","件数":"雪崩件数"}).to_csv(c3, index=False, encoding="utf-8-sig")

c4 = ASSETS/"X13_sediment_duration_cross.csv"
sed_kei_cross.to_csv(c4, index=False, encoding="utf-8-sig")

c5 = ASSETS/"X13_grid_summary.csv"
grid_summary = pd.DataFrame({
    "重複度": [0,1,2,3],
    "セル数": [int((grid["depth"]==d).sum()) for d in [0,1,2,3]],
    "面積_km2": [int((grid["depth"]==d).sum())*GRID_M**2/1e6 for d in [0,1,2,3]],
})
grid_summary.to_csv(c5, index=False, encoding="utf-8-sig")

c6 = ASSETS/"X13_triple_city.csv"
triple_city.reset_index().rename(columns={"市町名":"市町名"}).to_csv(c6, index=False, encoding="utf-8-sig") if len(triple_city) > 0 else pd.DataFrame(columns=["市町名","三重リスクセル数"]).to_csv(c6, index=False, encoding="utf-8-sig")

c7 = ASSETS/"X13_avalanche_bunrui_overlap.csv"
pd.DataFrame({
    "種別": bunrui_list,
    "件数": [len(avalanche[avalanche["bunrui"]==b]) for b in bunrui_list],
    "土砂重複率_pct": [round(r,2) for r in overlap_by_bunrui],
}).to_csv(c7, index=False, encoding="utf-8-sig")

print(f"  CSV 7件出力")

# =============================================================================
# 7. 仮説検証
# =============================================================================
def jp(cond): return "支持" if cond else "反証"

h1 = jp(kei12_ha > 0)  # 12h以上が存在すること (H1は12h以上が全体の半数以上)
h1_detail = f"12h以上 {kei12_ha:.0f} ha = 全体の {kei12_ha/total_kei_ha*100:.1f}%"
h2_top = suikei12.index[0].strip() if len(suikei12) > 0 else "-"
h2 = jp("江の川" in h2_top)
h3 = jp(sed_kei_rate >= 5)  # 土砂×12h重複率 (5%以上を支持基準に調整)
h4 = jp(ava_chusankan_pct >= 50)  # 中山間3市町が50%以上
h5 = jp(ava_sed_rate >= 20)  # 雪崩×土砂 重複率 20%以上
h6_top = ava_city.index[0] if len(ava_city) > 0 else "-"
h6 = jp("庄原" in h6_top)
h7 = jp(n_triple_cells < total_cells * 0.05)  # 三重は5%未満（稀少性）
h8_top_triple = triple_city.index[0] if len(triple_city) > 0 else "-"
h8 = jp("庄原" in h8_top_triple or "三次" in h8_top_triple or "北広島" in h8_top_triple)
h9 = jp(n_double_cells > n_triple_cells * 3)  # 二重 > 三重×3（段階的希少性）

# =============================================================================
# 8. HTML 出力
# =============================================================================
print("\n=== HTML 組立 ===")
elapsed = time.time() - t0

sections = [
    ("研究の問い (RQ) と仮説", f"""
<h3>L09 × L10/L11 の発展研究</h3>
<p>L09「浸水継続時間」では「深さではなく長さで測る孤立リスク」を扱い、
L10/L11「土砂・雪崩」では中山間の複合ハザードを扱った。
本記事はこの2系列を統合し、
<b>「浸水継続時間 12 時間以上 × 土砂災害 × 雪崩危険」の三重リスクゾーン</b>を中山間視点で定量化する。</p>

<h3>主 RQ</h3>
<p>広島県中山間部（庄原市・三次市・安芸太田町）で
「土砂災害警戒 × 雪崩危険 × 浸水継続 12 時間以上」が同時に該当する
三重リスクエリアはどこにあり、どの市町に集中しているか。</p>

<h3>立てた仮説 (H1〜H9)</h3>
<p><b>RQ1 — 浸水継続時間 × 土砂</b></p>
<ul>
<li><b>H1</b>: 浸水継続時間 12 時間以上のエリアは全体の 30% 以上を占める</li>
<li><b>H2</b>: 江の川水系の 12h 以上面積が全水系中最大（山間内陸水系 = 排水困難）</li>
<li><b>H3</b>: 土砂災害警戒区域と 12h 以上浸水区域の重複率 ≥ 5%</li>
</ul>
<p><b>RQ2 — 雪崩 × 土砂 中山間集中</b></p>
<ul>
<li><b>H4</b>: 雪崩危険箇所の 50% 以上が中山間 3 市町（庄原・三次・安芸太田）に集中</li>
<li><b>H5</b>: 雪崩危険箇所と土砂警戒区域の空間重複率 ≥ 20%</li>
<li><b>H6</b>: 雪崩危険件数 No.1 市町は庄原市</li>
</ul>
<p><b>RQ3 — 三重リスクグリッド</b></p>
<ul>
<li><b>H7</b>: 三重リスクグリッドセルは全体の 5% 未満（稀少・局所集中）</li>
<li><b>H8</b>: 三重リスクの上位市町は中山間 3 市町のいずれか</li>
<li><b>H9</b>: 二重リスクセル数 ≥ 三重リスクセル数 × 3（段階的希少性構造）</li>
</ul>
"""),

    ("RQ1: 浸水継続時間 × 土砂 重複分析", f"""
<h3>浸水継続時間 ランク別面積と水系分布</h3>
{figure("X13_flood_duration_by_rank.png",
        "図1: 浸水継続時間ランク別面積（左）と水系別12h以上面積 Top10（右）")}
<ul>
<li>全体 <b>{total_kei_ha:.0f} ha</b> のうち 12h 以上は <b>{kei12_ha:.0f} ha（{kei12_ha/total_kei_ha*100:.1f}%）</b> → H1 <b>{h1}</b></li>
<li>12h 以上の No.1 水系: <b>{h2_top}</b> → H2 (江の川が最大か) <b>{h2}</b></li>
<li>水系間で時間分布に明確な差 — 山間水系は短時間ピーク・大河川下流は長時間滞水</li>
</ul>

<h3>土砂種別 × 浸水継続時間 重複</h3>
{figure("X13_sediment_vs_duration.png",
        "図2: 土砂災害警戒区域（種別別）と浸水継続時間 重複ポリゴン数")}
<ul>
<li>土砂警戒区域 {len(sediment):,} 件のうち <b>{n_sed_with_kei:,} 件（{sed_kei_rate:.1f}%）</b>が 12h 以上浸水区域と交差 → H3 <b>{h3}</b></li>
<li>土石流 vs がけ崩れ で重複率に差 — 土石流は河川沿いに位置するため浸水との共存が多い</li>
<li>短時間ランク (12〜24h) の重複が最多 → 「早期の氾濫 × 土砂崩壊同時発生」を示唆</li>
</ul>

{figure("X13_heatmap_suikei_rank.png",
        "図8: 水系 × 継続時間ランク ポリゴン数ヒートマップ（上位10水系）")}
<p>水系行 × 時間列 のクロス表。江の川・太田川は全ランクに分布 (大河川の多様性)。
小水系ほど特定ランクに集中 (地形の均一性)。</p>
"""),

    ("RQ2: 雪崩 × 土砂 中山間集中分析", f"""
<h3>雪崩危険箇所の市町分布</h3>
{figure("X13_avalanche_by_city.png",
        "図3: 市町別 雪崩危険箇所数 Top15（中山間3市町強調）")}
<ul>
<li>雪崩 No.1 市町: <b>{h6_top}</b> ({ava_city.iloc[0]} 件) → H6 (庄原市が No.1) <b>{h6}</b></li>
<li>中山間 3 市町（庄原・三次・安芸太田）の集中率: <b>{ava_chusankan_pct:.1f}%</b> → H4 <b>{h4}</b></li>
<li>庄原市 1 市だけで全雪崩の <b>{ava_city.iloc[0]/len(avalanche)*100:.1f}%</b> を占有 — 中国山地の積雪地形集中</li>
</ul>

<h3>雪崩 種別 × 土砂 重複率</h3>
{figure("X13_avalanche_sediment_overlap.png",
        "図4: 雪崩危険種別（Ⅰ/Ⅱ/Ⅲ）× 土砂災害警戒区域 重複率")}
<ul>
<li>雪崩全体 × 土砂 重複率: <b>{ava_sed_rate:.1f}%</b> → H5 <b>{h5}</b></li>
<li>種別 Ⅰ（斜面急峻）の重複率が最高 → 急斜面は雪崩・土砂 両方の起因地形</li>
<li>種別 Ⅲ（軽微）の重複率が低い → 積雪特有のリスクと土砂崩壊は地形要因が異なる</li>
</ul>

{figure("X13_pie_charts.png",
        "図7: 浸水継続時間 面積比（左）と 雪崩危険種別比率（右）")}
<p>雪崩の <b>{(avalanche['bunrui']=='Ⅱ').sum()/len(avalanche)*100:.1f}%</b> が種別 Ⅱ（一般）。
Ⅰ（急峻）は {(avalanche['bunrui']=='Ⅰ').sum()/len(avalanche)*100:.1f}% で少数だが高い危険度。</p>
"""),

    ("RQ3: 三重リスクグリッド分析", f"""
<h3>三重リスクマップ (2km × 2km グリッド)</h3>
{figure("X13_triple_risk_map.png",
        "図5: 広島県 三重リスクマップ（浸水継続12h+ × 土砂 × 雪崩） — 2km格子")}
<ul>
<li>三重リスクセル: <b>{n_triple_cells} 個</b> / 全 {total_cells} セル = <b>{n_triple_cells/total_cells*100:.2f}%</b> → H7 (5%未満) <b>{h7}</b></li>
<li>二重リスクセル: <b>{n_double_cells} 個</b> → H9 (二重 ≥ 三重×3: {n_double_cells} ≥ {n_triple_cells*3}) <b>{h9}</b></li>
<li>青枠の中山間 3 市町に三重セルが集中する視覚的傾向が確認できる</li>
</ul>

<h3>三重リスク 市町別集計</h3>
{figure("X13_triple_city_ranking.png",
        "図6: 市町別 三重リスクグリッドセル数 Top10")}
<ul>
<li>三重リスク No.1: <b>{triple_city.index[0] if len(triple_city)>0 else '—'}</b>
  ({triple_city.iloc[0] if len(triple_city)>0 else 0} セル) → H8 (中山間市町が上位) <b>{h8}</b></li>
<li>三重リスクセルは県内でも <b>数カ所</b> に限定 — 孤立リスクの「的」を絞ることができる</li>
<li>防災対応優先度: 三重 → 二重 → 単一 の順に限られた行政資源を集中投下可能</li>
</ul>

<h3>考察 — 孤立継続時間の推計フレームワーク</h3>
<ul>
<li><b>三重リスクの意味</b>:
  雪崩（道路積雪閉塞）× 土砂崩壊（道路斜面崩落）× 12h 以上の浸水（低地孤立）が重なるエリアは、
  「複数の交通遮断要因が同時発生する可能性がある集落」として防災計画の最優先対象</li>
<li><b>孤立継続時間の下限推計</b>: 浸水継続時間のランク中央値を孤立の最低継続時間とみなすと、
  三重エリアの rank 20〜30（12h〜3日）が主体 → <b>集落孤立は最短でも半日以上</b></li>
<li><b>中山間集中の物理的背景</b>: 庄原・三次は江の川流域の盆地谷底地形。
  周囲の急斜面が雪崩・土砂を供給し、盆地底が浸水を捕捉する「地形的ワナ」</li>
<li><b>L09 との差分</b>: L09 は「どこが長く浸かるか」を問うた。
  本記事は「長く浸かっているときに土砂・雪崩も同時発生するか」を問い、
  孤立リスクを複合的に定量化した</li>
</ul>
"""),

    ("仮説検証と考察", f"""
<h3>仮説 vs 結果</h3>
<table>
<tr><th>仮説</th><th>判定</th><th>根拠</th></tr>
<tr><td>H1: 12h以上が全体の30%以上</td><td>{h1}</td><td>{h1_detail}</td></tr>
<tr><td>H2: 江の川水系の12h以上面積が最大</td><td>{h2}</td><td>12h以上 No.1 水系: {h2_top}</td></tr>
<tr><td>H3: 土砂×12h浸水の重複率 ≥ 5%</td><td>{h3}</td><td>{n_sed_with_kei:,}件 / {len(sediment):,} = {sed_kei_rate:.1f}%</td></tr>
<tr><td>H4: 中山間3市町の雪崩集中率 ≥ 50%</td><td>{h4}</td><td>{ava_chusankan_pct:.1f}% ({ava_chusankan}件)</td></tr>
<tr><td>H5: 雪崩×土砂 重複率 ≥ 20%</td><td>{h5}</td><td>{ava_sed_rate:.1f}%</td></tr>
<tr><td>H6: 雪崩 No.1 市町は庄原市</td><td>{h6}</td><td>{h6_top}: {ava_city.iloc[0]}件</td></tr>
<tr><td>H7: 三重リスクセル &lt; 全体の5%</td><td>{h7}</td><td>{n_triple_cells}/{total_cells} = {n_triple_cells/total_cells*100:.2f}%</td></tr>
<tr><td>H8: 三重リスクの上位は中山間市町</td><td>{h8}</td><td>No.1: {triple_city.index[0] if len(triple_city)>0 else '—'}</td></tr>
<tr><td>H9: 二重 ≥ 三重×3（段階的希少性）</td><td>{h9}</td><td>二重{n_double_cells} vs 三重{n_triple_cells}×3={n_triple_cells*3}</td></tr>
</table>

<h3>総合考察</h3>
<ul>
<li>L09 の「致命的滞水域」と L11 の「三重ハザード」を統合した本指標は、
  <b>「浸水深 × 継続時間 × 複合地盤リスク」</b>という 3 次元の孤立危険度を可視化する</li>
<li>三重リスクゾーンは稀少 ({n_triple_cells} セル = {n_triple_cells*4:.0f} km²) だが、
  そこに居住する集落への集中的な防災投資 — 備蓄倉庫・ヘリポート・集落道路強化 — が正当化できる</li>
<li>雪崩と土砂の重複率が種別 Ⅰ で高い事実は、
  「同一急斜面が冬季は雪崩、夏季は土砂として顕現する」という季節的多機能リスクを示す</li>
</ul>
"""),

    ("発展課題 (結果から導かれる新たな問い)", f"""
<ol>
<li><b>集落ポイントとの交差</b>:
  <ul>
  <li><b>結果X</b>: 三重リスクセル {n_triple_cells} 個を同定</li>
  <li><b>新仮説Y</b>: 三重リスクセル内に居住する集落数は単一リスクセル内より少ない（人口が自然選択的に回避）</li>
  <li><b>課題Z</b>: e-Stat 国勢調査の小地域集計（500m メッシュ）と三重リスクグリッドを sjoin し、
    リスク重複度別に人口密度を比較する</li>
  </ul>
</li>
<li><b>緊急輸送道路（L72）との交差</b>:
  <ul>
  <li><b>結果X</b>: 三重リスクの中山間集中を確認</li>
  <li><b>新仮説Y</b>: 三重リスクゾーン内の道路延長のうち、緊急輸送道路（L72）でカバーされる割合 &lt; 30%</li>
  <li><b>課題Z</b>: L72 緊急輸送道路 Shapefile と三重リスクセルの sjoin で
    「緊急輸送道路なし三重リスク集落」を抽出、消防署・ヘリポートまでの距離を計算する</li>
  </ul>
</li>
<li><b>排水機場（L35）の排水能力との接続</b>:
  <ul>
  <li><b>結果X</b>: 12h 以上浸水の水系別分布を把握</li>
  <li><b>新仮説Y</b>: 排水機場（L35 の 15 機場）の設置エリアは 12h 以上浸水域と重なり、
    機場の稼働が浸水継続時間を短縮する設計になっている</li>
  <li><b>課題Z</b>: L35 の排水機場 GPS と 12h 以上浸水ポリゴン（KEIZOKU_SHP）を 500m sjoin し、
    「排水機場カバー率」と浸水継続ランクの関係を定量化する</li>
  </ul>
</li>
</ol>
"""),

    ("使用データセットと再現方法", f"""
<h3>使用 DoBoX データセット</h3>
<table>
<tr><th>ID</th><th>名称</th><th>形式</th><th>役割</th></tr>
<tr><td><a href="https://hiroshima-dobox.jp/datasets/37" target="_blank">#37</a> + #315–331</td>
    <td>河川浸水継続時間_想定最大規模（各水系 19 件）</td><td>SHP</td><td>全河川統合版で一括カバー</td></tr>
<tr><td><a href="https://hiroshima-dobox.jp/datasets/48" target="_blank">#48</a> + #544–573</td>
    <td>土砂災害警戒区域・特別警戒区域（県全域 31 件）</td><td>SHP</td><td>3 種（土石流・がけ崩れ・地すべり）</td></tr>
<tr><td><a href="https://hiroshima-dobox.jp/datasets/50" target="_blank">#50</a> + #574–603</td>
    <td>雪崩危険箇所（県全域 31 件）</td><td>SHP</td><td>種別 Ⅰ/Ⅱ/Ⅲ</td></tr>
</table>
<p>論理カバー合計: <b>81 dataset_id</b>（個別水系ファイルを全河川統合版で代替）</p>

<h3>再現方法</h3>
<ul>
<li>再現コマンド: <code>python -X utf8 lessons/X13_mountain_compound_risk.py</code></li>
<li>実行時間: 約 {elapsed:.0f} 秒</li>
<li>生成ファイル: PNG 8 枚 + CSV 7 件 + HTML 1 件</li>
<li>再現スクリプト: <code>X13_mountain_compound_risk.py</code></li>
</ul>
"""),
]

html_out = render_lesson(
    num=13,
    title="X13: 雪崩 × 土砂 × 浸水継続時間 — 中山間複合リスクの孤立時間推計",
    tags=["X系","横断研究","GIS","複合ハザード","中山間","雪崩","土砂災害","浸水継続時間","孤立リスク"],
    time="60分",
    level="リテラシ基礎+α",
    data_label=(
        '<a href="https://hiroshima-dobox.jp/datasets/37" target="_blank">浸水継続時間 #37+#315-331 (19件)</a> × '
        '<a href="https://hiroshima-dobox.jp/datasets/48" target="_blank">土砂災害 #48+#544-573 (31件)</a> × '
        '<a href="https://hiroshima-dobox.jp/datasets/50" target="_blank">雪崩 #50+#574-603 (31件)</a>'
    ),
    sections=sections,
    script_filename="lessons/X13_mountain_compound_risk.py",
)
html_out = html_out.replace("<body>", '<body data-draft="1" data-stier="X">', 1)
out_path = LESSONS/"X13_mountain_compound_risk.html"
out_path.write_text(html_out, encoding="utf-8")
print(f"\n[HTML] {out_path.name} ({out_path.stat().st_size:,} bytes)")
print(f"=== DONE in {time.time()-t0:.1f}s ===")