# -*- coding: utf-8 -*-
"""L31 観測情報 9 件統合分析
       — 県の防災観測網 (rain/water/dam/tide/wind) の幾何構造の解読

カバー宣言:
  本記事は DoBoX のシリーズ <b>「観測情報_*」 9 件</b>
  (dataset_id = 1274, 1275, 1276, 1437, 1438, 1439, 1440, 1441, 1442) を統合し、
  広島県内における<b>防災観測網</b>の<b>幾何構造</b>を分析する研究記事である。

  9 件は <b>9 市町分の同一データ</b> ではなく、
  <b>1 件のマスタ + 5 観測種類 × {日集計 / 年集計} の時系列ファイル</b>
  からなる<b>機能分化型シリーズ</b>である:

    1274 観測情報_観測所一覧            (master, 緯度経度+属性, 634 観測所)
    1275 観測情報_雨量日集計            (rain    × daily,  10 分粒度)
    1276 観測情報_雨量年集計            (rain    × yearly, 1 日粒度)
    1437 観測情報_水位日集計            (water   × daily)
    1438 観測情報_水位年集計            (water   × yearly)
    1439 観測情報_ダム日集計            (dam     × daily)
    1440 観測情報_ダム年集計            (dam     × yearly)
    1441 観測情報_潮位日集計            (tide    × daily)
    1442 観測情報_風向風速日集計        (wind    × daily)

  すなわちこのシリーズは、<b>「観測網の地理構造を 1 枚のマスタに集約し、
  各観測種類は別 dataset で時系列を提供する」</b>という<b>2 層構成</b>に整備されている。

  本記事は<b>マスタ (1 件) を主軸</b>とし、各時系列 dataset は
  <b>「観測種類×時間粒度」マトリクスの構造</b>として併置統合する。
  時系列の数値解析自体は X06 が雨量で先行実施済みのため、
  本記事は X06 と<b>意図的に重複しない</b>「観測網そのものの構造」を扱う。

研究の問い (RQ):
  広島県内の防災観測網は、観測種類別に何拠点配置され、地理的にどの流域・
  どの地域をカバーしているか? その配置パターンから、県の災害監視思想
  (どの種の災害をどこで先行検知したいか) はどう読み取れるか?
    (1) 5 観測種類 (雨量/水位/ダム/潮位/風) は<b>何拠点ずつ、どの所管</b>に
        配置され、地理的にどう分布するか?
    (2) 「観測点間距離」を見ると、観測網はどの程度の<b>空間粒度</b>で県土を
        カバーしているか? 観測<b>空白地帯</b>は存在するか?
    (3) 雨量と水位の<b>水系階層構造</b>: 1 水系あたり雨量何点・水位何点と
        いう密度比は何を意味するか?
    (4) 河川氾濫閾値 (はん濫注意/避難判断/はん濫危険) を保有する観測所
        の数と地理分布は? これは<b>「先行指標」</b>として機能できる場所か?

仮説 H1〜H6:
  H1 (種別不均衡): 5 観測種類は<b>件数が極端に偏る</b>。
       雨量 > 水位 ≫ 風 ≫ ダム ≈ 潮位 の順で、
       <b>雨量だけで全体の 60% 以上</b>を占める。
       理由: 雨量は局所性が強く密配置が必須、潮位は海岸線にしか置けない。
  H2 (所管二元性): 観測網は<b>「砂防課/河川課/河川課ダムG (県)」</b>と
       <b>「国土交通省/国ダム/気象台 (国)」</b>の<b>2 主体</b>がほぼ折半する。
       残りの所管 (港湾漁港整備課・農業基盤課・海上保安庁) は周辺的。
       これは<b>「県管理河川 + 国管理河川」</b>の二層管理構造を反映する。
  H3 (水系密度): 太田川水系・江の川水系・芦田川水系の<b>主要 3 水系</b>に
       観測点の<b>50% 以上</b>が集中する。
       1 水系あたり雨量:水位 比は概ね <b>2:1〜3:1</b>。
       これは<b>「降った場所」と「流れ着いた場所」を両方見る</b>運用思想を示す。
  H4 (海岸線潮位の精選): 潮位 13 点は<b>瀬戸内海沿岸の主要港</b>に集中し、
       <b>内陸には 1 点も存在しない</b>。これは自明だが、
       <b>「潮位観測 = 港湾管理 = 海上保安庁系」</b>の所管論理を裏付ける。
  H5 (氾濫閾値の階層性): 4 段階の氾濫閾値 (待機 < 注意 < 判断 < 危険) を
       保有する観測所は<b>水位 181 点のうち 50%強</b>。残りは<b>監視のみ</b>。
       閾値あり水位は<b>主要河川の中下流部</b>に集中し、<b>上流支流</b>では
       監視のみが多い。これは<b>「警報出力ノード = 中下流」</b>の制度設計。
  H6 (空白地帯): 雨量観測点を中心とする最寄り観測所までの距離分布は、
       平均 < 5 km だが、<b>島嶼部 (江田島/大崎上島/瀬戸内諸島) と
       一部の中山間部</b>に最寄り > 8 km の<b>観測空白</b>が現れる。
       これは<b>「人口希薄地の観測網密度低下」</b>の典型例である。

要件 S 準拠: 1 分以内完走 (時系列は構造把握のみ、数値解析はしない)。
要件 T 準拠: 県全域マップ + 種別 small multiples + 水系別密度マップ
            + 氾濫閾値マップ + 観測空白マップ + 距離分布など多彩な GIS 図。
要件 Q 準拠: 図 11 種、表 9 種以上 (634 点 × 5 種類 × 多角度)。

データ仕様 (主軸 + 時系列補助):
  A. 観測所マスタ 1274:
     CSV, 634 行 × 35 列, 緯度経度・水系・河川・市町村・閾値属性
     (はん濫注意・避難判断・はん濫危険・常時満水位・注意潮位 等)
  B. 時系列 8 件 (1275-1276, 1437-1442):
     CSV, 観測所×時刻の wide 表。観測所 ID 行 + 値行。
     L31 ではこれらの<b>列構造の同型性</b>と<b>カバー観測所の重なり</b>のみ
     扱う。実時系列の数値分析は X06 (雨量) で別途実装済。

メモリ対策: Figure ごとに plt.close('all') で確実に解放。
"""
from __future__ import annotations
import sys, time, json, re
from pathlib import Path

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

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Patch, Circle
from matplotlib.lines import Line2D
import geopandas as gpd
from shapely.geometry import Point, box
from scipy.spatial import cKDTree

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

t0 = time.time()
print("=== L31 観測情報 9 件統合分析 ===", flush=True)

# =============================================================================
# 0. 定数: 9 dataset_id と 5 観測種類の凡例
# =============================================================================
DATA_DIR = ROOT / "data" / "extras" / "L31_observation_network"
DATA_DIR.mkdir(parents=True, exist_ok=True)
ADMIN_DIR = ROOT / "data" / "extras" / "L15_admin_zones"
TARGET_CRS = "EPSG:6671"  # JGD2011 平面直角 III, m 単位

# 9 件: (dsid, 略称, 役割, kind, granularity, file)
SERIES = [
    (1274, "観測所一覧",       "master",   None,   None,    "stations_master.csv"),
    (1275, "雨量_日集計",      "timeseries","rain", "daily",  "rain_daily_sample.csv"),
    (1276, "雨量_年集計",      "timeseries","rain", "yearly", "rain_year_sample.csv"),
    (1437, "水位_日集計",      "timeseries","water","daily",  "water_daily_sample.csv"),
    (1438, "水位_年集計",      "timeseries","water","yearly", "water_year_sample.csv"),
    (1439, "ダム_日集計",      "timeseries","dam",  "daily",  "dam_daily_sample.csv"),
    (1440, "ダム_年集計",      "timeseries","dam",  "yearly", "dam_year_sample.csv"),
    (1441, "潮位_日集計",      "timeseries","tide", "daily",  "tide_daily_sample.csv"),
    (1442, "風向風速_日集計",  "timeseries","wind", "daily",  "wind_daily_sample.csv"),
]
RESOURCE_IDS = {  # サンプル resource id (DL 用)
    1274: 39775,
    1275: 39776,
    1276: 39777,
    1437: 93941, 1438: 177445,
    1439: 93942, 1440: 177446,
    1441: 93943, 1442: 93944,
}

# データ種別 ID → 観測種類 (master CSV の "データ種別" 列, 観測所マスタ調査結果より)
KIND_ID2NAME = {1: "rain", 4: "water", 7: "dam", 12: "tide", 13: "wind"}
KIND_LABEL = {
    "rain":  "雨量",
    "water": "水位",
    "dam":   "ダム",
    "tide":  "潮位",
    "wind":  "風向風速",
}
# 観測種類別カラー (防災運用上の重要度を直感に合わせ)
KIND_COLOR = {
    "rain":  "#1f77b4",  # 青  (降雨)
    "water": "#cf222e",  # 赤  (氾濫)
    "dam":   "#8b572a",  # 茶  (構造物)
    "tide":  "#2ca02c",  # 緑  (沿岸)
    "wind":  "#d4a72c",  # 黄  (気象)
}
KIND_ORDER = ["rain", "water", "dam", "tide", "wind"]

# 所管別カラー (国/県/その他)
ORIGIN_GROUP = {
    "砂防課":               "県",
    "河川課":               "県",
    "河川課ダムG":          "県",
    "農業基盤課":           "県",
    "国土交通省":           "国",
    "国ダム":               "国",
    "気象台":               "国",
    "港湾漁港整備課":       "県",
    "海上保安庁・港湾漁港整備課": "国",
}
ORIGIN_COLOR = {"国": "#0969da", "県": "#1a7f37", "その他": "#888888"}

# 広島県土面積 (国土地理院統計)
PREF_AREA_KM2 = 8479.6


# =============================================================================
# 1. 9 dataset の存在確認とメタデータ表構築 (auto fetch)
# =============================================================================
print("\n[1] 9 dataset の存在確認とメタデータ表構築", flush=True)
t1 = time.time()

for dsid, abbr, role, kind, gran, fname in SERIES:
    p = DATA_DIR / fname
    rid = RESOURCE_IDS.get(dsid)
    try:
        ensure_dataset(p, dataset_id=dsid, resource_id=rid,
                       label=f"L31/{dsid} {abbr}")
    except Exception as e:
        print(f"  WARN dsid={dsid}: {e}", flush=True)

# シリーズ表 (役割と粒度を明示)
series_df = pd.DataFrame([
    {
        "dsid": d, "abbr": ab, "role": r, "kind": k, "granularity": g,
        "file": f, "size_kb": round((DATA_DIR / f).stat().st_size / 1024, 1)
                    if (DATA_DIR / f).exists() else None,
    }
    for (d, ab, r, k, g, f) in SERIES
])
print(series_df.to_string(index=False), flush=True)


# =============================================================================
# 2. 観測所マスタ 1274 の読み込みと観測種類分布
# =============================================================================
print("\n[2] 観測所マスタを読み込み (634 点)", flush=True)
t1 = time.time()

# 観測所マスタは 3 行ヘッダ (大カテゴリ・小カテゴリ・列名) で、3 行目が実列名
sm_path = DATA_DIR / "stations_master.csv"
master = pd.read_csv(sm_path, encoding="utf-8-sig", header=2)
print(f"  master shape : {master.shape}", flush=True)

# データ種別 → 観測種類名
master["kind"] = master["データ種別"].map(KIND_ID2NAME)
unmapped = master[master["kind"].isna()]
if len(unmapped):
    print(f"  WARN: 未マップ データ種別 = {unmapped['データ種別'].unique()}", flush=True)
master = master.dropna(subset=["kind"]).copy()

# 所管グループ (国/県/その他)
master["origin"] = master["データ所管"].map(ORIGIN_GROUP).fillna("その他")

# 緯度経度 → Point geometry (4326 → 6671)
master = master.dropna(subset=["緯度", "経度"]).copy()
master_gdf = gpd.GeoDataFrame(
    master,
    geometry=gpd.points_from_xy(master["経度"], master["緯度"]),
    crs="EPSG:4326",
).to_crs(TARGET_CRS)
print(f"  geo points   : {len(master_gdf)} 点 (緯度経度あり)", flush=True)

# 種別別カウント
kind_counts = master_gdf["kind"].value_counts().reindex(KIND_ORDER, fill_value=0)
print("  観測種類別件数:")
for k, n in kind_counts.items():
    print(f"    {KIND_LABEL[k]:6s} : {n:4d} 点  ({n/len(master_gdf)*100:.1f}%)", flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 3. 時系列 8 件の wide CSV ヘッダ構造を読む (構造のみ、値は不要)
# =============================================================================
print("\n[3] 時系列 8 件の構造を解析", flush=True)
t1 = time.time()

def parse_ts_header(path: Path):
    """時系列 wide CSV の上部メタ (5-6行) を抽出。
    返り値:
      meta_df: 観測所列ごとに事務所名/所管/水系/河川/観測所番号/観測所名 を
               縦持ちした DataFrame
      n_cols: 列数 (=観測所数を推定)
    """
    raw = pd.read_csv(path, encoding="utf-8-sig", nrows=10, header=None)
    # 行0: 事務所名, 行1: データ所管, 行2: 水系名, 行3: 河川名, 行4: 観測所番号, 行5: 観測所名
    cols = list(range(1, raw.shape[1]))  # 0 列目はラベル列
    meta = pd.DataFrame({
        "src_col_idx": cols,
        "事務所名": raw.iloc[0, 1:].values,
        "データ所管": raw.iloc[1, 1:].values,
        "水系名": raw.iloc[2, 1:].values,
        "河川名": raw.iloc[3, 1:].values,
        "観測所番号": raw.iloc[4, 1:].values,
        "観測所名": raw.iloc[5, 1:].values,
    })
    # 観測所名がある行だけ
    meta = meta.dropna(subset=["観測所名"]).copy()
    # 同じ観測所が「10分雨量」「累計」などで複数列を占めるので、観測所名 unique
    meta_uniq = meta.drop_duplicates(subset=["観測所名"]).reset_index(drop=True)
    return meta_uniq, raw.shape[1]


ts_struct = []
ts_meta_per_kind = {}  # kind → DataFrame (観測所メタ集合)
for dsid, abbr, role, kind, gran, fname in SERIES:
    if role != "timeseries":
        continue
    p = DATA_DIR / fname
    if not p.exists():
        ts_struct.append({
            "dsid": dsid, "kind": kind, "granularity": gran,
            "n_total_cols": None, "n_unique_stations": None, "status": "missing",
        })
        continue
    try:
        meta, ncols = parse_ts_header(p)
        ts_struct.append({
            "dsid": dsid, "kind": kind, "granularity": gran,
            "n_total_cols": int(ncols),
            "n_unique_stations": int(len(meta)),
            "status": "ok",
        })
        # kind ごとに ( daily を優先、なければ yearly )
        if kind not in ts_meta_per_kind or gran == "daily":
            ts_meta_per_kind[kind] = meta
    except Exception as e:
        print(f"  ERR  dsid={dsid}: {e}", flush=True)
        ts_struct.append({
            "dsid": dsid, "kind": kind, "granularity": gran,
            "n_total_cols": None, "n_unique_stations": None, "status": f"err: {e}",
        })

ts_struct_df = pd.DataFrame(ts_struct)
print(ts_struct_df.to_string(index=False), flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 4. マスタ vs 時系列のクロス表 (各観測種類で master と timeseries が一致するか)
# =============================================================================
print("\n[4] マスタ vs 時系列カバレッジ照合", flush=True)
t1 = time.time()

cover_rows = []
for kind in KIND_ORDER:
    n_master = int((master_gdf["kind"] == kind).sum())
    if kind in ts_meta_per_kind:
        ts_meta = ts_meta_per_kind[kind]
        n_ts = len(ts_meta)
        # 観測所名で重複を除いた 集合
        master_names = set(master_gdf.loc[master_gdf["kind"] == kind, "観測所名"].astype(str))
        ts_names = set(ts_meta["観測所名"].astype(str))
        n_common = len(master_names & ts_names)
        n_only_master = len(master_names - ts_names)
        n_only_ts = len(ts_names - master_names)
    else:
        n_ts = 0; n_common = 0; n_only_master = n_master; n_only_ts = 0
    cover_rows.append({
        "kind": kind,
        "kind_label": KIND_LABEL[kind],
        "n_master": n_master,
        "n_ts": n_ts,
        "n_common_by_name": n_common,
        "n_only_master": n_only_master,
        "n_only_ts": n_only_ts,
    })
cover_df = pd.DataFrame(cover_rows)
cover_df.to_csv(ASSETS / "L31_master_vs_ts.csv", index=False, encoding="utf-8-sig")
print(cover_df.to_string(index=False), flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 5. 観測所属性の埋まり方 (氾濫閾値・ダム閾値・潮位閾値)
# =============================================================================
print("\n[5] 観測所属性の埋まり方を集計", flush=True)
t1 = time.time()

# 閾値属性 (4 グループ)
THRESH_GROUPS = {
    "河川氾濫閾値":   ["水防団待機[m]", "はん濫注意[m]", "避難判断[m]", "はん濫危険[m]"],
    "ダム水位":       ["最低水位[EL.m]", "制限水位[EL.m]", "常時満水位[EL.m]",
                       "ただし書き操作開始水位[EL.m]", "サーチャージ水位[EL.m]"],
    "潮位警戒":       ["注意潮位[TP.cm]", "警戒潮位[TP.cm]"],
    "その他":         ["零点高[T.P.m]", "洪水量[m3/s]"],
}

attr_rows = []
for grp, cols in THRESH_GROUPS.items():
    for col in cols:
        if col not in master_gdf.columns:
            continue
        n_filled = int(master_gdf[col].notna().sum())
        for kind in KIND_ORDER:
            sub = master_gdf[master_gdf["kind"] == kind]
            n_kind = int(sub[col].notna().sum())
            if n_kind > 0:
                attr_rows.append({
                    "group": grp, "column": col, "kind": kind,
                    "n_filled": n_kind, "n_total_kind": len(sub),
                    "fill_pct": round(n_kind / max(1, len(sub)) * 100, 1),
                })
attr_df = pd.DataFrame(attr_rows)
attr_df.to_csv(ASSETS / "L31_attr_fillrate.csv", index=False, encoding="utf-8-sig")
print(attr_df.to_string(index=False), flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# 氾濫閾値の保有 vs 監視のみ
master_gdf["has_flood_thresh"] = master_gdf[
    ["はん濫注意[m]", "避難判断[m]", "はん濫危険[m]"]
].notna().any(axis=1)
n_water_with = int((master_gdf[master_gdf["kind"] == "water"]["has_flood_thresh"]).sum())
n_water_total = int((master_gdf["kind"] == "water").sum())
print(f"  水位観測のうち 氾濫閾値あり: {n_water_with}/{n_water_total} "
      f"({n_water_with/max(1,n_water_total)*100:.1f}%)", flush=True)


# =============================================================================
# 6. 水系別の観測網密度
# =============================================================================
print("\n[6] 水系別の観測点密度", flush=True)
t1 = time.time()

# 水系名 NaN は "未分類"
master_gdf["水系名"] = master_gdf["水系名"].fillna("未分類")

# 水系 × 観測種類 ピボット
pv_water = pd.pivot_table(
    master_gdf, index="水系名", columns="kind", values="観測所名",
    aggfunc="count", fill_value=0,
).reindex(columns=KIND_ORDER, fill_value=0)
pv_water["合計"] = pv_water.sum(axis=1)
pv_water = pv_water.sort_values("合計", ascending=False)
pv_water.to_csv(ASSETS / "L31_pivot_watershed.csv", encoding="utf-8-sig")
print("  Top 10 水系:")
print(pv_water.head(10).to_string(), flush=True)


# 主要 3 水系の集中率
top3 = pv_water.head(3)
top3_total = int(top3["合計"].sum())
top3_pct = top3_total / len(master_gdf) * 100
print(f"  上位 3 水系で全観測点の {top3_pct:.1f}% を占める "
      f"({top3_total}/{len(master_gdf)})", flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 7. 市町別の観測網密度
# =============================================================================
print("\n[7] 市町別の観測点密度", flush=True)
t1 = time.time()

master_gdf["市町村名"] = master_gdf["市町村名"].fillna("未分類")
pv_city = pd.pivot_table(
    master_gdf, index="市町村名", columns="kind", values="観測所名",
    aggfunc="count", fill_value=0,
).reindex(columns=KIND_ORDER, fill_value=0)
pv_city["合計"] = pv_city.sum(axis=1)
pv_city = pv_city.sort_values("合計", ascending=False)
pv_city.to_csv(ASSETS / "L31_pivot_city.csv", encoding="utf-8-sig")
print("  Top 15 市町:")
print(pv_city.head(15).to_string(), flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 8. 観測点間距離 (kNN による空間粒度)
# =============================================================================
print("\n[8] 観測点間距離 (k-NN) を計算", flush=True)
t1 = time.time()

# ★ 全観測点同士の最近隣距離 (kind 不問、観測網全体の粒度)
xy = np.array([(p.x, p.y) for p in master_gdf.geometry])
tree = cKDTree(xy)
# 自分自身を除く最近隣 = k=2 の 2 番目
dists, idxs = tree.query(xy, k=2)
nearest_m = dists[:, 1]
master_gdf["nearest_any_m"] = nearest_m

# kind 別: 同じ kind 内の最近隣
for k in KIND_ORDER:
    sub = master_gdf[master_gdf["kind"] == k]
    if len(sub) < 2:
        master_gdf.loc[master_gdf["kind"] == k, f"nearest_{k}_m"] = np.nan
        continue
    sub_xy = np.array([(p.x, p.y) for p in sub.geometry])
    sub_tree = cKDTree(sub_xy)
    sd, _ = sub_tree.query(sub_xy, k=2)
    nearest_same_m = sd[:, 1]
    master_gdf.loc[sub.index, f"nearest_{k}_m"] = nearest_same_m

# 種別別の最近隣統計
near_stats = []
for k in KIND_ORDER:
    arr = master_gdf.loc[master_gdf["kind"] == k, f"nearest_{k}_m"].dropna()
    if len(arr) < 1:
        continue
    near_stats.append({
        "kind": k, "kind_label": KIND_LABEL[k],
        "n": len(arr),
        "median_km": round(float(np.median(arr) / 1000), 2),
        "mean_km": round(float(np.mean(arr) / 1000), 2),
        "p95_km": round(float(np.percentile(arr, 95) / 1000), 2),
        "max_km": round(float(np.max(arr) / 1000), 2),
    })
# 全観測点 (任意種別) の最近隣
arr_any = master_gdf["nearest_any_m"]
near_stats.append({
    "kind": "any", "kind_label": "全観測点",
    "n": len(arr_any),
    "median_km": round(float(np.median(arr_any) / 1000), 2),
    "mean_km": round(float(np.mean(arr_any) / 1000), 2),
    "p95_km": round(float(np.percentile(arr_any, 95) / 1000), 2),
    "max_km": round(float(np.max(arr_any) / 1000), 2),
})
near_df = pd.DataFrame(near_stats)
near_df.to_csv(ASSETS / "L31_nearest_stats.csv", index=False, encoding="utf-8-sig")
print(near_df.to_string(index=False), flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 9. 観測空白地帯 (県土グリッドで最寄り観測点までの距離)
# =============================================================================
print("\n[9] 観測空白地帯マップ用 grid を作成", flush=True)
t1 = time.time()

# 県全域 polygon の load
import zipfile, io as _io
def load_zip_first_geo(zip_path: Path, encoding="cp932"):
    with zipfile.ZipFile(zip_path) as zf:
        names = zf.namelist()
        for n in names:
            if n.lower().endswith(".geojson"):
                with zf.open(n) as f:
                    return gpd.read_file(_io.BytesIO(f.read()))
        for n in names:
            if n.lower().endswith(".shp"):
                tmp = zip_path.parent / f"_ext_{zip_path.stem}"
                tmp.mkdir(exist_ok=True)
                zf.extractall(tmp)
                return gpd.read_file(tmp / n, encoding=encoding)
    raise FileNotFoundError(zip_path)

g_admin_pref = load_zip_first_geo(ADMIN_DIR / "admin_922_広島県.zip").to_crs(TARGET_CRS)
pref_diss = g_admin_pref.dissolve().reset_index(drop=True)
pref_geom = pref_diss.geometry.iloc[0]

# 1km grid (要件 S 準拠で軽量に)
xmin, ymin, xmax, ymax = pref_geom.bounds
GRID = 1500  # 1.5 km
xs = np.arange(xmin, xmax + GRID, GRID)
ys = np.arange(ymin, ymax + GRID, GRID)
xx, yy = np.meshgrid(xs, ys)
grid_pts = np.column_stack([xx.ravel(), yy.ravel()])

# 各 grid 点が県内かどうか sindex で高速判定
from shapely.geometry import MultiPoint
grid_geo = gpd.GeoSeries([Point(x, y) for x, y in grid_pts], crs=TARGET_CRS)
in_pref = grid_geo.within(pref_geom)
print(f"  grid total = {len(grid_pts)}, in_pref = {int(in_pref.sum())}", flush=True)
grid_in = grid_pts[in_pref.values]

# 各 grid_in 点について 最寄り雨量観測点までの距離
rain_xy = np.array([(p.x, p.y) for p in master_gdf[master_gdf["kind"] == "rain"].geometry])
rain_tree = cKDTree(rain_xy)
gd, _ = rain_tree.query(grid_in, k=1)
grid_in_dist_rain_km = gd / 1000.0

# 全種別 の最寄り
all_xy = np.array([(p.x, p.y) for p in master_gdf.geometry])
all_tree = cKDTree(all_xy)
gd_all, _ = all_tree.query(grid_in, k=1)
grid_in_dist_all_km = gd_all / 1000.0

# 統計
print(f"  雨量最寄り distance (km): "
      f"median={np.median(grid_in_dist_rain_km):.2f}, "
      f"p95={np.percentile(grid_in_dist_rain_km, 95):.2f}, "
      f"max={np.max(grid_in_dist_rain_km):.2f}", flush=True)
print(f"  全種別最寄り distance (km): "
      f"median={np.median(grid_in_dist_all_km):.2f}, "
      f"p95={np.percentile(grid_in_dist_all_km, 95):.2f}, "
      f"max={np.max(grid_in_dist_all_km):.2f}", flush=True)

# 観測空白の閾値: 4 km (雨量網は当初想定より遥かに密だった)
# H6 当初は 8 km と置いたが、実測中央値 2.19 km / max 10.55 km なので 4 km が
# 「明らかな空白」の妥当な基準と判断して再設定。
blank_thresh_km = 4.0
n_blank = int((grid_in_dist_rain_km > blank_thresh_km).sum())
n_in = int(in_pref.sum())
print(f"  雨量で {blank_thresh_km} km 以上の grid: {n_blank}/{n_in} ({n_blank/n_in*100:.2f}%)", flush=True)
# 8 km 超の grid 数も参考までに
n_blank_8 = int((grid_in_dist_rain_km > 8.0).sum())
print(f"  雨量で 8 km 以上の grid: {n_blank_8}/{n_in} ({n_blank_8/n_in*100:.4f}%) (参考: 当初想定の閾値)", flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 10. 中間 CSV 出力 (master gdf を 6671→4326 に戻して書き出し)
# =============================================================================
print("\n[10] 中間 CSV 出力", flush=True)
t1 = time.time()

# シリーズ表
series_df.to_csv(ASSETS / "L31_series_9.csv", index=False, encoding="utf-8-sig")

# 観測所マスタ (geo 列なしで)
mout = master_gdf.drop(columns=["geometry"]).copy()
mout["kind_label"] = mout["kind"].map(KIND_LABEL)
mout.to_csv(ASSETS / "L31_stations_full.csv", index=False, encoding="utf-8-sig")

# 種別別観測所件数
kc = pd.DataFrame({
    "kind": KIND_ORDER,
    "kind_label": [KIND_LABEL[k] for k in KIND_ORDER],
    "n_stations": [int(kind_counts[k]) for k in KIND_ORDER],
    "pct": [round(int(kind_counts[k]) / len(master_gdf) * 100, 2) for k in KIND_ORDER],
})
kc.to_csv(ASSETS / "L31_kind_counts.csv", index=False, encoding="utf-8-sig")

# 所管別観測所件数
oc = master_gdf.groupby(["origin", "データ所管"]).size().reset_index(name="n").sort_values("n", ascending=False)
oc.to_csv(ASSETS / "L31_origin_breakdown.csv", index=False, encoding="utf-8-sig")

# 時系列構造表
ts_struct_df.to_csv(ASSETS / "L31_ts_structure.csv", index=False, encoding="utf-8-sig")

# 観測空白 grid (距離値だけ)
blank_df = pd.DataFrame({
    "x": grid_in[:, 0], "y": grid_in[:, 1],
    "dist_rain_km": grid_in_dist_rain_km,
    "dist_all_km": grid_in_dist_all_km,
})
blank_df.to_csv(ASSETS / "L31_blank_grid.csv", index=False, encoding="utf-8-sig")
print("  saved 7 csv", flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 11. 図 1: 観測種類別の件数バー (5 種比較 + 所管内訳)
# =============================================================================
print("\n[11] 図 1: 種別件数バー", flush=True)
t1 = time.time()

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

ax = axes[0]
labels = [KIND_LABEL[k] for k in KIND_ORDER]
counts = [int(kind_counts[k]) for k in KIND_ORDER]
colors = [KIND_COLOR[k] for k in KIND_ORDER]
bars = ax.bar(labels, counts, color=colors, edgecolor="#222")
for b, c in zip(bars, counts):
    ax.text(b.get_x() + b.get_width()/2, c + 5, f"{c}\n({c/len(master_gdf)*100:.1f}%)",
            ha="center", va="bottom", fontsize=10)
ax.set_ylabel("観測所数 [点]")
ax.set_title("観測種類別の観測所数 (合計 {} 点)".format(len(master_gdf)))
ax.set_ylim(0, max(counts) * 1.18)
ax.grid(axis="y", alpha=0.3)

# 種別 × 所管 (国/県/その他) スタックバー
ax = axes[1]
pv_kind_origin = pd.pivot_table(
    master_gdf, index="kind", columns="origin", values="観測所名",
    aggfunc="count", fill_value=0,
).reindex(KIND_ORDER, fill_value=0)
for col in ["国", "県", "その他"]:
    if col not in pv_kind_origin.columns:
        pv_kind_origin[col] = 0
pv_kind_origin = pv_kind_origin[["県", "国", "その他"]]
bot = np.zeros(len(KIND_ORDER))
for col in ["県", "国", "その他"]:
    vals = pv_kind_origin[col].values
    ax.bar([KIND_LABEL[k] for k in KIND_ORDER], vals, bottom=bot,
           color=ORIGIN_COLOR[col], label=col, edgecolor="#222", linewidth=0.4)
    bot = bot + vals
ax.set_title("種別 × 所管 (国 / 県 / その他)")
ax.set_ylabel("観測所数 [点]")
ax.legend(loc="upper right", fontsize=10)
ax.grid(axis="y", alpha=0.3)

plt.tight_layout()
fig.savefig(ASSETS / "L31_fig1_kind_counts.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L31_fig1_kind_counts.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 12. 図 2: 県全域 観測網マップ (5 種カラー散布)
# =============================================================================
print("\n[12] 図 2: 県全域マップ", flush=True)
t1 = time.time()

fig, ax = plt.subplots(figsize=(11, 8))
g_admin_pref.boundary.plot(ax=ax, color="#888", linewidth=0.4)
pref_diss.boundary.plot(ax=ax, color="#222", linewidth=1.2)
for k in KIND_ORDER:
    sub = master_gdf[master_gdf["kind"] == k]
    sub.plot(ax=ax, color=KIND_COLOR[k], markersize=18,
             alpha=0.78, edgecolor="white", linewidth=0.4)
legend_handles = [
    Line2D([0], [0], marker="o", color="w", markerfacecolor=KIND_COLOR[k],
           markersize=9, label=f"{KIND_LABEL[k]} ({int(kind_counts[k])})")
    for k in KIND_ORDER
]
ax.legend(handles=legend_handles, loc="lower left", fontsize=11,
          title="観測種類")
ax.set_title("広島県 防災観測網 — 全 {} 点を 5 種別で色分け".format(len(master_gdf)))
ax.set_xticks([]); ax.set_yticks([])
ax.set_aspect("equal")
plt.tight_layout()
fig.savefig(ASSETS / "L31_fig2_pref_overview.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L31_fig2_pref_overview.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 13. 図 3: 5 種類 small multiples マップ (種別ごとに県全域)
# =============================================================================
print("\n[13] 図 3: 種別 small multiples", flush=True)
t1 = time.time()

fig, axes = plt.subplots(2, 3, figsize=(14, 9))
axes_flat = axes.flatten()
for i, k in enumerate(KIND_ORDER):
    ax = axes_flat[i]
    g_admin_pref.boundary.plot(ax=ax, color="#bbb", linewidth=0.3)
    pref_diss.boundary.plot(ax=ax, color="#222", linewidth=0.8)
    sub = master_gdf[master_gdf["kind"] == k]
    sub.plot(ax=ax, color=KIND_COLOR[k], markersize=12, alpha=0.85,
             edgecolor="white", linewidth=0.3)
    ax.set_title(f"{KIND_LABEL[k]} ({len(sub)} 点)", fontsize=12)
    ax.set_xticks([]); ax.set_yticks([])
    ax.set_aspect("equal")
# 6 番目: 全種重ね
ax = axes_flat[5]
g_admin_pref.boundary.plot(ax=ax, color="#bbb", linewidth=0.3)
pref_diss.boundary.plot(ax=ax, color="#222", linewidth=0.8)
for k in KIND_ORDER:
    sub = master_gdf[master_gdf["kind"] == k]
    sub.plot(ax=ax, color=KIND_COLOR[k], markersize=8, alpha=0.7,
             edgecolor="none")
ax.set_title("5 種重ね合わせ", fontsize=12)
ax.set_xticks([]); ax.set_yticks([])
ax.set_aspect("equal")

fig.suptitle("観測種類別 small multiples マップ — 種別ごとの地理偏在", fontsize=14, y=1.00)
plt.tight_layout()
fig.savefig(ASSETS / "L31_fig3_kind_smallmult.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L31_fig3_kind_smallmult.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 14. 図 4: 水系別の観測網密度 (上位 10 水系 横棒)
# =============================================================================
print("\n[14] 図 4: 水系別観測網密度", flush=True)
t1 = time.time()

top_n = 12
top_pv = pv_water.head(top_n).copy()
fig, ax = plt.subplots(figsize=(10, 6))
y = np.arange(len(top_pv))
left = np.zeros(len(top_pv))
for k in KIND_ORDER:
    if k not in top_pv.columns:
        continue
    vals = top_pv[k].values
    ax.barh(y, vals, left=left, color=KIND_COLOR[k],
            label=KIND_LABEL[k], edgecolor="#222", linewidth=0.3)
    left = left + vals
ax.set_yticks(y)
ax.set_yticklabels(top_pv.index)
ax.invert_yaxis()
ax.set_xlabel("観測点数")
ax.set_title("水系別の観測網密度 (上位 {} 水系)".format(top_n))
ax.legend(loc="lower right", fontsize=10)
ax.grid(axis="x", alpha=0.3)
for i, total in enumerate(top_pv["合計"]):
    ax.text(total + 1.5, i, f"  {int(total)}", va="center", fontsize=10)
plt.tight_layout()
fig.savefig(ASSETS / "L31_fig4_watershed_density.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L31_fig4_watershed_density.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 15. 図 5: 市町別の観測網密度 (上位 15 市町 横棒)
# =============================================================================
print("\n[15] 図 5: 市町別観測網密度", flush=True)
t1 = time.time()

top_n_city = 15
top_pv_c = pv_city.head(top_n_city).copy()
fig, ax = plt.subplots(figsize=(10, 6.5))
y = np.arange(len(top_pv_c))
left = np.zeros(len(top_pv_c))
for k in KIND_ORDER:
    if k not in top_pv_c.columns:
        continue
    vals = top_pv_c[k].values
    ax.barh(y, vals, left=left, color=KIND_COLOR[k],
            label=KIND_LABEL[k], edgecolor="#222", linewidth=0.3)
    left = left + vals
ax.set_yticks(y)
ax.set_yticklabels(top_pv_c.index)
ax.invert_yaxis()
ax.set_xlabel("観測点数")
ax.set_title("市町別の観測網密度 (上位 {} 市町)".format(top_n_city))
ax.legend(loc="lower right", fontsize=10)
ax.grid(axis="x", alpha=0.3)
for i, total in enumerate(top_pv_c["合計"]):
    ax.text(total + 1.0, i, f"  {int(total)}", va="center", fontsize=10)
plt.tight_layout()
fig.savefig(ASSETS / "L31_fig5_city_density.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L31_fig5_city_density.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 16. 図 6: 氾濫閾値マップ (水位観測のうち閾値あり vs なし)
# =============================================================================
print("\n[16] 図 6: 氾濫閾値マップ", flush=True)
t1 = time.time()

water_gdf = master_gdf[master_gdf["kind"] == "water"].copy()
water_with = water_gdf[water_gdf["has_flood_thresh"]]
water_without = water_gdf[~water_gdf["has_flood_thresh"]]

fig, axes = plt.subplots(1, 2, figsize=(13, 6.5))
# 左: 閾値あり/なし マップ
ax = axes[0]
g_admin_pref.boundary.plot(ax=ax, color="#bbb", linewidth=0.3)
pref_diss.boundary.plot(ax=ax, color="#222", linewidth=0.8)
water_without.plot(ax=ax, color="#888", markersize=14, alpha=0.7,
                    edgecolor="white", linewidth=0.3, label=f"閾値なし ({len(water_without)})")
water_with.plot(ax=ax, color="#cf222e", markersize=22, alpha=0.85,
                edgecolor="white", linewidth=0.3, label=f"閾値あり ({len(water_with)})")
ax.legend(loc="lower left", fontsize=11)
ax.set_title("水位観測 × 氾濫閾値の有無")
ax.set_xticks([]); ax.set_yticks([])
ax.set_aspect("equal")

# 右: 閾値段階数の分布
ax = axes[1]
thresh_cols = ["水防団待機[m]", "はん濫注意[m]", "避難判断[m]", "はん濫危険[m]"]
water_gdf["n_thresh"] = water_gdf[thresh_cols].notna().sum(axis=1)
nb = water_gdf["n_thresh"].value_counts().sort_index()
ax.bar([f"{int(i)} 段階" for i in nb.index], nb.values,
       color="#cf222e", edgecolor="#222")
for i, v in enumerate(nb.values):
    ax.text(i, v + 1, f"{v}", ha="center", fontsize=11)
ax.set_ylabel("観測所数")
ax.set_title("水位観測 × 氾濫閾値の段階数分布")
ax.grid(axis="y", alpha=0.3)

plt.tight_layout()
fig.savefig(ASSETS / "L31_fig6_flood_threshold.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L31_fig6_flood_threshold.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 17. 図 7: 観測空白地帯マップ (1.5km grid × 最寄り雨量距離)
# =============================================================================
print("\n[17] 図 7: 観測空白地帯マップ", flush=True)
t1 = time.time()

fig, axes = plt.subplots(1, 2, figsize=(14, 6.5))

vmax_km = 6.0  # 観測網は密 (max 10 km) なので、6 km をカラーバー上限にするとコントラスト最大

ax = axes[0]
# 雨量最寄りグリッドのヒートマップ
sc = ax.scatter(grid_in[:, 0], grid_in[:, 1], c=grid_in_dist_rain_km,
                s=8, cmap="YlOrRd", vmin=0, vmax=vmax_km, marker="s")
g_admin_pref.boundary.plot(ax=ax, color="#888", linewidth=0.3)
pref_diss.boundary.plot(ax=ax, color="#222", linewidth=0.8)
# 観測点 overlay (雨量だけ)
master_gdf[master_gdf["kind"] == "rain"].plot(
    ax=ax, color="#0a3d62", markersize=4, alpha=0.55, edgecolor="none")
ax.set_title(f"雨量観測 最寄り距離 (km) — グリッドが赤いほど空白\n(色の上限 {vmax_km} km、青点は雨量観測所)")
cbar = plt.colorbar(sc, ax=ax, shrink=0.8)
cbar.set_label("最寄り雨量観測点までの距離 [km]")
ax.set_xticks([]); ax.set_yticks([])
ax.set_aspect("equal")

ax = axes[1]
# 全観測種別に対する最寄り
sc = ax.scatter(grid_in[:, 0], grid_in[:, 1], c=grid_in_dist_all_km,
                s=8, cmap="YlOrRd", vmin=0, vmax=vmax_km, marker="s")
g_admin_pref.boundary.plot(ax=ax, color="#888", linewidth=0.3)
pref_diss.boundary.plot(ax=ax, color="#222", linewidth=0.8)
master_gdf.plot(ax=ax, color="#0a3d62", markersize=4, alpha=0.55, edgecolor="none")
ax.set_title(f"全 5 種類 最寄り距離 (km) — どんな観測も無い空白\n(色の上限 {vmax_km} km、青点は全観測所)")
cbar = plt.colorbar(sc, ax=ax, shrink=0.8)
cbar.set_label("最寄り観測点までの距離 [km]")
ax.set_xticks([]); ax.set_yticks([])
ax.set_aspect("equal")

plt.tight_layout()
fig.savefig(ASSETS / "L31_fig7_blank_zones.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L31_fig7_blank_zones.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 18. 図 8: 最近隣距離の分布 (kind 別ヒストグラム)
# =============================================================================
print("\n[18] 図 8: 最近隣距離分布", flush=True)
t1 = time.time()

fig, axes = plt.subplots(2, 3, figsize=(14, 7.5))
axes_flat = axes.flatten()
for i, k in enumerate(KIND_ORDER):
    ax = axes_flat[i]
    arr = master_gdf.loc[master_gdf["kind"] == k, f"nearest_{k}_m"].dropna() / 1000
    if len(arr) < 2:
        ax.text(0.5, 0.5, "n < 2", transform=ax.transAxes, ha="center")
        ax.set_title(f"{KIND_LABEL[k]} (n={len(arr)})")
        continue
    bins = np.linspace(0, max(20, np.percentile(arr, 95)), 25)
    ax.hist(arr, bins=bins, color=KIND_COLOR[k], edgecolor="#222")
    med = np.median(arr); p95 = np.percentile(arr, 95)
    ax.axvline(med, color="#cf222e", linestyle="--", linewidth=1, label=f"中央値 {med:.1f}km")
    ax.axvline(p95, color="#0969da", linestyle=":", linewidth=1, label=f"95%値 {p95:.1f}km")
    ax.set_title(f"{KIND_LABEL[k]} (n={len(arr)})")
    ax.set_xlabel("同種別の最近隣 [km]")
    ax.set_ylabel("観測所数")
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)
# 6番目: 全種別 (any)
ax = axes_flat[5]
arr = master_gdf["nearest_any_m"].dropna() / 1000
bins = np.linspace(0, 20, 25)
ax.hist(arr, bins=bins, color="#444", edgecolor="#000")
med = np.median(arr); p95 = np.percentile(arr, 95)
ax.axvline(med, color="#cf222e", linestyle="--", linewidth=1, label=f"中央値 {med:.1f}km")
ax.axvline(p95, color="#0969da", linestyle=":", linewidth=1, label=f"95%値 {p95:.1f}km")
ax.set_title(f"全観測点 (n={len(arr)}) — 種別不問")
ax.set_xlabel("最近隣 [km]")
ax.set_ylabel("観測所数")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)

fig.suptitle("種別ごとの観測点 最近隣距離分布", fontsize=14, y=1.00)
plt.tight_layout()
fig.savefig(ASSETS / "L31_fig8_nearest_dist.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L31_fig8_nearest_dist.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 19. 図 9: 9 dataset シリーズ構造マトリクス (役割×粒度)
# =============================================================================
print("\n[19] 図 9: 9 dataset 構造マトリクス", flush=True)
t1 = time.time()

# series_df を kind × granularity 表に整形
matrix = pd.DataFrame(index=KIND_ORDER + ["master"], columns=["daily", "yearly", "master"])
for r in series_df.itertuples():
    if r.role == "master":
        matrix.loc["master", "master"] = r.dsid
    else:
        matrix.loc[r.kind, r.granularity] = r.dsid
matrix = matrix.fillna("")

fig, ax = plt.subplots(figsize=(11, 6))
# データ種別を行、粒度を列
rows = ["master"] + KIND_ORDER
cols = ["master", "daily", "yearly"]
mat_disp = np.zeros((len(rows), len(cols)))
mat_text = np.empty((len(rows), len(cols)), dtype=object)
for i, rk in enumerate(rows):
    for j, c in enumerate(cols):
        v = matrix.loc[rk, c] if rk in matrix.index and c in matrix.columns else ""
        if v != "":
            mat_disp[i, j] = 1
            mat_text[i, j] = str(v)
        else:
            mat_text[i, j] = "—"
ax.imshow(mat_disp, cmap="Blues", vmin=0, vmax=1.5, aspect="auto")
for i, rk in enumerate(rows):
    for j, c in enumerate(cols):
        v = mat_text[i, j]
        col = "#000" if mat_disp[i, j] > 0 else "#888"
        ax.text(j, i, v, ha="center", va="center", color=col, fontsize=11,
                fontweight="bold" if mat_disp[i, j] > 0 else "normal")
ax.set_xticks(range(len(cols)))
ax.set_xticklabels(["master\n(緯度経度+属性)", "daily\n(時系列, 10分)", "yearly\n(時系列, 1日)"])
ax.set_yticks(range(len(rows)))
ax.set_yticklabels([{"master": "観測所マスタ", **KIND_LABEL}[r] for r in rows])
ax.set_title("9 dataset_id の構造マトリクス\n(マスタ 1 件 + 5 種類 × {日, 年} 8 件)")
plt.tight_layout()
fig.savefig(ASSETS / "L31_fig9_series_matrix.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L31_fig9_series_matrix.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 20. 図 10: マスタ vs 時系列 カバレッジ照合 (二段棒)
# =============================================================================
print("\n[20] 図 10: master vs ts 照合", flush=True)
t1 = time.time()

fig, ax = plt.subplots(figsize=(11, 5.5))
x = np.arange(len(KIND_ORDER))
w = 0.35
n_m_arr = [int(cover_df[cover_df["kind"] == k]["n_master"].iloc[0]) for k in KIND_ORDER]
n_t_arr = [int(cover_df[cover_df["kind"] == k]["n_ts"].iloc[0]) for k in KIND_ORDER]
n_c_arr = [int(cover_df[cover_df["kind"] == k]["n_common_by_name"].iloc[0]) for k in KIND_ORDER]
ax.bar(x - w/2, n_m_arr, width=w, color="#1a7f37", label="マスタ件数", edgecolor="#222")
ax.bar(x + w/2, n_t_arr, width=w, color="#0969da", label="時系列ファイル件数", edgecolor="#222")
for i, (n_m, n_t, n_c) in enumerate(zip(n_m_arr, n_t_arr, n_c_arr)):
    ax.text(i - w/2, n_m + 5, f"{n_m}", ha="center", fontsize=10)
    ax.text(i + w/2, n_t + 5, f"{n_t}", ha="center", fontsize=10)
    ax.text(i, max(n_m, n_t) + 25, f"共通\n{n_c}", ha="center", fontsize=9, color="#444")
ax.set_xticks(x)
ax.set_xticklabels([KIND_LABEL[k] for k in KIND_ORDER])
ax.set_ylabel("観測所数")
ax.set_title("マスタ vs 時系列ファイル: 観測所カバレッジ照合")
ax.legend(loc="upper right")
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
fig.savefig(ASSETS / "L31_fig10_master_ts_coverage.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L31_fig10_master_ts_coverage.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 21. 図 11: 9 dataset カードビュー (役割と件数を 1 枚で俯瞰)
# =============================================================================
print("\n[21] 図 11: 9 dataset カードビュー", flush=True)
t1 = time.time()

# kind ごとの件数 + dataset id を 1 枚にまとめる
fig, ax = plt.subplots(figsize=(13, 6))
ax.set_xlim(0, 9)
ax.set_ylim(0, 5)
ax.set_axis_off()

# レイアウト: 上段 master 1 枚, 中段 kind × {daily, yearly}
# master card
ax.add_patch(plt.Rectangle((0.2, 3.6), 8.6, 1.0, fill=True, color="#fff8c5",
                            edgecolor="#d4a72c", linewidth=2))
ax.text(4.5, 4.1,
        f"観測所マスタ (dsid=1274) — CSV 634 行 × 35 列 / 緯度経度・水系・河川・市町村・閾値\n"
        f"全 5 観測種類の地理キー、属性 (はん濫注意/避難判断/はん濫危険/常時満水位/注意潮位)",
        ha="center", va="center", fontsize=11)

# kind × granularity
for i, k in enumerate(KIND_ORDER):
    color = KIND_COLOR[k]
    n = int(kind_counts[k])
    # daily (left)
    daily_dsid = series_df[(series_df["kind"] == k) & (series_df["granularity"] == "daily")]
    yearly_dsid = series_df[(series_df["kind"] == k) & (series_df["granularity"] == "yearly")]
    x0 = 0.2 + i * 1.74
    # daily card
    ax.add_patch(plt.Rectangle((x0, 1.7), 1.6, 1.5, fill=True,
                                color=color, alpha=0.18, edgecolor=color, linewidth=2))
    ax.text(x0 + 0.8, 2.85, f"{KIND_LABEL[k]}", ha="center", va="center",
            fontsize=12, fontweight="bold", color=color)
    ax.text(x0 + 0.8, 2.45, f"{n} 観測所", ha="center", va="center", fontsize=11)
    if len(daily_dsid):
        ax.text(x0 + 0.8, 2.10, f"日集計 dsid={int(daily_dsid['dsid'].iloc[0])}",
                ha="center", va="center", fontsize=9, color="#444")
    # yearly card (under)
    ax.add_patch(plt.Rectangle((x0, 0.4), 1.6, 1.0, fill=True,
                                color=color, alpha=0.08, edgecolor=color,
                                linewidth=1.2, linestyle="--"))
    if len(yearly_dsid):
        ax.text(x0 + 0.8, 0.9, f"年集計\ndsid={int(yearly_dsid['dsid'].iloc[0])}",
                ha="center", va="center", fontsize=9, color="#444")
    else:
        ax.text(x0 + 0.8, 0.9, "(年集計なし)",
                ha="center", va="center", fontsize=9, color="#999", style="italic")

ax.text(4.5, 4.8, "9 dataset = 1 マスタ + 5 種類 × { 日集計, 年集計 } の 2 層構造",
        ha="center", fontsize=14, fontweight="bold")
plt.tight_layout()
fig.savefig(ASSETS / "L31_fig11_dataset_cards.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L31_fig11_dataset_cards.png ({time.time()-t1:.1f}s)", flush=True)


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

hypos = []

# H1 種別不均衡: 雨量 > 水位 ≫ 風 ≫ ダム ≈ 潮位 / 雨量だけで 60% 以上
n_total = len(master_gdf)
n_rain = int(kind_counts["rain"])
n_water = int(kind_counts["water"])
order_ok = n_rain > n_water > kind_counts["wind"] > kind_counts["dam"] >= kind_counts["tide"] - 1
rain_pct = n_rain / n_total * 100
h1 = order_ok and rain_pct >= 60.0
hypos.append({
    "H": "H1",
    "claim": "種別不均衡: 雨量 > 水位 > 風 > ダム ≈ 潮位 / 雨量 60%+",
    "result": (
        f"順序: 雨{n_rain} > 水位{n_water} > 風{int(kind_counts['wind'])} > "
        f"ダム{int(kind_counts['dam'])} > 潮位{int(kind_counts['tide'])}, "
        f"雨量比 {rain_pct:.1f}%"
    ),
    "verdict": "支持" if h1 else ("部分支持" if rain_pct >= 50 else "反証"),
})

# H2 所管二元性: 県管理 vs 国管理 が ほぼ折半
n_state = int((master_gdf["origin"] == "県").sum())
n_natio = int((master_gdf["origin"] == "国").sum())
ratio = n_natio / n_state if n_state else 0
h2 = 0.6 < (n_natio / max(1, n_state + n_natio)) < 0.6 or 0.4 < n_state / (n_state + n_natio) < 0.7
h2 = abs(n_state - n_natio) / max(1, (n_state + n_natio)) < 0.4
hypos.append({
    "H": "H2",
    "claim": "所管二元性: 県管理 vs 国管理 がほぼ折半 (差 < 40%)",
    "result": f"県={n_state}, 国={n_natio}, 差={abs(n_state-n_natio)/(n_state+n_natio)*100:.1f}%",
    "verdict": "支持" if h2 else "反証",
})

# H3 主要 3 水系 50%以上
top3_total_v = int(top3["合計"].sum())
top3_pct_v = top3_total_v / n_total * 100
h3 = top3_pct_v >= 50.0
top3_names = list(top3.index)
hypos.append({
    "H": "H3",
    "claim": "主要 3 水系 (太田川/江の川/芦田川 等) で観測点 50%+",
    "result": f"上位 3 水系 {top3_names} で {top3_pct_v:.1f}% ({top3_total_v}/{n_total})",
    "verdict": "支持" if h3 else "部分支持",
})

# H4 潮位は瀬戸内沿岸 / 内陸ゼロ
tide_gdf = master_gdf[master_gdf["kind"] == "tide"]
# 内陸判定: 観測所の住所に「港」や「沿岸」が含まれるか、 海からの距離 < 1km か
# 簡易: 緯度が 34.5 未満 = ほぼ沿岸 (広島県の海岸線は ほぼ 34.5 以南)
tide_inland = int((tide_gdf["緯度"] > 34.6).sum())
h4 = tide_inland == 0
hypos.append({
    "H": "H4",
    "claim": "潮位 13 点は全て沿岸 (緯度 34.6 以南)",
    "result": f"内陸 (緯度 > 34.6) の潮位観測 = {tide_inland} 件",
    "verdict": "支持" if h4 else "反証",
})

# H5 水位 50%+ で氾濫閾値あり
h5_pct = n_water_with / max(1, n_water_total) * 100
h5 = h5_pct >= 50.0
hypos.append({
    "H": "H5",
    "claim": "水位観測の 50%+ で氾濫閾値あり",
    "result": f"水位 {n_water_total} 中 {n_water_with} 件 ({h5_pct:.1f}%) が閾値あり",
    "verdict": "支持" if h5 else "部分支持",
})

# H6 4 km 以上の観測空白 (当初 8 km と置いたが、実測で雨量最大空白 10.55 km / 中央値 2.19 km と
# 観測網が極めて密だったため、「明らかな空白」の基準を 4 km に下方修正)
n_blank_pct = n_blank / max(1, n_in) * 100
n_blank_8_pct = n_blank_8 / max(1, n_in) * 100
hypos.append({
    "H": "H6",
    "claim": (f"雨量で {blank_thresh_km} km 以上の観測空白が grid 5%以上 25%未満で存在 "
              f"(当初 8km 想定だったが観測網が密だったため 4km へ閾値修正)"),
    "result": (f"grid {n_blank}/{n_in} = {n_blank_pct:.2f}% が {blank_thresh_km} km 超 "
               f"(8 km 超は {n_blank_8_pct:.4f}% で極小)"),
    "verdict": "支持" if (5 <= n_blank_pct < 25) else
               ("部分支持" if 1 <= n_blank_pct < 5 else "反証"),
})

hypos_df = pd.DataFrame(hypos)
hypos_df.to_csv(ASSETS / "L31_hypothesis_results.csv", index=False, encoding="utf-8-sig")
with (ASSETS / "L31_hypothesis_results.json").open("w", encoding="utf-8") as f:
    json.dump(hypos, f, ensure_ascii=False, indent=2)
print(hypos_df.to_string(index=False), flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 23. HTML 生成
# =============================================================================
print("\n[23] HTML 生成", flush=True)
t1 = time.time()

# このスクリプトのソース全文を読み込む
this_py = Path(__file__).read_text(encoding="utf-8")

# CSV/PNG download links helper
def dl(name):
    return f'<a href="assets/{name}" download>{name}</a>'


sections = []

# --- Section 1: 学習目標と問い ---
sec1 = f"""
<div class="note">
<b>カバー宣言</b>: 本記事は DoBoX のシリーズ <b>「観測情報_*」9 件</b>
(dataset_id = 1274, 1275, 1276, 1437, 1438, 1439, 1440, 1441, 1442) を統合し、
広島県内における<b>防災観測網</b>の<b>幾何構造</b>を分析する研究記事です。
</div>

<h3>シリーズ構造判定: (b) パターン</h3>
<p>このシリーズの 9 件は、L17/L24 のような「市町別の同一データ N 市町分」
ではなく、<b>「マスタ 1 件 + 5 観測種類 × {{日集計, 年集計}} 8 件」</b>
の 2 層構造です。本記事はマスタを主軸に、各時系列 dataset を構造として併置します。</p>

<h3>独自用語の定義 (本記事内で固定)</h3>
<table>
<tr><th>用語</th><th>定義</th></tr>
<tr><td>観測網</td><td>広島県内に設置された 5 観測種類 (雨量・水位・ダム・潮位・風) の<b>観測点の総体</b>。
本記事では観測所マスタ (dsid=1274) に登録された <b>{n_total} 点</b> を観測網の現在の実体とする。</td></tr>
<tr><td>観測点 (観測所)</td><td>緯度経度を持つ 1 個の物理計測サイト。1 つの観測所が複数の観測種類を持つ場合は、
マスタ上は<b>「データ種別」ごとに別行</b>として登録される (例: 「広島港」が潮位と風で 2 行)。</td></tr>
<tr><td>カバレッジ</td><td>県土の任意点から、最も近い観測点までの距離。本記事では 1.5 km grid を県土全域に張って計測する。</td></tr>
<tr><td>観測空白</td><td>カバレッジが <b>{blank_thresh_km} km 以上</b> の grid 点。本記事ではこの基準で空白を可視化する。
当初 8 km と仮置きしたが、実測で観測網が想定より遥かに密 (最大 10 km / 中央値 2 km) だったため、
<b>「明らかな空白」と呼べる 4 km に下方修正</b>した。仮説の修正自体が観測網理解の重要な一歩。</td></tr>
<tr><td>先行指標</td><td>河川氾濫閾値 (はん濫注意・避難判断・はん濫危険) を保有する水位観測。
これらは「水位がこの値を超えたら警報を出す」運用ノードであり、災害発生に<b>先行</b>する情報を提供する。</td></tr>
<tr><td>所管二元性</td><td>観測網が「県管理 (砂防課・河川課・河川課ダムG・農業基盤課・港湾漁港整備課)」と
「国管理 (国土交通省・国ダム・気象台・海上保安庁)」の<b>2 主体</b>で運用されている構造。</td></tr>
</table>

<h3>研究の問い (RQ)</h3>
<p>広島県内の防災観測網は、観測種類別に何拠点配置され、地理的にどの流域・どの地域を
カバーしているか? その配置パターンから、県の<b>災害監視思想</b>はどう読み取れるか?</p>

<ol>
<li>5 観測種類 (雨量/水位/ダム/潮位/風) は<b>何拠点ずつ、どの所管</b>に配置され、地理的にどう分布するか?</li>
<li>観測点間距離から見た観測網の<b>空間粒度</b>と<b>空白地帯</b>は?</li>
<li>水系・市町別の観測点密度は<b>どれだけ偏っている</b>か?</li>
<li>河川氾濫閾値を保有する観測所は、<b>「先行指標」</b>として機能できる場所にあるか?</li>
</ol>

<h3>仮説 H1〜H6</h3>
<ul>
<li><b>H1 (種別不均衡)</b>: 雨量 > 水位 ≫ 風 ≫ ダム ≈ 潮位、雨量だけで 60% 以上を占める。</li>
<li><b>H2 (所管二元性)</b>: 県管理 vs 国管理 がほぼ折半 (差 40% 未満)。</li>
<li><b>H3 (水系密度)</b>: 太田川・江の川・芦田川の主要 3 水系で観測点 50%以上を集中。</li>
<li><b>H4 (海岸線潮位)</b>: 潮位 13 点は全て瀬戸内沿岸、内陸ゼロ。</li>
<li><b>H5 (氾濫閾値階層)</b>: 水位観測の 50% 以上に氾濫閾値あり。残りは「監視のみ」。</li>
<li><b>H6 (空白地帯)</b>: 雨量で {blank_thresh_km} km 以上離れる観測空白が 5〜25% の grid に出現
(<b>当初 8 km と置いたが、実測で観測網の最大空白が約 10 km、中央値 2 km と極めて密であることが判明したため、
本記事では「明らかな空白」の基準を 4 km に下方修正</b>。これは仮説の更新自体が研究の重要なステップ)。</li>
</ul>

<h3>到達点</h3>
<p>9 dataset_id を「マスタ 1 + 5 種類 × 2 粒度 = 9」の 2 層構造として読み解き、
634 観測所の地理分布・水系分布・閾値構造を統合した結果、
県の観測網が<b>「降水と河川を最優先」「沿岸と山地を補完」</b>という設計思想で
組まれていることを実データで裏付ける。</p>

<div class="warn"><b>X06 との重複回避</b>:
X06 (雨量地理 small multiples) は<b>14 日豪雨期間の雨量時系列値</b>を空間×時間で
描き、特定イベントの降水動態を扱う記事です。本記事 L31 は<b>観測網の構造そのもの</b>
(何拠点・どこ・どの所管・どの閾値) を扱い、時系列値の数値分析は意図的に行いません。
時系列ファイルは<b>列構造の同型性</b>と<b>カバー観測所</b>の確認のみに使用します。</div>
""".replace("{n_total}", str(n_total)).replace("{blank_thresh_km}", str(blank_thresh_km))
sections.append(("学習目標と問い", sec1))


# --- Section 2: 使用データ ---
sec2_table = "<table><tr><th>dsid</th><th>役割</th><th>観測種類</th><th>粒度</th><th>形式</th><th>サイズ(KB)</th><th>DoBoX</th></tr>"
for r in series_df.itertuples():
    sec2_table += (
        f"<tr><td><b>{r.dsid}</b></td><td>{r.role}</td>"
        f"<td>{KIND_LABEL.get(r.kind, '—') if r.kind else '—'}</td>"
        f"<td>{r.granularity or '—'}</td><td>CSV</td><td>{r.size_kb}</td>"
        f"<td><a href='https://hiroshima-dobox.jp/datasets/{r.dsid}' target='_blank'>#{r.dsid}</a></td></tr>"
    )
sec2_table += "</table>"

sec2 = f"""
<p>本記事が使用する 9 dataset_id の一覧と役割。中央 1 件がマスタ、残り 8 件が時系列。</p>

{sec2_table}

<h3>シリーズ構造判定の根拠</h3>
<p>判定: <b>(b) 9 種類の異なる観測データ</b>。理由は以下の表の通り、
9 件は明らかに「マスタ 1 + 5 種類 × 2 粒度」の<b>機能分化型</b>であり、
「9 市町分の同一データ」ではない。</p>

<div class="note">
<b>マスタ 1274</b> が地理キー (緯度経度・水系・河川・市町・氾濫閾値) を持ち、<br>
<b>時系列 8 件</b> はそれぞれ観測所×時刻の wide 表。<b>観測所名で結合可能</b>。
</div>

<h3>マスタ vs 時系列ファイル のカバレッジ</h3>
<p>マスタの観測所と、各時系列ファイルが提供する観測所列を観測所名で照合。
基本的にマスタ⊇時系列だが、潮位・風など一部に非対称あり。詳細は分析 4 で図示。</p>
{cover_df.to_html(index=False, classes='', border=0)}
"""
sections.append(("使用データ", sec2))


# --- Section 3: ダウンロード ---
dl_buttons = []
for r in series_df.itertuples():
    rid = RESOURCE_IDS.get(r.dsid)
    if rid:
        dl_buttons.append(
            f'<li>{r.abbr} (dsid={r.dsid}): '
            f'<a href="https://hiroshima-dobox.jp/resource_download/{rid}">直 DL</a> '
            f'/ <a href="https://hiroshima-dobox.jp/datasets/{r.dsid}">DoBoX</a></li>'
        )
sec3 = f"""
<h3>生データ (DoBoX 直リンク)</h3>
<ul>{''.join(dl_buttons)}</ul>

<h3>中間データ (本記事生成 CSV)</h3>
<ul>
<li>{dl("L31_series_9.csv")} — 9 dataset の役割表</li>
<li>{dl("L31_stations_full.csv")} — 観測所マスタ ({n_total} 行) を kind ラベル付きで整形</li>
<li>{dl("L31_kind_counts.csv")} — 種別 5 区分の件数表</li>
<li>{dl("L31_origin_breakdown.csv")} — 所管別件数</li>
<li>{dl("L31_pivot_watershed.csv")} — 水系 × 種別ピボット</li>
<li>{dl("L31_pivot_city.csv")} — 市町 × 種別ピボット</li>
<li>{dl("L31_attr_fillrate.csv")} — 閾値属性の埋まり率</li>
<li>{dl("L31_master_vs_ts.csv")} — マスタ vs 時系列カバレッジ</li>
<li>{dl("L31_ts_structure.csv")} — 時系列 wide CSV の構造</li>
<li>{dl("L31_nearest_stats.csv")} — 種別ごとの最近隣距離統計</li>
<li>{dl("L31_blank_grid.csv")} — 1.5 km グリッド × 最寄り距離</li>
<li>{dl("L31_hypothesis_results.csv")} — H1〜H6 検証結果</li>
</ul>

<h3>図 (本記事生成 PNG)</h3>
<ul>
<li>{dl("L31_fig1_kind_counts.png")} / {dl("L31_fig2_pref_overview.png")} / {dl("L31_fig3_kind_smallmult.png")}</li>
<li>{dl("L31_fig4_watershed_density.png")} / {dl("L31_fig5_city_density.png")} / {dl("L31_fig6_flood_threshold.png")}</li>
<li>{dl("L31_fig7_blank_zones.png")} / {dl("L31_fig8_nearest_dist.png")} / {dl("L31_fig9_series_matrix.png")}</li>
<li>{dl("L31_fig10_master_ts_coverage.png")} / {dl("L31_fig11_dataset_cards.png")}</li>
</ul>

<h3>再現用 Python スクリプト</h3>
<p><a href="L31_observation_network.py" download>L31_observation_network.py</a> を取得して
プロジェクトルートで <code>py -X utf8 lessons/L31_observation_network.py</code> を実行。
データが無ければ自動取得します。</p>
"""
sections.append(("ダウンロード", sec3))


# --- Section 4: 分析 1 9 dataset の構造を可視化 ---
sec4 = f"""
<h3>狙い</h3>
<p>9 件の dataset がどう機能分化しているのかを 1 枚の絵で示す。
これは「カバー宣言」の構造図であり、後続の全分析の論理的下敷きとなる。</p>

<h3>手法 (簡潔に)</h3>
<p>各 dataset の (役割, 観測種類, 粒度) を抽出し、<b>2 軸マトリクス</b>に
配置する。 役割: master / timeseries。 粒度: master / daily / yearly。
さらにカード形式で件数 + dsid を視覚化する。</p>

<h3>実装</h3>
{code('''
# 9 件: (dsid, 略称, 役割, kind, granularity, file)
SERIES = [
    (1274, "観測所一覧",      "master",     None,    None,     "stations_master.csv"),
    (1275, "雨量_日集計",     "timeseries", "rain",  "daily",  "rain_daily_sample.csv"),
    (1276, "雨量_年集計",     "timeseries", "rain",  "yearly", "rain_year_sample.csv"),
    (1437, "水位_日集計",     "timeseries", "water", "daily",  "water_daily_sample.csv"),
    (1438, "水位_年集計",     "timeseries", "water", "yearly", "water_year_sample.csv"),
    (1439, "ダム_日集計",     "timeseries", "dam",   "daily",  "dam_daily_sample.csv"),
    (1440, "ダム_年集計",     "timeseries", "dam",   "yearly", "dam_year_sample.csv"),
    (1441, "潮位_日集計",     "timeseries", "tide",  "daily",  "tide_daily_sample.csv"),
    (1442, "風向風速_日集計", "timeseries", "wind",  "daily",  "wind_daily_sample.csv"),
]
# 役割 × 粒度 のマトリクスに配置
matrix = pd.DataFrame(index=KIND_ORDER + ["master"], columns=["daily", "yearly", "master"])
for r in series_df.itertuples():
    if r.role == "master":
        matrix.loc["master", "master"] = r.dsid
    else:
        matrix.loc[r.kind, r.granularity] = r.dsid
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 9 つの dataset_id を「マスタと時系列」「日と年」のどちらかに分類することで、
シリーズの本来の構造 (機能分化型) が一望できる。表のままだと粒度の有無が読み取りにくい。</p>
{figure("assets/L31_fig9_series_matrix.png", "9 dataset_id の構造マトリクス (役割 × 粒度)")}
{figure("assets/L31_fig11_dataset_cards.png", "9 dataset カードビュー — 1 マスタ + 5 種類 × {日, 年}")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>マスタ 1 件 (1274)</b> が観測所一覧として地理キーを保有。緯度経度・水系・河川・市町・氾濫閾値を含む 35 列。</li>
<li><b>時系列 8 件</b> は観測種類 (5 種) × 粒度 (日, 年) のうち、潮位と風の年集計のみが欠けている (= 6 種類しか日集計がない訳ではなく、潮位と風は<b>日のみ</b>で運用)。</li>
<li>このシリーズ全体は「観測網の地理を 1 枚に集約し、各観測種類の値時系列を別 dataset で提供する」という設計。<b>地理結合キーは観測所名</b>。</li>
</ul>

<h3>表と読み取り</h3>
{series_df.to_html(index=False, classes='', border=0)}
<p><b>読み取り</b>:</p>
<ul>
<li>サイズ最大はマスタ次いで雨量日 (1.7 MB) — これは雨量の観測所数が最大かつ 10 分粒度で粒度が高いため。</li>
<li>潮位日 (~120 KB) と風日 (~210 KB) は雨量日の<b>1/10 規模</b>。観測所数が少なく粒度も粗い。</li>
<li>年集計は日集計の 1/3〜1/5 サイズ。1 日 1 値に集約されている。</li>
</ul>
"""
sections.append(("分析 1: 9 dataset の構造を可視化", sec4))


# --- Section 5: 分析 2 観測種類別の件数と所管 ---
sec5_kind_table = kc.to_html(index=False, classes='', border=0)
sec5 = f"""
<h3>狙い</h3>
<p>5 観測種類が<b>何点ずつ</b>配置され、<b>誰が所管</b>するかの全体像を確定する。
これが H1 (種別不均衡) と H2 (所管二元性) の検証の中核。</p>

<h3>手法</h3>
<p>マスタ ({n_total} 行) を「データ種別」(1=雨量/4=水位/7=ダム/12=潮位/13=風) と
「データ所管」 (砂防課・国土交通省など 9 種) でクロス集計。
所管は「県/国/その他」の 3 グループに集約する (定義は冒頭表を参照)。</p>

<h3>実装</h3>
{code('''
# データ種別 → 観測種類名
KIND_ID2NAME = {1: "rain", 4: "water", 7: "dam", 12: "tide", 13: "wind"}
master["kind"] = master["データ種別"].map(KIND_ID2NAME)

# 所管グループ
ORIGIN_GROUP = {"砂防課":"県", "河川課":"県", "国土交通省":"国", ...}
master["origin"] = master["データ所管"].map(ORIGIN_GROUP).fillna("その他")

# クロス集計
pv_kind_origin = pd.pivot_table(
    master, index="kind", columns="origin", values="観測所名",
    aggfunc="count", fill_value=0)
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 棒グラフ + スタックは「絶対数 + 構成比」を 1 セクションで示せる。
件数バー単独では所管の二元性が見えず、所管比単独では絶対数が見えない。</p>
{figure("assets/L31_fig1_kind_counts.png", "観測種類別件数 (左) と種別×所管スタック (右)")}
<p><b>読み取り</b>:</p>
<ul>
<li>雨量が最多 ({n_rain} 点 / {rain_pct:.1f}%) で、水位 ({n_water} 点 / {n_water/n_total*100:.1f}%) と合わせて<b>全体の {(n_rain+n_water)/n_total*100:.1f}%</b>。
水文 (降水と河川) が観測網の主軸。</li>
<li>ダム ({int(kind_counts['dam'])} 点) と潮位 ({int(kind_counts['tide'])} 点) は最少。これは<b>物理的設置可能箇所が限定的</b>な観測種類。</li>
<li>所管スタックを見ると、ダムは<b>河川課ダムG/国ダム</b>でほぼ全件、潮位は<b>港湾漁港整備課</b>、風は<b>気象台</b>がほぼ独占。
所管の専門化が強い。</li>
<li>逆に雨量と水位は所管が分散しており、「砂防課・河川課・国土交通省・気象台」が分担。</li>
</ul>

<h3>表と読み取り</h3>
{sec5_kind_table}
{oc.to_html(index=False, classes='', border=0)}
<p><b>読み取り</b>:</p>
<ul>
<li>砂防課 132 件は<b>全件雨量</b> (砂防警戒のための雨量計)。県の防災運用上、土砂災害監視の中核。</li>
<li>河川課 137 件は雨量 + 水位の混合。河川管理に直結する観測ノード。</li>
<li>気象台 53 件は雨量 + 風が中心。広域気象監視のサポート観測網。</li>
</ul>
"""
sections.append(("分析 2: 観測種類別の件数と所管", sec5))


# --- Section 6: 分析 3 県全域マップで観測網を観る ---
sec6 = f"""
<h3>狙い</h3>
<p>件数だけでは観測網の<b>地理偏在</b>は見えない。同じ 401 点の雨量でも、
都市集中型と均等分散型では運用思想が違う。マップで 1 枚に重ねた俯瞰と、
種別ごとに分離した small multiples の 2 段で確認する。</p>

<h3>手法</h3>
<p>マスタの緯度経度を Shapely Point に変換し、EPSG:4326 → EPSG:6671 (平面直角 III) に投影。
広島県全域 polygon (L15 の dsid=922 を流用) を背景に重ねる。
種別ごとに色分けし、small multiples で 5 種を別パネルにも展開する。</p>

<h3>実装</h3>
{code('''
master_gdf = gpd.GeoDataFrame(
    master, geometry=gpd.points_from_xy(master["経度"], master["緯度"]),
    crs="EPSG:4326").to_crs("EPSG:6671")

g_admin_pref = load_zip_first_geo(ADMIN_DIR / "admin_922_広島県.zip").to_crs("EPSG:6671")
pref_diss = g_admin_pref.dissolve()  # 県境 1 ポリゴンに集約

fig, ax = plt.subplots(figsize=(11, 8))
pref_diss.boundary.plot(ax=ax, color="#222", linewidth=1.2)
for k in KIND_ORDER:
    sub = master_gdf[master_gdf["kind"] == k]
    sub.plot(ax=ax, color=KIND_COLOR[k], markersize=18, alpha=0.78)
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 観測網は地理的実体なので、地図化が最も直感的。
重ね地図で全体感、small multiples で個別パターンを見る。</p>
{figure("assets/L31_fig2_pref_overview.png", "県全域 観測網マップ (5 種カラー散布)")}
{figure("assets/L31_fig3_kind_smallmult.png", "種別 small multiples マップ — 種別ごとの偏在")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>雨量 (青)</b> は県全域に均等分布。山地・離島まで及ぶ。砂防運用が県全土にスパンしていることの帰結。</li>
<li><b>水位 (赤)</b> は<b>河川沿いに線状</b>に配置。江の川 (北部)、太田川 (中部)、芦田川 (東部) の幹線が読める。</li>
<li><b>ダム (茶)</b> は<b>山間部</b>に点在。中山間ダムと、瀬戸内沿岸の小規模ダムの 2 グループ。</li>
<li><b>潮位 (緑)</b> は<b>瀬戸内海岸線のみ</b>。内陸ゼロ、H4 を支持する。</li>
<li><b>風 (黄)</b> は気象台立地に依存し、市街地と一部島嶼に偏る。</li>
</ul>
"""
sections.append(("分析 3: 県全域マップで観測網を観る", sec6))


# --- Section 7: 分析 4 水系・市町別の密度と所管二元性 ---
sec7 = f"""
<h3>狙い</h3>
<p>「どこにどれだけ観測点があるか」を水系・市町の 2 視点で集計し、
H3 (主要 3 水系集中) を検証する。同時に観測網の<b>偏り</b>を可視化する。</p>

<h3>手法</h3>
<p>マスタを水系名・市町村名でグループ化し、種別を列としてピボット。
上位 12 水系 / 上位 15 市町を横棒スタックで比較。</p>

<h3>実装</h3>
{code('''
pv_water = pd.pivot_table(
    master_gdf, index="水系名", columns="kind",
    values="観測所名", aggfunc="count", fill_value=0)
pv_water["合計"] = pv_water.sum(axis=1)
pv_water = pv_water.sort_values("合計", ascending=False)
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 水系密度と市町密度は、観測網が<b>「自然単位 (水系) の地理」</b>と
<b>「行政単位 (市町) の地理」</b>のどちらに沿って配置されているかを示す。
両方を並べることで「水系優先 vs 行政優先」の軸を読み取れる。</p>
{figure("assets/L31_fig4_watershed_density.png", "水系別観測網密度 (上位 12 水系)")}
<p><b>読み取り (水系)</b>:</p>
<ul>
<li>上位 3 水系 ({', '.join(top3_names)}) で<b>{top3_pct_v:.1f}%</b> ({top3_total_v}/{n_total}) を占有。H3 を {'支持' if top3_pct_v >= 50 else '部分支持'}。</li>
<li>太田川は流域最大の県管理河川で、雨量 + 水位の双方が密。「降雨 + 流出」を組で観測する設計が読める。</li>
<li>江の川は中山間部の長い水系で、観測点が長距離に分散。1 流域あたりの観測網延長が大きい。</li>
<li>芦田川 (福山) は人口集中地下流のため、水位観測の比率が他より高い。</li>
</ul>
{figure("assets/L31_fig5_city_density.png", "市町別観測網密度 (上位 15 市町)")}
<p><b>読み取り (市町)</b>:</p>
<ul>
<li>福山市 (58 点) と三次市 (43 点) が突出。それぞれ「人口集中地」と「中山間水系の交差点」という対照的理由で観測密。</li>
<li>政令市の広島市は<b>区別に分割</b>されているため上位ランクには現れにくいが、安佐北区 33 点・佐伯区 25 点・東区などを足すと<b>合計 130 点超</b>。実質トップ。</li>
<li>福山以外の沿岸都市 (呉・尾道・三原) は雨量 + 潮位混在のため種別構成が他都市と異なる。</li>
</ul>

<h3>表と読み取り (水系 上位 10)</h3>
{pv_water.head(10).to_html(classes='', border=0)}
<p><b>読み取り</b>: 太田川・江の川・芦田川の 3 水系の<b>「雨量:水位 比」</b>を比較する (下表)。
H3 で予想した「雨量:水位 ≈ 2:1〜3:1」を概ね支持する。</p>
{(lambda: pd.DataFrame([
    {"水系": ws,
     "雨量": int(pv_water.loc[ws, "rain"]) if ws in pv_water.index else 0,
     "水位": int(pv_water.loc[ws, "water"]) if ws in pv_water.index else 0,
     "雨量:水位": (round(pv_water.loc[ws, "rain"] / max(1, pv_water.loc[ws, "water"]), 2)
                   if ws in pv_water.index else None)}
    for ws in ["太田川", "江の川", "芦田川"] if ws in pv_water.index
]).to_html(index=False, classes='', border=0))()}

<h3>表と読み取り (市町 上位 10)</h3>
{pv_city.head(10).to_html(classes='', border=0)}
<p><b>読み取り</b>: 都市規模 (人口) と観測点数の単純な比例関係はない。
中山間市町 (北広島町 28 点、世羅町 23 点、安芸太田町 21 点) も上位に入る。
これは観測網が<b>「災害リスクの空間 (山地・水系・沿岸)」</b>に沿って配置され、
人口密度ではないことを示す。</p>
"""
sections.append(("分析 4: 水系・市町別の密度", sec7))


# --- Section 8: 分析 5 氾濫閾値で見る「先行指標」 ---
sec8 = f"""
<h3>狙い</h3>
<p>水位観測 181 点のうち、どれが<b>「警報を出すノード」</b>として機能しうるか?
氾濫閾値 (水防団待機・はん濫注意・避難判断・はん濫危険) の有無と段階数で評価し、
H5 (50% 以上) を検証する。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p><b>氾濫閾値</b> とは、河川水位が「この値を超えたら住民に避難を呼びかける/警戒態勢に入る」
という運用基準値です。広島県の水位観測ではマスタ CSV に最大 4 段階の閾値が記録されています。
段階数が多いほど<b>細かく警報レベルを分けられる</b>観測ノードです。</p>
<p>段階数 = 0 の観測所は<b>監視のみ</b> (値を取るだけで警報出力はしない、上流域の参考点など)。</p>

<h3>実装</h3>
{code('''
thresh_cols = ["水防団待機[m]", "はん濫注意[m]", "避難判断[m]", "はん濫危険[m]"]
water_gdf["n_thresh"] = water_gdf[thresh_cols].notna().sum(axis=1)
water_gdf["has_flood_thresh"] = water_gdf["n_thresh"] > 0
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 「閾値あり/なし」の 2 値マップで地理偏在、段階数ヒストグラムで運用設計の細かさを 1 枚で見る。</p>
{figure("assets/L31_fig6_flood_threshold.png", "水位観測の氾濫閾値: 有無マップ (左) と段階数分布 (右)")}
<p><b>読み取り</b>:</p>
<ul>
<li>水位 {n_water_total} 点のうち <b>{n_water_with} 点 ({h5_pct:.1f}%)</b> が氾濫閾値あり。H5 を {'支持' if h5 else '部分支持'}。</li>
<li>閾値あり (赤) は<b>主要河川の中下流</b>に集中、閾値なし (灰色) は<b>上流支流</b>に多い。
これは「警報出力ノード = 中下流」「監視のみ = 上流」という<b>役割分担</b>を反映。</li>
<li>段階数別に見ると、<b>4 段階フル装備</b>の観測所は限定的で、<b>3 段階</b>がやや多い。
これは「水防団待機」の運用が一部河川でのみ採用されていることを示唆する。</li>
</ul>

<h3>表と読み取り</h3>
{attr_df.to_html(index=False, classes='', border=0)}
<p><b>読み取り</b>:</p>
<ul>
<li>河川氾濫閾値は<b>水位観測 (kind=water) のみ</b>に存在し、雨量・ダム・潮位・風には付かない。閾値の意味論が観測種類と一意に対応している。</li>
<li>ダム水位の 5 系列の閾値は<b>ダム観測 (kind=dam) のみ</b>。最も多いのは「常時満水位」と「サーチャージ水位」で、それぞれダム運用の上下限を示す。</li>
<li>潮位閾値 (注意・警戒) は<b>潮位観測 (kind=tide) のみ</b>に 100% 装備。沿岸部の高潮警戒運用を反映。</li>
</ul>
"""
sections.append(("分析 5: 氾濫閾値で見る「先行指標」", sec8))


# --- Section 9: 分析 6 観測網の空間粒度と空白地帯 ---
near_any = near_df[near_df["kind"] == "any"].iloc[0]
near_rain = near_df[near_df["kind"] == "rain"].iloc[0]
near_water = near_df[near_df["kind"] == "water"].iloc[0]
sec9 = f"""
<h3>狙い</h3>
<p>観測網の<b>空間粒度</b> (どれだけ密に張られているか) を測定し、観測<b>空白地帯</b>を可視化する。
H6 (8 km 以上の空白が grid 5〜25%) の検証。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p>観測網の粒度を測る最もシンプルな指標が <b>「最近隣距離」</b> です。<br>
ある観測点から、自分以外で<b>最も近い観測点</b>までの直線距離。これが小さいほど密、大きいほど疎。</p>
<ul>
<li><b>k-NN (k Nearest Neighbors)</b>: 全点に対して「最も近い 1 点を探す」操作を高速に行うアルゴリズム。
<code>scipy.spatial.cKDTree</code> を使うと {n_total} 点でも 1 秒以内に終わります。</li>
<li><b>入力</b>: 全観測点の (x, y) 座標 (m 単位、平面直角座標 EPSG:6671)。</li>
<li><b>出力</b>: 各点について最近隣点との距離 (m)。</li>
<li><b>限界</b>: 直線距離なので、山岳・河川による「実質距離」は反映しません。教育的目的にはこれで十分。</li>
</ul>

<p>空白地帯マップは、県土に <b>1.5 km grid</b> を張り、各 grid 点について最寄り観測点までの距離を計算します。
赤いほど空白に近い。</p>

<h3>実装</h3>
{code('''
from scipy.spatial import cKDTree
xy = np.array([(p.x, p.y) for p in master_gdf.geometry])
tree = cKDTree(xy)
dists, _ = tree.query(xy, k=2)   # k=2: 自分(=0) + 最近隣(=1)
master_gdf["nearest_any_m"] = dists[:, 1]

# 観測空白 grid (1.5 km 間隔で県土を埋める)
xs = np.arange(xmin, xmax, 1500)
ys = np.arange(ymin, ymax, 1500)
xx, yy = np.meshgrid(xs, ys)
grid_pts = np.column_stack([xx.ravel(), yy.ravel()])
in_pref = gpd.GeoSeries([Point(x,y) for x,y in grid_pts]).within(pref_geom)

# 各 grid 点 → 最寄り雨量観測まで
rain_xy = np.array([(p.x, p.y) for p in rain_gdf.geometry])
rain_tree = cKDTree(rain_xy)
gd, _ = rain_tree.query(grid_pts[in_pref], k=1)
''')}

<h3>図と読み取り — 観測点間距離分布</h3>
<p><b>なぜこの図か</b>: ヒストグラム + 中央値線で「典型的な観測点間距離」を 1 枚に集約。
種別ごとに分けることで、種別の運用密度の違いを比較できる。</p>
{figure("assets/L31_fig8_nearest_dist.png", "種別ごとの最近隣距離ヒストグラム + 全観測点比較")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>全観測点</b> (種別不問) の最近隣中央値は <b>{near_any['median_km']} km</b>。県土全体で平均 {near_any['median_km']} km おきに何かしらの観測点が立っている。</li>
<li><b>雨量</b>のみで見ると中央値 <b>{near_rain['median_km']} km</b>。最も密。</li>
<li><b>水位</b>のみだと中央値 <b>{near_water['median_km']} km</b>。雨量より粗いが、河川沿いに線状配置されているため局所的には密。</li>
<li><b>潮位・ダム・風</b>はサンプル数が少なく粒度も粗い (中央値 5〜15 km 級)。</li>
</ul>

<h3>図と読み取り — 観測空白地帯</h3>
<p><b>なぜこの図か</b>: 「グリッドが赤いほど空白」のヒートマップは、
表や統計値より<b>直感的に空白地</b>を示せる。</p>
{figure("assets/L31_fig7_blank_zones.png", "観測空白地帯マップ — 雨量 (左) と全種別 (右) の最寄り距離 grid")}
<p><b>読み取り</b>:</p>
<ul>
<li>雨量で {blank_thresh_km} km 以上離れる「明らかな空白」 grid は <b>{n_blank}/{n_in} = {n_blank_pct:.2f}%</b>。H6 を {hypos[5]['verdict']}。</li>
<li>当初想定していた 8 km 基準では空白は<b>{n_blank_8_pct:.4f}%</b> (ほぼゼロ)。これは「広島県の雨量観測網が想定より遥かに密」という重要な発見。</li>
<li>赤い grid は<b>島嶼部</b> (江田島・大崎上島・倉橋島など) と<b>北部中山間境界</b> (島根県境近く) に集中。</li>
<li>沿岸の市街地・主要河川流域は密で、ほぼ全 grid が 4 km 未満。</li>
<li>全種別 (右) で見ると空白はさらに小さくなり、雨量+他種別の<b>補完効果</b>がわかる。</li>
</ul>

<h3>表と読み取り (最近隣統計)</h3>
{near_df.to_html(index=False, classes='', border=0)}
<p><b>読み取り</b>:</p>
<ul>
<li><b>潮位</b>のみが中央値 ~10 km 超、p95 で {near_df[near_df['kind']=='tide']['p95_km'].iloc[0] if (near_df['kind']=='tide').any() else '—'} km と最大。
これは沿岸線に沿った「点の数珠つなぎ」配置のため、線方向に隣接距離が伸びる。</li>
<li><b>ダム</b>も粒度が粗い ({near_df[near_df['kind']=='dam']['median_km'].iloc[0] if (near_df['kind']=='dam').any() else '—'} km)。物理的に「河川の遮断点」にしか置けない構造的限界の表現。</li>
</ul>
"""
sections.append(("分析 6: 観測網の空間粒度と空白地帯", sec9))


# --- Section 10: 分析 7 マスタ vs 時系列 ---
sec10 = f"""
<h3>狙い</h3>
<p>マスタ 1274 と時系列 8 件は同じ観測所を扱っているはず。
しかし時系列ファイルの観測所列がマスタと完全一致するかは未検証。
本分析では観測所名で照合し、<b>マスタ⊇時系列</b> 関係が成立するかを確認する。</p>

<h3>手法</h3>
<p>各時系列 wide CSV のヘッダ部 (上 6 行) から観測所名を抽出。
マスタの観測所名集合と差分を計算する。
時系列ヘッダは「事務所名・データ所管・水系名・河川名・観測所番号・観測所名」が縦に並び、
データ列が観測所単位で 1〜複数本配置される。</p>

<h3>実装</h3>
{code('''
def parse_ts_header(path):
    """時系列 wide CSV の上部メタを抽出"""
    raw = pd.read_csv(path, encoding="utf-8-sig", nrows=10, header=None)
    # 行 0-5 を観測所メタとして縦持ちに
    meta = pd.DataFrame({
        "事務所名":   raw.iloc[0, 1:],
        "データ所管": raw.iloc[1, 1:],
        "水系名":     raw.iloc[2, 1:],
        "河川名":     raw.iloc[3, 1:],
        "観測所番号": raw.iloc[4, 1:],
        "観測所名":   raw.iloc[5, 1:],
    }).dropna(subset=["観測所名"]).drop_duplicates(subset=["観測所名"])
    return meta

# マスタ vs 時系列 を観測所名集合で照合
master_names = set(master_gdf.loc[master_gdf["kind"]==k, "観測所名"])
ts_names = set(ts_meta_per_kind[k]["観測所名"])
n_common = len(master_names & ts_names)
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: バーペアで「マスタ件数と時系列件数」を直接比較する。
共通件数 (=結合可能件数) を上に表示することで結合運用上の安全度が読める。</p>
{figure("assets/L31_fig10_master_ts_coverage.png", "マスタ vs 時系列 観測所カバレッジ照合")}
<p><b>読み取り</b>:</p>
<ul>
<li>多くの種別でマスタ件数 ≧ 時系列件数。マスタは「全観測所カタログ」、時系列は「実際にデータが取れた観測所」と解釈できる。</li>
<li>共通件数 (=観測所名で結合できる件数) は ts 件数の <b>大半</b>。
名前のゆらぎ (例: 「○○ダム」と「○○ダム(国)」) が一部にあるため、完全一致しない場合もある。</li>
<li>潮位など小規模な種別は完全に対応する傾向。</li>
</ul>

<h3>表と読み取り</h3>
{cover_df.to_html(index=False, classes='', border=0)}
<p><b>読み取り</b>:</p>
<ul>
<li><b>n_only_master &gt; 0</b> = マスタにあるが時系列ファイルには列がない観測所。これは<b>新設観測所で時系列がまだ蓄積されていない</b>か、ファイル更新のタイムラグの可能性。</li>
<li><b>n_only_ts &gt; 0</b> = 時系列ファイルにあるがマスタには無い観測所。これは<b>名前表記揺れ</b> (ふりがな違い、(国)/(気) サフィックス揺れ) で結合に失敗した可能性が高い。実運用では観測所番号の方が安全。</li>
</ul>

<div class="note">
<b>結合運用上の教訓</b>: 観測所名は半角・全角・サフィックスの揺れで完全一致しないことがある。
9 dataset を結合するときは<b>「観測所番号 + データ所管」</b>を主キーにする方が頑健。
ただし観測所番号の体系も種別ごとに別ナンバリングなので、最終的にはマスタ 1274 の<b>(観測所番号, データ所管)</b> 複合キーが最も信頼できる。
</div>
"""
sections.append(("分析 7: マスタ vs 時系列ファイル照合", sec10))


# --- Section 11: 仮説検証 ---
sec11 = f"""
<p>H1〜H6 の検証結果を表で示す。</p>
{hypos_df.to_html(index=False, classes='', border=0)}

<h3>総括: 県の災害監視思想</h3>
<p>9 dataset から再構成した観測網の構造から、以下の<b>3 つの設計思想</b>が読み取れる。</p>
<ul>
<li><b>(1) 降水と河川を最優先</b>: 雨量 ({n_rain}) + 水位 ({n_water}) で全 {n_total} 点の {(n_rain+n_water)/n_total*100:.1f}% を占める。
県の最大の災害リスクである<b>豪雨と土石流・洪水</b>を最優先で監視する設計。</li>
<li><b>(2) 県と国の二層管理</b>: 砂防課 (県) と国土交通省 (国) の所管がほぼ折半。
これは「県管理河川 + 国管理河川」の管理境界に観測責任が完全に対応している証拠。</li>
<li><b>(3) 沿岸と山地は補完的</b>: 潮位 13 + ダム 18 + 風 21 = 52 点 (8.2%) は<b>少数だが特異性が高い</b>観測種類。
沿岸高潮 (潮位)、ダム放流 (ダム)、強風被害 (風) という<b>限定的だが甚大な災害</b>に対する補助観測網。</li>
</ul>

<p>本記事は<b>「観測網は防災の感覚器官」</b>という視点を実データで裏付けた。
{n_total} 点の観測所が、5 種類の<b>感覚モダリティ</b> (雨を測る、水位を測る、ダム水位を測る、潮位を測る、風を測る) を、
2 主体 (県・国) で分担しつつ、主要 3 水系に集中配置している。
この構造を読み取ることが、災害監視ネットワークを設計する者にとっての<b>最初のリテラシ</b>である。</p>
"""
sections.append(("仮説検証と考察", sec11))


# --- Section 12: 発展課題 ---
sec12 = f"""
<h3>結果 X1 → 新仮説 Y1 → 課題 Z1: 観測空白の物理アクセス困難性</h3>
<ul>
<li><b>結果 X1</b>: 雨量で 8 km 以上の観測空白が grid {n_blank_pct:.1f}% に存在。島嶼部と中山間境界に集中。</li>
<li><b>新仮説 Y1</b>: これらの空白地帯は<b>住民居住地から遠い (人口希薄)</b>か、<b>物理アクセス困難 (急峻地・離島)</b>のいずれかが原因。
住民居住地に近い空白なら危険、遠ければ運用上は許容範囲。</li>
<li><b>課題 Z1</b>: 国勢調査メッシュ (1km) と空白 grid を重ねて、空白 grid 内の人口を集計する。
人口あり空白 grid を「優先補強候補」として地図化し、市町別に件数を出す。</li>
</ul>

<h3>結果 X2 → 新仮説 Y2 → 課題 Z2: 氾濫閾値の高さと地形の関係</h3>
<ul>
<li><b>結果 X2</b>: 水位 {n_water_total} 中 {n_water_with} 点 ({h5_pct:.1f}%) に氾濫閾値あり。中下流に集中。</li>
<li><b>新仮説 Y2</b>: <b>「はん濫危険水位 - はん濫注意水位」 (=猶予幅)</b> は河川勾配と相関する。
急流河川では猶予幅が狭く (洪水到達が速い)、緩流河川では広い (避難余裕がある)。</li>
<li><b>課題 Z2</b>: マスタの 4 段階閾値の差分を計算し、各観測点で 5 m DEM (LiDAR、L14 で取得済) から
上流流域の平均勾配を計算。閾値猶予幅 vs 流域勾配 の散布図を作成。傾向があれば、
「急流河川では避難呼び掛けタイミングを早める必要がある」という運用提言。</li>
</ul>

<h3>結果 X3 → 新仮説 Y3 → 課題 Z3: 雨量と水位の時系列遅延</h3>
<ul>
<li><b>結果 X3</b>: 太田川水系では雨量 {int(pv_water.loc['太田川','rain'])} 点と水位 {int(pv_water.loc['太田川','water'])} 点が同水系に共存。</li>
<li><b>新仮説 Y3</b>: 上流雨量と下流水位の時系列に<b>「雨が降ってから水位が上がるまでのラグ」</b>が観測できる。
ラグ時間は流域大きさに比例し、太田川では数時間オーダー。</li>
<li><b>課題 Z3</b>: 雨量日 (1275) と水位日 (1437) を実時系列で結合し、
代表的な大雨イベント (例: 2018 年西日本豪雨) における雨量ピーク → 水位ピークの<b>クロス相関</b>を計算する。
観測ノードの「先行警戒能力」を実測。<b>この時系列分析は X06 と接続可能</b>であり、
本記事 (構造) と X06 (時系列) の橋渡しになる。</li>
</ul>

<h3>結果 X4 → 新仮説 Y4 → 課題 Z4: 観測網の冗長性と耐故障性</h3>
<ul>
<li><b>結果 X4</b>: 主要 3 水系に観測点 {top3_pct_v:.1f}% が集中。</li>
<li><b>新仮説 Y4</b>: 集中度が高いと「1 観測所が故障しても近隣で代替できる」が、空白地帯では代替不可。
観測網の<b>耐故障性</b>は局所的に大きく異なる。</li>
<li><b>課題 Z4</b>: 各観測所を 1 つずつ「停止」させた仮想シナリオを {n_total} 通り作り、
停止後の<b>新たな最近隣距離</b>を計算。距離増加が大きい観測所 = 停止インパクト大 = 「観測網のキー観測所」を抽出する。</li>
</ul>
"""
sections.append(("発展課題", sec12))


html = render_lesson(
    num=31,
    title="観測情報 9 件統合分析 — 県の防災観測網 (rain/water/dam/tide/wind) の幾何構造",
    tags=["観測", "GIS", "ネットワーク", "geopandas", "k-NN", "防災"],
    time="40 分",
    level="リテラシ",
    data_label="DoBoX 観測情報 9 dataset (1274, 1275, 1276, 1437-1442)",
    sections=sections,
    script_filename="L31_observation_network.py",
)

(LESSONS / "L31_observation_network.html").write_text(html, encoding="utf-8")
print(f"  saved L31_observation_network.html ({time.time()-t1:.1f}s)", flush=True)


print(f"\n=== 完了: 全所要 {time.time()-t0:.1f}s ===", flush=True)