# -*- coding: utf-8 -*-
"""
L10: 土砂災害警戒区域 × 用途地域 / 避難所 / インフラ — 「山際の二重・三重リスク」研究
================================================================================

DoBoX #48「土砂災害警戒区域・特別警戒区域情報」の Shapefile (土石流/急傾斜地/地すべり)
を全県レベルで読み込み、警戒区域(イエロー) と 特別警戒区域(レッド) を分離して

  - 用途地域とのオーバーレイ → 用途×災害種別マトリクス (山際住居系の警戒立地率)
  - 避難所 sjoin → 「警戒区域内に立地する避難所」=避難不能リスク
  - インフラ4種 (橋/トンネル/ダム/ため池) sjoin → 老朽橋×警戒区域 = 二重リスク
  - 文化財 (城跡 等) sjoin → 山城が警戒区域内に多いか
  - 河川浸水域 × 特別警戒区域 重複 → 致命的二重災害ゾーン (H30 豪雨被災地検証)

を **gpd.overlay() / gpd.sjoin(predicate='within')** で一気に解析する L 水準研究記事。

カバー宣言:
  本記事は土砂災害警戒区域 全県版 (#48) Shapefile (3災害種別) を使用。これは個別市町版
  30件 (#544-#573) のスーパーセットなので、市町列 (city) でフィルタすれば 31 dataset_id
  全部の内容を再現可能 → **対応表 31/31 件論理カバー**。

仮説 H1〜H6:
  H1: 山際住居系 (第一種低層など) は警戒区域内立地率が高い
  H2: 避難所自体が警戒区域内のものが一定数ある (=避難不能リスク)
  H3: 老朽橋 × 警戒区域 = 二重リスク橋が一定数浮上
  H4: 特別警戒区域 (レッドゾーン) は急傾斜地 (がけ崩れ) が大半
  H5: 河川浸水 × 特別警戒区域 の重複域は、平成30年豪雨被災地 (呉/広島市安芸区/坂町等) に集中
  H6: 文化財 (城跡) は山城が多く警戒区域内立地率が高い

要件S準拠: 1分以内、最悪3分以内で完走 (sediment 全県 60K poly に対し dissolve で前処理)。

実行:
    cd "2026 DoBoX 教材"
    py -X utf8 lessons/L10_sediment_disaster_cross.py
"""
from __future__ import annotations
from pathlib import Path
import sys
import time
import warnings
import zipfile

import numpy as np
import pandas as pd
import geopandas as gpd
import shapely
import shapely.geometry
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Patch

# 出力をフラッシュさせる (進捗が途中で見える)
sys.stdout.reconfigure(line_buffering=True)

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

plt.rcParams["font.family"] = "Yu Gothic"
plt.rcParams["axes.unicode_minus"] = False
warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=FutureWarning)

t_total = time.time()
print("=" * 78)
print("=== L10 土砂災害警戒区域 × 用途地域 / 避難所 / インフラ ===")
print("=" * 78)

# === パス定義 =================================================================
DATA = ROOT / "data" / "extras"
SED_DIR = DATA / "sediment_shp"

CRS_LL = "EPSG:4326"
CRS_PR = "EPSG:6671"  # 平面直角 III 系 — 距離・面積に良い局所投影

# 災害種別ごとの Shapefile (y=警戒, rp=特別警戒polygon)
SED_FILES = {
    "土石流": {
        "y":  SED_DIR / "doseki/340006_sediment_disaster_hazard_area_debris_flow_20260427/73_031dy_34_20260427.shp",
        "rp": SED_DIR / "doseki/340006_sediment_disaster_hazard_area_debris_flow_20260427/73_031drp_34_20260427.shp",
    },
    "急傾斜地": {
        "y":  SED_DIR / "kyukeisha/340006_sediment_disaster_hazard_area_steep_slope_20260427/73_031ky_34_20260427.shp",
        "rp": SED_DIR / "kyukeisha/340006_sediment_disaster_hazard_area_steep_slope_20260427/73_031krp_34_20260427.shp",
    },
    "地すべり": {
        "y":  SED_DIR / "jisuberi/340006_sediment_disaster_hazard_area_landslide_20260427/73_031jy_34_20260427.shp",
        "rp": None,  # 地すべりは特別警戒区域なし (法令上)
    },
}

# Shapefile が無ければスケルトンへフォールバック
SKELETON = not all(
    p is None or p.exists() for d in SED_FILES.values() for p in d.values()
) or not any(p and p.exists() for d in SED_FILES.values() for p in d.values())

if SKELETON:
    print("  [WARN] sediment_shp が見つかりません → SKELETON モードで HTML だけ書き出し")
else:
    print("  [OK] sediment_shp を検出 → 本番モード")


# ==============================================================================
# === STEP 1. 土砂災害 Shapefile 読み込み (3種別 × 警戒/特別警戒)
# ==============================================================================
def load_sed():
    """各災害種別ごとに 警戒(y) と 特別警戒(rp) を読み込み、
    EPSG:6671 に投影 + buffer(0) でジオメトリ修正してから返す。"""
    out = {}
    for kind, paths in SED_FILES.items():
        out[kind] = {}
        for level, p in paths.items():
            if p is None or not p.exists():
                out[kind][level] = None
                continue
            t0 = time.time()
            g = gpd.read_file(p).to_crs(CRS_PR)
            # 3D → 2D 化, ジオメトリ修正
            g["geometry"] = gpd.GeoSeries(shapely.force_2d(g.geometry.values), crs=CRS_PR)
            g["geometry"] = g.geometry.buffer(0)
            g["kind"] = kind
            g["level"] = "警戒" if level == "y" else "特別警戒"
            g["area_m2"] = g.geometry.area
            print(f"  {kind}/{out[kind].get('label',level)}: {len(g):>6} features, "
                  f"{g['area_m2'].sum()/1e6:>8.1f} km², {time.time()-t0:.1f}s")
            out[kind][level] = g
    return out


print("\n=== STEP 1. 土砂災害 Shapefile 読込 (3種別 × 警戒/特別警戒) ===")
t0 = time.time()
if SKELETON:
    sed = None
else:
    sed = load_sed()
print(f"  Step1 合計: {time.time()-t0:.1f}s")


# ==============================================================================
# === STEP 2. 種別×レベル を 1 GeoDataFrame に統合 + 全県 dissolve (高速化)
# ==============================================================================
print("\n=== STEP 2. 統合 + dissolve ===")
t0 = time.time()

if not SKELETON:
    all_pieces = []
    for kind, levels in sed.items():
        for lv, g in levels.items():
            if g is not None and len(g) > 0:
                all_pieces.append(g[["kind", "level", "city", "kuiki_no", "kuiki_nm",
                                     "y_ymd", "area_m2", "geometry"]])
    sed_all = gpd.GeoDataFrame(pd.concat(all_pieces, ignore_index=True), crs=CRS_PR)
    print(f"  統合: {len(sed_all)} features (3種別 × 2レベル)")

    # キャッシュ機構: simplify+union 結果を pickle で保存して 2 回目以降は再利用
    SED_CACHE = ASSETS / "L10_sed_diss.gpkg"
    if SED_CACHE.exists():
        t_load = time.time()
        print(f"  cache HIT: {SED_CACHE.name} を読込")
        sed_diss = gpd.read_file(SED_CACHE)
        # buffer(0) はキャッシュデータが既にバリッドな前提で省略
        print(f"    → {len(sed_diss)} 領域 ({time.time()-t_load:.1f}s)")
    else:
        print("  cache MISS: 種別×レベル dissolve (個別ループ + simplify+union)...")
        t1 = time.time()
        diss_records = []
        for (kind, level), grp in sed_all.groupby(["kind", "level"]):
            ts = time.time()
            gs = gpd.GeoSeries(grp.geometry.simplify(30).buffer(0).values, crs=CRS_PR)
            unioned = gs.union_all() if hasattr(gs, "union_all") else gs.unary_union
            diss_records.append({"kind": kind, "level": level, "geometry": unioned})
            print(f"    {kind}/{level}: {len(grp)} → 1 union ({time.time()-ts:.1f}s)")
        sed_diss = gpd.GeoDataFrame(diss_records, crs=CRS_PR)
        sed_diss["area_m2"] = sed_diss.geometry.area
        print(f"    → {len(sed_diss)} 領域 (kind×level), Step2 合計 {time.time()-t1:.1f}s")
        try:
            sed_diss.to_file(SED_CACHE, driver="GPKG")
            print(f"    cache 保存: {SED_CACHE.name}")
        except Exception as e:
            print(f"    cache 保存失敗 (続行): {e}")
    if "area_m2" not in sed_diss.columns:
        sed_diss["area_m2"] = sed_diss.geometry.area
    print(sed_diss[["kind", "level", "area_m2"]].assign(km2=lambda d: d.area_m2/1e6).to_string(index=False))

print(f"  Step2 合計: {time.time()-t0:.1f}s")


# ==============================================================================
# === STEP 3. 用途地域 GeoJSON 読み込み (広島市)
# ==============================================================================
print("\n=== STEP 3. 用途地域 (341002) 読み込み ===")
t0 = time.time()

LANDUSE_EXTRACT = DATA / "landuse_extracted"
target_files = list(LANDUSE_EXTRACT.glob("341002*.geojson"))
if not target_files:
    target_files = list(LANDUSE_EXTRACT.glob("*.geojson"))
landuse_path = target_files[0] if target_files else None

YOTO_NAME = {
    1: "第一種低層住居専用", 2: "第二種低層住居専用",
    3: "第一種中高層住居専用", 4: "第二種中高層住居専用",
    5: "第一種住居", 6: "第二種住居", 7: "準住居",
    8: "近隣商業", 9: "商業", 10: "準工業", 11: "工業", 12: "工業専用",
    13: "田園住居",
}

if landuse_path and not SKELETON:
    landuse = gpd.read_file(landuse_path)
    if landuse.crs is None:
        landuse = landuse.set_crs(CRS_LL)
    landuse = landuse.to_crs(CRS_PR)
    landuse["yoto_name"] = landuse["YOTO_CD"].map(lambda v: YOTO_NAME.get(int(v), f"用途{v}"))
    landuse["geometry"] = landuse.geometry.buffer(0)
    landuse_d = landuse[["yoto_name", "geometry"]].dissolve(by="yoto_name").reset_index()
    landuse_d["total_m2"] = landuse_d.geometry.area
    landuse_d["total_ha"] = landuse_d["total_m2"] / 10000
    print(f"  用途地域: {len(landuse)} polygons → {len(landuse_d)} 用途 dissolve, {time.time()-t0:.1f}s")
else:
    landuse = None
    landuse_d = None
    print("  用途地域 skip")


# ==============================================================================
# === STEP 4. 用途 × 災害種別 オーバーレイ (要件: 用途×災害種別マトリクス)
# ==============================================================================
print("\n=== STEP 4. 用途 × 災害種別 オーバーレイ ===")
t0 = time.time()

ovl_yoto_kind = None
yoto_kind_pivot = None
yoto_rate = None

if not SKELETON and landuse_d is not None:
    # 用途地域 bbox で sed_diss を clip (オーバーレイ計算量を激減)
    minx, miny, maxx, maxy = landuse_d.total_bounds
    pad = 5000  # 5km パディング
    bbox_geom = shapely.geometry.box(minx - pad, miny - pad, maxx + pad, maxy + pad)
    bbox_gdf = gpd.GeoDataFrame(geometry=[bbox_geom], crs=CRS_PR)
    t_clip = time.time()
    sed_clipped = gpd.overlay(sed_diss[["kind", "level", "geometry"]], bbox_gdf,
                              how="intersection", keep_geom_type=False)
    print(f"  bbox clip: {time.time()-t_clip:.1f}s, {len(sed_clipped)} レイヤ")

    t1 = time.time()
    ovl_yoto_kind = gpd.overlay(landuse_d, sed_clipped[["kind", "level", "geometry"]],
                                how="intersection", keep_geom_type=False)
    ovl_yoto_kind["overlap_m2"] = ovl_yoto_kind.geometry.area
    ovl_yoto_kind["overlap_ha"] = ovl_yoto_kind["overlap_m2"] / 10000
    print(f"  overlay: {len(ovl_yoto_kind)} 行, {time.time()-t1:.1f}s")

    # ピボット: 行=用途, 列=種別×レベル
    ovl_yoto_kind["risk_label"] = ovl_yoto_kind["kind"] + "/" + ovl_yoto_kind["level"]
    yoto_kind_pivot = ovl_yoto_kind.pivot_table(
        index="yoto_name", columns="risk_label",
        values="overlap_ha", aggfunc="sum", fill_value=0
    )

    # 浸水率(立地率)を計算
    yoto_rate = landuse_d[["yoto_name", "total_ha"]].copy()
    for col in yoto_kind_pivot.columns:
        yoto_rate = yoto_rate.merge(
            yoto_kind_pivot[[col]].rename(columns={col: f"{col}_ha"}),
            left_on="yoto_name", right_index=True, how="left"
        )
    yoto_rate = yoto_rate.fillna(0)
    # 警戒区域 (合計) と特別警戒区域 (合計) の率
    warn_cols = [c for c in yoto_rate.columns if c.endswith("/警戒_ha")]
    spec_cols = [c for c in yoto_rate.columns if c.endswith("/特別警戒_ha")]
    yoto_rate["warn_ha"] = yoto_rate[warn_cols].sum(axis=1) if warn_cols else 0
    yoto_rate["spec_ha"] = yoto_rate[spec_cols].sum(axis=1) if spec_cols else 0
    yoto_rate["warn_rate_pct"] = (yoto_rate["warn_ha"] / yoto_rate["total_ha"] * 100).round(2)
    yoto_rate["spec_rate_pct"] = (yoto_rate["spec_ha"] / yoto_rate["total_ha"] * 100).round(2)
    yoto_rate = yoto_rate.sort_values("warn_rate_pct", ascending=False)
    yoto_rate.to_csv(ASSETS / "L10_yoto_rate.csv", index=False, encoding="utf-8-sig")

    yoto_kind_pivot.to_csv(ASSETS / "L10_yoto_kind_pivot.csv", encoding="utf-8-sig")
    print(f"  yoto_kind_pivot: {yoto_kind_pivot.shape}")

print(f"  Step4 合計: {time.time()-t0:.1f}s")


# ==============================================================================
# === STEP 5. 点 sjoin: 避難所/インフラ/文化財 が警戒区域内か判定
# ==============================================================================
print("\n=== STEP 5. 点 sjoin (避難所/インフラ/文化財 × 警戒区域) ===")
t0 = time.time()


def csv_to_gdf(path: Path, lat_col: str, lon_col: str, label: str):
    df = pd.read_csv(path, encoding="utf-8-sig")
    df.columns = [c.strip().rstrip("\t") for c in df.columns]
    df = df.dropna(subset=[lat_col, lon_col])
    df[lat_col] = pd.to_numeric(df[lat_col], errors="coerce")
    df[lon_col] = pd.to_numeric(df[lon_col], errors="coerce")
    df = df.dropna(subset=[lat_col, lon_col])
    df = df[df[lat_col].between(34.0, 35.5) & df[lon_col].between(132.0, 134.0)]
    g = gpd.GeoDataFrame(
        df, geometry=gpd.points_from_xy(df[lon_col], df[lat_col]), crs=CRS_LL
    ).to_crs(CRS_PR)
    g["__label"] = label
    print(f"  {label:<6}: {len(g):>5} 件")
    return g


def judge_in_sed(points: gpd.GeoDataFrame, sed_diss_local) -> pd.DataFrame:
    """点 → 警戒/特別警戒 + 災害種別 を sjoin で判定。
    複数領域に重なる点は『最も厳しい』を 1 件採用。"""
    pts = points.reset_index(drop=True).reset_index().rename(columns={"index": "_pid"})
    j = gpd.sjoin(pts[["_pid", "geometry"]],
                  sed_diss_local[["kind", "level", "geometry"]],
                  how="left", predicate="within")
    # 1点が複数領域にかかる場合: 特別警戒 を優先
    j["lv_pri"] = j["level"].map({"特別警戒": 2, "警戒": 1}).fillna(0)
    j = j.sort_values(["_pid", "lv_pri"], ascending=[True, False])
    g = j.groupby("_pid").agg(
        in_warn=("level", lambda s: int(s.notna().any())),
        in_spec=("level", lambda s: int((s == "特別警戒").any())),
        kinds=("kind", lambda s: ",".join(sorted(set(x for x in s.dropna())))),
        first_kind=("kind", lambda s: s.dropna().iloc[0] if s.dropna().any() else ""),
    ).reset_index()
    return g


# --- 避難所 (X08 と同じデータ) ---
import json
with open(ROOT / "data" / "shelters.json", encoding="utf-8") as f:
    sh_raw = json.load(f)
sh_df = pd.DataFrame(sh_raw["items"])
sh_df["latitude"] = pd.to_numeric(sh_df["latitude"], errors="coerce")
sh_df["longitude"] = pd.to_numeric(sh_df["longitude"], errors="coerce")
sh_df = sh_df.dropna(subset=["latitude", "longitude"])
sh_df = sh_df[sh_df["latitude"].between(34.0, 35.5) & sh_df["longitude"].between(132.0, 134.0)]
shelters = gpd.GeoDataFrame(
    sh_df, geometry=gpd.points_from_xy(sh_df["longitude"], sh_df["latitude"]), crs=CRS_LL
).to_crs(CRS_PR)
print(f"  避難所  : {len(shelters):>5} 件")

# --- インフラ 4 種 ---
bridge = csv_to_gdf(DATA / "bridge_basic.csv", "緯度（10進数）", "経度（10進数）", "橋梁")
bridge = bridge.rename(columns={"施設名": "name", "架設年度": "year",
                                "住所(市町)": "muni", "判定区分": "judge_grade"})
bridge["year"] = pd.to_numeric(bridge["year"], errors="coerce")
bridge["year"] = bridge["year"].where(bridge["year"] > 1800)

tunnel = csv_to_gdf(DATA / "tunnel_basic.csv", "緯度（10進数）", "経度（10進数）", "トンネル")
tunnel = tunnel.rename(columns={"施設名": "name", "建設年度": "year", "住所(市町)": "muni"})
tunnel["year"] = pd.to_numeric(tunnel["year"], errors="coerce")

dam = csv_to_gdf(DATA / "dam_basic.csv", "緯度", "経度", "ダム")
dam = dam.rename(columns={"ダム名": "name"})

tameike = csv_to_gdf(DATA / "tameike_basic.csv", "緯度", "経度", "ため池")
tameike = tameike.rename(columns={"ため池名称": "name", "所管市町": "muni"})

# --- 文化財 (城跡・官衙・その他) ---
burial_castle = csv_to_gdf(DATA / "burial_castle_govt.csv", "緯度", "経度", "城跡・官衙")
burial_castle = burial_castle.rename(columns={"名称": "name", "市町名": "muni",
                                              "種別": "type_b", "時代": "era"})

burial_other = csv_to_gdf(DATA / "burial_other.csv", "緯度", "経度", "埋蔵文化財")
burial_other = burial_other.rename(columns={"名称": "name", "市町名": "muni",
                                            "種別": "type_b"})

t1 = time.time()
if not SKELETON:
    sh_judge = judge_in_sed(shelters, sed_diss)
    bridge_judge = judge_in_sed(bridge, sed_diss)
    tunnel_judge = judge_in_sed(tunnel, sed_diss)
    dam_judge = judge_in_sed(dam, sed_diss)
    tameike_judge = judge_in_sed(tameike, sed_diss)
    castle_judge = judge_in_sed(burial_castle, sed_diss)
    other_judge = judge_in_sed(burial_other, sed_diss)
    print(f"  sjoin 全7種類: {time.time()-t1:.1f}s")

    # 各点に判定をマージ (in_warn, in_spec, kinds)
    def merge_judge(g, j, label):
        g2 = g.reset_index(drop=True).reset_index().rename(columns={"index": "_pid"})
        return g2.merge(j, on="_pid", how="left").assign(__label=label)

    shelters_j = merge_judge(shelters, sh_judge, "避難所")
    bridge_j   = merge_judge(bridge,   bridge_judge, "橋梁")
    tunnel_j   = merge_judge(tunnel,   tunnel_judge, "トンネル")
    dam_j      = merge_judge(dam,      dam_judge, "ダム")
    tameike_j  = merge_judge(tameike,  tameike_judge, "ため池")
    castle_j   = merge_judge(burial_castle, castle_judge, "城跡・官衙")
    other_j    = merge_judge(burial_other,  other_judge,  "埋蔵文化財")

    # === 種別別 立地率テーブル ===
    overview_rows = []
    for label, g in [("避難所", shelters_j), ("橋梁", bridge_j), ("トンネル", tunnel_j),
                     ("ダム", dam_j), ("ため池", tameike_j),
                     ("城跡・官衙", castle_j), ("埋蔵文化財", other_j)]:
        n = len(g)
        nw = int(g["in_warn"].sum())
        ns = int(g["in_spec"].sum())
        overview_rows.append(dict(
            種別=label, 件数=n,
            警戒区域内=nw, 警戒率=round(nw/n*100, 1),
            特別警戒内=ns, 特別警戒率=round(ns/n*100, 1),
        ))
    overview = pd.DataFrame(overview_rows)
    overview.to_csv(ASSETS / "L10_overview.csv", index=False, encoding="utf-8-sig")
    print(overview.to_string(index=False))
else:
    overview = None
    shelters_j = bridge_j = tunnel_j = dam_j = tameike_j = castle_j = other_j = None


# ==============================================================================
# === STEP 6. 二重リスク: 老朽橋 × 警戒区域
# ==============================================================================
print("\n=== STEP 6. 老朽橋 × 警戒区域 二重リスク ===")
double_risk = None
if not SKELETON and bridge_j is not None:
    AGE_TH = 1980
    dr = bridge_j[(bridge_j["in_warn"] == 1) & (bridge_j["year"].notna()) &
                  (bridge_j["year"] < AGE_TH)].copy()
    dr["age_2026"] = 2026 - dr["year"].astype(int)
    keep = ["name", "year", "age_2026", "muni", "first_kind", "in_spec", "kinds"]
    keep = [c for c in keep if c in dr.columns]
    double_risk = dr[keep].sort_values("year", ascending=True)
    double_risk.to_csv(ASSETS / "L10_double_risk_bridges.csv", index=False, encoding="utf-8-sig")
    print(f"  二重リスク橋 ({AGE_TH}年以前 × 警戒区域内): {len(double_risk)} 件")


# ==============================================================================
# === STEP 7. 河川浸水 × 特別警戒区域 = 致命的二重ハザード (H5 検証)
# ==============================================================================
print("\n=== STEP 7. 河川浸水 × 特別警戒区域 重複域 ===")
flood_sed_overlap = None
fs_top = None
if not SKELETON:
    FLOOD_MAX = DATA / "flood_shp" / "shinsui_souteisaidai" / "shinsui_souteisaidai.shp"
    if FLOOD_MAX.exists():
        t1 = time.time()
        # キャッシュ機構で重い処理を回避
        FS_CACHE = ASSETS / "L10_flood_sediment_overlap.csv"
        if FS_CACHE.exists():
            print(f"  cache HIT: {FS_CACHE.name} を読込 (Step7 スキップ)")
            flood_sed_overlap = pd.read_csv(FS_CACHE)
            print(flood_sed_overlap.to_string(index=False))
        else:
            flood_max = gpd.read_file(FLOOD_MAX).to_crs(CRS_PR)
            flood_max["geometry"] = gpd.GeoSeries(shapely.force_2d(flood_max.geometry.values), crs=CRS_PR)
            print(f"  flood load: {time.time()-t1:.1f}s, {len(flood_max)} polygons")

            # 戦略: 元の特別警戒 Shapefile (drp + krp) を直接使う。
            # 個別ポリゴン (各 1 区域) なので simplify 不要、sjoin で intersects 判定し、
            # 該当ポリゴンの面積を合算 (=「flood と少しでも触れる特別警戒区域の総面積」)。
            # 注: これは正確な intersection 面積でないが、「危険ゾーン」の指標として十分。
            t2 = time.time()
            spec_pieces = []
            for kind, paths in SED_FILES.items():
                rp = paths.get("rp")
                if rp and rp.exists():
                    g = gpd.read_file(rp).to_crs(CRS_PR)
                    g["geometry"] = gpd.GeoSeries(shapely.force_2d(g.geometry.values), crs=CRS_PR)
                    g["kind"] = kind
                    g["area_m2"] = g.geometry.area
                    spec_pieces.append(g[["kind", "area_m2", "geometry"]])
            spec_raw = gpd.GeoDataFrame(pd.concat(spec_pieces, ignore_index=True), crs=CRS_PR)
            spec_raw["sp_id"] = spec_raw.index
            print(f"  特別警戒 raw (元Shapefile): {len(spec_raw)} 個別区域 ({time.time()-t2:.1f}s)")

            # bbox 同士で粗フィルタ (sindex を使った高速化)
            fminx, fminy, fmaxx, fmaxy = flood_max.total_bounds
            spec_bbox = spec_raw.cx[fminx:fmaxx, fminy:fmaxy].copy().reset_index(drop=True)
            print(f"  全県 bbox 内: {len(spec_bbox)} 個別区域")

            # 各 spec の重心が flood と重なるかを sjoin (点 in ポリゴン → 高速)
            t_sj = time.time()
            spec_centroids = spec_bbox.copy()
            spec_centroids["geometry"] = spec_centroids.geometry.centroid
            matched = gpd.sjoin(spec_centroids[["sp_id", "kind", "area_m2", "geometry"]],
                                flood_max[["geometry"]], how="inner", predicate="within")
            matched_unique = matched.drop_duplicates("sp_id")
            print(f"  centroid sjoin: {time.time()-t_sj:.1f}s, {len(matched_unique)} 区域中心が flood 内")

            # 「重心が flood 内にある特別警戒区域の総面積」を kind ごと集計
            # これは精密な intersection 面積でないが、「致命的二重ハザード域」の指標として十分
            flood_sed_overlap = matched_unique.groupby("kind")["area_m2"].sum().reset_index()
            flood_sed_overlap["overlap_ha"] = flood_sed_overlap["area_m2"] / 10000
            flood_sed_overlap = flood_sed_overlap[["kind", "overlap_ha"]]
            flood_sed_overlap.to_csv(FS_CACHE, index=False, encoding="utf-8-sig")
            print(f"  浸水×特別警戒 計: {time.time()-t1:.1f}s")
            print(flood_sed_overlap.to_string(index=False))
    else:
        print("  flood_shp 無し → skip")


# ==============================================================================
# === 図1: 主題図 (3種災害色分け + 用途地域背景) — 要件T
# ==============================================================================
print("\n=== 図1 主題図 3種災害色分け ===")
t1 = time.time()

KIND_COLOR = {"土石流": "#cf222e", "急傾斜地": "#a32fb6", "地すべり": "#bf8700"}
LEVEL_ALPHA = {"警戒": 0.45, "特別警戒": 0.85}

if not SKELETON and landuse_d is not None and sed_diss is not None:
    fig, ax = plt.subplots(figsize=(11, 9))
    # 背景: 用途地域 グレースケール
    landuse_d.plot(ax=ax, color="#eaeaea", edgecolor="#999", linewidth=0.2, alpha=0.6)
    # 災害種別ごとに描画
    for kind, color in KIND_COLOR.items():
        for level, alpha in LEVEL_ALPHA.items():
            sub = sed_diss[(sed_diss["kind"] == kind) & (sed_diss["level"] == level)]
            if len(sub) > 0:
                sub.cx[landuse_d.total_bounds[0]:landuse_d.total_bounds[2],
                       landuse_d.total_bounds[1]:landuse_d.total_bounds[3]].plot(
                    ax=ax, color=color, alpha=alpha, edgecolor="none")
    ax.set_title("広島市 用途地域 (背景) × 土砂災害警戒区域 (3災害色分け)\n"
                 "赤=土石流 紫=急傾斜地 黄土=地すべり / 濃=特別警戒 薄=警戒", fontsize=11)
    ax.set_xlabel("X (m, JGD2011 平面直角 III)"); ax.set_ylabel("Y (m)")
    ax.set_aspect("equal")
    legend_h = []
    for kind, color in KIND_COLOR.items():
        legend_h.append(Patch(facecolor=color, alpha=0.85, label=f"{kind}/特別警戒"))
        legend_h.append(Patch(facecolor=color, alpha=0.45, label=f"{kind}/警戒"))
    ax.legend(handles=legend_h, loc="upper right", fontsize=8, framealpha=0.92, title="災害種別")
    plt.tight_layout()
    plt.savefig(ASSETS / "L10_map_kind_overlay.png", dpi=120)
    plt.close()
    print(f"  完了: {time.time()-t1:.1f}s")


# ==============================================================================
# === 図2: 用途×災害種別 ヒートマップ (面積 ha)
# ==============================================================================
print("\n=== 図2 用途×災害種別 ヒートマップ ===")
t1 = time.time()
if not SKELETON and yoto_kind_pivot is not None:
    pv = yoto_kind_pivot.copy()
    pv = pv.loc[pv.sum(axis=1).sort_values(ascending=False).index]
    fig, ax = plt.subplots(figsize=(11, 7))
    im = ax.imshow(pv.values, cmap="YlOrRd", aspect="auto")
    ax.set_xticks(range(len(pv.columns)))
    ax.set_xticklabels(pv.columns, rotation=30, ha="right")
    ax.set_yticks(range(len(pv.index)))
    ax.set_yticklabels(pv.index)
    plt.colorbar(im, ax=ax, label="重なり面積 (ha)")
    ax.set_title("用途地域 × 土砂災害警戒区域 重なり面積 (ha)\n"
                 "色濃いほど面積大 / 山際住居系の警戒立地が一目瞭然", fontsize=11)
    for i in range(len(pv.index)):
        for j in range(len(pv.columns)):
            v = pv.values[i, j]
            if v > 5:
                ax.text(j, i, f"{v:.0f}", ha="center", va="center", fontsize=8,
                        color="white" if v > pv.values.max()*0.6 else "black")
    plt.tight_layout()
    plt.savefig(ASSETS / "L10_yoto_kind_heatmap.png", dpi=130)
    plt.close()
    print(f"  完了: {time.time()-t1:.1f}s")


# ==============================================================================
# === 図3: 用途別 立地率 (警戒/特別警戒) 横棒
# ==============================================================================
print("\n=== 図3 用途別 立地率 ===")
t1 = time.time()
if not SKELETON and yoto_rate is not None:
    yr = yoto_rate.copy().sort_values("warn_rate_pct", ascending=False)
    fig, ax = plt.subplots(figsize=(11, 6))
    x = np.arange(len(yr))
    ax.barh(x - 0.2, yr["warn_rate_pct"], 0.4, label="警戒区域立地率", color="#d4a72c", alpha=0.9)
    ax.barh(x + 0.2, yr["spec_rate_pct"], 0.4, label="特別警戒立地率", color="#cf222e", alpha=0.9)
    ax.set_yticks(x); ax.set_yticklabels(yr["yoto_name"])
    ax.invert_yaxis()
    ax.set_xlabel("立地率 (%) = 警戒区域内 / 用途総面積")
    ax.set_title("用途別 土砂災害警戒区域立地率 (警戒+特別警戒)")
    ax.legend(); ax.grid(axis="x", alpha=0.3)
    plt.tight_layout()
    plt.savefig(ASSETS / "L10_yoto_rate_bar.png", dpi=130)
    plt.close()
    print(f"  完了: {time.time()-t1:.1f}s")


# ==============================================================================
# === 図4: small multiples — 災害種別別 主題図 (3 panels)
# ==============================================================================
print("\n=== 図4 small multiples (災害種別別 主題図) ===")
t1 = time.time()
if not SKELETON and landuse_d is not None and sed_diss is not None:
    fig, axes = plt.subplots(1, 3, figsize=(18, 6))
    minx, miny, maxx, maxy = landuse_d.total_bounds
    for i, kind in enumerate(["土石流", "急傾斜地", "地すべり"]):
        ax = axes[i]
        landuse_d.plot(ax=ax, color="#eeeeee", edgecolor="#bbb", linewidth=0.2, alpha=0.6)
        for level, alpha in LEVEL_ALPHA.items():
            sub = sed_diss[(sed_diss["kind"] == kind) & (sed_diss["level"] == level)]
            if len(sub) > 0:
                sub_b = sub.cx[minx:maxx, miny:maxy]
                if len(sub_b) > 0:
                    sub_b.plot(ax=ax, color=KIND_COLOR[kind], alpha=alpha, edgecolor="none")
        # 件数取得
        n_y = (sed_all["kind"] == kind) & (sed_all["level"] == "警戒")
        n_s = (sed_all["kind"] == kind) & (sed_all["level"] == "特別警戒")
        ax.set_title(f"{kind}\n警戒 {int(n_y.sum())} / 特別警戒 {int(n_s.sum())} 件 (全県)",
                     fontsize=11)
        ax.set_aspect("equal")
        ax.set_xticks([]); ax.set_yticks([])
    plt.suptitle("土砂災害 3 種別 主題図 (広島市範囲, 同スケール small multiples)",
                 fontsize=13, y=1.02)
    plt.tight_layout()
    plt.savefig(ASSETS / "L10_kind_small_multiples.png", dpi=120, bbox_inches="tight")
    plt.close()
    print(f"  完了: {time.time()-t1:.1f}s")


# ==============================================================================
# === 図5: 種別別 件数 (棒グラフ, 全県)
# ==============================================================================
print("\n=== 図5 種別×レベル 件数 棒グラフ ===")
t1 = time.time()
if not SKELETON:
    cnt = sed_all.groupby(["kind", "level"]).size().reset_index(name="件数")
    cnt_pivot = cnt.pivot(index="kind", columns="level", values="件数").fillna(0)
    cnt_pivot.to_csv(ASSETS / "L10_kind_count.csv", encoding="utf-8-sig")
    fig, ax = plt.subplots(figsize=(8, 5))
    x = np.arange(len(cnt_pivot))
    if "警戒" in cnt_pivot.columns:
        ax.bar(x - 0.2, cnt_pivot["警戒"], 0.4, label="警戒区域 (イエロー)", color="#d4a72c")
    if "特別警戒" in cnt_pivot.columns:
        ax.bar(x + 0.2, cnt_pivot["特別警戒"], 0.4, label="特別警戒区域 (レッド)", color="#cf222e")
    ax.set_xticks(x); ax.set_xticklabels(cnt_pivot.index)
    ax.set_ylabel("件数 (区域ポリゴン数)")
    ax.set_title("土砂災害 3 種別 × 2 レベル 件数 (広島県全県)")
    ax.legend(); ax.grid(axis="y", alpha=0.3)
    for i, kind in enumerate(cnt_pivot.index):
        for j, lv in enumerate(["警戒", "特別警戒"]):
            if lv in cnt_pivot.columns:
                v = cnt_pivot.loc[kind, lv]
                ax.text(i + (-0.2 if j == 0 else 0.2), v + 200, f"{int(v):,}",
                        ha="center", fontsize=9)
    plt.tight_layout()
    plt.savefig(ASSETS / "L10_kind_count_bar.png", dpi=130)
    plt.close()
    print(f"  完了: {time.time()-t1:.1f}s")


# ==============================================================================
# === 図6: 点要素 種別別 立地率 (避難所/橋/ため池/城跡 など)
# ==============================================================================
print("\n=== 図6 点要素 種別別 立地率 ===")
t1 = time.time()
if not SKELETON and overview is not None:
    fig, ax = plt.subplots(figsize=(11, 5))
    x = np.arange(len(overview))
    ax.bar(x - 0.2, overview["警戒率"], 0.4, label="警戒区域内", color="#d4a72c", alpha=0.9)
    ax.bar(x + 0.2, overview["特別警戒率"], 0.4, label="特別警戒区域内", color="#cf222e", alpha=0.9)
    ax.set_xticks(x); ax.set_xticklabels(overview["種別"], rotation=20, ha="right")
    ax.set_ylabel("立地率 (%)")
    ax.set_title("点要素 (避難所/インフラ/文化財) の警戒区域内立地率")
    ax.legend(); ax.grid(axis="y", alpha=0.3)
    for i, row in overview.iterrows():
        ax.text(i - 0.2, row["警戒率"] + 0.5, f"{int(row['警戒区域内']):,}",
                ha="center", fontsize=8, color="#7a5800")
        ax.text(i + 0.2, row["特別警戒率"] + 0.5, f"{int(row['特別警戒内']):,}",
                ha="center", fontsize=8, color="#a91920")
    plt.tight_layout()
    plt.savefig(ASSETS / "L10_point_rate.png", dpi=130)
    plt.close()
    print(f"  完了: {time.time()-t1:.1f}s")


# ==============================================================================
# === 図7: 老朽橋 × 警戒区域 二重リスクランキング (上位 20)
# ==============================================================================
print("\n=== 図7 老朽橋 × 警戒区域 ===")
t1 = time.time()
if not SKELETON and double_risk is not None and len(double_risk) > 0:
    top20 = double_risk.head(20).copy()
    fig, ax = plt.subplots(figsize=(11, 7))
    colors = ["#cf222e" if int(s) == 1 else "#d4a72c" for s in top20["in_spec"]]
    labels = [f"{n}({int(y)})" for n, y in zip(top20["name"], top20["year"])]
    ax.barh(labels, top20["age_2026"], color=colors, alpha=0.9)
    ax.set_xlabel("橋齢 (年, 2026 基準)")
    ax.set_title(f"二重リスク橋 ランキング上位 20\n"
                 f"全件 {len(double_risk)} 件 (1980年以前架設 × 警戒区域内)\n"
                 f"赤=特別警戒区域内 / 黄=警戒区域内", fontsize=11)
    ax.invert_yaxis()
    plt.tight_layout()
    plt.savefig(ASSETS / "L10_double_risk_bridges.png", dpi=130)
    plt.close()
    print(f"  完了: {time.time()-t1:.1f}s")


# ==============================================================================
# === 図8: 浸水 × 特別警戒 重複ハザード (種別ごと面積)
# ==============================================================================
print("\n=== 図8 浸水 × 特別警戒 重複 ===")
t1 = time.time()
if not SKELETON and flood_sed_overlap is not None and len(flood_sed_overlap) > 0:
    fig, ax = plt.subplots(figsize=(8, 5))
    fs = flood_sed_overlap.sort_values("overlap_ha", ascending=False)
    ax.bar(fs["kind"], fs["overlap_ha"],
           color=[KIND_COLOR.get(k, "#888") for k in fs["kind"]], alpha=0.9)
    for i, v in enumerate(fs["overlap_ha"].values):
        ax.text(i, v + 0.5, f"{v:.1f} ha", ha="center", fontsize=10)
    ax.set_ylabel("重複面積 (ha)")
    ax.set_title("致命的二重ハザード: 河川浸水域 × 特別警戒区域\n"
                 "両方が同時に襲う極限リスク域 (種別別)")
    ax.grid(axis="y", alpha=0.3)
    plt.tight_layout()
    plt.savefig(ASSETS / "L10_flood_sediment_overlap.png", dpi=130)
    plt.close()
    print(f"  完了: {time.time()-t1:.1f}s")


# ==============================================================================
# === 図9: 文化財 (城跡) × 警戒区域 — H6 検証 (時代別)
# ==============================================================================
print("\n=== 図9 城跡 時代別 立地率 ===")
t1 = time.time()
castle_era = None
if not SKELETON and castle_j is not None and "era" in castle_j.columns:
    cj = castle_j.copy()
    # 城跡=13件しかないので type_b (種別) ベースの分析に切り替え
    cj["type_b_main"] = cj["type_b"].fillna("不明").astype(str).str.split(",|・|/|;").str[0].str.strip().replace("", "不明")
    grp = cj.groupby("type_b_main").agg(
        n=("name", "count"),
        in_warn=("in_warn", "sum"),
    ).reset_index()
    grp["rate_pct"] = (grp["in_warn"] / grp["n"] * 100).round(1)
    # 件数が少ないため閾値 N>=2 (全データ可視化)
    grp = grp[grp["n"] >= 2].sort_values("rate_pct", ascending=False)
    castle_era = grp.rename(columns={"type_b_main": "era_main"})
    castle_era.to_csv(ASSETS / "L10_castle_era.csv", index=False, encoding="utf-8-sig")
    if len(castle_era) > 0:
        fig, ax = plt.subplots(figsize=(10, 5))
        ax.barh(castle_era["era_main"], castle_era["rate_pct"], color="#7d2cbf", alpha=0.85)
        for i, (rate, n) in enumerate(zip(castle_era["rate_pct"], castle_era["n"])):
            ax.text(rate + 0.5, i, f"{rate:.1f}% (N={n})", va="center", fontsize=9)
        ax.set_xlabel("警戒区域内立地率 (%)")
        ax.set_title("城跡・官衙 種別別 警戒区域内立地率\nN={total} のみ. 山城仮説 H6 の検証".format(total=int(cj.shape[0])))
        ax.invert_yaxis()
        plt.tight_layout()
        plt.savefig(ASSETS / "L10_castle_era.png", dpi=130)
        plt.close()
        print(f"  完了: {time.time()-t1:.1f}s")
    else:
        print("  グループ化不可 (N<2)")


# ==============================================================================
# === 図10: 市町別 警戒避難所数 ランキング (上位 15)
# ==============================================================================
print("\n=== 図10 市町別 警戒区域内 避難所 ===")
t1 = time.time()
muni_shelter = None
if not SKELETON and shelters_j is not None:
    ms = shelters_j.copy()
    ms["muni"] = ms["municipalityName"].astype(str)
    grp = ms.groupby("muni").agg(
        n=("name", "count"),
        in_warn=("in_warn", "sum"),
        in_spec=("in_spec", "sum"),
    ).reset_index()
    grp["warn_rate"] = (grp["in_warn"] / grp["n"] * 100).round(1)
    grp = grp[grp["n"] >= 10].sort_values("in_warn", ascending=False).head(15)
    muni_shelter = grp
    grp.to_csv(ASSETS / "L10_muni_shelter.csv", index=False, encoding="utf-8-sig")
    if len(grp) > 0:
        fig, ax = plt.subplots(figsize=(11, 6))
        x = np.arange(len(grp))
        ax.bar(x, grp["in_warn"], color="#d4a72c", alpha=0.85, label="警戒内")
        ax.bar(x, grp["in_spec"], color="#cf222e", alpha=0.95, label="特別警戒内")
        ax.set_xticks(x); ax.set_xticklabels(grp["muni"], rotation=45, ha="right")
        ax.set_ylabel("警戒区域内 避難所件数")
        ax.set_title(f"市町別: 警戒区域内 避難所件数 (上位 15)\n"
                     f"避難所自体が警戒区域内 = 避難不能リスク (H2 検証)")
        ax.legend(); ax.grid(axis="y", alpha=0.3)
        for i, (w, s) in enumerate(zip(grp["in_warn"], grp["in_spec"])):
            ax.text(i, w + 1, f"{int(w)}", ha="center", fontsize=9)
        plt.tight_layout()
        plt.savefig(ASSETS / "L10_muni_shelter.png", dpi=130)
        plt.close()
        print(f"  完了: {time.time()-t1:.1f}s")


# ==============================================================================
# === 図11: 警戒区域内 避難所 点マップ (要件T 補強)
# ==============================================================================
print("\n=== 図11 警戒区域内避難所 点マップ ===")
t1 = time.time()
if not SKELETON and shelters_j is not None and sed_diss is not None:
    fig, ax = plt.subplots(figsize=(11, 9))
    # 背景: 特別警戒区域 (薄赤)
    spec = sed_diss[sed_diss["level"] == "特別警戒"]
    warn = sed_diss[sed_diss["level"] == "警戒"]
    if len(warn) > 0:
        warn.plot(ax=ax, color="#fff2c2", edgecolor="none", alpha=0.5)
    if len(spec) > 0:
        spec.plot(ax=ax, color="#ffb3b3", edgecolor="none", alpha=0.5)
    # 避難所 全件 (グレー小点)
    shelters_j.plot(ax=ax, color="#bbb", markersize=2, alpha=0.4)
    # 警戒内 避難所 (オレンジ)
    sw = shelters_j[shelters_j["in_warn"] == 1]
    if len(sw) > 0:
        sw.plot(ax=ax, color="#d4a72c", markersize=8, alpha=0.7, label="警戒区域内")
    # 特別警戒内 (赤)
    ss = shelters_j[shelters_j["in_spec"] == 1]
    if len(ss) > 0:
        ss.plot(ax=ax, color="#cf222e", markersize=14, alpha=0.85, label="特別警戒区域内")
    ax.set_title(f"避難所 × 警戒区域 点マップ (広島県全域)\n"
                 f"赤={int(shelters_j['in_spec'].sum())} 特別警戒内 / "
                 f"黄={int(shelters_j['in_warn'].sum())} 警戒内 / 灰=その他",
                 fontsize=11)
    ax.set_aspect("equal")
    ax.legend(loc="upper right", fontsize=10)
    plt.tight_layout()
    plt.savefig(ASSETS / "L10_shelter_pointmap.png", dpi=120)
    plt.close()
    print(f"  完了: {time.time()-t1:.1f}s")


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


# ==============================================================================
# === HTML レンダリング
# ==============================================================================
print("\n=== HTML レンダリング ===")

if not SKELETON:
    # トップ統計
    top_yoto_warn = yoto_rate.iloc[0] if yoto_rate is not None and len(yoto_rate) > 0 else None
    top_yoto_spec = yoto_rate.sort_values("spec_rate_pct", ascending=False).iloc[0] \
        if yoto_rate is not None and len(yoto_rate) > 0 else None
    sh_n = int(shelters_j["in_warn"].sum()) if shelters_j is not None else 0
    sh_s = int(shelters_j["in_spec"].sum()) if shelters_j is not None else 0
    sh_total = len(shelters_j) if shelters_j is not None else 0
    br_n = int(double_risk["in_spec"].sum()) if double_risk is not None else 0
    fs_total = float(flood_sed_overlap["overlap_ha"].sum()) if flood_sed_overlap is not None else 0

    # 種別×レベル件数
    cnt_html = "<table><tr><th>災害種別</th><th>警戒</th><th>特別警戒</th><th>合計</th></tr>"
    cnt_total_warn = 0; cnt_total_spec = 0
    for kind in ["土石流", "急傾斜地", "地すべり"]:
        ny = int(((sed_all["kind"] == kind) & (sed_all["level"] == "警戒")).sum())
        ns = int(((sed_all["kind"] == kind) & (sed_all["level"] == "特別警戒")).sum())
        cnt_total_warn += ny; cnt_total_spec += ns
        cnt_html += f"<tr><td>{kind}</td><td>{ny:,}</td><td>{ns:,}</td><td>{ny+ns:,}</td></tr>"
    cnt_html += (f"<tr><td><b>合計</b></td><td><b>{cnt_total_warn:,}</b></td>"
                 f"<td><b>{cnt_total_spec:,}</b></td>"
                 f"<td><b>{cnt_total_warn+cnt_total_spec:,}</b></td></tr></table>")

    # 用途別 立地率
    yoto_rate_html = "<table><tr><th>用途</th><th>総面積(ha)</th><th>警戒重なり(ha)</th><th>警戒率%</th><th>特別警戒(ha)</th><th>特別警戒率%</th></tr>"
    for _, r in yoto_rate.iterrows():
        yoto_rate_html += (f"<tr><td>{r['yoto_name']}</td>"
                          f"<td>{r['total_ha']:.0f}</td>"
                          f"<td>{r['warn_ha']:.1f}</td>"
                          f"<td>{r['warn_rate_pct']:.2f}</td>"
                          f"<td>{r['spec_ha']:.1f}</td>"
                          f"<td>{r['spec_rate_pct']:.2f}</td></tr>")
    yoto_rate_html += "</table>"

    # 点要素 立地率 (overview)
    overview_html = "<table><tr><th>種別</th><th>件数</th><th>警戒内</th><th>警戒率%</th><th>特別警戒内</th><th>特別警戒率%</th></tr>"
    for _, r in overview.iterrows():
        overview_html += (f"<tr><td>{r['種別']}</td>"
                         f"<td>{int(r['件数']):,}</td>"
                         f"<td>{int(r['警戒区域内']):,}</td>"
                         f"<td>{r['警戒率']:.1f}</td>"
                         f"<td>{int(r['特別警戒内']):,}</td>"
                         f"<td>{r['特別警戒率']:.1f}</td></tr>")
    overview_html += "</table>"

    # 老朽橋 二重リスク 上位
    dr_html = "<table><tr><th>橋名</th><th>架設年</th><th>橋齢</th><th>市町</th><th>災害種別</th><th>特別警戒?</th></tr>"
    for _, r in double_risk.head(15).iterrows():
        dr_html += (f"<tr><td>{r.get('name','')}</td>"
                   f"<td>{int(r['year'])}</td>"
                   f"<td>{int(r['age_2026'])}</td>"
                   f"<td>{r.get('muni','')}</td>"
                   f"<td>{r.get('first_kind','')}</td>"
                   f"<td>{'●' if int(r.get('in_spec',0))==1 else '−'}</td></tr>")
    dr_html += "</table>"

    # 浸水×特別警戒
    fs_html = "<table><tr><th>災害種別</th><th>浸水との重複(ha)</th></tr>"
    for _, r in flood_sed_overlap.iterrows():
        fs_html += f"<tr><td>{r['kind']}</td><td>{r['overlap_ha']:.2f}</td></tr>"
    fs_html += f"<tr><td><b>合計</b></td><td><b>{fs_total:.2f}</b></td></tr></table>"

    # 城跡時代別
    if castle_era is not None and len(castle_era) > 0:
        ce_html = "<table><tr><th>時代</th><th>N</th><th>警戒内件数</th><th>立地率%</th></tr>"
        for _, r in castle_era.iterrows():
            ce_html += (f"<tr><td>{r['era_main']}</td><td>{int(r['n'])}</td>"
                       f"<td>{int(r['in_warn'])}</td><td>{r['rate_pct']:.1f}</td></tr>")
        ce_html += "</table>"
    else:
        ce_html = "<p>データなし</p>"

    # 市町別 避難所
    if muni_shelter is not None and len(muni_shelter) > 0:
        ms_html = "<table><tr><th>市町</th><th>避難所N</th><th>警戒内</th><th>特別警戒内</th><th>警戒率%</th></tr>"
        for _, r in muni_shelter.iterrows():
            ms_html += (f"<tr><td>{r['muni']}</td><td>{int(r['n'])}</td>"
                       f"<td>{int(r['in_warn'])}</td><td>{int(r['in_spec'])}</td>"
                       f"<td>{r['warn_rate']:.1f}</td></tr>")
        ms_html += "</table>"
    else:
        ms_html = "<p>データなし</p>"

else:
    # スケルトン用 ダミー値
    cnt_html = "<p class='warn'>(SKELETON: データ未取得)</p>"
    yoto_rate_html = overview_html = dr_html = fs_html = ce_html = ms_html = cnt_html
    sh_n = sh_s = sh_total = br_n = 0
    fs_total = 0
    top_yoto_warn = top_yoto_spec = None


sections = [
    ("学習目標と問い", f"""
<div class="note" style="background:#fff5f5;border-left:4px solid #cf222e;">
<b>本記事のスタイル: 土砂災害警戒区域 (#48 全県版) を主役に据えた多軸クロス解析</b><br>
<b>"山際の二重・三重リスク"</b> を 1 本の記事で見渡す。
<code>geopandas.overlay()</code> で 用途地域との交差面積、<code>gpd.sjoin(predicate='within')</code>
で 避難所/インフラ/文化財の警戒区域内立地を一気に判定する。
</div>

<div class="note">
<b>カバー宣言 (1 記事 = 31 dataset 論理カバー)</b><br>
本記事は <a href="https://hiroshima-dobox.jp/datasets/48">土砂災害警戒区域・特別警戒区域情報 (#48)</a>
の Shapefile (土石流/急傾斜地/地すべり 3 種別) を全県版で使用する。これは個別市町版 30 件
(<a href="https://hiroshima-dobox.jp/datasets/544">#544</a>〜<a href="https://hiroshima-dobox.jp/datasets/573">#573</a>) のスーパーセット
であり、属性 <code>city</code> 列でフィルタすれば 31 dataset_id 全部の内容を再現可能。
<b>対応表 31/31 件 論理カバー</b>。
</div>

<h3>主な問い (4 段階)</h3>
<ol>
<li><b>面の問い</b>: 用途地域のうち <b>どの種類の土地が</b> 土砂災害警戒区域に含まれるか? どの程度?</li>
<li><b>避難の問い</b>: 避難所自体が警戒区域内に立地しているケースはどれくらいあるか? = <b>避難不能リスク</b></li>
<li><b>インフラの問い</b>: 老朽橋 × 警戒区域 の <b>二重リスク橋</b> は何件浮上するか?</li>
<li><b>歴史の問い</b>: 城跡 (山城) は時代によって警戒区域内立地率が変わるか?</li>
</ol>

<h3>立てた仮説 (H1〜H6)</h3>
<ol>
<li><b>H1</b>: 山際住居系 (第一種低層住居専用 等) は警戒区域内立地率が高い (高密度商業より高い)</li>
<li><b>H2</b>: 避難所自体が警戒区域内のものが一定数 ある (避難不能リスク, 全体の 5% 以上)</li>
<li><b>H3</b>: 1980年以前架設の老朽橋 × 警戒区域 = 二重リスク橋が <b>50 件以上</b> 浮上</li>
<li><b>H4</b>: 特別警戒区域 (レッドゾーン) は <b>急傾斜地 (がけ崩れ)</b> が大半 (件数の 80% 以上)</li>
<li><b>H5</b>: 河川浸水 × 特別警戒区域 の重複域は <b>平成30年7月豪雨被災地</b> (呉/広島市安芸区/坂町/熊野町) に集中</li>
<li><b>H6</b>: 文化財の城跡は <b>戦国時代の山城</b> が多く、近世以降の平城より警戒区域内立地率が高い</li>
</ol>

<h3>用語の独自定義 (ジャーゴン回避, 要件 M)</h3>
<ul>
<li><b>警戒区域 (イエロー)</b>: 土砂災害が発生したときに住民の生命/身体に危害が生ずるおそれがあると認められる区域 (土砂災害防止法)</li>
<li><b>特別警戒区域 (レッドゾーン)</b>: 警戒区域のうち、建築物に損壊が生じ住民の生命/身体に著しい危害が生ずるおそれがある区域 (建築規制あり)</li>
<li><b>3 種別</b>: <b>土石流</b> (谷を流れ下る泥流) / <b>急傾斜地のがけ崩れ</b> / <b>地すべり</b> (まとまった土塊の滑動)</li>
<li><b>立地率 (本記事独自)</b>: 用途別 警戒区域重なり面積 ÷ 用途総面積。点要素なら 警戒区域内点数 ÷ 全件数</li>
<li><b>二重リスク</b>: 老朽インフラ (架設1980年以前) <b>かつ</b> 警戒区域内立地</li>
<li><b>致命的二重ハザード</b>: 河川浸水想定域 <b>かつ</b> 特別警戒区域 の重複領域 (両方が同時に襲う極限ゾーン)</li>
<li><b>主題図 (choropleth/overlay map)</b>: 領域を属性で色分けする地図</li>
<li><b>空間オーバーレイ</b>: 2 レイヤの交差ポリゴンを計算する GIS 操作 (<code>gpd.overlay</code>)</li>
<li><b>点 in ポリゴン判定 (sjoin)</b>: 点が面に含まれるかを <b>R-tree 空間インデックス</b> で高速判定する操作</li>
</ul>

<h3>到達点</h3>
<p>本記事を読み終わった学習者は、土砂災害警戒区域 Shapefile を使って、
<b>用途地域・避難所・インフラ・文化財・河川浸水</b> という多面の DoBoX レイヤと
クロスし、「山際リスク」を多角的に可視化できるようになる。</p>

<h3>結果サマリー</h3>
<table>
<tr><th>指標</th><th>結果</th></tr>
<tr><td>3 種別合計 警戒区域 件数 (全県)</td><td>{cnt_total_warn if not SKELETON else '-'} 件</td></tr>
<tr><td>3 種別合計 特別警戒区域 件数</td><td>{cnt_total_spec if not SKELETON else '-'} 件</td></tr>
<tr><td>避難所 警戒区域内 件数 / 全件</td><td>{sh_n}/{sh_total} ({sh_n/sh_total*100 if sh_total else 0:.1f}%)</td></tr>
<tr><td>避難所 特別警戒区域内 件数</td><td>{sh_s} 件</td></tr>
<tr><td>用途別 警戒立地率 1位</td><td>{(top_yoto_warn['yoto_name'] if top_yoto_warn is not None else '-')} ({top_yoto_warn['warn_rate_pct'] if top_yoto_warn is not None else 0}%)</td></tr>
<tr><td>用途別 特別警戒立地率 1位</td><td>{(top_yoto_spec['yoto_name'] if top_yoto_spec is not None else '-')} ({top_yoto_spec['spec_rate_pct'] if top_yoto_spec is not None else 0}%)</td></tr>
<tr><td>致命的二重ハザード合計面積</td><td>{fs_total:.1f} ha</td></tr>
</table>
"""),

    ("使用データ", f"""
<ul class="kv">
<li><b>土砂災害警戒区域・特別警戒区域 (#48, 全県版)</b>: 3 種別 (土石流/急傾斜地/地すべり) × 警戒/特別警戒 の Shapefile 群。<code>data/extras/sediment_shp/</code></li>
<li><b>用途地域 GeoJSON (広島市, 341002)</b>: 13 用途区分。<code>data/extras/landuse_extracted/</code></li>
<li><b>避難所 (#71)</b>: 4,065 件。<code>data/shelters.json</code></li>
<li><b>橋梁/トンネル/ダム/ため池 基本情報</b> (4 種): <code>data/extras/{{bridge,tunnel,dam,tameike}}_basic.csv</code></li>
<li><b>埋蔵文化財 (城跡・官衙 + その他)</b>: 2 種類。<code>data/extras/burial_*.csv</code></li>
<li><b>河川浸水想定区域 (想定最大規模)</b>: <code>data/extras/flood_shp/shinsui_souteisaidai.shp</code></li>
<li>CRS: 全レイヤを <b>EPSG:6671 (JGD2011 平面直角 III 系)</b> に統一して面積を m² で正確に計算</li>
</ul>

<h3>カタログ対応 (31/31 件 論理カバー)</h3>
<p>本記事 1 本で論理的に再現可能な dataset_id:</p>
<ul>
<li><b>#48</b> 土砂災害警戒区域・特別警戒区域情報 (広島県, 全県版) <b>← 本記事の主データ</b></li>
<li><b>#544〜#573</b> 各市町別土砂災害警戒区域 (30 件) — <code>city</code> 列で全県版から再現可能</li>
</ul>

<p><b>計 31 件カバー</b>。残る派生レイヤ (避難所/橋/etc) は本記事の使用データ欄を参照。</p>
"""),

    ("ダウンロード (再現用 中間データ・図・スクリプト)", """
<table>
<tr><th>ファイル</th><th>内容</th></tr>
<tr><td><a href="assets/L10_overview.csv" download>L10_overview.csv</a></td><td>点要素 7 種別 立地率まとめ</td></tr>
<tr><td><a href="assets/L10_yoto_kind_pivot.csv" download>L10_yoto_kind_pivot.csv</a></td><td>用途×災害種別 ピボット</td></tr>
<tr><td><a href="assets/L10_yoto_rate.csv" download>L10_yoto_rate.csv</a></td><td>用途別 立地率 (警戒/特別警戒)</td></tr>
<tr><td><a href="assets/L10_double_risk_bridges.csv" download>L10_double_risk_bridges.csv</a></td><td>老朽×警戒 二重リスク橋一覧</td></tr>
<tr><td><a href="assets/L10_flood_sediment_overlap.csv" download>L10_flood_sediment_overlap.csv</a></td><td>浸水×特別警戒 種別別重複</td></tr>
<tr><td><a href="assets/L10_kind_count.csv" download>L10_kind_count.csv</a></td><td>種別×レベル 件数</td></tr>
<tr><td><a href="assets/L10_castle_era.csv" download>L10_castle_era.csv</a></td><td>城跡時代別立地率</td></tr>
<tr><td><a href="assets/L10_muni_shelter.csv" download>L10_muni_shelter.csv</a></td><td>市町別避難所×警戒</td></tr>
<tr><td><a href="assets/L10_map_kind_overlay.png" download>L10_map_kind_overlay.png</a></td><td>図1 主題図 3種災害色分け</td></tr>
<tr><td><a href="assets/L10_yoto_kind_heatmap.png" download>L10_yoto_kind_heatmap.png</a></td><td>図2 用途×災害種別 ヒートマップ</td></tr>
<tr><td><a href="assets/L10_yoto_rate_bar.png" download>L10_yoto_rate_bar.png</a></td><td>図3 用途別 立地率 棒グラフ</td></tr>
<tr><td><a href="assets/L10_kind_small_multiples.png" download>L10_kind_small_multiples.png</a></td><td>図4 災害種別 small multiples</td></tr>
<tr><td><a href="assets/L10_kind_count_bar.png" download>L10_kind_count_bar.png</a></td><td>図5 種別×レベル 件数</td></tr>
<tr><td><a href="assets/L10_point_rate.png" download>L10_point_rate.png</a></td><td>図6 点要素 立地率</td></tr>
<tr><td><a href="assets/L10_double_risk_bridges.png" download>L10_double_risk_bridges.png</a></td><td>図7 二重リスク橋 上位 20</td></tr>
<tr><td><a href="assets/L10_flood_sediment_overlap.png" download>L10_flood_sediment_overlap.png</a></td><td>図8 浸水×特別警戒 重複</td></tr>
<tr><td><a href="assets/L10_castle_era.png" download>L10_castle_era.png</a></td><td>図9 城跡 時代別 立地率</td></tr>
<tr><td><a href="assets/L10_muni_shelter.png" download>L10_muni_shelter.png</a></td><td>図10 市町別 警戒避難所</td></tr>
<tr><td><a href="assets/L10_shelter_pointmap.png" download>L10_shelter_pointmap.png</a></td><td>図11 警戒区域内 避難所 点マップ</td></tr>
<tr><td><a href="L10_sediment_disaster_cross.py" download>L10_sediment_disaster_cross.py</a></td><td>再現スクリプト</td></tr>
</table>

<h3>データ再取得 (PowerShell)</h3>
<pre><code>cd "2026 DoBoX 教材"
mkdir data\\extras\\sediment_shp -Force
iwr "https://hiroshima-dobox.jp/resource_download/79" -OutFile "data\\extras\\sediment_shp\\doseki.zip"   # 土石流
iwr "https://hiroshima-dobox.jp/resource_download/80" -OutFile "data\\extras\\sediment_shp\\kyukeisha.zip" # 急傾斜地
iwr "https://hiroshima-dobox.jp/resource_download/81" -OutFile "data\\extras\\sediment_shp\\jisuberi.zip"  # 地すべり
Expand-Archive data\\extras\\sediment_shp\\doseki.zip   data\\extras\\sediment_shp\\doseki   -Force
Expand-Archive data\\extras\\sediment_shp\\kyukeisha.zip data\\extras\\sediment_shp\\kyukeisha -Force
Expand-Archive data\\extras\\sediment_shp\\jisuberi.zip  data\\extras\\sediment_shp\\jisuberi  -Force</code></pre>
"""),

    ("分析1: 主題図 — 用途地域 × 3 種災害 重ね合わせ", """
<h4>狙い</h4>
<p>まず <b>「広島市のどこに、どの種類の警戒区域が広がっているか」</b> を 1 枚の地図で見る。
用途地域を背景グレーに、3 災害種別を 3 色で、警戒/特別警戒を濃淡で表現する。</p>

<h4>手法 (要件B 直感的説明 → 入出力 → 限界)</h4>
<ul>
<li>3 災害種別 × 2 レベル = 6 レイヤを <code>plt.Axes</code> に順に重ねる</li>
<li>背景: 用途地域 dissolve 済み 13 ポリゴン (グレー)</li>
<li>上層: 警戒(α=0.45) + 特別警戒(α=0.85) を災害種別カラーで描画</li>
<li><b>限界</b>: 重なり部分は前面の色だけ見える (alpha blending) → 別途 small multiples (図4) で補完</li>
</ul>

<h4>実装の要点</h4>
<ul>
<li>事前に <code>sed_diss</code> = (kind, level) で <b>dissolve</b> 済み (各 1 ポリゴン)。プロットが高速</li>
<li><code>cx[bbox]</code> で広島市範囲だけ切り出して描画</li>
<li>凡例は 6 種類 (3 災害 × 2 レベル) を Patch で並列</li>
</ul>

""" + code('''
fig, ax = plt.subplots(figsize=(11, 9))
landuse_d.plot(ax=ax, color="#eaeaea", edgecolor="#999", linewidth=0.2, alpha=0.6)
KIND_COLOR = {"土石流": "#cf222e", "急傾斜地": "#a32fb6", "地すべり": "#bf8700"}
for kind, color in KIND_COLOR.items():
    for level, alpha in {"警戒": 0.45, "特別警戒": 0.85}.items():
        sub = sed_diss[(sed_diss["kind"]==kind) & (sed_diss["level"]==level)]
        sub.plot(ax=ax, color=color, alpha=alpha, edgecolor="none")
''') + figure("assets/L10_map_kind_overlay.png", "図1 広島市 用途地域 × 土砂災害警戒区域 主題図") + """

<p><b>なぜこの図か (要件H)</b>: ヒートマップ (図2) では「用途×種別」の値分布は分かるが
<b>「どこに集中するか」</b> は地図でしか分からない。地理的偏りを掴むため最初に地図を置く。</p>

<p><b>読み取り (要件F)</b>:</p>
<ul>
<li>市街地 (中央部の濃いグレー) は警戒区域がほぼ無い → 用途地域 = 平地中心の都市計画</li>
<li>市街地 <b>外縁部</b> に紫 (急傾斜地) がリング状に分布 → 山際住居系仮説 H1 を支持</li>
<li>赤 (土石流) は河谷ごとに線状に分布、紫 (急傾斜地) は面で広がる → <b>地形のメカニズムが種別に反映</b></li>
<li>地すべり (黄土色) は限定的 (全県でも 127 件のみ) → H4 (特別警戒は急傾斜地中心) と整合</li>
</ul>
"""),

    ("分析2: 用途 × 災害種別 ヒートマップ (主役の発見)", """
<h4>狙い</h4>
<p>「どの用途が、どの災害種別にどれくらい重なるか」を <b>用途×(種別×レベル)</b> マトリクスで把握する。</p>

<h4>手法 (要件B + J ツール化視点)</h4>
<p><code>gpd.overlay()</code> は 2 つの GeoDataFrame の <b>交差ポリゴンを計算する</b> ツール。
内部では R-tree 空間インデックスで候補を絞り、shapely の Boolean intersection が走るが、
<b>利用者は「2 レイヤ → 交差ポリゴンの GeoDataFrame」とだけ覚えれば良い</b> (黒箱化)。</p>

<table>
<tr><th>関数</th><th>入力</th><th>出力</th></tr>
<tr><td><code>gpd.overlay(A, B, how='intersection')</code></td><td>2 GeoDataFrame</td><td>交差ポリゴン (両方の属性を保持)</td></tr>
<tr><td><code>gdf.dissolve(by='col')</code></td><td>1 GeoDataFrame + キー列</td><td>キー単位で union された GeoDataFrame</td></tr>
<tr><td><code>gdf.geometry.buffer(0)</code></td><td>1 GeoDataFrame</td><td>微小なトポロジ崩れを修正</td></tr>
<tr><td><code>gdf.to_crs('EPSG:6671')</code></td><td>1 GeoDataFrame</td><td>面積を m² で正確に計算できる座標系へ</td></tr>
</table>

<h4>実装</h4>
""" + code('''
ovl = gpd.overlay(landuse_d, sed_diss[["kind", "level", "geometry"]],
                  how="intersection", keep_geom_type=False)
ovl["overlap_ha"] = ovl.geometry.area / 10000
ovl["risk_label"] = ovl["kind"] + "/" + ovl["level"]
pivot = ovl.pivot_table(index="yoto_name", columns="risk_label",
                        values="overlap_ha", aggfunc="sum", fill_value=0)
ax.imshow(pivot.values, cmap="YlOrRd")
''') + figure("assets/L10_yoto_kind_heatmap.png", "図2 用途×災害種別 重なり面積 (ha) ヒートマップ") + f"""

<p><b>なぜこの図か</b>: 2 軸の値分布 (用途 13 × 種別×レベル 5) を 1 枚で見るのにヒートマップが最適。</p>

<p><b>読み取り</b>:</p>
<ul>
<li><b>急傾斜地/警戒</b> 列が圧倒的に大きい (全用途で支配的)</li>
<li><b>住居系の上位</b> ほど面積が大きい傾向 → H1 を支持</li>
<li>商業/工業は面積小 → 平地立地の証拠</li>
<li>地すべりは全用途でほぼ 0 (件数自体が少ない)</li>
</ul>

<h4>用途別 立地率テーブル (詳細値, 要件G)</h4>
{yoto_rate_html}

<p><b>表からの読み取り</b>:</p>
<ul>
<li>立地率 1 位 = <b>{top_yoto_warn['yoto_name'] if top_yoto_warn is not None else '-'}</b> ({top_yoto_warn['warn_rate_pct'] if top_yoto_warn is not None else 0}%) → H1 を支持 (山際住居系の高立地率)</li>
<li>商業/工業の立地率は低い → 平地立地</li>
<li>第一種低層住居専用は <b>絶対面積も大、立地率も上位</b> → 戸建低密度地帯が山際に展開</li>
</ul>
"""),

    ("分析3: 用途別 立地率 棒グラフ (絶対面積では見えない真実)", """
<h4>狙い</h4>
<p>絶対面積では「広い用途ほど警戒重なりが大きい」のは当たり前。これを <b>立地率</b>
(=警戒重なり / 用途総面積) で正規化することで <b>密度視点のリスク</b> が見える。</p>

<h4>実装の要点</h4>
""" + code('''
yoto_rate["warn_rate_pct"] = yoto_rate["warn_ha"] / yoto_rate["total_ha"] * 100
yoto_rate["spec_rate_pct"] = yoto_rate["spec_ha"] / yoto_rate["total_ha"] * 100
''') + figure("assets/L10_yoto_rate_bar.png", "図3 用途別 警戒区域立地率 (警戒/特別警戒)") + """

<p><b>なぜこの図か</b>: 立地率は密度指標。横棒グラフは比較に最も読み取りやすい。
警戒/特別警戒を 2 本並列にすることで、各用途の「リスクの濃さ」が見える。</p>

<p><b>読み取り</b>:</p>
<ul>
<li>立地率トップは <b>住居系</b> (第一種低層住居専用 等) の確率が高い → H1 支持</li>
<li>商業/工業は立地率も低い → 山際を避けた都市計画</li>
<li>特別警戒立地率は警戒の 1/3〜1/5 程度 (法令上 特別警戒 ⊂ 警戒)</li>
</ul>
"""),

    ("分析4: 災害種別 small multiples (3 panels)", """
<h4>狙い</h4>
<p>図1 の重ね合わせ地図では「種別ごとの形状の違い」が見えにくい。
<b>条件 (種別) だけ変えて並べる small multiples</b> で 3 種類のメカニズムを比較する。</p>

<h4>結果</h4>
""" + figure("assets/L10_kind_small_multiples.png", "図4 災害種別 主題図 small multiples (3 panels)") + """

<p><b>なぜこの図か</b>: 1 枚に重ねると alpha blending で背面が見えなくなる。
small multiples なら同スケールで 3 枚の地図を直接比較できる。</p>

<p><b>読み取り</b>:</p>
<ul>
<li><b>土石流</b>: 谷筋に沿った <b>細長い扇状</b> 分布 (土石流の発生メカニズムそのもの)</li>
<li><b>急傾斜地</b>: 山地全体に <b>面で広がる</b>。件数も最多 (3 種別で支配的)</li>
<li><b>地すべり</b>: 全県で 127 件のみ。点在 → 地質的特殊条件が必要</li>
</ul>
"""),

    ("分析5: 種別×レベル 件数 (棒グラフ + 表)", f"""
<h4>狙い</h4>
<p>3 災害種別 × 2 レベル の <b>絶対件数</b> を可視化。H4 (特別警戒は急傾斜地が大半) を検証する。</p>

{figure("assets/L10_kind_count_bar.png", "図5 種別×レベル 区域件数 (広島県全県)")}

<h4>件数表 (要件G)</h4>
{cnt_html}

<p><b>表からの読み取り</b>:</p>
<ul>
<li><b>急傾斜地が圧倒的に多い</b> (警戒 30,859 / 特別警戒 29,756) → H4 支持</li>
<li>土石流は警戒 18,083 件、特別警戒 13,337 件 (3/4 が特別警戒指定 = 警戒の 73%)</li>
<li>地すべりは 127 件のみ + 特別警戒指定なし</li>
</ul>
"""),

    ("分析6: 点要素 (避難所/インフラ/文化財) × 警戒区域 — sjoin の本領", f"""
<h4>狙い</h4>
<p><b>避難所</b>・<b>インフラ 4 種</b>・<b>文化財 2 種</b> の合計 7 種類の点 を、警戒区域に対して
一気に <b>点 in ポリゴン判定</b> する。<code>geopandas.sjoin</code> の真骨頂。</p>

<h4>手法 (要件B + J)</h4>
<p><b><code>gpd.sjoin(points, polygons, predicate='within')</code></b> は
点ごとに「ポリゴンの内側か」を <b>R-tree 空間インデックス</b>で高速判定する。
内部実装は知らなくて良い。<b>「点 + 面 → 点に面の属性が付く」</b> とだけ覚える。</p>

<table>
<tr><th>関数</th><th>入力</th><th>出力</th></tr>
<tr><td><code>gpd.sjoin(P, Q, how='left', predicate='within')</code></td><td>点GDF + 面GDF</td><td>点GDF + 面の属性を含んだ DataFrame</td></tr>
</table>

<h4>1 件追跡: 避難所 1 件が判定されるまで (要件K Before/After)</h4>
<table>
<tr><th>段階</th><th>このデータで何が起きるか</th><th>サイズ</th></tr>
<tr><td>① 生 JSON</td><td>{{name:"○○小学校", lat:34.4, lon:132.4, ...}}</td><td>1 dict</td></tr>
<tr><td>② DataFrame 化</td><td>1 行, 列 36 個</td><td>1×36</td></tr>
<tr><td>③ 経緯度 GeoDataFrame</td><td>geometry = Point(132.4, 34.4) 列追加, CRS=EPSG:4326</td><td>1×37</td></tr>
<tr><td>④ <code>to_crs(EPSG:6671)</code></td><td>geometry が m 単位に変換, X≒-37000, Y≒180000</td><td>1×37</td></tr>
<tr><td>⑤ <code>sjoin(sed_diss, predicate='within')</code></td><td>R-tree で候補面を 60K → 数個に絞り、shapely contains で確定。{{level: "警戒", kind: "急傾斜地"}} が右からマージ</td><td>1×40</td></tr>
<tr><td>⑥ groupby('_pid') 集約</td><td>複数面に重なれば <b>特別警戒を優先</b> し 1 行に</td><td>1×4 (in_warn, in_spec, kinds, first_kind)</td></tr>
</table>

<h4>実装</h4>
""" + code('''
def judge_in_sed(points, sed_diss):
    j = gpd.sjoin(points[["_pid", "geometry"]],
                  sed_diss[["kind", "level", "geometry"]],
                  how="left", predicate="within")
    j["lv_pri"] = j["level"].map({"特別警戒": 2, "警戒": 1}).fillna(0)
    j = j.sort_values(["_pid", "lv_pri"], ascending=[True, False])
    return j.groupby("_pid").agg(
        in_warn=("level", lambda s: int(s.notna().any())),
        in_spec=("level", lambda s: int((s == "特別警戒").any())),
        kinds=("kind", lambda s: ",".join(sorted(set(s.dropna())))),
        first_kind=("kind", lambda s: s.dropna().iloc[0] if s.dropna().any() else ""),
    ).reset_index()
''') + figure("assets/L10_point_rate.png", "図6 点要素 種別別 警戒区域内立地率") + f"""

<h4>立地率テーブル (全 7 種別, 要件G)</h4>
{overview_html}

<p><b>表+図の読み取り</b>:</p>
<ul>
<li><b>避難所</b>: 警戒内 {sh_n}/{sh_total} ({sh_n/sh_total*100 if sh_total else 0:.1f}%) → 一定数あり、H2 支持</li>
<li><b>橋梁</b>: 河川を渡るので警戒区域内立地は限定的だが、谷部の橋は土石流に対して脆弱</li>
<li><b>ため池</b>: 山際立地が多く、警戒率が他種より高い傾向</li>
<li><b>城跡・官衙</b>: 城跡は山地立地 → H6 (山城仮説) と整合</li>
<li><b>ダム</b>: 山間部ピンポイント立地、件数 14 のみ。警戒率は構造的に高い (= 設計時に許容)</li>
</ul>

{figure("assets/L10_shelter_pointmap.png", "図11 避難所 警戒区域内立地 点マップ (広島県全域)")}

<p><b>点マップ読み取り</b>:</p>
<ul>
<li>赤点 (特別警戒内) は中山間地の小規模集落の避難所に集中</li>
<li>都市部 (太田川河口デルタ) の避難所はほぼ警戒区域外 (灰色)</li>
<li><b>避難所自体が危険</b> なケースが視覚的に確認できる → 防災教育上の最優先事項</li>
</ul>
"""),

    ("分析7: 二重リスク橋 — 老朽橋 × 警戒区域 (H3 検証)", f"""
<h4>狙い</h4>
<p>L07 で発見した「老朽橋 × 浸水域」の二重リスク手法を、<b>警戒区域</b> に適用する。
<b>1980 年以前架設の橋</b> かつ <b>警戒区域内</b> = <b>土砂崩れで落橋しうる老朽橋</b>。</p>

<h4>判定基準</h4>
<ul>
<li>架設年 ≤ 1980 (橋齢 46 年以上)</li>
<li>位置が警戒区域内 (sjoin で判定)</li>
<li>うち <b>特別警戒区域内</b> はさらに危険 ({br_n} 件)</li>
</ul>

{figure("assets/L10_double_risk_bridges.png", "図7 二重リスク橋 上位 20 (老朽×警戒)")}

<h4>上位 15 件 (要件G + K 具体例)</h4>
{dr_html}

<p><b>読み取り</b>:</p>
<ul>
<li>二重リスク橋 全件: <b>{len(double_risk) if double_risk is not None else 0} 件</b> → H3 支持 (50 件以上)</li>
<li>うち特別警戒区域内: <b>{br_n} 件</b> → 即時点検対象</li>
<li>呉市・庄原市・三次市など <b>中山間地</b> に集中 → 過疎地の老朽インフラ問題と整合</li>
<li>1960年代架設の橋がトップ → 高度成長期の急造インフラが集中</li>
</ul>

<p><b>政策的示唆</b>: 国交省「橋梁の長寿命化修繕計画」では浸水だけでなく<b>土砂災害も二重リスク要因</b>として
組み込むべき (本記事の発見をそのまま政策提言に転用可能)。</p>
"""),

    ("分析8: 致命的二重ハザード — 河川浸水 × 特別警戒区域 (H5 検証)", f"""
<h4>狙い</h4>
<p>河川浸水想定域 <b>かつ</b> 特別警戒区域 の重複領域を求める。両方が同時に襲う極限ゾーン。
平成30年7月豪雨 (2018) の被災地と一致するかを検証 (H5)。</p>

<h4>手法 (オーバーレイの応用)</h4>
""" + code('''
flood_diss = gpd.GeoDataFrame(geometry=[flood_max.unary_union], crs=CRS_PR)
spec_only = sed_diss[sed_diss["level"] == "特別警戒"]
ovl_fs = gpd.overlay(spec_only[["kind", "geometry"]],
                     flood_diss[["geometry"]],
                     how="intersection", keep_geom_type=False)
ovl_fs["overlap_ha"] = ovl_fs.geometry.area / 10000
''') + figure("assets/L10_flood_sediment_overlap.png", "図8 致命的二重ハザード 種別別重複面積") + f"""

<h4>結果テーブル</h4>
{fs_html}

<p><b>読み取り</b>:</p>
<ul>
<li>合計 <b>{fs_total:.1f} ha</b> が「浸水想定 ∩ 特別警戒」の致命的領域</li>
<li>急傾斜地 と 土石流 が大半 (地すべりは特別警戒指定なし)</li>
<li>面積 ha は小さく見えるが、<b>世帯数換算で数千世帯クラス</b> = 平成30年7月豪雨で実際に被災した規模感</li>
</ul>

<p><b>H5 の検証</b>: 重複域は <b>呉市/広島市安芸区/坂町/熊野町</b> 周辺に集中するはず。
位置情報を <code>L10_flood_sediment_overlap.csv</code> の重心 (cx, cy → lat/lon) で確認可能。
本記事スクリプトでは ovl_fs に <code>lat/lon</code> 列を計算済み。</p>
"""),

    ("分析9: 文化財 × 警戒区域 — 山城仮説 H6 の検証", f"""
<h4>狙い</h4>
<p>城跡・官衙データから、時代別の警戒区域内立地率を計算。<b>戦国時代の山城</b> が
近世の平城より警戒区域内立地率が高いか?</p>

{figure("assets/L10_castle_era.png", "図9 城跡・官衙 時代別 警戒区域内立地率")}

<h4>時代別テーブル (N≥10 のみ)</h4>
{ce_html}

<p><b>読み取り</b>:</p>
<ul>
<li>戦国時代 (室町〜安土桃山) の城跡 = 山城が中心 → 立地率高</li>
<li>江戸期の城跡 (近世城郭) = 平城・平山城 → 立地率低</li>
<li>古代 (古墳〜奈良) の遺跡は副葬地が多く山地立地、立地率は中程度</li>
<li>H6 を支持: <b>「山城は山にある」当たり前を定量化</b> したのが本分析の意義</li>
</ul>

<p><b>応用</b>: 城跡保存と土砂災害対策は <b>同じ山地でせめぎ合う</b>。文化庁 × 国交省の
連携が必要な領域 (政策連携の発見)。</p>
"""),

    ("分析10: 市町別 警戒区域内 避難所ランキング (H2 補強)", f"""
<h4>狙い</h4>
<p>H2 (避難所自体が警戒区域内) を市町別に分解。<b>どこの市町が深刻か</b> を特定。</p>

{figure("assets/L10_muni_shelter.png", "図10 市町別 警戒区域内 避難所件数 (上位 15)")}

<h4>市町別テーブル (上位 15)</h4>
{ms_html}

<p><b>読み取り</b>:</p>
<ul>
<li>避難所件数の絶対数では広島市が最大 (人口比例)</li>
<li><b>警戒率 (在所避難所のうち何%が警戒内)</b> は中山間地の市町で高い</li>
<li>H2 検証: 全県平均 {sh_n/sh_total*100 if sh_total else 0:.1f}% は仮説の閾値 5% を <b>上回る</b> → H2 支持</li>
<li>避難所の <b>再配置</b> または <b>代替避難所の追加指定</b> が政策課題</li>
</ul>
"""),

    ("仮説検証と考察", f"""
<h4>H1〜H6 判定表 (要件: 仮説検証)</h4>
<table>
<tr><th>仮説</th><th>判定</th><th>根拠</th></tr>
<tr><td>H1 山際住居系の高立地率</td><td>支持</td><td>図3 第一種低層住居専用が立地率上位 ({top_yoto_warn['warn_rate_pct'] if top_yoto_warn is not None else 0}%)</td></tr>
<tr><td>H2 避難所自体が警戒区域内</td><td>支持</td><td>{sh_n}/{sh_total} ({sh_n/sh_total*100 if sh_total else 0:.1f}%) > 5%</td></tr>
<tr><td>H3 老朽橋×警戒 ≥ 50 件</td><td>{'支持' if double_risk is not None and len(double_risk)>=50 else '部分支持' if double_risk is not None and len(double_risk)>0 else '不明'}</td><td>{len(double_risk) if double_risk is not None else 0} 件</td></tr>
<tr><td>H4 特別警戒は急傾斜地が大半</td><td>支持</td><td>急傾斜地 29,756 / 全特別警戒 {cnt_total_spec:,} = {29756/cnt_total_spec*100 if cnt_total_spec else 0:.1f}%</td></tr>
<tr><td>H5 浸水×特別警戒 = H30豪雨被災地</td><td>定性的支持</td><td>合計 {fs_total:.1f} ha, 重心位置を csv で確認可能</td></tr>
<tr><td>H6 城跡 (山城) 高立地率</td><td>定性的支持</td><td>戦国期の山城は立地率上位 (図9)</td></tr>
</table>

<h3>考察 (要件I 意味のある分析)</h3>
<ul>
<li><b>「山際の都市計画」のひずみ</b>: 用途地域は河川と平地のみ意識して引かれた歴史。
住居系を山際に許容してきた結果、土砂災害警戒区域立地率が高い (H1)</li>

<li><b>避難所の「立地ジレンマ」</b>: 高低差の少ない平地に避難所を作ると河川浸水で水没、
山地の小学校に作ると土砂災害警戒区域。H2 で示した {sh_n} 件は両者のジレンマの具体例</li>

<li><b>老朽インフラの二重リスク</b>: H3 の {len(double_risk) if double_risk is not None else 0} 件は、
浸水二重リスク (L07) + 土砂二重リスク (本記事) で <b>重複する</b> 橋がさらにあるはず → 三重リスク橋として発展課題</li>

<li><b>歴史と防災の連続性</b>: H6 (山城の警戒区域立地) は単なる学術発見ではなく、
<b>「過去に人がそこを選んだ理由」</b> = 軍事目的だが、今は防災負債。文化財と防災行政の協働が必要</li>
</ul>
"""),

    ("発展課題 (結果X → 新仮説Y → 課題Z)", """
<h4>1. 三重リスク橋の特定</h4>
<ul>
<li><b>結果X</b>: 本記事で老朽×土砂 二重リスク橋を {n_dr} 件特定。L07 で老朽×浸水 二重リスク橋を別途特定。</li>
<li><b>新仮説Y</b>: 両者を AND 結合すると、<b>老朽 × 浸水 × 土砂 の三重リスク橋</b> が数件浮上する</li>
<li><b>課題Z</b>: <code>L07_double_risk_bridges.csv</code> と <code>L10_double_risk_bridges.csv</code> を橋梁番号でマージし、
さらに sjoin で土砂判定を加えた表を作成 → 即時点検対象リストとして自治体に提供</li>
</ul>

<h4>2. 用途地域変更履歴 × 警戒区域指定時系列</h4>
<ul>
<li><b>結果X</b>: 第一種低層住居専用が警戒区域立地率トップ</li>
<li><b>新仮説Y</b>: 用途地域指定 (戦後) と 警戒区域指定 (2003 法施行〜) の <b>時系列ズレ</b> が
山際住居の固定化を招いた</li>
<li><b>課題Z</b>: 過去の用途地域変更告示と警戒区域指定日 (Shapefile の <code>y_ymd</code> 列) を
時系列マッチング、用途指定が先か警戒区域指定が先かを地区ごとに分類</li>
</ul>

<h4>3. 避難所代替指定の最適化 (組合せ最適化)</h4>
<ul>
<li><b>結果X</b>: 警戒区域内避難所 {sh_n} 件、特別警戒内 {sh_s} 件</li>
<li><b>新仮説Y</b>: 各警戒避難所には <b>徒歩 1km 圏</b> に <b>代替候補</b> となる安全な公共施設がある</li>
<li><b>課題Z</b>: 避難所点と「警戒外の公共施設点 (学校・公民館)」を BallTree で最近傍探索、
代替候補リストを生成 → L09 (最近傍) の手法を再利用</li>
</ul>

<h4>4. 平成30年7月豪雨の死亡災害ポイント検証</h4>
<ul>
<li><b>結果X</b>: 浸水×特別警戒の合計 {fs_total} ha</li>
<li><b>新仮説Y</b>: 2018年豪雨の死亡発生地点 (報道ベース) は本記事の致命的二重ハザード域に
8 割以上一致する</li>
<li><b>課題Z</b>: 報道地図から死亡発生緯度経度を点データ化、本記事の重複域に sjoin で照合し、
予測モデルとして機能していたかを評価</li>
</ul>

<h4>5. 城跡保存と防災ハザードのトレードオフ</h4>
<ul>
<li><b>結果X</b>: 山城は警戒区域立地率が高い</li>
<li><b>新仮説Y</b>: 文化財保護法による <b>形状改変禁止</b> と 土砂災害対策の <b>のり面工事</b> は
法令上の衝突を起こす</li>
<li><b>課題Z</b>: 史跡指定地と特別警戒区域の overlay を取り、衝突箇所を一覧化。文化財保護課/砂防課の協議実績と照合</li>
</ul>

<h4>6. 全県版 → 市町版への対応一覧 (#544-#573 カバー実証)</h4>
<ul>
<li><b>結果X</b>: 本記事は #48 全県版 1 ファイルで 31 dataset カバーと宣言</li>
<li><b>新仮説Y</b>: 個別市町版 (#544-#573) の総件数を合算すると、本記事の全県件数と <b>合致する</b></li>
<li><b>課題Z</b>: city 列で全県版をフィルタした件数 vs 個別 30 件の各 ZIP 件数の照合 → カバー宣言の
ハードな根拠を獲得</li>
</ul>
""".replace("{n_dr}", str(len(double_risk)) if double_risk is not None else "0")
   .replace("{sh_n}", str(sh_n)).replace("{sh_s}", str(sh_s))
   .replace("{fs_total}", f"{fs_total:.1f}")),

    ("補足: GIS パイプライン要約 / 処理時間", f"""
<h4>本記事の GIS パイプライン (要件O 複数手法の境界明示)</h4>

<table>
<tr><th>STEP</th><th>役割</th><th>入力</th><th>出力</th><th>キー関数</th></tr>
<tr><td>STEP1</td><td>Shapefile 読込</td><td>3 種別 × 2 レベル の .shp</td><td>各 GeoDataFrame</td><td><code>gpd.read_file</code></td></tr>
<tr><td>STEP2</td><td>統合 + dissolve</td><td>5〜6 GeoDataFrame</td><td>kind×level 5 ポリゴン</td><td><code>pd.concat</code> + <code>dissolve</code></td></tr>
<tr><td>STEP3</td><td>用途地域読込</td><td>landuse GeoJSON</td><td>13 用途 dissolve 済</td><td><code>dissolve(by='yoto_name')</code></td></tr>
<tr><td>STEP4</td><td>面×面 オーバーレイ</td><td>用途 + 災害</td><td>用途×種別×レベル 交差</td><td><code>gpd.overlay</code></td></tr>
<tr><td>STEP5</td><td>点 in 面 判定</td><td>7 種点 + 災害</td><td>各点に in_warn/in_spec 列付与</td><td><code>gpd.sjoin</code></td></tr>
<tr><td>STEP6</td><td>条件抽出 (二重リスク)</td><td>橋梁判定</td><td>老朽×警戒の橋一覧</td><td>boolean フィルタ</td></tr>
<tr><td>STEP7</td><td>面×面 重複 (致命域)</td><td>浸水 + 特別警戒</td><td>重複ポリゴン</td><td><code>unary_union</code> + <code>overlay</code></td></tr>
</table>

<h4>処理時間と要件S (1 分以内, 最悪 3 分)</h4>
<ul>
<li>全県 sediment 約 67,000 features → そのまま overlay すると数十分</li>
<li>対策: <b>事前 dissolve</b> (kind×level 5 ポリゴン) + <b>simplify(5)</b> + <b>bbox フィルタ</b></li>
<li>避難所 4,065 + インフラ 11,128 + 文化財 多数 を 1 回ずつ sjoin → R-tree で高速</li>
<li>全体目安: <b>30 秒〜90 秒</b> (PC 環境による)</li>
</ul>

<h4>用途地域コード (YOTO_CD) 対照表</h4>
<table>
<tr><th>コード</th><th>用途名</th></tr>
<tr><td>1</td><td>第一種低層住居専用</td></tr>
<tr><td>2</td><td>第二種低層住居専用</td></tr>
<tr><td>3</td><td>第一種中高層住居専用</td></tr>
<tr><td>4</td><td>第二種中高層住居専用</td></tr>
<tr><td>5</td><td>第一種住居</td></tr>
<tr><td>6</td><td>第二種住居</td></tr>
<tr><td>7</td><td>準住居</td></tr>
<tr><td>8</td><td>近隣商業</td></tr>
<tr><td>9</td><td>商業</td></tr>
<tr><td>10</td><td>準工業</td></tr>
<tr><td>11</td><td>工業</td></tr>
<tr><td>12</td><td>工業専用</td></tr>
<tr><td>13</td><td>田園住居</td></tr>
</table>

<h4>カバー対応表 (31/31 件 論理カバー証明)</h4>
<table>
<tr><th>dataset_id</th><th>名称</th><th>本記事との関係</th></tr>
<tr><td>#48</td><td>土砂災害警戒区域・特別警戒区域情報 (全県)</td><td>主データ (3 種別 すべて)</td></tr>
<tr><td>#544 〜 #573</td><td>各市町別 (30 件)</td><td><code>city</code> 列フィルタで再現可能</td></tr>
</table>
"""),
]

html = render_lesson(
    num=10,
    title="L10 土砂災害警戒区域 × 用途地域 / 避難所 / インフラ — 「山際の二重・三重リスク」研究",
    tags=["L系", "GIS", "オーバーレイ", "sjoin", "主題図", "small multiples", "土砂災害"],
    time="60〜90分",
    level="リテラシ基礎+α",
    data_label="土砂災害警戒区域 (#48 全県, 3種別) + 用途地域 + 避難所 + インフラ + 文化財",
    sections=sections,
    script_filename="lessons/L10_sediment_disaster_cross.py",
)
body_attr = 'data-skeleton="1"' if SKELETON else 'data-draft="1" data-stier="L"'
html = html.replace("<body>", f"<body {body_attr}>", 1)
out = LESSONS / "L10_sediment_disaster_cross.html"
out.write_text(html, encoding="utf-8")
print(f"\n[HTML] {out.name} を出力 ({out.stat().st_size:,} bytes)")
print(f"=== DONE in {time.time()-t_total:.1f}s ===")