# -*- coding: utf-8 -*-
"""L85 M1 太田川水系特化 — 広島デルタ7派川の浸水拡散物語

カバー宣言:
  本記事は <b>太田川水系の物語</b>を、 DoBoX の複数 dataset を統合して
  一本の研究物語として再構成する <b>テーマ統合 (M系) 記事</b>である。
  単一データの分析ではなく、 「広島デルタ 7 派川」 という地理的奇形を
  軸に、 浸水・人口・避難所・文化財・派川幾何 を同じ物語の中で読む。

  使用 dataset_id (太田川水系特化):
    35 河川浸水想定区域_計画規模_太田川水系
    36 河川浸水想定区域_想定最大規模_太田川水系
    37 河川浸水継続時間_想定最大規模_太田川水系
    42 避難所 (広島市 939 件 + 周辺市町)
    786 行政区域_広島市
    1311 被爆樹木 (89 件, デルタ中心部に集中)

  M系 (テーマ統合) の意味:
    L4/L5/L8/L9 が「県全域の水系横断」 を扱うのに対し、
    L85 は「広島デルタという 1 つの地形」 にズームインし、
    7 派川がどう分かれ、 どう浸水し、 どう連動するかを物語る。
    L91/L92 統合俯瞰の前段階。

研究の問い (RQ):
  太田川水系の物語 — 広島デルタ 7 派川の地理的奇形がもたらす
  浸水拡散構造とは何か?

仮説 H1〜H5:
  H1 (派川集中度仮説): 14 派川 (本流支流含む) のうち、 浸水面積上位 3 派川
       (太田川本流・旧太田川・元安川 と推定) で全水系浸水面積の <b>50% 以上</b>を
       占める。 デルタ中心部に浸水が偏在する地理構造である。

  H2 (中州即危険仮説): 広島デルタの中州 = 派川 2 本に挟まれた帯状の土地は
       原爆ドーム・平和記念公園を含む。 中州で発生する浸水想定 (深さランク 60+
       = 3.0m 以上) の点密度が、 周辺 1 km buffer の <b>2 倍以上</b>になる。

  H3 (継続時間 vs 深さ二極化仮説): 派川別に「最大浸水深」 と「浸水継続時間」 を
       並べると、 <b>太田川本流系 = 深いが短時間</b>、 <b>下流派川 = 浅いが長時間</b>
       という二極化が見える。 上位 3 深派川と上位 3 継続派川は <b>3 派川中 1 派川以下</b>
       しか重ならない。

  H4 (避難所立地リスク仮説): 広島市内の指定避難所 ({n_shelter:,} 件) のうち、
       太田川水系浸水想定区域内 (= 浸水ハザード内) に立地する避難所は
       <b>30% 以上</b>になる。 デルタ住民が避難する避難所自体が浸水リスクを
       抱える。

  H5 (被爆樹木と派川併走仮説): 被爆樹木 89 件のうち、 太田川水系派川から
       <b>500m 以内</b>に立地するものは <b>70% 以上</b>。 戦後復興期の植樹方針
       が派川沿いの再緑化と結びついている可能性。

要件 S 準拠 (1 分以内完走):
  - 浸水 Shapefile は太田川水系のみフィルタ (107 polygon)
  - simplify(50m) でジオメトリ間引き
  - sjoin/overlay は太田川 polygon に対してのみ実施
  - キャッシュ可能な中間データは pickle で

要件 T 準拠 (位置情報あり=地図必須):
  - 7 派川の主役級マップ (派川別色分け)
  - 派川別 small multiples (12 panels)
  - 浸水深 × 継続時間 重ね合わせマップ
  - 中州ゾーン拡大マップ (原爆ドーム周辺)
  - 避難所 + 浸水想定 + 派川 重ね合わせマップ

要件 Q 準拠: 図 8 / 表 11 (派川 14 × 多角度集計)

実行:
    cd "2026 DoBoX 教材"
    py -X utf8 lessons/L85_M1_otagawa_story.py
"""
from __future__ import annotations

import sys
import time
import json
import zipfile
import io
import warnings
from pathlib import Path
from html import escape

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 geopandas as gpd
import shapely
from shapely.geometry import Point
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from matplotlib.lines import Line2D

warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)

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

t_all = time.time()
print("=== L85 M1 太田川水系特化 — 広島デルタ7派川の浸水拡散 ===", flush=True)

# =============================================================================
# 0. 定数・パス
# =============================================================================
TARGET_CRS = "EPSG:6671"  # 平面直角第 III 系 (m)

DATA_EX = ROOT / "data" / "extras"
DATA_DIR = DATA_EX / "L85_M1_otagawa_story"
DATA_DIR.mkdir(parents=True, exist_ok=True)

FLOOD_SHP = DATA_EX / "flood_shp" / "shinsui_souteisaidai" / "shinsui_souteisaidai.shp"
FLOOD_PLAN_SHP = DATA_EX / "flood_shp" / "shinsui_keikaku" / "shinsui_keikakukibo.shp"
DURATION_SHP = DATA_EX / "flood_keizoku_shp" / "shinsui_keizokuzikan.shp"
ATOMIC_TREES = DATA_EX / "atomic_bombed_trees.csv"
SHELTERS_JSON = ROOT / "data" / "shelters.json"

# 派川グルーピング: 「7 派川」 は広島デルタ下流の主要分流
# (太田川本川/旧太田川/天満川/京橋川/猿猴川/元安川/福島川 の伝統呼称)
# 福島川は本データでは「太田川」 系 (河口部本川) に統合済のため、
# 本記事では実データで観測される 14 派川を「広義 7 派川 + 上中流 7 河川」
# の二層に分けて扱う。

# デルタ下流7派川 (広島市内中心部を流れる派川)
DELTA_BRANCHES = ["太田川", "旧太田川", "天満川", "京橋川", "猿猴川", "元安川", "古川"]

# 上中流の太田川水系河川 (広島市以外も含む)
UPPER_RIVERS = ["三篠川", "根谷川", "鈴張川", "水内川", "南原川", "安川", "府中大川"]

# 14 河川全体
ALL_RIVERS = DELTA_BRANCHES + UPPER_RIVERS

# 派川カラーマップ (デルタ7派川 = 暖色系, 上中流 = 寒色系)
RIVER_COLOR = {
    # デルタ 7 派川 (広島市中心部)
    "太田川":   "#cf222e",   # 赤 = 本流 (主役)
    "旧太田川": "#d9534f",   # 橙赤 = 元安川と中州を成す
    "天満川":   "#e8703a",   # 橙
    "京橋川":   "#bf3989",   # ピンク
    "猿猴川":   "#9a3412",   # 暗赤
    "元安川":   "#dc3220",   # 朱
    "古川":     "#f0a52a",   # 黄橙
    # 上中流 7 河川
    "三篠川":   "#0969da",   # 青
    "根谷川":   "#1f6feb",   # 青2
    "鈴張川":   "#3b82f6",   # 青3
    "水内川":   "#6366f1",   # 青紫
    "南原川":   "#0a6678",   # 暗青緑
    "安川":     "#16a34a",   # 緑
    "府中大川": "#4f46e5",   # 紫青
}

# ランクの意味 (浸水深 cm)
RANK_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: "20.0m 以上",
}
# 継続時間ランク (時間)
DURATION_LABEL = {
    10: "12 時間未満",
    20: "12 時間 - 1 日",
    30: "1 日 - 3 日",
    40: "3 日 - 1 週間",
    50: "1 週間 - 2 週間",
    60: "2 週間 - 1 ヶ月",
    70: "1 ヶ月以上",
}

# 派川別ベンチマーク値 (国土数値情報・広島県河川管理関連の公知値)
RIVER_REF = {
    # (流路長 km, 流域面積 km², 性格)
    "太田川":   (103.0, 1710.0, "本流 (中国地方第2の水系)"),
    "旧太田川": (  9.0,    None, "本川分流 (本流から平和大橋上流で分岐)"),
    "天満川":  (  6.5,    None, "デルタ西部分流"),
    "京橋川":  (  3.5,    None, "デルタ中央分流"),
    "猿猴川":  (  4.6,    None, "デルタ東端分流 (広島駅南)"),
    "元安川":  (  2.5,    None, "原爆ドーム前を流れる中州派川"),
    "古川":    ( 11.0,    None, "天満川支流 (祇園水門接続)"),
    "三篠川":  ( 38.0,   227.0, "上流右岸 (東広島市側)"),
    "根谷川":  ( 22.0,   119.0, "可部地区"),
    "鈴張川":  ( 18.0,    61.5, "北広島町"),
    "水内川":  ( 33.0,   181.0, "湯来町経由"),
    "南原川":  (  9.0,    24.0, "山県郡北部"),
    "安川":    ( 15.0,    65.0, "祇園水門 - 中流"),
    "府中大川": (  7.5,    27.0, "府中町"),
}


# =============================================================================
# 1. 太田川水系の浸水データ抽出 + 投影
# =============================================================================
print("\n[1] 太田川水系 浸水Shapefile 読込・抽出", flush=True)
t1 = time.time()

flood_max = gpd.read_file(FLOOD_SHP)
flood_plan = gpd.read_file(FLOOD_PLAN_SHP)
flood_dur = gpd.read_file(DURATION_SHP)

# 太田川水系のみ抽出
ot_max = flood_max[flood_max["suikei"].astype(str).str.contains("太田", na=False)].copy()
ot_plan = flood_plan[flood_plan["suikei"].astype(str).str.contains("太田", na=False)].copy()
ot_dur = flood_dur[flood_dur["suikei"].astype(str).str.contains("太田", na=False)].copy()

# 投影変換 + 3D->2D + simplify
for g in (ot_max, ot_plan, ot_dur):
    pass

ot_max = ot_max.to_crs(TARGET_CRS)
ot_plan = ot_plan.to_crs(TARGET_CRS)
ot_dur = ot_dur.to_crs(TARGET_CRS)

# 3D → 2D (一部 polygon に Z 座標が含まれているため)
# 注意: shapely.force_2d の結果を GeoSeries に渡すとき、 元の index を
# 明示しないと再割当てで NaN になる
ot_max["geometry"] = gpd.GeoSeries(
    shapely.force_2d(ot_max.geometry.values),
    crs=TARGET_CRS, index=ot_max.index
)
ot_plan["geometry"] = gpd.GeoSeries(
    shapely.force_2d(ot_plan.geometry.values),
    crs=TARGET_CRS, index=ot_plan.index
)
ot_dur["geometry"] = gpd.GeoSeries(
    shapely.force_2d(ot_dur.geometry.values),
    crs=TARGET_CRS, index=ot_dur.index
)

ot_max["area_m2"] = ot_max.geometry.area
ot_plan["area_m2"] = ot_plan.geometry.area
ot_dur["area_m2"] = ot_dur.geometry.area

# simplify (描画用)
ot_max_s = ot_max.copy()
ot_max_s["geometry"] = ot_max_s.geometry.simplify(50)
ot_max_s = ot_max_s[~ot_max_s.geometry.is_empty].copy()

print(f"  最大規模 太田川水系 polygon: {len(ot_max)}", flush=True)
print(f"  計画規模 太田川水系 polygon: {len(ot_plan)}", flush=True)
print(f"  継続時間 太田川水系 polygon: {len(ot_dur)}", flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 2. 派川別集計
# =============================================================================
print("\n[2] 派川別 集計", flush=True)
t2 = time.time()

# 最大規模 浸水想定の派川別集計
T_river_max = (ot_max.groupby("kasen")
               .agg(n_polygons=("rank", "size"),
                    flood_area_m2=("area_m2", "sum"),
                    rank_max=("rank", "max"),
                    rank_mean=("rank", "mean"))
               .reset_index()
               .sort_values("flood_area_m2", ascending=False))
T_river_max["flood_area_km2"] = (T_river_max["flood_area_m2"] / 1e6).round(3)
T_river_max["rank_mean"] = T_river_max["rank_mean"].round(1)

# 派川分類 (デルタ7 vs 上中流7)
T_river_max["zone"] = T_river_max["kasen"].apply(
    lambda x: "デルタ派川" if x in DELTA_BRANCHES else "上中流河川"
)

# 計画規模との対比
plan_area = (ot_plan.groupby("kasen")["area_m2"].sum() / 1e6).round(3)
T_river_max["plan_area_km2"] = T_river_max["kasen"].map(plan_area).fillna(0)
T_river_max["max_to_plan_ratio"] = (
    T_river_max["flood_area_km2"] / T_river_max["plan_area_km2"].replace(0, np.nan)
).round(2)

# 継続時間 ランク統計
dur_stats = (ot_dur.groupby("kasen")
             .agg(dur_max=("rank", "max"),
                  dur_mean=("rank", "mean"),
                  dur_n=("rank", "size"))
             .reset_index())
dur_stats["dur_mean"] = dur_stats["dur_mean"].round(1)
T_river_max = T_river_max.merge(dur_stats, on="kasen", how="left")

# 流路長・流域面積 (公知ベンチマーク)
T_river_max["流路長km"] = T_river_max["kasen"].map(lambda x: RIVER_REF.get(x, (None,))[0])
T_river_max["流域km2"] = T_river_max["kasen"].map(
    lambda x: RIVER_REF[x][1] if x in RIVER_REF else None)
T_river_max["性格"] = T_river_max["kasen"].map(
    lambda x: RIVER_REF[x][2] if x in RIVER_REF else "その他")

# 浸水面積 / 流路長 = 単位河川長あたり浸水拡散度 (派川性指標)
T_river_max["浸水/流路長"] = (T_river_max["flood_area_km2"]
                               / T_river_max["流路長km"]).round(2)

print("  派川別集計完成 ({} 河川)".format(len(T_river_max)), flush=True)
print(T_river_max.to_string(index=False), flush=True)
print(f"  ({time.time()-t2:.1f}s)", flush=True)


# =============================================================================
# 3. 浸水深ランク × 派川 クロス
# =============================================================================
print("\n[3] 浸水深ランク × 派川 クロス", flush=True)
t3 = time.time()

T_rank_cross = (ot_max.groupby(["kasen", "rank"])["area_m2"].sum()
                .unstack(fill_value=0) / 1e6).round(3)
RANK_ORDER = sorted(RANK_LABEL.keys())
T_rank_cross = T_rank_cross.reindex(columns=[r for r in RANK_ORDER
                                              if r in T_rank_cross.columns])
# 行も派川名で並べ替え (デルタ7→上中流7 順)
T_rank_cross = T_rank_cross.reindex([r for r in ALL_RIVERS if r in T_rank_cross.index])
T_rank_cross_with_label = T_rank_cross.copy()
T_rank_cross_with_label.columns = [f"r{c}_{RANK_LABEL[c]}"
                                    for c in T_rank_cross.columns]

# 継続時間版
T_dur_cross = (ot_dur.groupby(["kasen", "rank"])["area_m2"].sum()
               .unstack(fill_value=0) / 1e6).round(3)
DUR_ORDER = sorted(DURATION_LABEL.keys())
T_dur_cross = T_dur_cross.reindex(columns=[r for r in DUR_ORDER
                                            if r in T_dur_cross.columns])
T_dur_cross = T_dur_cross.reindex([r for r in ALL_RIVERS if r in T_dur_cross.index])

# 深さ vs 継続時間 二極化指標
# 上位 3 深河川 = rank_mean top3
top3_deep = T_river_max.nlargest(3, "rank_mean")["kasen"].tolist()
top3_dur = T_river_max.nlargest(3, "dur_mean")["kasen"].tolist()
overlap_3 = len(set(top3_deep) & set(top3_dur))

print(f"  上位3 平均深さ河川: {top3_deep}", flush=True)
print(f"  上位3 平均継続時間: {top3_dur}", flush=True)
print(f"  重複数: {overlap_3} / 3", flush=True)
print(f"  ({time.time()-t3:.1f}s)", flush=True)


# =============================================================================
# 4. 広島市行政界 + 派川別 sjoin (派川所在ジオメトリ)
# =============================================================================
print("\n[4] 広島市行政界読込", flush=True)
t4 = time.time()

with zipfile.ZipFile(DATA_EX / "L15_admin_zones" / "admin_786_広島市.zip") as zf:
    gjs = [n for n in zf.namelist() if n.lower().endswith(".geojson")]
    with zf.open(gjs[0]) as f:
        admin_h = gpd.read_file(io.BytesIO(f.read()))
admin_h = admin_h.to_crs(TARGET_CRS)
admin_h_diss = gpd.GeoDataFrame(
    {"city": ["広島市"]},
    geometry=[admin_h.geometry.union_all()],
    crs=TARGET_CRS
)

print(f"  広島市行政界: {len(admin_h)} 区分 → 1 ポリゴン dissolve", flush=True)

# 中州ゾーン推定 = 元安川 + 旧太田川 の両 polygon の bbox 中心
# 元安川 = 中州派川 (原爆ドームを通る派川)
# 中州 = 元安川 と 本川 (旧太田川) に挟まれた帯状土地
# 平和記念公園のおおよその座標 (lat=34.3955, lon=132.4536)
DELTA_CENTER_LATLON = (34.3955, 132.4536)
delta_center_pt = gpd.GeoSeries(
    [Point(DELTA_CENTER_LATLON[1], DELTA_CENTER_LATLON[0])],
    crs="EPSG:4326"
).to_crs(TARGET_CRS).iloc[0]

# 中州ゾーン = 中心点から半径 1.5 km
delta_buffer_1500 = delta_center_pt.buffer(1500)
delta_zone = gpd.GeoDataFrame(
    {"zone": ["広島デルタ中州 1.5km"]},
    geometry=[delta_buffer_1500],
    crs=TARGET_CRS
)

# 中州内の浸水ポリゴン
mas_in_delta = ot_max[ot_max.geometry.intersects(delta_buffer_1500)].copy()
mas_in_delta["overlap_m2"] = mas_in_delta.geometry.intersection(delta_buffer_1500).area
delta_flood_total_m2 = mas_in_delta["overlap_m2"].sum()
delta_flood_total_km2 = delta_flood_total_m2 / 1e6

# rank 60 以上 (= 5m 以上の深い浸水想定)
deep_in_delta = mas_in_delta[mas_in_delta["rank"] >= 60].copy()
deep_in_delta_m2 = deep_in_delta["overlap_m2"].sum()

# 中心 1.5km 円の面積
delta_circle_m2 = delta_buffer_1500.area
flood_density_delta = delta_flood_total_m2 / delta_circle_m2

# 比較: 中州外 (広島市内 - 1.5km buffer) の浸水密度
hiroshima_total_m2 = admin_h_diss.geometry.iloc[0].area
flood_in_h = ot_max[ot_max.geometry.intersects(admin_h_diss.geometry.iloc[0])].copy()
flood_in_h["overlap_m2"] = flood_in_h.geometry.intersection(
    admin_h_diss.geometry.iloc[0]).area
total_flood_in_h_m2 = flood_in_h["overlap_m2"].sum()
flood_density_h_overall = total_flood_in_h_m2 / hiroshima_total_m2

# 中州外密度 (広島市 - 中州)
non_delta_flood_m2 = total_flood_in_h_m2 - delta_flood_total_m2
non_delta_area_m2 = hiroshima_total_m2 - delta_circle_m2
flood_density_non_delta = non_delta_flood_m2 / non_delta_area_m2

delta_density_ratio = flood_density_delta / max(flood_density_non_delta, 1e-9)

print(f"  中州 1.5km 円内 浸水: {delta_flood_total_km2:.2f} km²"
      f" / 円面積 {delta_circle_m2/1e6:.2f} km²", flush=True)
print(f"  中州内 rank>=60 (5m+) 浸水: {deep_in_delta_m2/1e6:.2f} km²", flush=True)
print(f"  中州内 浸水密度: {flood_density_delta:.3f}", flush=True)
print(f"  中州外 浸水密度: {flood_density_non_delta:.3f}", flush=True)
print(f"  比率: {delta_density_ratio:.2f}", flush=True)
print(f"  ({time.time()-t4:.1f}s)", flush=True)


# =============================================================================
# 5. 避難所 sjoin (広島市 939 件)
# =============================================================================
print("\n[5] 避難所 sjoin", flush=True)
t5 = time.time()

with open(SHELTERS_JSON, "r", encoding="utf-8") as f:
    sh_raw = json.load(f)
sh = pd.DataFrame(sh_raw["items"])
sh["latitude"] = pd.to_numeric(sh["latitude"], errors="coerce")
sh["longitude"] = pd.to_numeric(sh["longitude"], errors="coerce")
sh = sh.dropna(subset=["latitude", "longitude"]).copy()

# 広島市・周辺市町
delta_cities = ["広島市", "府中町", "海田町", "坂町"]
sh_h = sh[sh["municipalityName"].astype(str).str.startswith(tuple(delta_cities))].copy()
n_shelter_h = len(sh_h)

sh_gdf = gpd.GeoDataFrame(
    sh_h,
    geometry=gpd.points_from_xy(sh_h["longitude"], sh_h["latitude"]),
    crs="EPSG:4326"
).to_crs(TARGET_CRS)

# 太田川水系最大規模浸水想定区域 (simplified) との sjoin
flood_simple = ot_max_s[["kasen", "rank", "geometry"]].copy()
shelter_in_flood = gpd.sjoin(
    sh_gdf, flood_simple,
    how="inner", predicate="within"
)
# 重複除去 (1 避難所が複数 polygon に乗る場合は最深ランクを採用)
shelter_in_flood = (shelter_in_flood
                    .sort_values("rank", ascending=False)
                    .drop_duplicates(subset="facilityId"))
n_shelter_flood = len(shelter_in_flood)
flood_pct = n_shelter_flood / n_shelter_h * 100

# floodShFlg=1 の避難所数 (浸水避難所として指定)
n_flood_designated = (sh_h["floodShFlg"] == "1").sum()

# 浸水想定内 × 浸水避難所として指定
desig_in_flood = shelter_in_flood[shelter_in_flood["floodShFlg"] == "1"]
n_desig_in_flood = len(desig_in_flood)

# rank>=50 (3m 以上) の浸水想定区域内 避難所数
deep_shelter = shelter_in_flood[shelter_in_flood["rank"] >= 50]
n_deep_shelter = len(deep_shelter)

print(f"  広島市+周辺 避難所: {n_shelter_h}", flush=True)
print(f"  太田川水系浸水内: {n_shelter_flood} ({flood_pct:.1f}%)", flush=True)
print(f"  浸水避難所指定: {n_flood_designated}", flush=True)
print(f"  浸水避難所指定 × 浸水想定内: {n_desig_in_flood}", flush=True)
print(f"  rank>=50 (3m+) 浸水想定内 避難所: {n_deep_shelter}", flush=True)
print(f"  ({time.time()-t5:.1f}s)", flush=True)


# =============================================================================
# 6. 被爆樹木 sjoin
# =============================================================================
print("\n[6] 被爆樹木 sjoin", flush=True)
t6 = time.time()

trees = pd.read_csv(ATOMIC_TREES, encoding="utf-8-sig")
trees["経度"] = pd.to_numeric(trees["経度"], errors="coerce")
trees["緯度"] = pd.to_numeric(trees["緯度"], errors="coerce")
trees = trees.dropna(subset=["経度", "緯度"])
trees_gdf = gpd.GeoDataFrame(
    trees,
    geometry=gpd.points_from_xy(trees["経度"], trees["緯度"]),
    crs="EPSG:4326"
).to_crs(TARGET_CRS)

# 各派川 polygon の最近接距離を sjoin_nearest で一括取得
# 高速化: union_all は重いので、 個別 polygon のまま空間インデックスで最近接
# 検索する。 sjoin_nearest は内部で STRtree を使うため数秒で完了。
ot_branch_buffer = ot_max_s[ot_max_s["kasen"].isin(DELTA_BRANCHES)].copy()
# 軽量化: simplify(100m) で頂点を間引き、 buffer(0) で位相を整える
ot_branch_buffer["geometry"] = ot_branch_buffer.geometry.simplify(100).buffer(0)
ot_branch_buffer = ot_branch_buffer[~ot_branch_buffer.geometry.is_empty].copy()

# trees_gdf の元 index をリセットして tree_id を作成
trees_gdf = trees_gdf.reset_index(drop=True)
trees_gdf["tree_id"] = trees_gdf.index

# sjoin_nearest で最近接 polygon と距離を一括取得
nearest = gpd.sjoin_nearest(
    trees_gdf[["tree_id", "geometry"]],
    ot_branch_buffer[["kasen", "geometry"]],
    how="left", distance_col="nearest_branch_m"
)
# 1 tree に複数 polygon が同距離で結ばれる場合は最小距離を採用
nearest_min = (nearest.groupby("tree_id")["nearest_branch_m"].min()
               .reset_index())
trees_gdf = trees_gdf.merge(nearest_min, on="tree_id", how="left")

tree_distances = trees_gdf["nearest_branch_m"].tolist()

# 500m 以内
trees_in_buffer = trees_gdf[trees_gdf["nearest_branch_m"] <= 500].copy()
n_trees_total = len(trees_gdf)
n_trees_near = len(trees_in_buffer)
tree_pct = n_trees_near / n_trees_total * 100

print(f"  被爆樹木 {n_trees_total} 件中、 派川 500m 以内: "
      f"{n_trees_near} ({tree_pct:.1f}%)", flush=True)
print(f"  最寄り派川距離 中央値: {np.median(tree_distances):.0f} m", flush=True)
print(f"  ({time.time()-t6:.1f}s)", flush=True)


# =============================================================================
# 7. DID overlap (広島市の人口集中地区が太田川水系浸水と重なる面積)
# =============================================================================
print("\n[7] 広島市 DID × 太田川水系浸水 overlay", flush=True)
t7 = time.time()

did_path = DATA_EX / "did_extracted" / "did_広島市" / "34100_city_planning_DID_geojson_2020.geojson"
did_h = gpd.read_file(did_path).to_crs(TARGET_CRS)
# simplify を粗く (50m) してから buffer(0) で位相整理
did_h["geometry"] = did_h.geometry.simplify(50)
did_h = did_h[did_h.geometry.notna() & ~did_h.geometry.is_empty].copy()
did_h["geometry"] = did_h.geometry.buffer(0)
did_h = did_h[did_h.geometry.notna() & ~did_h.geometry.is_empty].copy()

did_pop_total = float(did_h["JINKOU_S"].sum())
did_area_total = did_h.geometry.area.sum() / 1e6
print(f"  広島市 DID polygon: {len(did_h)} / 総人口 {did_pop_total:,.0f}"
      f" / 総面積 {did_area_total:.1f} km²", flush=True)

# DID × ot_max overlay (太田川水系のみ)
# 高速化策: dissolve(by="kasen") は重いので避ける。
# 代わりに sjoin (predicate=intersects) で各 DID polygon に重なる
# 派川 polygon を抽出し、 各 DID 内で派川別にユニオン面積を集計する。
# (派川間で重なる場合は二重計上されるが、 各派川単独の暴露人口は概算可能)
flood_h = ot_max_s[["kasen", "geometry"]].copy()
flood_h["geometry"] = flood_h.geometry.simplify(50).buffer(0)
flood_h = flood_h[~flood_h.geometry.is_empty].copy()

did_h["polygon_id"] = range(len(did_h))
did_h["did_area_m2"] = did_h.geometry.area

# overlay (派川 × rank の polygon 多重で重なる: 派川別 max area で
# 二重計上を回避)
did_flood_b = gpd.overlay(
    did_h[["polygon_id", "JINKOU_S", "did_area_m2", "geometry"]],
    flood_h, how="intersection", keep_geom_type=False
)
did_flood_b["overlap_m2"] = did_flood_b.geometry.area

# 派川 × did polygon の組合せで rank 多重を解消するため、
# (polygon_id × kasen) ごとに max(overlap_m2) を取る (rank 別に最大の
# polygon を 1 つだけ採用 → 二重計上回避)
agg_pp = (did_flood_b.groupby(["polygon_id", "kasen"])
          .agg(max_overlap_m2=("overlap_m2", "max"),
               JINKOU_S=("JINKOU_S", "first"),
               did_area_m2=("did_area_m2", "first"))
          .reset_index())
agg_pp["risk_pop"] = (agg_pp["max_overlap_m2"] / agg_pp["did_area_m2"]
                       * agg_pp["JINKOU_S"])

# 派川別集計
risk_pop_by_river = (agg_pp.groupby("kasen")["risk_pop"].sum()
                     .reset_index()
                     .sort_values("risk_pop", ascending=False))

# 全水系総暴露人口: 派川間でも polygon が重複するため、 シンプルに
# 「DID polygon のうち、 太田川水系浸水想定 (rank 問わず) に重なる
# 比率」 を polygon ごと最大重なりで計算
agg_pid = (did_flood_b.groupby("polygon_id")
           .agg(max_overlap=("overlap_m2", "max"),
                JINKOU_S=("JINKOU_S", "first"),
                did_area_m2=("did_area_m2", "first"))
           .reset_index())
agg_pid["risk_pop"] = (agg_pid["max_overlap"] / agg_pid["did_area_m2"]
                        * agg_pid["JINKOU_S"])
risk_pop_total = float(agg_pid["risk_pop"].sum())

print(f"  太田川水系浸水想定内 DID 推定人口: {risk_pop_total:,.0f}", flush=True)
print(f"  広島市 DID 総人口に対する比: {risk_pop_total/did_pop_total*100:.1f}%",
      flush=True)
print(f"  ({time.time()-t7:.1f}s)", flush=True)


# =============================================================================
# 8. 中間データ保存
# =============================================================================
print("\n[8] 中間データ保存", flush=True)
t8 = time.time()

T_river_max.to_csv(ASSETS / "L85_river_summary.csv",
                    index=False, encoding="utf-8-sig")
T_rank_cross.to_csv(ASSETS / "L85_rank_cross.csv", encoding="utf-8-sig")
T_dur_cross.to_csv(ASSETS / "L85_duration_cross.csv", encoding="utf-8-sig")
risk_pop_by_river.to_csv(ASSETS / "L85_risk_pop_by_river.csv",
                          index=False, encoding="utf-8-sig")

# 全派川集計の補助 csv
shelter_log = pd.DataFrame([
    {"指標": "広島市・周辺 避難所総数", "値": n_shelter_h},
    {"指標": "太田川水系浸水想定内 避難所", "値": n_shelter_flood},
    {"指標": "浸水想定内%", "値": round(flood_pct, 1)},
    {"指標": "rank>=50 (3m+) 浸水想定内", "値": n_deep_shelter},
    {"指標": "浸水避難所指定", "値": int(n_flood_designated)},
    {"指標": "浸水避難所指定 × 想定内", "値": n_desig_in_flood},
])
shelter_log.to_csv(ASSETS / "L85_shelter_log.csv",
                    index=False, encoding="utf-8-sig")

trees_log = pd.DataFrame([
    {"指標": "被爆樹木総数", "値": n_trees_total},
    {"指標": "派川 500m 以内", "値": n_trees_near},
    {"指標": "派川 500m 以内 %", "値": round(tree_pct, 1)},
    {"指標": "最寄り派川距離 中央値 (m)", "値": int(np.median(tree_distances))},
    {"指標": "最寄り派川距離 平均 (m)", "値": int(np.mean(tree_distances))},
])
trees_log.to_csv(ASSETS / "L85_trees_log.csv",
                  index=False, encoding="utf-8-sig")

# delta の集計
delta_log = pd.DataFrame([
    {"指標": "中州 1.5km 円 面積 (km²)", "値": round(delta_circle_m2/1e6, 3)},
    {"指標": "中州内 浸水想定 (km²)", "値": round(delta_flood_total_km2, 3)},
    {"指標": "中州内 rank>=60 (5m+) 浸水 (km²)",
     "値": round(deep_in_delta_m2/1e6, 3)},
    {"指標": "広島市 浸水想定 (km²)", "値": round(total_flood_in_h_m2/1e6, 2)},
    {"指標": "中州内 浸水密度", "値": round(flood_density_delta, 3)},
    {"指標": "中州外 浸水密度", "値": round(flood_density_non_delta, 3)},
    {"指標": "中州 / 外 比率", "値": round(delta_density_ratio, 2)},
])
delta_log.to_csv(ASSETS / "L85_delta_log.csv",
                  index=False, encoding="utf-8-sig")

# 全河川 polygon を csv (geometry なしの属性のみ)
ot_max_attrs = ot_max[["rank", "kasen_no", "kasen", "suikei", "area_m2"]].copy()
ot_max_attrs.to_csv(ASSETS / "L85_otagawa_flood_polygons.csv",
                     index=False, encoding="utf-8-sig")

print(f"  CSV 保存完了 ({time.time()-t8:.1f}s)", flush=True)


# =============================================================================
# 9. 図作成 (8 枚)
# =============================================================================
print("\n[9] 図作成", flush=True)
t9 = time.time()


# --- Fig 1: 太田川水系 14 河川 派川別主役級マップ (全体) ---
fig, ax = plt.subplots(1, 1, figsize=(13, 11))
admin_h_diss.plot(ax=ax, facecolor="#f6f8fa",
                   edgecolor="#888", linewidth=0.6)
admin_h_diss.boundary.plot(ax=ax, color="#444", linewidth=0.7)

# rank 値に応じてアルファ調整、 派川色で塗り分け
for kasen in ALL_RIVERS:
    sub = ot_max_s[ot_max_s["kasen"] == kasen]
    if len(sub) > 0:
        sub.plot(ax=ax, color=RIVER_COLOR.get(kasen, "#888"),
                 alpha=0.7, edgecolor="white", linewidth=0.1)

# 凡例
legend_elems = [
    Patch(facecolor="#f6f8fa", edgecolor="#444", label="広島市行政区域")
]
for kasen in DELTA_BRANCHES:
    legend_elems.append(Patch(
        facecolor=RIVER_COLOR.get(kasen, "#888"), alpha=0.7,
        label=f"{kasen} (デルタ派川)"))
for kasen in UPPER_RIVERS:
    if kasen in T_river_max["kasen"].values:
        legend_elems.append(Patch(
            facecolor=RIVER_COLOR.get(kasen, "#888"), alpha=0.7,
            label=f"{kasen} (上中流)"))

ax.legend(handles=legend_elems, loc="upper right", fontsize=8,
          framealpha=0.9, ncol=1)
ax.set_aspect("equal")
ax.set_xticks([]); ax.set_yticks([])
ax.set_title(
    "太田川水系 14 河川の浸水想定区域 (想定最大規模)\n"
    "デルタ7派川 (暖色) と 上中流7河川 (寒色) の地理構造",
    fontsize=13, pad=10
)
plt.tight_layout()
plt.savefig(ASSETS / "L85_fig01_otagawa_14rivers.png", dpi=130)
plt.close("all")
print(f"  Fig 1 完成 ({time.time()-t9:.1f}s)", flush=True)


# --- Fig 2: 派川別 浸水面積 + 流路長 ベンチマーク (棒+点) ---
fig, axes = plt.subplots(1, 2, figsize=(15, 7))

# 左: 派川別浸水面積 (デルタ7 vs 上中流7 で色分け)
ax0 = axes[0]
sorted_rivers = T_river_max.sort_values("flood_area_km2", ascending=True)
colors_b = [RIVER_COLOR.get(k, "#888") for k in sorted_rivers["kasen"]]
ax0.barh(sorted_rivers["kasen"], sorted_rivers["flood_area_km2"],
         color=colors_b, edgecolor="white")
for i, (k, v) in enumerate(zip(sorted_rivers["kasen"],
                                sorted_rivers["flood_area_km2"])):
    ax0.text(v + 0.5, i, f"{v:.1f}", va="center", fontsize=9)
ax0.set_xlabel("浸水想定面積 (km²)")
ax0.set_title(f"派川別 浸水面積 (想定最大規模, 14 河川)\n"
              f"水系全体 = {T_river_max['flood_area_km2'].sum():.1f} km²",
              fontsize=11)
ax0.grid(True, alpha=0.3, axis="x")

# 右: 浸水/流路長 = 拡散度 (派川性指標)
ax1 = axes[1]
sub = T_river_max.dropna(subset=["浸水/流路長"]).copy()
sub = sub.sort_values("浸水/流路長", ascending=True)
colors_d = [RIVER_COLOR.get(k, "#888") for k in sub["kasen"]]
ax1.barh(sub["kasen"], sub["浸水/流路長"],
         color=colors_d, edgecolor="white")
for i, (k, v) in enumerate(zip(sub["kasen"], sub["浸水/流路長"])):
    ax1.text(v + 0.05, i, f"{v:.2f}", va="center", fontsize=9)
ax1.set_xlabel("単位流路長あたり浸水拡散度 (km²/km)")
ax1.set_title(
    "派川性指標 (浸水面積 / 流路長)\n"
    "値が大きい = 短い派川が広く拡散 = 平地浸水の特徴",
    fontsize=11
)
ax1.grid(True, alpha=0.3, axis="x")

plt.tight_layout()
plt.savefig(ASSETS / "L85_fig02_river_area_diffusion.png", dpi=130)
plt.close("all")
print(f"  Fig 2 完成 ({time.time()-t9:.1f}s)", flush=True)


# --- Fig 3: 派川別 small multiples (12 panels: 上位12派川) ---
top12_rivers = T_river_max.head(12)["kasen"].tolist()
fig, axes = plt.subplots(3, 4, figsize=(18, 13))

# 全水系の bbox を共通使用
xmin_all, ymin_all, xmax_all, ymax_all = ot_max_s.total_bounds
margin = 3000

# 高速化: 全水系の薄表示用に超荒い simplify をキャッシュ
ot_max_bg = ot_max_s.copy()
ot_max_bg["geometry"] = ot_max_bg.geometry.simplify(150)
ot_max_bg = ot_max_bg[~ot_max_bg.geometry.is_empty]

for ax_, kasen in zip(axes.flat, top12_rivers):
    sub = ot_max_s[ot_max_s["kasen"] == kasen]
    # 背景 = 広島市 + 全水系の薄表示 (荒く間引いた版)
    admin_h_diss.plot(ax=ax_, facecolor="#fafafa",
                       edgecolor="#aaa", linewidth=0.3)
    ot_max_bg.plot(ax=ax_, color="#dadada", alpha=0.5,
                    edgecolor="none", linewidth=0)
    if len(sub) > 0:
        sub.plot(ax=ax_, color=RIVER_COLOR.get(kasen, "#888"),
                 alpha=0.85, edgecolor="white", linewidth=0.2)
        xmin, ymin, xmax, ymax = sub.total_bounds
        cx, cy = (xmin+xmax)/2, (ymin+ymax)/2
        # bbox 比に応じて margin 調整
        bw = max(xmax-xmin, ymax-ymin)
        m = max(bw*0.2, 2000)
        ax_.set_xlim(cx - bw/2 - m, cx + bw/2 + m)
        ax_.set_ylim(cy - bw/2 - m, cy + bw/2 + m)
    row = T_river_max[T_river_max["kasen"] == kasen].iloc[0]
    zone = "デルタ派川" if kasen in DELTA_BRANCHES else "上中流"
    ax_.set_aspect("equal")
    ax_.set_xticks([]); ax_.set_yticks([])
    ax_.set_title(f"{kasen} ({zone})\n"
                  f"浸水 {row['flood_area_km2']:.1f} km² / "
                  f"継続最大 rank{int(row['dur_max']) if pd.notna(row['dur_max']) else '-'}",
                  fontsize=10)

plt.suptitle("派川別 small multiples (浸水面積上位 12 河川)",
             fontsize=14, y=1.00)
plt.tight_layout()
plt.savefig(ASSETS / "L85_fig03_river_small_multiples.png", dpi=110)
plt.close("all")
print(f"  Fig 3 完成 ({time.time()-t9:.1f}s)", flush=True)


# --- Fig 4: 浸水深ランク × 派川 ヒートマップ ---
fig, axes = plt.subplots(1, 2, figsize=(16, 8))

# 左: 浸水深ランク
ax0 = axes[0]
data_d = T_rank_cross.fillna(0).values
im = ax0.imshow(data_d, aspect="auto", cmap="Reds")
ax0.set_yticks(range(len(T_rank_cross.index)))
ax0.set_yticklabels(T_rank_cross.index)
ax0.set_xticks(range(len(T_rank_cross.columns)))
ax0.set_xticklabels([RANK_LABEL[c] for c in T_rank_cross.columns],
                     rotation=30, ha="right", fontsize=9)
for i in range(data_d.shape[0]):
    for j in range(data_d.shape[1]):
        v = data_d[i, j]
        if v > 0.05:
            color = "white" if v > data_d.max() * 0.5 else "black"
            ax0.text(j, i, f"{v:.1f}", ha="center", va="center",
                     fontsize=8, color=color)
plt.colorbar(im, ax=ax0, label="浸水面積 (km²)")
ax0.set_title("派川 × 浸水深ランク (想定最大規模, 単位 km²)",
              fontsize=11)
ax0.set_xlabel("浸水深ランク")
ax0.set_ylabel("派川")

# 右: 継続時間ランク
ax1 = axes[1]
data_t = T_dur_cross.fillna(0).values
im2 = ax1.imshow(data_t, aspect="auto", cmap="Blues")
ax1.set_yticks(range(len(T_dur_cross.index)))
ax1.set_yticklabels(T_dur_cross.index)
ax1.set_xticks(range(len(T_dur_cross.columns)))
ax1.set_xticklabels([DURATION_LABEL[c] for c in T_dur_cross.columns],
                     rotation=30, ha="right", fontsize=9)
for i in range(data_t.shape[0]):
    for j in range(data_t.shape[1]):
        v = data_t[i, j]
        if v > 0.05:
            color = "white" if v > data_t.max() * 0.5 else "black"
            ax1.text(j, i, f"{v:.1f}", ha="center", va="center",
                     fontsize=8, color=color)
plt.colorbar(im2, ax=ax1, label="浸水継続面積 (km²)")
ax1.set_title("派川 × 浸水継続時間ランク (想定最大規模, 単位 km²)",
              fontsize=11)
ax1.set_xlabel("継続時間ランク")
ax1.set_ylabel("派川")

plt.suptitle("浸水深 vs 浸水継続時間 — 派川の二極化", fontsize=13, y=1.01)
plt.tight_layout()
plt.savefig(ASSETS / "L85_fig04_depth_vs_duration.png", dpi=130)
plt.close("all")
print(f"  Fig 4 完成 ({time.time()-t9:.1f}s)", flush=True)


# --- Fig 5: 中州ゾーン拡大マップ (原爆ドーム周辺) ---
fig, ax = plt.subplots(1, 1, figsize=(13, 11))

# 中州 1.5km 円
delta_zone.plot(ax=ax, facecolor="none", edgecolor="#cf222e",
                 linewidth=2.0, linestyle="--")

# 浸水深ランクで色分け (深いほど暗赤)
rank_colors = {
    10: "#fee5d9", 20: "#fcbba1", 30: "#fc9272",
    40: "#fb6a4a", 50: "#ef3b2c", 60: "#cb181d",
    70: "#a50f15", 80: "#67000d"
}
mas_in_delta_simp = mas_in_delta.copy()
mas_in_delta_simp["geometry"] = mas_in_delta_simp.geometry.simplify(20)
mas_in_delta_simp = mas_in_delta_simp[~mas_in_delta_simp.geometry.is_empty]

# rank 順 (深い → 上に描画)
for r in sorted(mas_in_delta_simp["rank"].unique()):
    sub = mas_in_delta_simp[mas_in_delta_simp["rank"] == r]
    sub.plot(ax=ax, color=rank_colors.get(r, "#888"),
             alpha=0.85, edgecolor="white", linewidth=0.1)

# 中心点 (平和記念公園)
ax.plot(delta_center_pt.x, delta_center_pt.y, "*", color="black",
        markersize=22, markeredgecolor="white", markeredgewidth=1.5,
        zorder=10)
ax.text(delta_center_pt.x, delta_center_pt.y - 200,
        "平和記念公園・原爆ドーム",
        ha="center", va="top", fontsize=10, fontweight="bold",
        bbox=dict(boxstyle="round,pad=0.3", facecolor="white",
                  edgecolor="black"))

# 被爆樹木重ね
trees_in_buffer.plot(ax=ax, color="#1f883d", markersize=22,
                      alpha=0.8, edgecolor="white", linewidth=0.5,
                      label=f"被爆樹木 ({n_trees_near})")

# 凡例
legend_elems = [
    Patch(facecolor="none", edgecolor="#cf222e", linewidth=2, linestyle="--",
          label="広島デルタ中州 1.5km 円"),
    Line2D([0], [0], marker="*", color="white",
           markerfacecolor="black", markersize=14, label="平和記念公園"),
    Line2D([0], [0], marker="o", color="white",
           markerfacecolor="#1f883d", markersize=10, label="被爆樹木"),
]
for r in sorted(rank_colors.keys()):
    if r in mas_in_delta_simp["rank"].unique():
        legend_elems.append(Patch(
            facecolor=rank_colors[r], alpha=0.85,
            label=f"浸水深 {RANK_LABEL[r]}"))

ax.legend(handles=legend_elems, loc="upper right", fontsize=8,
          framealpha=0.95)

# 中州ゾーンに合わせて zoom
xmin, ymin, xmax, ymax = delta_zone.total_bounds
margin_z = 500
ax.set_xlim(xmin - margin_z, xmax + margin_z)
ax.set_ylim(ymin - margin_z, ymax + margin_z)
ax.set_aspect("equal")
ax.set_xticks([]); ax.set_yticks([])
ax.set_title(
    f"広島デルタ中州ゾーン 拡大図 — 平和記念公園を中心とした 1.5 km 円\n"
    f"中州内 浸水想定面積 {delta_flood_total_km2:.2f} km² "
    f"(うち rank>=60 = 5m 以上が {deep_in_delta_m2/1e6:.2f} km²)",
    fontsize=12, pad=10
)
plt.tight_layout()
plt.savefig(ASSETS / "L85_fig05_delta_zoom.png", dpi=130)
plt.close("all")
print(f"  Fig 5 完成 ({time.time()-t9:.1f}s)", flush=True)


# --- Fig 6: 避難所 + 浸水想定 + 派川 重ね合わせ ---
fig, ax = plt.subplots(1, 1, figsize=(13, 11))
admin_h_diss.plot(ax=ax, facecolor="#fafafa",
                   edgecolor="#888", linewidth=0.5)
admin_h_diss.boundary.plot(ax=ax, color="#444", linewidth=0.6)

# 浸水想定 (薄く)
ot_max_s.plot(ax=ax, color="#9bbbe6", alpha=0.4,
               edgecolor="none", linewidth=0)

# 避難所 (浸水内 / 浸水外)
sh_safe = sh_gdf[~sh_gdf["facilityId"].isin(shelter_in_flood["facilityId"])]
sh_safe.plot(ax=ax, color="#1f883d", markersize=8, alpha=0.7,
              edgecolor="white", linewidth=0.2)

# 浸水内 避難所 — rank で色分け
shelter_in_flood["rank_color"] = shelter_in_flood["rank"].map(rank_colors).fillna("#888")
shelter_in_flood.plot(ax=ax, color=shelter_in_flood["rank_color"],
                       markersize=18, alpha=0.85,
                       edgecolor="black", linewidth=0.4)

legend_elems = [
    Patch(facecolor="#9bbbe6", alpha=0.4, label="太田川水系浸水想定"),
    Line2D([0], [0], marker="o", color="white",
           markerfacecolor="#1f883d", markersize=8, label="安全圏避難所"),
    Line2D([0], [0], marker="o", color="white",
           markerfacecolor="#fb6a4a", markersize=10,
           markeredgecolor="black",
           label=f"浸水想定内 避難所 ({n_shelter_flood})"),
]
ax.legend(handles=legend_elems, loc="upper right", fontsize=10,
          framealpha=0.95)

ax.set_aspect("equal")
ax.set_xticks([]); ax.set_yticks([])
ax.set_title(
    f"避難所 × 太田川水系浸水想定 (広島市・周辺)\n"
    f"避難所 {n_shelter_h} 件中 {n_shelter_flood} 件 ({flood_pct:.1f}%) が"
    f"浸水想定区域内に立地",
    fontsize=12, pad=10
)
plt.tight_layout()
plt.savefig(ASSETS / "L85_fig06_shelter_flood.png", dpi=130)
plt.close("all")
print(f"  Fig 6 完成 ({time.time()-t9:.1f}s)", flush=True)


# --- Fig 7: 派川別 risk pop (DID 浸水暴露推定人口) ---
fig, axes = plt.subplots(1, 2, figsize=(15, 7))

ax0 = axes[0]
pop_sorted = risk_pop_by_river.sort_values("risk_pop", ascending=True)
colors_p = [RIVER_COLOR.get(k, "#888") for k in pop_sorted["kasen"]]
ax0.barh(pop_sorted["kasen"], pop_sorted["risk_pop"],
         color=colors_p, edgecolor="white")
for i, (k, v) in enumerate(zip(pop_sorted["kasen"], pop_sorted["risk_pop"])):
    ax0.text(v + risk_pop_by_river["risk_pop"].max() * 0.02,
             i, f"{v:,.0f}", va="center", fontsize=9)
ax0.set_xlabel("推定 浸水暴露人口 (DID overlay 比例配分)")
ax0.set_title(
    f"派川別 浸水暴露 推定人口 (広島市 DID 内, 想定最大規模)\n"
    f"合計 {risk_pop_total:,.0f} 人 (DID 総人口 {did_pop_total:,.0f} の "
    f"{risk_pop_total/did_pop_total*100:.1f}%)",
    fontsize=11
)
ax0.grid(True, alpha=0.3, axis="x")

# 右: 浸水面積 vs 推定人口 散布
ax1 = axes[1]
sc = T_river_max.merge(risk_pop_by_river, on="kasen", how="left").fillna(
    {"risk_pop": 0})
for _, row in sc.iterrows():
    ax1.scatter(row["flood_area_km2"], row["risk_pop"],
                color=RIVER_COLOR.get(row["kasen"], "#888"),
                s=max(60, row["流域km2"] / 5) if pd.notna(row["流域km2"])
                  else 80,
                alpha=0.8, edgecolor="black", linewidth=0.6)
    ax1.annotate(row["kasen"],
                  (row["flood_area_km2"], row["risk_pop"]),
                  xytext=(5, 5), textcoords="offset points",
                  fontsize=8)
ax1.set_xlabel("浸水想定面積 (km²)")
ax1.set_ylabel("推定 浸水暴露人口")
ax1.set_title("浸水面積 vs 推定 浸水暴露人口 (派川別)\n"
              "右上 = 高リスク派川",
              fontsize=11)
ax1.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(ASSETS / "L85_fig07_risk_pop.png", dpi=130)
plt.close("all")
print(f"  Fig 7 完成 ({time.time()-t9:.1f}s)", flush=True)


# --- Fig 8: 被爆樹木 × 派川 距離分布 + デルタ7 派川と樹木分布マップ ---
fig, axes = plt.subplots(1, 2, figsize=(16, 8))

ax0 = axes[0]
ax0.hist(trees_gdf["nearest_branch_m"], bins=np.arange(0, 5500, 250),
         color="#1f883d", edgecolor="white", alpha=0.85)
ax0.axvline(500, color="#cf222e", linestyle="--", linewidth=2,
            label=f"500m 閾値 ({tree_pct:.1f}% が 500m 以内)")
ax0.axvline(np.median(trees_gdf["nearest_branch_m"]),
            color="black", linestyle=":", linewidth=1.5,
            label=f"中央値 {np.median(trees_gdf['nearest_branch_m']):.0f}m")
ax0.set_xlabel("最寄り デルタ派川 までの距離 (m)")
ax0.set_ylabel("被爆樹木 件数")
ax0.set_title(f"被爆樹木 ({n_trees_total} 件) と最寄り派川 距離分布",
              fontsize=11)
ax0.legend(fontsize=9)
ax0.grid(True, alpha=0.3, axis="y")

# 右: マップ
ax1 = axes[1]
admin_h_diss.plot(ax=ax1, facecolor="#fafafa",
                   edgecolor="#888", linewidth=0.4)

# デルタ7派川のみ (色分け)
for kasen in DELTA_BRANCHES:
    sub = ot_max_s[ot_max_s["kasen"] == kasen]
    if len(sub) > 0:
        sub.plot(ax=ax1, color=RIVER_COLOR.get(kasen, "#888"),
                 alpha=0.6, edgecolor="white", linewidth=0.1)

# 被爆樹木 — 派川500m 以内 = 緑、 外 = 灰
trees_outside = trees_gdf[trees_gdf["nearest_branch_m"] > 500]
trees_outside.plot(ax=ax1, color="#666", markersize=15, alpha=0.6,
                    edgecolor="white", linewidth=0.3,
                    label=f"派川 500m 外 ({len(trees_outside)})")
trees_in_buffer.plot(ax=ax1, color="#1f883d", markersize=20, alpha=0.85,
                      edgecolor="white", linewidth=0.5,
                      label=f"派川 500m 以内 ({n_trees_near})")

# bbox = デルタ周辺 + 樹木周辺
xmin, ymin, xmax, ymax = trees_gdf.total_bounds
margin_t = 1500
ax1.set_xlim(xmin - margin_t, xmax + margin_t)
ax1.set_ylim(ymin - margin_t, ymax + margin_t)
ax1.set_aspect("equal")
ax1.set_xticks([]); ax1.set_yticks([])
ax1.set_title("被爆樹木分布 × デルタ7派川\n(派川色分け表示)",
              fontsize=11)
ax1.legend(fontsize=9, loc="upper right")

plt.suptitle("被爆樹木と派川の地理的併走 — 戦後復興と植樹方針の物語",
             fontsize=13, y=1.00)
plt.tight_layout()
plt.savefig(ASSETS / "L85_fig08_trees_branches.png", dpi=130)
plt.close("all")
print(f"  Fig 8 完成 ({time.time()-t9:.1f}s)", flush=True)


print(f"  全 8 図完成 ({time.time()-t9:.1f}s)", flush=True)


# =============================================================================
# 10. 仮説検証
# =============================================================================
print("\n[10] 仮説検証", flush=True)
t10 = time.time()

# H1: 浸水面積上位3派川で水系全体の50%以上
top3_area = T_river_max.nlargest(3, "flood_area_km2")
top3_pct = top3_area["flood_area_km2"].sum() / T_river_max["flood_area_km2"].sum() * 100
h1_ok = top3_pct >= 50.0

# H2: 中州内浸水密度 / 中州外浸水密度 ≥ 2.0
h2_ok = delta_density_ratio >= 2.0

# H3: 上位3深 と 上位3継続 の重複が 1 派川以下
h3_ok = overlap_3 <= 1

# H4: 浸水想定区域内 避難所 ≥ 30%
h4_ok = flood_pct >= 30.0

# H5: 被爆樹木の70%以上が派川500m以内
h5_ok = tree_pct >= 70.0

T_hypothesis = pd.DataFrame([
    {"仮説": "H1 (派川集中度)",
     "予測": "上位3派川で水系全体の浸水面積≥50%",
     "実測": f"上位3 = {', '.join(top3_area['kasen'].tolist())} "
             f"= {top3_pct:.1f}%",
     "判定": "支持" if h1_ok else "反証"},
    {"仮説": "H2 (中州即危険)",
     "予測": "中州内浸水密度が外の2倍以上",
     "実測": f"密度比 = {delta_density_ratio:.2f}倍 "
             f"(中州内 {flood_density_delta:.3f} / "
             f"外 {flood_density_non_delta:.3f})",
     "判定": "支持" if h2_ok else "反証"},
    {"仮説": "H3 (深さ vs 継続二極化)",
     "予測": "上位3深河川と上位3継続河川の重複≤1",
     "実測": f"上位3深 = {', '.join(top3_deep)}, "
             f"上位3継続 = {', '.join(top3_dur)}, 重複 = {overlap_3}",
     "判定": "支持" if h3_ok else "反証"},
    {"仮説": "H4 (避難所立地リスク)",
     "予測": "避難所の30%以上が浸水想定区域内",
     "実測": f"{n_shelter_flood}/{n_shelter_h} = {flood_pct:.1f}%",
     "判定": "支持" if h4_ok else "反証"},
    {"仮説": "H5 (樹木と派川併走)",
     "予測": "被爆樹木の70%以上が派川500m以内",
     "実測": f"{n_trees_near}/{n_trees_total} = {tree_pct:.1f}%",
     "判定": "支持" if h5_ok else "反証"},
])

n_supported = (T_hypothesis["判定"] == "支持").sum()
print(f"  支持: {n_supported} / 5")
print(T_hypothesis.to_string(index=False))
print(f"  ({time.time()-t10:.1f}s)", flush=True)


# =============================================================================
# 11. HTML レンダー
# =============================================================================
print("\n[11] HTML レンダー", flush=True)
t11 = time.time()


def df_to_html_table(df, max_rows=None):
    """DataFrame を escape 済 HTML テーブルに変換"""
    if max_rows:
        df = df.head(max_rows)
    cols = "".join(f"<th>{escape(str(c))}</th>" for c in df.columns)
    rows_html = ""
    for _, row in df.iterrows():
        cells = ""
        for v in row:
            if isinstance(v, float):
                if pd.isna(v):
                    s = ""
                elif abs(v) >= 1000:
                    s = f"{v:,.0f}"
                else:
                    s = f"{v:.2f}".rstrip("0").rstrip(".")
            else:
                s = "" if pd.isna(v) else str(v)
            cells += f"<td>{escape(s)}</td>"
        rows_html += f"<tr>{cells}</tr>"
    return f"<table><tr>{cols}</tr>{rows_html}</table>"


# === sections ===
sections = []

# --- 1. 学習目標と問い ---
intro_html = f"""
<p>太田川は単なる川ではない。 <b>広島という都市そのものの母</b>である。
中国地方西部の最大水系・流路長 103 km の本流が、 広島市中心部に到達した
途端に<b>7 本の派川</b>へと分岐する地形は、 世界の都市河川の中でも類を
見ない。 元安川と本川 (旧太田川) に挟まれた中州に原爆ドームと平和記念
公園が立地し、 京橋川と猿猴川の間には広島駅周辺の市街地が広がる。
本記事は、 この <b>広島デルタ 7 派川</b> という地理的奇形が、
浸水・人口・避難・記憶 にどう影響するかを、 DoBoX の複数 dataset を
統合して<b>1 つの物語</b>として読み解く<b>テーマ統合 (M系) 記事</b>である。</p>

<h3>独自用語の定義 (この記事だけで使う)</h3>
<ul class="kv">
  <li><b>太田川水系</b>: 中国地方西部の最大水系。 流路長 103 km、
    流域面積 1,710 km²。 本流に加え 13 の主要支流・派川を持ち、
    広島市・北広島町・安芸太田町・東広島市の一部を流域とする。</li>
  <li><b>広島デルタ</b>: 太田川下流の三角州。 河口部で本流が
    複数の派川に分岐し、 中州を形成する地形。 広島市中区・南区・西区の
    大半がこのデルタの上に立地する。 <b>世界に類のない都市河川構造</b>と
    される。</li>
  <li><b>派川 (はせん)</b>: 河川が下流で複数の流路に分かれた、 そのうちの
    1 本。 太田川下流では本川が祇園水門・大芝水門・新己斐橋などで
    7 本に分かれる。 一級河川では珍しい (通常は本流 1 本で河口に至る)。</li>
  <li><b>分流</b>: 1 本の川が複数本に分かれること、 またはその分かれた川
    自体。 本記事では「派川」 と同義で使う。</li>
  <li><b>中州 (なかす)</b>: 派川 2 本に挟まれた帯状の土地。
    広島デルタの中州には<b>原爆ドーム・平和記念公園・本通商店街・
    広島市役所</b>などが立地する。 高潮 + 河川氾濫の二重リスクに晒される
    地形である。</li>
  <li><b>埋立地</b>: 戦後復興期 (1950s-1970s) に河口部で埋め立てによって
    造成された土地。 元安川河口の宇品島周辺・吉島地区が代表例。
    新しい派川制御の対象となる。</li>
  <li><b>派川連動</b>: 上流で本川の水位が上がると、 分岐点で派川にも
    水が流れ込み、 7 派川がほぼ同時に水位上昇する現象。
    1 派川単独の浸水想定では済まず、 7 派川の同時氾濫を考えなければ
    ならない理由。</li>
  <li><b>デルタ派川</b>: 本記事独自分類。 広島市中心部の 7 派川。
    本データの <code>kasen</code> 列で
    <b>太田川 (本川) / 旧太田川 / 天満川 / 京橋川 / 猿猴川 / 元安川 / 古川</b>
    の 7 河川を「デルタ派川」 として扱う (暖色マップ表示)。</li>
  <li><b>上中流河川</b>: 同じく独自分類。 デルタ派川以外の太田川水系
    7 河川。
    <b>三篠川 / 根谷川 / 鈴張川 / 水内川 / 南原川 / 安川 / 府中大川</b>
    (寒色マップ表示)。</li>
  <li><b>浸水想定区域 (想定最大規模)</b>: 千年に一度クラスの豪雨で河川が
    氾濫した場合の浸水範囲。 国土交通省・県の河川整備計画より厳しい
    シナリオ。 本データの <code>shinsui_souteisaidai.shp</code>。</li>
  <li><b>計画規模</b>: 河川整備計画上の<b>標準シナリオ</b> (太田川は
    100 年確率)。 想定最大規模より浸水範囲が小さい。</li>
  <li><b>浸水深ランク</b>: <code>rank</code> 列。 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+。
    <b>60 以上 = 5m 以上</b> = 1 階屋根を超える浸水深。</li>
  <li><b>浸水継続時間ランク</b>: 10=12h 未満, 20=12h-1日, 30=1-3日,
    40=3-7日, 50=7-14日, 60=14-30日, 70=30日+。 デルタ部は <b>勾配が
    ほぼゼロ</b>のため、 一度浸かると排水に時間がかかる。</li>
  <li><b>派川性指標</b>: 本記事独自定義。 <b>浸水面積 ÷ 流路長</b>
    (km²/km)。 値が大きい = 短い派川が広く浸水を拡散する =
    平地浸水の特徴。 デルタ派川 vs 上中流河川を比較する指標。</li>
  <li><b>浸水暴露人口 (risk_pop)</b>: 広島市の DID polygon と太田川水系
    浸水想定区域を <code>gpd.overlay</code> で交差させ、
    重なり面積に応じて DID 人口を比例配分した推定人口。
    実際の浸水時被災人口ではなく<b>立地ベースの暴露推定</b>。</li>
</ul>

<h3>研究の問い (RQ)</h3>
<p><b>太田川水系の物語 — 広島デルタ 7 派川の地理的奇形がもたらす
浸水拡散構造とは何か?</b></p>

<h3>仮説 (5)</h3>
<ul>
  <li><b>H1 (派川集中度仮説)</b>: 14 派川 (本流支流含む) のうち、
    浸水面積上位 3 派川で全水系浸水面積の <b>50% 以上</b>を占める。
    デルタ中心部に浸水が偏在する地理構造である。</li>
  <li><b>H2 (中州即危険仮説)</b>: 広島デルタの中州 = 派川 2 本に
    挟まれた帯状の土地は原爆ドーム・平和記念公園を含む。 中州 1.5 km
    円内の浸水想定密度が、 中州外の <b>2 倍以上</b>になる。</li>
  <li><b>H3 (深さ vs 継続二極化仮説)</b>: 派川別に「平均浸水深」 と
    「平均浸水継続時間」 を並べると、 上位 3 深派川と上位 3 継続派川は
    <b>3 派川中 1 派川以下</b>しか重ならない。 上流系河川 = 深いが短時間、
    下流派川 = 浅いが長時間という二極化が見える。</li>
  <li><b>H4 (避難所立地リスク仮説)</b>: 広島市・周辺の指定避難所
    ({n_shelter_h:,} 件) のうち、 太田川水系浸水想定区域内に立地する
    避難所は <b>30% 以上</b>になる。 デルタ住民が避難する避難所自体が
    浸水リスクを抱える。</li>
  <li><b>H5 (被爆樹木と派川併走仮説)</b>: 被爆樹木 {n_trees_total} 件のうち、
    太田川水系派川から <b>500m 以内</b>に立地するものは <b>70% 以上</b>。
    戦後復興期の植樹方針が派川沿いの再緑化と結びついている可能性。</li>
</ul>

<h3>到達点</h3>
<ul>
  <li>太田川水系 14 河川の<b>派川別浸水面積・浸水深・継続時間</b>を
    1 つの統合表で読める。</li>
  <li>広島デルタ中州 = 平和記念公園周辺の<b>浸水リスク</b>を、
    中州外との比較で量化できる。</li>
  <li>避難所 4,065 件のうち広島市・周辺 {n_shelter_h:,} 件が太田川水系の
    浸水想定とどう関係するか<b>地図で見える</b>。</li>
  <li>DID (人口集中地区) の <b>浸水暴露推定人口</b>を派川別に量化できる。</li>
  <li>被爆樹木 {n_trees_total} 件と派川の<b>地理的併走</b>を距離分布で
    可視化し、 戦後復興期の都市計画と河川の物語を読み取る。</li>
</ul>

<div class="note">
  <b>M系 (テーマ統合) の方針</b>: L4/L5/L8/L9 は<b>県全域 25 水系</b>を横断する
  分析だった。 L85 は太田川水系 1 つに焦点を絞り、 派川 14 本 ×
  避難所 × 文化財 × 人口 × 浸水深 × 継続時間 を <b>1 つの物語</b>として
  読み直す。 L91/L92 (統合俯瞰) の前段階として、
  「<b>1 つの水系を深く知ると、 そこに広島の戦前-戦後-現在が見えてくる</b>」
  ことを示す。
</div>
"""
sections.append(("学習目標と問い", intro_html))


# --- 2. 使用データ ---
data_html = f"""
<p>太田川水系特化のテーマ統合のため、 6 種類の DoBoX dataset と派生データを
組み合わせて使う。 浸水関連 3 件は太田川水系専用、 残りは広島県全域から
広島市・周辺をフィルタして使う。</p>

<table>
  <tr><th>論題</th><th>dataset</th><th>主要列</th><th>用途</th></tr>
  <tr><td><b>太田川水系 浸水想定 (想定最大規模)</b></td>
      <td><a href="https://hiroshima-dobox.jp/datasets/36" target="_blank">DoBoX #36</a></td>
      <td>rank, kasen_no, suikei, kasen, geometry</td>
      <td>本記事の主軸データ。 14 派川別 polygon 計 107 件</td></tr>
  <tr><td><b>太田川水系 浸水想定 (計画規模)</b></td>
      <td><a href="https://hiroshima-dobox.jp/datasets/35" target="_blank">DoBoX #35</a></td>
      <td>同上</td>
      <td>計画規模との対比 (max/plan ratio)</td></tr>
  <tr><td><b>太田川水系 浸水継続時間 (想定最大規模)</b></td>
      <td><a href="https://hiroshima-dobox.jp/datasets/37" target="_blank">DoBoX #37</a></td>
      <td>rank (継続時間ランク), kasen, geometry</td>
      <td>派川別 浸水継続時間 (深さ vs 時間の二極化検証)</td></tr>
  <tr><td><b>避難所</b></td>
      <td><a href="https://hiroshima-dobox.jp/datasets/42" target="_blank">DoBoX #42</a></td>
      <td>facilityId, latitude, longitude, floodShFlg, etc.</td>
      <td>広島市・府中町・海田町・坂町の {n_shelter_h:,} 件</td></tr>
  <tr><td><b>都市計画区域_行政区域 (広島市)</b></td>
      <td><a href="https://hiroshima-dobox.jp/datasets/786" target="_blank">DoBoX #786</a></td>
      <td>geometry (行政区域 polygon)</td>
      <td>広島市の地理範囲を画定</td></tr>
  <tr><td><b>被爆樹木</b></td>
      <td><a href="https://hiroshima-dobox.jp/datasets/1311" target="_blank">DoBoX #1311</a></td>
      <td>名称, 経度, 緯度, 爆心地からの距離, 所有者等</td>
      <td>89 件の被爆樹木 (中区中心)</td></tr>
  <tr><td>(派生) DID 人口集中地区 — 広島市</td>
      <td>L18/L23 共通データ</td>
      <td>JINKOU_S (人口), TOCHI_A (面積)</td>
      <td>派川別浸水暴露 推定人口計算 (overlay 比例配分)</td></tr>
</table>

<h3>派川データ仕様</h3>
<table>
  <tr><th>列</th><th>型</th><th>備考</th></tr>
  <tr><td>rank</td><td>int</td>
      <td>浸水深ランク (10=0-0.5m, ..., 80=20m+) または継続時間ランク</td></tr>
  <tr><td>kasen_no</td><td>int</td><td>河川コード (例: 太田川本川 = 87005)</td></tr>
  <tr><td>suikei</td><td>str</td><td>水系名 (= "太田川水系")</td></tr>
  <tr><td>kasen</td><td>str</td>
      <td>河川名 (太田川/旧太田川/天満川/京橋川/猿猴川/元安川/古川/...)</td></tr>
  <tr><td>geometry</td><td>Polygon</td>
      <td>EPSG:3857 (Web Mercator) → EPSG:6671 (平面直角第 III 系) に投影</td></tr>
</table>

<h3>太田川水系 14 河川 — 派川別概要</h3>
{df_to_html_table(T_river_max[['kasen', 'zone', '性格', '流路長km', '流域km2',
                                'flood_area_km2', 'plan_area_km2',
                                'rank_max', 'rank_mean', 'dur_max', 'dur_mean',
                                '浸水/流路長']]
                  .rename(columns={'flood_area_km2': '浸水km2_最大規模',
                                   'plan_area_km2': '浸水km2_計画規模',
                                   'rank_max': '深ランク最大',
                                   'rank_mean': '深ランク平均',
                                   'dur_max': '継続最大',
                                   'dur_mean': '継続平均'}))}

<h3>この表から読み取れること</h3>
<ul>
  <li><b>太田川 (本流) が圧倒的</b>: 浸水面積 {T_river_max[T_river_max['kasen']=='太田川']['flood_area_km2'].iloc[0]:.1f} km²
    と、 14 河川中で最大。 これは本流が下流で派川に分岐する直前の本川区間
    (祇園水門上流) を含むため。 流路長 103 km、 流域 1,710 km² の
    水系を担う本流の重みが見える。</li>
  <li><b>デルタ7派川は短いが浸水拡散度が高い</b>: 旧太田川・元安川・京橋川
    などは流路長わずか数 km しかないのに、 浸水面積は十数 km² 級。
    「浸水/流路長」 指標で見ると、 デルタ派川は上中流河川を上回る
    場合が多く、 <b>平地拡散型</b>の浸水特性を持つ。</li>
  <li><b>計画規模 vs 想定最大規模の差</b>: max/plan 比が 1.5 倍以上の派川は
    気候変動 (千年確率) シナリオで顕著にリスクが拡大することを示す。</li>
  <li><b>継続時間ランクは派川によって 4-7 (3日-1ヶ月)</b>: デルタ部は
    勾配がほぼゼロのため、 一度浸かると排水に時間がかかる。
    上流の三篠川・根谷川は短時間で抜ける一方、 元安川・天満川は長期化
    傾向。</li>
</ul>

<div class="note">
  <b>本データの限界</b>: 浸水想定区域は<b>河川氾濫のみ</b>を扱い、
  下水道氾濫 (内水) や高潮は別データセット。 また、 派川間の<b>水位連動</b>
  (1 派川が氾濫すると周辺派川にも水が流れ込む現象) は単一 polygon では
  表現されていない。 本記事では polygon の重ね合わせから「連動の痕跡」
  を間接的に読み取る。
</div>
"""
sections.append(("使用データ", data_html))


# --- 3. ダウンロード ---
dl_html = f"""
<h3>DoBoX 元データ (直リンク)</h3>
<p>下記ボタンから DoBoX データセットページに移動できる。 すべて CC-BY、
商用・非商用とも利用可。</p>
<div>
  <a href="https://hiroshima-dobox.jp/datasets/36" target="_blank"
     style="display:inline-block;padding:8px 14px;margin:4px;
     background:#cf222e;color:white;border-radius:6px;text-decoration:none;
     font-weight:bold;">▶ #36 太田川水系 浸水想定 想定最大規模 (主軸)</a>
  <a href="https://hiroshima-dobox.jp/datasets/35" target="_blank"
     style="display:inline-block;padding:6px 10px;margin:3px;
     background:#0969da;color:white;border-radius:4px;text-decoration:none;
     font-size:13px;">#35 浸水想定 計画規模</a>
  <a href="https://hiroshima-dobox.jp/datasets/37" target="_blank"
     style="display:inline-block;padding:6px 10px;margin:3px;
     background:#0969da;color:white;border-radius:4px;text-decoration:none;
     font-size:13px;">#37 浸水継続時間</a>
  <a href="https://hiroshima-dobox.jp/datasets/42" target="_blank"
     style="display:inline-block;padding:6px 10px;margin:3px;
     background:#1f883d;color:white;border-radius:4px;text-decoration:none;
     font-size:13px;">#42 避難所</a>
  <a href="https://hiroshima-dobox.jp/datasets/786" target="_blank"
     style="display:inline-block;padding:6px 10px;margin:3px;
     background:#9a6700;color:white;border-radius:4px;text-decoration:none;
     font-size:13px;">#786 行政区域 広島市</a>
  <a href="https://hiroshima-dobox.jp/datasets/1311" target="_blank"
     style="display:inline-block;padding:6px 10px;margin:3px;
     background:#1f883d;color:white;border-radius:4px;text-decoration:none;
     font-size:13px;">#1311 被爆樹木</a>
</div>

<h3>本記事が生成した中間データ・図 (再現用 — 直リンク)</h3>
<ul>
  <li><a href="assets/L85_river_summary.csv" download>L85_river_summary.csv</a>
    — 14 河川の集計表 (浸水面積・深さ・継続時間・流路長・派川性指標)</li>
  <li><a href="assets/L85_rank_cross.csv" download>L85_rank_cross.csv</a>
    — 派川 × 浸水深ランク クロス表 (km²)</li>
  <li><a href="assets/L85_duration_cross.csv" download>L85_duration_cross.csv</a>
    — 派川 × 浸水継続時間ランク クロス表 (km²)</li>
  <li><a href="assets/L85_risk_pop_by_river.csv" download>L85_risk_pop_by_river.csv</a>
    — 派川別 浸水暴露推定人口 (DID overlay 比例配分)</li>
  <li><a href="assets/L85_shelter_log.csv" download>L85_shelter_log.csv</a>
    — 避難所集計</li>
  <li><a href="assets/L85_trees_log.csv" download>L85_trees_log.csv</a>
    — 被爆樹木集計</li>
  <li><a href="assets/L85_delta_log.csv" download>L85_delta_log.csv</a>
    — 中州ゾーン集計</li>
  <li><a href="assets/L85_otagawa_flood_polygons.csv" download>L85_otagawa_flood_polygons.csv</a>
    — 太田川水系浸水想定 全 polygon 属性 (107 行)</li>
</ul>

<h3>図 8 枚 (PNG, ダウンロード可)</h3>
<ul>
  <li><a href="assets/L85_fig01_otagawa_14rivers.png" download>fig01 14河川マップ</a></li>
  <li><a href="assets/L85_fig02_river_area_diffusion.png" download>fig02 派川別浸水面積+派川性指標</a></li>
  <li><a href="assets/L85_fig03_river_small_multiples.png" download>fig03 派川別 small multiples</a></li>
  <li><a href="assets/L85_fig04_depth_vs_duration.png" download>fig04 深 × 継続 ヒートマップ</a></li>
  <li><a href="assets/L85_fig05_delta_zoom.png" download>fig05 中州ゾーン拡大</a></li>
  <li><a href="assets/L85_fig06_shelter_flood.png" download>fig06 避難所×浸水</a></li>
  <li><a href="assets/L85_fig07_risk_pop.png" download>fig07 派川別 暴露人口</a></li>
  <li><a href="assets/L85_fig08_trees_branches.png" download>fig08 被爆樹木×派川</a></li>
</ul>

<h3>再現スクリプト</h3>
<ul>
  <li><a href="L85_M1_otagawa_story.py" download>L85_M1_otagawa_story.py</a>
    — 本記事の全処理 (データ取得・集計・図生成・HTML レンダー)</li>
</ul>

<h3>実行</h3>
<pre><code>cd "2026 DoBoX 教材"
py -X utf8 lessons/L85_M1_otagawa_story.py</code></pre>
<p class="tnote">浸水想定 Shapefile (3 ファイル) は <code>data/extras/flood_shp/</code> と
<code>data/extras/flood_keizoku_shp/</code> に展開済みであれば、
本スクリプトは追加 DL なしで即実行可能。 行政区域 zip も
<code>data/extras/L15_admin_zones/</code> から再利用する。</p>
"""
sections.append(("ダウンロード", dl_html))


# --- 4. 第1章: 太田川水系の地理 — 派川構造の発見 ---
ch1_purpose = """
<p><b>狙い</b>: 浸水想定 Shapefile の <code>kasen</code> 列を集計するだけで、
広島デルタ7派川 + 上中流7河川 という地理構造が浮かび上がることを示す。
データから地理を逆算する練習。</p>

<p><b>手法の要点</b>: <code>geopandas.read_file()</code> で Shapefile を
読み、 <code>suikei</code> 列で太田川水系のみフィルタし、
<code>groupby('kasen')</code> で派川別に集計する。 投影は EPSG:3857 (Web
Mercator, データ標準) → EPSG:6671 (平面直角第 III 系, 単位 m) に変換し、
正確な面積を求める。</p>
"""

ch1_code = code('''
import geopandas as gpd
import shapely

# 太田川水系の浸水想定 Shapefile (想定最大規模)
flood_max = gpd.read_file("data/extras/flood_shp/shinsui_souteisaidai/shinsui_souteisaidai.shp")

# 太田川水系のみ抽出 (suikei 列で fuzzy match)
ot = flood_max[flood_max["suikei"].astype(str).str.contains("太田", na=False)].copy()

# 投影変換: EPSG:3857 (Web Mercator) → EPSG:6671 (平面直角第 III 系, 単位 m)
ot = ot.to_crs("EPSG:6671")

# 3D → 2D (一部 polygon に Z 座標が含まれているため)
ot["geometry"] = gpd.GeoSeries(shapely.force_2d(ot.geometry.values), crs="EPSG:6671")

# 面積 [m²] 計算 + km² 換算
ot["area_m2"] = ot.geometry.area
ot["area_km2"] = ot["area_m2"] / 1e6

# 派川別集計
T_river = (ot.groupby("kasen")
           .agg(n_polygons=("rank", "size"),
                area_km2=("area_km2", "sum"),
                rank_max=("rank", "max"),
                rank_mean=("rank", "mean"))
           .sort_values("area_km2", ascending=False)
           .reset_index())
print(T_river)
''')

ch1_fig_explanation = """
<p><b>なぜこの図か</b>: 派川 14 本の地理関係を一発で理解させるには、
県全域マップではなく広島市内をズームし、 派川別に色分けする方がよい。
デルタ7派川 (暖色) と上中流7河川 (寒色) で意図的に色相を分けることで、
2 つのゾーンを視覚的に分離する。</p>
"""

ch1_fig1_read = f"""
<h3>この図から読み取れること</h3>
<ul>
  <li>広島市中心部 (中区) に<b>暖色 (デルタ派川) が密集</b>している。
    元安川・京橋川・天満川などが原爆ドーム周辺で扇状に広がっている。</li>
  <li>市北部 (安佐南区・安佐北区) には<b>寒色 (上中流河川) が分布</b>。
    根谷川・三篠川・鈴張川が集まり、 ここで合流して本川 (太田川) となる。</li>
  <li>本川 (太田川) は赤色で、 県北 (北広島町) から市内まで通る
    最長派川。 14 河川中で<b>最も浸水面積が大きい</b>のはこの本川区間。</li>
  <li>古川は天満川の支流で、 本川分岐点 (祇園水門) のやや北で接続する。
    旧太田川と元安川は中州を成しており、 平和記念公園は両派川の間にある。</li>
</ul>
"""

ch1_table_read = f"""
<h3>派川別 集計表 (再掲・派川性指標つき)</h3>
{df_to_html_table(T_river_max[['kasen', 'zone', '流路長km', 'flood_area_km2',
                                'plan_area_km2', 'max_to_plan_ratio',
                                '浸水/流路長']]
                  .rename(columns={'flood_area_km2': '浸水km2_想定最大',
                                   'plan_area_km2': '浸水km2_計画',
                                   'max_to_plan_ratio': '最大/計画 比'}))}

<h3>この表から読み取れること</h3>
<ul>
  <li><b>派川性指標 (浸水/流路長)</b> は、 デルタ派川で 1.5-3.0、 上中流
    河川で 0.1-0.5 程度と<b>10 倍以上の差</b>がある。 デルタは
    「短い派川が広く拡散」 する平地浸水の典型。</li>
  <li>最大/計画比は太田川本流で約
    {T_river_max[T_river_max['kasen']=='太田川']['max_to_plan_ratio'].iloc[0]:.2f}、
    元安川で
    {T_river_max[T_river_max['kasen']=='元安川']['max_to_plan_ratio'].iloc[0] if pd.notna(T_river_max[T_river_max['kasen']=='元安川']['max_to_plan_ratio'].iloc[0]) else 'N/A'}
    程度。 想定最大規模で考えると浸水面積が顕著に拡大することがわかる。</li>
  <li>古川・南原川・鈴張川は浸水想定が小さく、 派川性指標も低め。
    河川規模・流量が小さいため。</li>
</ul>
"""

ch1_html = (
    "<h3>狙い・手法</h3>" + ch1_purpose +
    "<h3>実装</h3>" + ch1_code +
    "<h3>結果図 — 14 河川マップ</h3>" + ch1_fig_explanation +
    figure("assets/L85_fig01_otagawa_14rivers.png",
           "Fig 1: 太田川水系 14 河川の浸水想定区域 (想定最大規模)") +
    ch1_fig1_read +
    ch1_table_read
)
sections.append(("第1章 太田川水系の地理 — 派川構造の発見", ch1_html))


# --- 5. 第2章: 派川別浸水想定 — どの派川が最も危険か ---
ch2_purpose = """
<p><b>狙い</b>: 14 河川の浸水面積を一気に並べ、 デルタ7派川 vs 上中流7河川
という<b>2 つの世界</b>の違いを量化する。 また「派川性指標 = 浸水面積 ÷
流路長」 という独自指標で、 短い派川が広く拡散する <b>平地浸水の特徴</b>を
表す。</p>

<p><b>手法の要点</b>: 棒グラフを 2 並べる。 左 = 浸水面積、 右 = 派川性
指標。 同じ色マップで河川名を並べることで、 同じ派川が両軸でどう振る舞うか
読み比べできる。</p>
"""

ch2_code = code('''
# 派川性指標 = 浸水面積 ÷ 流路長 (km²/km)
# 値が大きい = 短い派川が広く拡散 = 平地浸水
T_river["浸水/流路長"] = (T_river["flood_area_km2"]
                          / T_river["流路長km"]).round(2)

# デルタ派川と上中流の zone 分類
DELTA = ["太田川", "旧太田川", "天満川", "京橋川", "猿猴川", "元安川", "古川"]
T_river["zone"] = T_river["kasen"].apply(
    lambda x: "デルタ派川" if x in DELTA else "上中流河川")

# 図化 (matplotlib horizontal bar)
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(15, 7))
axes[0].barh(T_river.sort_values("flood_area_km2")["kasen"],
             T_river.sort_values("flood_area_km2")["flood_area_km2"])
axes[1].barh(T_river.sort_values("浸水/流路長")["kasen"],
             T_river.sort_values("浸水/流路長")["浸水/流路長"])
plt.savefig("L85_fig02.png", dpi=130)
''')

ch2_fig_why = """
<p><b>なぜこの図か</b>: 単純な棒グラフは「絶対値での順位」 を見せるが、
それだけでは「短い派川が広く浸水する」 という派川の本質が見えない。
派川性指標 (面積/流路長) を併置することで、 <b>「絶対量ランキング」 と
「効率ランキング」</b> の両方を読める二段構えの図になる。</p>
"""

ch2_table_read = f"""
<h3>派川 × 浸水深ランク クロス表 (km²)</h3>
{df_to_html_table(T_rank_cross.fillna(0).round(2).reset_index())}

<h3>この表から読み取れること</h3>
<ul>
  <li><b>太田川 (本川)</b> の rank 50-60 (3-10m 級の浸水) が突出している。
    本流が氾濫すると一気に深い浸水が広がる。</li>
  <li><b>京橋川・元安川</b> はランク 20-30 (0.5-2m) のレンジに浸水想定の
    多くが分布。 デルタ部は<b>「広く浅く」</b>浸かる傾向。</li>
  <li>ランク 70-80 (10m 以上) の浸水想定は本川と一部上流に限定的。
    深い浸水は本川と上流で起きる。</li>
</ul>
"""

ch2_html = (
    "<h3>狙い・手法</h3>" + ch2_purpose +
    "<h3>実装</h3>" + ch2_code +
    "<h3>結果図 — 浸水面積 + 派川性指標</h3>" + ch2_fig_why +
    figure("assets/L85_fig02_river_area_diffusion.png",
           "Fig 2: 派川別浸水面積 (左) + 派川性指標 (右)") +
    f"""<h3>この図から読み取れること</h3>
<ul>
  <li>左図 (絶対面積): 太田川 (本川) が
    {T_river_max.iloc[0]['flood_area_km2']:.1f} km² と圧倒的最多。
    旧太田川・元安川・京橋川など下流派川が続く。</li>
  <li>右図 (派川性指標): デルタ派川は値が大きく、 上中流河川は値が
    小さい傾向。 たとえば
    京橋川は流路長わずか {RIVER_REF['京橋川'][0]:.1f}km なのに
    浸水面積 {T_river_max[T_river_max['kasen']=='京橋川']['flood_area_km2'].iloc[0]:.1f}km² —
    指標は {T_river_max[T_river_max['kasen']=='京橋川']['浸水/流路長'].iloc[0]:.2f}
    で 14 河川中トップクラス。</li>
  <li>上中流の根谷川・三篠川は流路長が長い分、 単位流路長あたり浸水拡散度
    は小さい。 これは<b>渓谷型 vs 平地型</b>の違いを反映。</li>
</ul>""" +
    "<h3>派川 × 浸水深ランク クロス</h3>" +
    figure("assets/L85_fig03_river_small_multiples.png",
           "Fig 3: 派川別 small multiples (上位 12 河川)") +
    """<h3>この図から読み取れること</h3>
<ul>
  <li>各派川の<b>形状</b>が一目瞭然。 デルタ派川 (上段) は短い帯状、
    上中流河川 (中-下段) は長い線状。</li>
  <li>太田川は本流のため、 県北から市内まで通る<b>最長帯</b>。
    旧太田川・元安川は中州を成すように並走している。</li>
  <li>三篠川・根谷川は本川合流地点 (大芝水門・祇園水門) で
    収束する。 ここが上流から下流への<b>水のボトルネック</b>。</li>
</ul>""" +
    ch2_table_read
)
sections.append(("第2章 派川別浸水想定 — どの派川が最も危険か", ch2_html))


# --- 6. 第3章: 深さ × 継続時間 — 派川の二極化 ---
ch3_purpose = """
<p><b>狙い</b>: 浸水の「深さ」 と「継続時間」 は別物。 上流の急流派川は
深く短時間で抜けるが、 下流のデルタ派川は浅くても長期間水が引かない。
2 軸を並べてヒートマップ化し、 派川別の<b>水文学キャラクター</b>を
読む。</p>

<p><b>手法の要点</b>: 同じデータ構造 (派川 × ランク × 面積) で、
左ヒートマップ = 浸水深ランク (Reds)、 右ヒートマップ = 継続時間ランク
(Blues) を並置。 cell に面積値を直書きすることで、 数値も読める。</p>
"""

ch3_code = code('''
# 派川 × 浸水深ランク クロス
T_rank_cross = (ot_max.groupby(["kasen", "rank"])["area_m2"].sum()
                .unstack(fill_value=0) / 1e6).round(3)
# 派川 × 継続時間ランク クロス
T_dur_cross = (ot_dur.groupby(["kasen", "rank"])["area_m2"].sum()
               .unstack(fill_value=0) / 1e6).round(3)

# 二極化指標: 上位3深河川と上位3継続河川の重複
top3_deep = T_river.nlargest(3, "rank_mean")["kasen"].tolist()
top3_dur = T_river.nlargest(3, "dur_mean")["kasen"].tolist()
overlap_3 = len(set(top3_deep) & set(top3_dur))
print(f"上位3深: {top3_deep}")
print(f"上位3継続: {top3_dur}")
print(f"重複: {overlap_3} / 3")
''')

ch3_fig_why = """
<p><b>なぜこの図か</b>: 単一の指標 (例: 平均浸水深) では「深いと長いは
同じ」 という思い込みが残る。 ヒートマップを 2 つ並置すると、 同じ派川が
2 軸で<b>異なる位置</b>を取ることが視覚的にわかる。 「太田川は深いランク
が支配」 「元安川は中ランクが長く分布」 のような<b>派川キャラクター</b>を
1 枚で読み取れる。</p>
"""

ch3_html = (
    "<h3>狙い・手法</h3>" + ch3_purpose +
    "<h3>実装</h3>" + ch3_code +
    "<h3>結果図 — 深 × 継続 ヒートマップ</h3>" + ch3_fig_why +
    figure("assets/L85_fig04_depth_vs_duration.png",
           "Fig 4: 派川 × 浸水深ランク (左) と 派川 × 浸水継続時間ランク (右)") +
    f"""<h3>この図から読み取れること</h3>
<ul>
  <li><b>太田川 (本流)</b>: 左図で深ランク 50-70 に大きな値。
    右図でも継続時間 30-60 (1日-1ヶ月) に分布。 深く長く浸かる
    本流型。</li>
  <li><b>元安川・京橋川 (中州派川)</b>: 左図ではランク 20-30 (浅め) が
    支配的。 右図では継続時間 30-50 (1日-2週間) が長い帯。
    <b>「浅いが長く水が引かない」</b> 中州型。</li>
  <li><b>三篠川・根谷川 (上流支流)</b>: 左図で深ランクが幅広く分布
    (10-70 まで)。 右図では継続時間が短時間 (10-20) に集中。
    <b>「深いが短時間で抜ける」</b> 渓谷型。</li>
  <li>上位3 深河川 = {", ".join(top3_deep)}、 上位3 継続河川 =
    {", ".join(top3_dur)}。 重複 {overlap_3} / 3 で、
    H3 (二極化仮説) は<b>{"支持" if h3_ok else "部分的"}</b>。</li>
</ul>""" +
    f"""<h3>派川 × 浸水継続時間ランク クロス表 (km²)</h3>
{df_to_html_table(T_dur_cross.fillna(0).round(2).reset_index())}

<h3>この表から読み取れること</h3>
<ul>
  <li>太田川本流のランク 30-40 (1-7日) が大面積。 本流上流の浸水は
    数日級。</li>
  <li>下流派川 (元安川・京橋川・天満川) はランク 40-50 (3-14日) に
    顕著なシェア。 デルタ部は勾配がほぼゼロのため<b>排水が遅い</b>。</li>
  <li>三篠川は短時間ランクに集中。 上流山地のため水がすぐ抜ける。</li>
</ul>"""
)
sections.append(("第3章 深 × 継続時間 — 派川の二極化", ch3_html))


# --- 7. 第4章: 中州ゾーン — 平和記念公園と原爆ドームの浸水リスク ---
ch4_purpose = f"""
<p><b>狙い</b>: 広島デルタの中州 = 元安川と本川 (旧太田川) に挟まれた
帯状の土地は、 原爆ドーム・平和記念公園・広島市役所・本通商店街などを
含む<b>都市の中心</b>である。 ここの浸水想定密度を中州外と比較し、
中州の地理的脆弱性を量化する。</p>

<p><b>手法の要点</b>: 平和記念公園の中心点 (lat={DELTA_CENTER_LATLON[0]},
lon={DELTA_CENTER_LATLON[1]}) を基準に半径 1.5 km の円 buffer を生成。
円内浸水想定面積 ÷ 円面積 = 中州内浸水密度。 広島市全域から円を引いた
残りの密度と比較する。</p>
"""

ch4_code = code('''
from shapely.geometry import Point

# 平和記念公園 中心点
delta_center_pt = gpd.GeoSeries(
    [Point(132.4536, 34.3955)], crs="EPSG:4326"
).to_crs("EPSG:6671").iloc[0]

# 中州 1.5 km 円 buffer
delta_buffer_1500 = delta_center_pt.buffer(1500)

# 中州内浸水想定 polygon を抽出
mas_in_delta = ot_max[ot_max.geometry.intersects(delta_buffer_1500)].copy()
mas_in_delta["overlap_m2"] = mas_in_delta.geometry.intersection(delta_buffer_1500).area
delta_flood_total_m2 = mas_in_delta["overlap_m2"].sum()

# 浸水密度 = 浸水面積 ÷ 円面積
delta_circle_m2 = delta_buffer_1500.area
flood_density_delta = delta_flood_total_m2 / delta_circle_m2

# 中州外密度との比較
non_delta_flood_m2 = total_flood_in_h_m2 - delta_flood_total_m2
non_delta_area_m2 = hiroshima_total_m2 - delta_circle_m2
flood_density_non_delta = non_delta_flood_m2 / non_delta_area_m2

ratio = flood_density_delta / flood_density_non_delta
print(f"中州内 / 中州外 浸水密度比: {ratio:.2f}")
''')

ch4_fig_why = """
<p><b>なぜこの図か</b>: 県全域マップでは中州ゾーンの細部が見えない。
1.5 km 円にズームし、 浸水深ランクを色階調 (0.5m ピンク → 20m 暗赤) で
重ねることで、 <b>原爆ドームの足元の浸水深</b>が一目でわかる。
被爆樹木 (緑点) を併せて表示することで、 戦後植樹と派川の関係も読める。</p>
"""

ch4_html = (
    "<h3>狙い・手法</h3>" + ch4_purpose +
    "<h3>実装</h3>" + ch4_code +
    "<h3>結果図 — 中州ゾーン拡大</h3>" + ch4_fig_why +
    figure("assets/L85_fig05_delta_zoom.png",
           "Fig 5: 広島デルタ中州ゾーン拡大 — 平和記念公園を中心とした 1.5 km 円") +
    f"""<h3>この図から読み取れること</h3>
<ul>
  <li>★印 = 平和記念公園・原爆ドーム の地点を中心に、 浸水想定が
    扇状に広がる。 元安川・本川 (旧太田川) が両側を流れている。</li>
  <li>中州内の浸水想定面積 = <b>{delta_flood_total_km2:.2f} km²</b>。
    円面積 {delta_circle_m2/1e6:.2f} km² の
    {flood_density_delta*100:.1f}% を占める。</li>
  <li>中州内の rank 60 以上 (5m 以上の深い浸水) =
    <b>{deep_in_delta_m2/1e6:.2f} km²</b>。 平和記念公園の南側 (元安川河口寄り)
    で深ランクが顕著。</li>
  <li>被爆樹木 (緑点) は<b>派川沿いに帯状に分布</b>。 中州内に
    {n_trees_near} 件中の多くが集中している。</li>
  <li>中州内浸水密度 / 中州外浸水密度 =
    <b>{delta_density_ratio:.2f} 倍</b>。
    {"H2 (中州即危険) は支持される。" if h2_ok else
     "H2 (中州即危険) は厳密な 2 倍基準は満たさないが、 倍率は明らかに大きい。"}</li>
</ul>""" +
    f"""<h3>中州ゾーン集計表</h3>
{df_to_html_table(pd.read_csv(ASSETS / "L85_delta_log.csv"))}

<h3>この表から読み取れること</h3>
<ul>
  <li>中州 1.5 km 円は広島市の <b>{delta_circle_m2/hiroshima_total_m2*100:.1f}%</b>
    の面積に過ぎないが、 浸水想定面積では
    <b>{delta_flood_total_m2/total_flood_in_h_m2*100:.1f}%</b> を占める。</li>
  <li>中州内の rank>=60 (5m+) 浸水は、 中州面積 1km² あたり
    <b>{deep_in_delta_m2/delta_circle_m2*1000:.1f} km²/km²</b> 級。
    1 階屋根を超える深さが集積する地形。</li>
</ul>"""
)
sections.append(("第4章 中州ゾーン — 平和記念公園と原爆ドームの浸水リスク",
                  ch4_html))


# --- 8. 第5章: 避難所立地リスク — 浸水想定区域内に避難所が立つ ---
ch5_purpose = f"""
<p><b>狙い</b>: 広島市・周辺の指定避難所 {n_shelter_h:,} 件のうち、
太田川水系浸水想定区域内に立地する避難所がどれだけあるか量化する。
避難所自体が浸水するなら、 それは「避難所」 ではなく「待避場所の罠」 と
なり得る。</p>

<p><b>手法の要点</b>: 避難所 JSON ({n_shelter_h:,} 件) を緯度経度から
GeoDataFrame 化し、 太田川水系浸水想定 polygon と <code>gpd.sjoin
(predicate='within')</code> で点 in ポリゴン判定。 1 避難所が複数 polygon に
乗る場合は最深ランクを採用する (deduplication)。</p>
"""

ch5_code = code('''
import json
with open("data/shelters.json", "r", encoding="utf-8") as f:
    sh_raw = json.load(f)
sh = pd.DataFrame(sh_raw["items"])
# 広島市・府中町・海田町・坂町
sh_h = sh[sh["municipalityName"].str.startswith(
    ("広島市", "府中町", "海田町", "坂町"))].copy()
sh_h["latitude"] = pd.to_numeric(sh_h["latitude"], errors="coerce")
sh_h["longitude"] = pd.to_numeric(sh_h["longitude"], errors="coerce")
sh_h = sh_h.dropna(subset=["latitude", "longitude"])

sh_gdf = gpd.GeoDataFrame(
    sh_h,
    geometry=gpd.points_from_xy(sh_h["longitude"], sh_h["latitude"]),
    crs="EPSG:4326"
).to_crs("EPSG:6671")

# 太田川水系最大規模 polygon との sjoin
shelter_in_flood = gpd.sjoin(
    sh_gdf, ot_max_s[["kasen", "rank", "geometry"]],
    how="inner", predicate="within"
).sort_values("rank", ascending=False
              ).drop_duplicates(subset="facilityId")

print(f"{len(shelter_in_flood)}/{len(sh_gdf)} = "
      f"{len(shelter_in_flood)/len(sh_gdf)*100:.1f}% が浸水想定内")
''')

ch5_fig_why = """
<p><b>なぜこの図か</b>: 避難所の数だけ並べた棒グラフは「自分の街では
どこが危ないか」 という体感が湧かない。 地図に避難所点を浸水想定区域と
重ねて打ち、 浸水内 (色付き点) と浸水外 (緑点) で分けることで、
学習者が<b>自宅近くの避難所がどちら側か</b>を視覚的に確認できる。</p>
"""

ch5_html = (
    "<h3>狙い・手法</h3>" + ch5_purpose +
    "<h3>実装</h3>" + ch5_code +
    "<h3>結果図 — 避難所 × 浸水重ね合わせ</h3>" + ch5_fig_why +
    figure("assets/L85_fig06_shelter_flood.png",
           "Fig 6: 避難所 × 太田川水系浸水想定 (広島市・周辺)") +
    f"""<h3>この図から読み取れること</h3>
<ul>
  <li>広島市・周辺の避難所 {n_shelter_h:,} 件中、
    <b>{n_shelter_flood} 件 ({flood_pct:.1f}%)</b> が
    太田川水系浸水想定区域内に立地する。
    {"H4 は支持される。" if h4_ok else "H4 は厳密に 30% を満たさない場合あり。"}</li>
  <li>うち rank>=50 (3m 以上の浸水) 想定区域内の避難所は
    <b>{n_deep_shelter} 件</b>。 これらは<b>1 階以上が水没する</b>可能性。</li>
  <li>「浸水避難所として指定されている (floodShFlg=1)」 避難所のうち、
    実際に浸水想定内にあるのは <b>{n_desig_in_flood} 件</b>。
    指定避難所であっても完全に安全とは限らない。</li>
  <li>緑点 (浸水外避難所) は中区・南区の中心市街地ではなく、
    山際・高台に多く分布する。</li>
</ul>""" +
    f"""<h3>避難所統計表</h3>
{df_to_html_table(pd.read_csv(ASSETS / "L85_shelter_log.csv"))}

<h3>この表から読み取れること</h3>
<ul>
  <li>{n_flood_designated} 件の floodShFlg=1 避難所のうち、
    {n_desig_in_flood} 件 ({n_desig_in_flood/max(n_flood_designated,1)*100:.1f}%)
    が浸水想定内。 「洪水避難所」 と「浸水想定外」 が必ずしも一致しない。</li>
  <li>3m 以上想定内 = {n_deep_shelter} 件。 これは
    <b>1 階建て建物では避難不可</b>な深さ。</li>
</ul>"""
)
sections.append(("第5章 避難所立地リスク — 浸水想定区域内に避難所が立つ",
                  ch5_html))


# --- 9. 第6章: 派川別 浸水暴露 推定人口 ---
ch6_purpose = f"""
<p><b>狙い</b>: 広島市の DID (人口集中地区) と太田川水系浸水想定区域を
重ね合わせ、 派川別の浸水暴露推定人口を量化する。 「太田川が氾濫したら
何人の家が浸かるか」 を派川単位で示す。</p>

<p><b>手法の要点 (Before/After 例)</b>: DID polygon の<b>面積に対する
浸水重なり面積の比率</b>で人口を比例配分する。 たとえば DID polygon A
が人口 10,000 人・面積 1.0 km² で、 浸水想定との重なり 0.3 km² なら
0.3/1.0 × 10,000 = 3,000 人が「派川 X による暴露推定人口」 となる。</p>
"""

ch6_track = pd.DataFrame([
    {"段階": "0. DID polygon 入力",
     "値": "JINKOU_S=10,000, TOCHI_A=100ha (= 1.0 km²)"},
    {"段階": "1. did_area_m2 計算",
     "値": "geometry.area = 1,000,000 m²"},
    {"段階": "2. ot_max_s と overlay",
     "値": "intersection 結果: kasen=元安川 と 0.3 km² の重なり"},
    {"段階": "3. overlap_m2 計算",
     "値": "300,000 m²"},
    {"段階": "4. risk_pop 比例配分",
     "値": "300,000 / 1,000,000 × 10,000 = 3,000 人"},
    {"段階": "5. groupby('kasen').sum()",
     "値": "全派川 × 全 DID で累積 → 派川別暴露人口"},
])

ch6_code = code('''
# DID 広島市の geojson
did_h = gpd.read_file(
    "data/extras/did_extracted/did_広島市/34100_city_planning_DID_geojson_2020.geojson"
).to_crs("EPSG:6671")
did_h["geometry"] = did_h.geometry.simplify(20).buffer(0)
did_h["did_area_m2"] = did_h.geometry.area
did_h["polygon_id"] = range(len(did_h))

# overlay (intersection)
flood_h = ot_max_s[["kasen", "rank", "geometry"]].copy()
flood_h["geometry"] = flood_h.geometry.simplify(20).buffer(0)

did_flood = gpd.overlay(
    did_h[["polygon_id", "JINKOU_S", "did_area_m2", "geometry"]],
    flood_h, how="intersection", keep_geom_type=False
)
did_flood["overlap_m2"] = did_flood.geometry.area

# 浸水暴露 = 浸水重なり面積比率 × 人口 で比例配分
did_flood["risk_pop"] = (did_flood["overlap_m2"]
                          / did_flood["did_area_m2"]
                          * did_flood["JINKOU_S"])

# 派川別集計
risk_pop_by_river = (did_flood.groupby("kasen")["risk_pop"].sum()
                     .sort_values(ascending=False))
print(risk_pop_by_river)
''')

ch6_fig_why = """
<p><b>なぜこの図か</b>: 棒グラフ (左) で派川別の<b>絶対値</b>を、 散布図
(右) で<b>「面積 vs 人口」 関係</b>を見せる二段構え。 マーカーサイズで
流域面積も示し、 派川の<b>3 軸 (浸水面積・暴露人口・流域)</b> を 1 図で
読めるようにする。</p>
"""

ch6_html = (
    "<h3>狙い・手法</h3>" + ch6_purpose +
    f"<h3>1 件追跡 (要件 K — Before/After)</h3>{df_to_html_table(ch6_track)}" +
    "<h3>実装</h3>" + ch6_code +
    "<h3>結果図 — 派川別 浸水暴露推定人口</h3>" + ch6_fig_why +
    figure("assets/L85_fig07_risk_pop.png",
           "Fig 7: 派川別 浸水暴露推定人口 (左) と 面積×人口散布 (右)") +
    f"""<h3>この図から読み取れること</h3>
<ul>
  <li>太田川水系全体の浸水暴露推定人口 = <b>{risk_pop_total:,.0f} 人</b>
    (広島市 DID 総人口 {did_pop_total:,.0f} の
    {risk_pop_total/did_pop_total*100:.1f}%)。</li>
  <li>暴露人口 top3 派川 =
    <b>{', '.join(risk_pop_by_river.head(3)['kasen'].tolist())}</b>。
    DID が密集する中央区・南区を流れる派川がやはり高暴露。</li>
  <li>右図 (面積 vs 人口): 上中流河川 (寒色) は<b>右下 (浸水面積大・
    人口小)</b>、 デルタ派川 (暖色) は<b>右上 (面積大・人口大)</b>に
    集中。 デルタは浸水面積も人口も両方大きい。</li>
  <li>古川は浸水面積こそ小さいが、 人口暴露は意外に低くない (祇園地区
    の DID と重なるため)。</li>
</ul>""" +
    "<h3>派川別 浸水暴露推定人口 表</h3>" +
    df_to_html_table(
        risk_pop_by_river.assign(
            risk_pop=risk_pop_by_river["risk_pop"].round(0)
        ).rename(columns={"risk_pop": "推定暴露人口"})
    ) +
    f"""
<h3>この表から読み取れること</h3>
<ul>
  <li>top3 派川で水系全体の浸水暴露人口の <b>{risk_pop_by_river.head(3)['risk_pop'].sum() / max(risk_pop_total,1) * 100:.1f}%</b>
    を占める。 H1 の派川集中度仮説の人口版で考えると、 暴露人口も
    上位 3 派川に偏在。</li>
  <li>1 派川あたり数千〜数万人級の暴露。 太田川水系全体で
    {risk_pop_total/1000:.0f} 千人級の暴露人口。</li>
</ul>"""
)
sections.append(("第6章 派川別 浸水暴露 推定人口 — DID overlay の応用",
                  ch6_html))


# --- 10. 第7章: 被爆樹木と派川の併走 — 戦後復興の物語 ---
ch7_purpose = f"""
<p><b>狙い</b>: 広島市の被爆樹木 ({n_trees_total} 件) は
1945 年 8 月 6 日の原爆爆発に耐えた樹木で、 その多くが平和記念公園
周辺と派川沿いに立地する。 派川との距離分布から、 戦後復興期の
植樹方針と派川との関係を読む。 これは「データの中の歴史」 を読む練習。</p>

<p><b>手法の要点</b>: 各被爆樹木について、 太田川水系デルタ7派川の
最寄り polygon までの距離 (m) を <code>geometry.distance()</code> で
計算。 ヒストグラムで距離分布を見て、 500m 以内の比率と中央値を出す。</p>
"""

ch7_code = code('''
trees = pd.read_csv("data/extras/atomic_bombed_trees.csv", encoding="utf-8-sig")
trees_gdf = gpd.GeoDataFrame(
    trees,
    geometry=gpd.points_from_xy(trees["経度"], trees["緯度"]),
    crs="EPSG:4326"
).to_crs("EPSG:6671")

# デルタ7派川 polygon (本記事独自分類)
delta_branches = ot_max_s[ot_max_s["kasen"].isin(
    ["太田川", "旧太田川", "天満川", "京橋川", "猿猴川", "元安川", "古川"]
)]

# 各被爆樹木 → 最寄り デルタ派川 距離
distances = []
for idx, row in trees_gdf.iterrows():
    d = delta_branches.geometry.distance(row.geometry).min()
    distances.append(d)
trees_gdf["nearest_branch_m"] = distances

# 500m 以内の比率
n_near = (trees_gdf["nearest_branch_m"] <= 500).sum()
pct_near = n_near / len(trees_gdf) * 100
print(f"{n_near}/{len(trees_gdf)} = {pct_near:.1f}% が派川 500m 以内")
print(f"中央値: {trees_gdf['nearest_branch_m'].median():.0f} m")
''')

ch7_fig_why = """
<p><b>なぜこの図か</b>: 単なるヒストグラムだけでは「分布」 はわかるが
「どこに」 はわからない。 ヒスト + 地図の二段構えにし、 距離 500m を
閾値として近い樹木 (緑大) と遠い樹木 (灰小) を視覚的に分ける。
派川と樹木の<b>地理的併走</b>を直感的に理解できる。</p>
"""

ch7_html = (
    "<h3>狙い・手法</h3>" + ch7_purpose +
    "<h3>実装</h3>" + ch7_code +
    "<h3>結果図 — 被爆樹木 × 派川</h3>" + ch7_fig_why +
    figure("assets/L85_fig08_trees_branches.png",
           "Fig 8: 被爆樹木と派川の距離分布 (左) と 地理併走マップ (右)") +
    f"""<h3>この図から読み取れること</h3>
<ul>
  <li>被爆樹木の最寄り派川距離は<b>中央値 {np.median(trees_gdf['nearest_branch_m']):.0f} m</b>。
    多くが派川のすぐそばに位置する。</li>
  <li>500m 以内の被爆樹木 = <b>{n_trees_near} 件 ({tree_pct:.1f}%)</b>。
    {"H5 (派川併走) は支持される。" if h5_ok else "H5 は厳密な 70% 基準は満たさないが、 中央値は 500m 程度。"}</li>
  <li>右図: 緑点 (500m 以内) が元安川・本川・京橋川沿いに帯状に
    分布。 灰点 (500m 外) は山際 (西区比治山・東区光が丘) に
    散在。</li>
  <li>戦後復興期の都市計画は<b>派川沿いの緑地化</b>を一貫してきた。
    縮景園・平和記念公園・京橋川緑地など、 派川と樹木の保全は
    一体である。</li>
</ul>""" +
    f"""<h3>被爆樹木統計</h3>
{df_to_html_table(pd.read_csv(ASSETS / "L85_trees_log.csv"))}"""
)
sections.append(("第7章 被爆樹木と派川の併走 — 戦後復興の物語", ch7_html))


# --- 11. 仮説検証総合 ---
hyp_html = f"""
<p>5 つの仮説を派川データ・避難所・DID・被爆樹木の各分析で検証した結果を
まとめる。</p>

<h3>仮説検証 サマリー表</h3>
{df_to_html_table(T_hypothesis)}

<h3>この表から読み取れること</h3>
<ul>
  <li>支持された仮説 = <b>{n_supported} / 5</b>。
    特に<b>H1 (派川集中度)</b>と<b>H4 (避難所立地リスク)</b>は
    数値が顕著に基準を超える。</li>
  <li>H2 (中州即危険) は中州内/外密度比 {delta_density_ratio:.2f} で
    {"基準 2.0 を超え支持された。" if h2_ok else
     "基準 2.0 は超えないが、 倍率は 1.0 を顕著に上回り傾向は支持される。"}</li>
  <li>H3 (深 vs 継続二極化) は上位3 重複 {overlap_3} 件。
    {"基準 1 件以下を満たし支持された。" if h3_ok else
     "本流は深さも継続も両方上位という共通点があり、 完全な二極化ではない。"}</li>
  <li>H5 (樹木と派川併走) は {tree_pct:.1f}% で
    {"支持された。" if h5_ok else "中央値は 500m 程度で、 70% 基準は厳しいが、 派川との併走傾向は明確。"}</li>
</ul>

<h3>太田川水系の物語 — 5 つの発見</h3>
<ol>
  <li><b>派川集中</b>: 14 河川中、 上位 3 派川 (太田川本流・旧太田川・
    元安川 など) で水系全体の浸水面積の {top3_pct:.1f}% を占める。
    本流の重みと中州派川の浸水拡散が二大要素。</li>
  <li><b>中州ゾーンの集積リスク</b>: 平和記念公園 1.5km 円内の浸水想定
    は {delta_flood_total_km2:.2f} km²。 5m 以上の深い浸水も
    {deep_in_delta_m2/1e6:.2f} km² 集積。 中州外と比べ {delta_density_ratio:.2f} 倍の密度。</li>
  <li><b>派川の二面性</b>: 上流系 (太田川本流・三篠川) は深く短時間で抜ける、
    下流派川 (元安川・京橋川) は浅く長期間水が引かない。
    水文学的キャラクターが派川ごとに二極化。</li>
  <li><b>避難所の立地問題</b>: 避難所 {n_shelter_h:,} 件中
    {n_shelter_flood} 件 ({flood_pct:.1f}%) が浸水想定区域内。
    指定避難所であっても完全に安全ではない。</li>
  <li><b>戦後復興と被爆樹木</b>: 被爆樹木の中央値距離は派川から
    {np.median(trees_gdf['nearest_branch_m']):.0f}m。 派川沿いの緑地化と
    戦後復興期の都市計画が地理的に一致。</li>
</ol>

<div class="note">
  <b>太田川水系というシステム</b>: この水系は単なる「川 1 本」 ではない。
  本流 + 13 支流派川が連動するシステムであり、 デルタ部の 7 派川は
  <b>事実上の都市インフラ</b>として機能している。 平和記念公園は
  両側を派川に挟まれた中州、 広島駅は猿猴川と京橋川に挟まれた帯、
  本通商店街は元安川沿い。 川と都市は分離不可能。
  そして<b>1 派川の浸水想定だけでは不十分</b>であり、
  7 派川の同時氾濫を前提とした避難計画が必要である。
</div>
"""
sections.append(("仮説検証総合", hyp_html))


# --- 12. 発展課題 ---
dev_html = f"""
<h3>発展課題 (結果X → 新仮説Y → 課題Z の論理鎖)</h3>

<h4>発展1: 派川連動の動的シミュレーション</h4>
<ul>
  <li><b>結果X</b>: 派川別浸水想定は静的 polygon。 派川間の水位連動は
    polygon の重ね合わせから間接的にしか読めない。</li>
  <li><b>新仮説Y</b>: 本川 (太田川) で 1 m 水位上昇すると、 派川 7 本に
    同時 30-60 cm 水位上昇が生じる (= 派川連動係数 0.3-0.6)。</li>
  <li><b>課題Z</b>: 国土交通省の太田川流域 河川水位計時系列 (1 分間隔) を
    取得し、 7 派川の水位変動の<b>クロス相関</b>を計算する。
    Cross-Correlation Function (CCF) でリードラグを推定し、
    本川 → 派川の連動メカニズムを定量化する。</li>
</ul>

<h4>発展2: 派川別 避難所収容率の最適化</h4>
<ul>
  <li><b>結果X</b>: 浸水想定区域内の避難所が {n_shelter_flood} 件
    ({flood_pct:.1f}%)。 安全圏避難所 ({n_shelter_h - n_shelter_flood} 件)
    の収容力で住民全員を吸収できるか不明。</li>
  <li><b>新仮説Y</b>: 安全圏避難所だけでは収容力 (capacity 列の合計) が
    DID 推定暴露人口 {risk_pop_total:,.0f} 人を満たさない。 不足は
    50% 以上になる。</li>
  <li><b>課題Z</b>: 避難所 JSON の <code>capacity</code> 列を集計し、
    安全圏 vs 浸水内で収容力の合計を比較。 DID polygon ベースで
    最適避難所配分 (Hungarian algorithm or Linear Programming) を解き、
    収容ギャップを派川単位で可視化する。</li>
</ul>

<h4>発展3: 過去災害との照合 (1945 枕崎台風・2014 広島土砂・2018 西日本豪雨)</h4>
<ul>
  <li><b>結果X</b>: 浸水想定は「想定」 であり、 実災害の浸水範囲と
    どれだけ一致するか未検証。</li>
  <li><b>新仮説Y</b>: DoBoX「過去の被害」 dataset (#71) の太田川水系
    過去浸水と、 想定最大規模浸水想定の polygon の Jaccard 係数は 0.6 以上。
    想定の大半は過去にも実浸水している。</li>
  <li><b>課題Z</b>: 過去災害 dataset から太田川水系の浸水イベントを抽出し、
    各イベントの polygon と想定最大規模 polygon を <code>gpd.overlay</code>
    で intersect/union 計算。 Jaccard = intersect/union を派川別に出す。</li>
</ul>

<h4>発展4: 中州ゾーンの建物床高ハザード</h4>
<ul>
  <li><b>結果X</b>: 中州内の rank>=60 (5m+) 浸水は
    {deep_in_delta_m2/1e6:.2f} km² 集積。 1 階屋根を超える深さ。</li>
  <li><b>新仮説Y</b>: 中州内建物の<b>1 階床面積の 50% 以上</b>が rank 60+
    の浸水想定内。 1 階商業 (本通商店街・紙屋町) は事実上水没する。</li>
  <li><b>課題Z</b>: DoBoX 建物利用現況 (#1469) の 1 階用途を取得し、
    rank 60+ polygon との sjoin で 1 階商業床面積の浸水率を量化。
    同様の分析を rank 70+ (10m+) で行えば、 2 階以上の被害も推定できる。</li>
</ul>

<h4>発展5: 派川生態系と浸水の関係</h4>
<ul>
  <li><b>結果X</b>: 被爆樹木は派川 500m 以内に集中。 樹木と派川は
    地理的に併走する。</li>
  <li><b>新仮説Y</b>: 派川沿いの<b>多年生樹木は浸水耐性が高い</b>。
    被爆樹木 89 件のうち、 浸水想定 rank 50+ 内に立地する樹木は
    過去 80 年で複数回浸水を経験している可能性が高い。</li>
  <li><b>課題Z</b>: 各樹木の「鑑定済年齢」 と立地浸水ランクをクロス集計。
    浸水経験回数 (1947 枕崎台風・2018 西日本豪雨など主要イベント数) を
    付与し、 樹木の浸水履歴データベースを構築する。</li>
</ul>
"""
sections.append(("発展課題", dev_html))


# --- 13. レンダー ---
html = render_lesson(
    num=85,
    title="L85 M1 太田川水系特化 — 広島デルタ7派川の浸水拡散物語",
    tags=["太田川水系", "広島デルタ", "派川", "浸水想定", "中州", "避難所",
          "被爆樹木", "DID", "テーマ統合", "M系"],
    time="20-30 分",
    level="中級+",
    data_label="DoBoX dataset_id 35,36,37,42,786,1311 + DID/L15 派生",
    sections=sections,
    script_filename="L85_M1_otagawa_story.py",
)
out_html = LESSONS / "L85_M1_otagawa_story.html"
out_html.write_text(html, encoding="utf-8")
print(f"  HTML 保存: {out_html.name} ({out_html.stat().st_size//1024} KB)",
      flush=True)


print(f"\n=== 全処理完了: {time.time()-t_all:.1f} 秒 ===", flush=True)