"""L28 都市機能誘導区域 9 件 統合分析
       — 広島県内 8 LIP 市の「拠点配置」と L19 居住誘導との入れ子構造の解明

カバー宣言:
  本記事は DoBoX のシリーズ
    都市計画区域情報_区域データ_*_都市機能誘導区域 (LIP 策定 8 市 + 県全域 = 9 件)
  を統合し、広島県の<b>「コンパクトシティ＋ネットワーク」の核 = 都市機能誘導区域</b>
  の幾何学的実装を分析する研究記事である。

  9 dataset_id (L19 居住誘導 dsid と +1 ずれの隣接 dsid):
    広島市 (796), 呉市 (806), 竹原市 (813), 三原市 (823), 福山市 (839),
    府中市 (849), 東広島市 (877), 廿日市市 (887), 広島県全域 (934)

研究の問い (RQ):
  L19 居住誘導区域シリーズ (jukyo_*) と L28 都市機能誘導区域シリーズ (toshikinou_*) は、
  各市町で <b>+1 ずれの隣接 dataset_id</b> として別々に提供されている。
  両シリーズはどう異なるのか？
  L28「都市機能誘導区域」シリーズは、
  名前通り <b>都市機能誘導区域 (KUIKI=3) のみ</b>を含むのか、
  それとも L19 と類似の混合データなのか？
  そして 8 LIP 市の<b>都市機能誘導区域</b>は、
  立地適正化計画の理論モデル「<b>居住誘導 ⊃ 都市機能誘導</b>の二重入れ子」を
  幾何学的に満たしているか？

仮説 H1〜H5:
  H1 (シリーズ実態): L28 シリーズはファイル名と裏腹に、
       中身は <b>L19 居住誘導 GeoJSON のサブセット</b>である。
       具体的には、L19 jukyo の KUIKI_CD=1 (居住誘導) ポリゴンの中で、
       <b>都市機能誘導区域 (KUIKI=3) に近接する部分</b>だけを抜き出した
       <b>狭い切り出し版</b>になっている。
       (期待: 8 市平均で L28 KUIKI=1 面積 / L19 KUIKI=1 面積 = 30〜50%)
  H2 (包含関係): L28 KUIKI=1 ポリゴンは <b>L19 KUIKI=1 ポリゴンに完全に含まれる</b>。
       両 GeoJSON で同じ KUIKI=1 ポリゴンの差分は <b>0.1 km² 未満</b>であり、
       L28 が L19 のサブセットであることが幾何的に証明される。
  H3 (拠点ネットワーク型の検証 — 反証されると重要発見): L28 KUIKI=1 (近接居住誘導) と
       L28 KUIKI=3 (都市機能誘導) は<b>幾何的に重複する</b>と理論モデルは予測。
       検証指標: 拠点 → 居住誘導 最短距離 < 1 km、包含率 > 50% を期待。
       <b>もし反証されれば</b>「広島県の LIP は<b>分離型</b> ─ 居住誘導域と
       都市機能誘導が km 単位で離れて配置されている」という重要発見になる。
       L19 居住誘導全体 (jukyo_*) との距離も計測し、より広い範囲との関係も検証する。
  H4 (8 市の都市機能誘導の規模差): 都市機能誘導区域 (KUIKI=3) の総面積は、
       8 市で <b>1〜2 km²/市</b> の小規模な拠点として指定されている。
       人口の多い市ほど多数の拠点を持つ ─ <b>多核拠点ネットワーク型</b>。
       (期待: 政令市・中核市は 2 拠点以上、人口 50万人以下の市は 1 拠点)
  H5 (整合性): 8 市別 8 ファイルの合計と県全域 ds=934 (1150 ポリゴン) は
       <b>件数・面積ともに完全一致</b>する。CITY_CD 集合も同一。

要件 S 準拠: 1 分以内完走 (1150 ポリゴン × 各市少数なので極めて軽量)。
要件 T 準拠: 主題図、市町別 small multiples、都市機能誘導 vs 居住誘導
              重ね地図、拠点間距離 choropleth など複数 GIS 図。
要件 Q 準拠: 図 9 枚以上、表 9 枚以上。

データ仕様 (6 列, L19 と完全同型):
  FID         int     ポリゴン番号
  TOKEI_CD    int     統計区分 (全 1)
  CITY_CD     int     市区町村コード
  KUIKI_CD    int     <b>5 値 (0/1/2/3/4)</b> — L19 と同じ意味体系
  RITTEKI_CD  int     立地適正化フラグ (全 1)
  geometry    Polygon EPSG:6671

L19 との関係 (本記事の主結論):
  - L19 シリーズ = 居住誘導区域に焦点 (jukyo_*.zip)
  - L28 シリーズ = 都市機能誘導区域に焦点 (toshikinou_*.zip)
  - 両シリーズは同じテーブル構造を持ち、KUIKI_CD で類型を識別
  - 実データ確認の結果、<b>L28 は L19 居住誘導の都市機能誘導近接サブセット</b>
    (廿日市市で L28 KUIKI=1 ⊆ L19 KUIKI=1 が 99.95% 包含 で実証済)

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

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

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
import geopandas as gpd

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

t0 = time.time()
print("=== L28 都市機能誘導区域 9 件 統合分析 ===", flush=True)

# =============================================================================
# 0. 定数: 9 dataset_id (L19 居住誘導と +1 ずれの隣接 dsid)
# =============================================================================
DATA_DIR = ROOT / "data" / "extras" / "L28_urban_function_induction"
DATA_DIR_L19 = ROOT / "data" / "extras" / "L19_residential_induction"
ADMIN_DIR = ROOT / "data" / "extras" / "L15_admin_zones"
TARGET_CRS = "EPSG:6671"  # JGD2011 平面直角 III, 広島県, m 単位

# (L28 dsid, L19 dsid, 市町名, タイプ, 行政_dsid)
LIP_CITY_DEFS = [
    (796, 795, "広島市",   "政令市",       786),
    (806, 805, "呉市",     "中核市",       797),
    (813, 812, "竹原市",   "市",           807),
    (823, 822, "三原市",   "市",           814),
    (839, 838, "福山市",   "中核市",       832),
    (849, 848, "府中市",   "市",           840),
    (877, 876, "東広島市", "施行時特例市", 868),
    (887, 886, "廿日市市", "市",           878),
]
KEN_DSID_L28 = 934
KEN_DSID_L19 = 933

# KUIKI_CD 凡例 (L19 と完全に同じ意味体系)
KUIKI_INFO = {
    0: ("その他/非誘導",       "#cccccc", "誘導指定外の参考"),
    1: ("居住誘導区域",        "#1f78b4", "都計法 81 条 居住誘導 (面)"),
    2: ("居住誘導 境界線",     "#a6cee3", "面ではなく線 (細片)"),
    3: ("都市機能誘導区域",    "#e31a1c", "都計法 81 条 都市機能誘導 (面)"),
    4: ("都市機能誘導 境界線", "#fb9a99", "面ではなく線"),
}

# 都市タイプカラー
CTYPE_COLOR = {
    "政令市": "#cf222e", "中核市": "#bf3989",
    "施行時特例市": "#a371f7", "市": "#0969da", "町": "#1f883d",
}

# 公的統計参考値 (L19 と同値)
CITY_REF = {
    "広島市":   {"area_km2": 906.7, "pop_k": 1189},
    "呉市":     {"area_km2": 352.8, "pop_k":  210},
    "竹原市":   {"area_km2": 118.2, "pop_k":   24},
    "三原市":   {"area_km2": 471.6, "pop_k":   90},
    "福山市":   {"area_km2": 518.1, "pop_k":  459},
    "府中市":   {"area_km2": 195.8, "pop_k":   37},
    "東広島市": {"area_km2": 635.3, "pop_k":  198},
    "廿日市市": {"area_km2": 489.5, "pop_k":  117},
}
LIP_CITY_ORDER = [d[2] for d in LIP_CITY_DEFS]


# =============================================================================
# 1. 9 GeoJSON 統合読み込み + L19 比較データ + 8 市の行政区域
# =============================================================================
print("\n[1] L28 / L19 GeoJSON 読み込み", flush=True)
t1 = time.time()


def load_geojson_zip(zip_path: Path) -> gpd.GeoDataFrame:
    """ZIP 内の単一 .geojson を BytesIO 経由で読む (一時ファイル不要)"""
    with zipfile.ZipFile(zip_path) as zf:
        gjs = [n for n in zf.namelist() if n.lower().endswith(".geojson")]
        if not gjs:
            raise FileNotFoundError(f"no .geojson in {zip_path.name}")
        with zf.open(gjs[0]) as f:
            return gpd.read_file(io.BytesIO(f.read()))


# L28: 8 市別 GeoJSON
l28_frames = []
for ds28, ds19, name, ctype, _adm_ds in LIP_CITY_DEFS:
    g = load_geojson_zip(DATA_DIR / f"toshikinou_{ds28}_{name}.zip")
    g["source_city"] = name
    g["source_dsid_l28"] = ds28
    g["ctype"] = ctype
    l28_frames.append(g)

l28 = gpd.GeoDataFrame(pd.concat(l28_frames, ignore_index=True),
                       geometry="geometry", crs=l28_frames[0].crs)
l28 = l28.to_crs(TARGET_CRS)
l28["KUIKI_NAME"] = l28["KUIKI_CD"].map(lambda c: KUIKI_INFO[c][0])
l28["poly_area_km2"] = l28.geometry.area / 1e6
l28["poly_perim_km"] = l28.geometry.length / 1e3
l28["is_area"] = l28["KUIKI_CD"].isin([1, 3])
l28["is_line"] = l28["KUIKI_CD"].isin([2, 4])

# L19: 同 8 市別 GeoJSON (比較用)
l19_frames = []
for _ds28, ds19, name, ctype, _adm_ds in LIP_CITY_DEFS:
    g = load_geojson_zip(DATA_DIR_L19 / f"jukyo_{ds19}_{name}.zip")
    g["source_city"] = name
    g["source_dsid_l19"] = ds19
    g["ctype"] = ctype
    l19_frames.append(g)
l19 = gpd.GeoDataFrame(pd.concat(l19_frames, ignore_index=True),
                       geometry="geometry", crs=l19_frames[0].crs)
l19 = l19.to_crs(TARGET_CRS)
l19["poly_area_km2"] = l19.geometry.area / 1e6

# 県全域 2 件 (L28 ds=934, L19 ds=933)
ken28 = load_geojson_zip(DATA_DIR / f"toshikinou_{KEN_DSID_L28}_広島県.zip").to_crs(TARGET_CRS)
ken28["poly_area_km2"] = ken28.geometry.area / 1e6
ken19 = load_geojson_zip(DATA_DIR_L19 / f"jukyo_{KEN_DSID_L19}_広島県.zip").to_crs(TARGET_CRS)
ken19["poly_area_km2"] = ken19.geometry.area / 1e6

# 行政区域 (L15 共有, 8 LIP 市のみ)
admin_frames = []
for _, _, name, _, adm_ds in LIP_CITY_DEFS:
    g = load_geojson_zip(ADMIN_DIR / f"admin_{adm_ds}_{name}.zip")
    g["city"] = name
    admin_frames.append(g)
admin_all = gpd.GeoDataFrame(pd.concat(admin_frames, ignore_index=True),
                             geometry="geometry", crs=admin_frames[0].crs)
admin_all = admin_all.to_crs(TARGET_CRS)
admin_diss = admin_all.dissolve(by="city").reset_index()
admin_diss["admin_area_km2"] = admin_diss.geometry.area / 1e6

print(f"  L28 統合: {len(l28)} polys (8 市別)", flush=True)
print(f"  L19 統合: {len(l19)} polys (8 市別, 比較用)", flush=True)
print(f"  L28 県全域 ds={KEN_DSID_L28}: {len(ken28)} polys", flush=True)
print(f"  L19 県全域 ds={KEN_DSID_L19}: {len(ken19)} polys", flush=True)
print(f"  行政区域 8 市: {admin_diss['admin_area_km2'].sum():.0f} km²"
      f" ({time.time()-t1:.2f}s)", flush=True)


# =============================================================================
# 2. KUIKI_CD 構造の解読: L28 vs L19 で同じ KUIKI_CD でも面積が違う
# =============================================================================
print("\n[2] KUIKI_CD 構造解読", flush=True)
t1 = time.time()


def n_parts(g):
    if g is None or g.is_empty:
        return 0
    if g.geom_type == "MultiPolygon":
        return len(list(g.geoms))
    return 1


# L28 / L19 の市町×KUIKI 集計を別々に作る
def kuiki_agg(gdf, src_label):
    rows = []
    for name in LIP_CITY_ORDER:
        sub = gdf[gdf["source_city"] == name]
        for k in sorted(sub["KUIKI_CD"].unique()):
            sub_k = sub[sub["KUIKI_CD"] == k]
            rows.append({
                "src": src_label,
                "city": name,
                "KUIKI_CD": int(k),
                "KUIKI_NAME": KUIKI_INFO[int(k)][0],
                "n_polys": int(len(sub_k)),
                "area_km2": float(sub_k["poly_area_km2"].sum()),
            })
    return pd.DataFrame(rows)


l28_agg = kuiki_agg(l28, "L28 (toshikinou_*)")
l19_agg = kuiki_agg(l19, "L19 (jukyo_*)")
two_agg = pd.concat([l28_agg, l19_agg], ignore_index=True)

# L28 / L19 同 KUIKI_CD 横並び
side_by_side = []
for name in LIP_CITY_ORDER:
    for k in [0, 1, 2, 3, 4]:
        a28 = l28_agg[(l28_agg["city"] == name) &
                      (l28_agg["KUIKI_CD"] == k)]["area_km2"].sum()
        n28 = l28_agg[(l28_agg["city"] == name) &
                      (l28_agg["KUIKI_CD"] == k)]["n_polys"].sum()
        a19 = l19_agg[(l19_agg["city"] == name) &
                      (l19_agg["KUIKI_CD"] == k)]["area_km2"].sum()
        n19 = l19_agg[(l19_agg["city"] == name) &
                      (l19_agg["KUIKI_CD"] == k)]["n_polys"].sum()
        side_by_side.append({
            "city": name,
            "KUIKI_CD": k,
            "KUIKI_NAME": KUIKI_INFO[k][0],
            "L28_n": int(n28),
            "L19_n": int(n19),
            "L28_area_km2": round(float(a28), 4),
            "L19_area_km2": round(float(a19), 4),
            "L28_div_L19_pct": round(a28 / a19 * 100, 2) if a19 > 0 else np.nan,
        })
side_by_side_df = pd.DataFrame(side_by_side)
print(f"  L28/L19 KUIKI 比較表 完成 ({time.time()-t1:.2f}s)", flush=True)


# =============================================================================
# 3. 入れ子検証: L28 KUIKI=1 ⊆ L19 KUIKI=1 の幾何包含
# =============================================================================
print("\n[3] 入れ子検証 (L28 ⊆ L19 K1)", flush=True)
t1 = time.time()

nest_rows = []
for ds28, ds19, name, ctype, _ in LIP_CITY_DEFS:
    g28_1 = l28[(l28["source_city"] == name) & (l28["KUIKI_CD"] == 1)]
    g19_1 = l19[(l19["source_city"] == name) & (l19["KUIKI_CD"] == 1)]
    if len(g28_1) == 0 or len(g19_1) == 0:
        nest_rows.append({"city": name, "L28_K1_km2": 0.0, "L19_K1_km2": 0.0,
                          "L28_in_L19_pct": np.nan,
                          "L28_outside_L19_km2": 0.0})
        continue
    g28_diss = g28_1.dissolve().geometry.iloc[0]
    g19_diss = g19_1.dissolve().geometry.iloc[0]
    a28 = g28_diss.area / 1e6
    a19 = g19_diss.area / 1e6
    inter = g28_diss.intersection(g19_diss).area / 1e6
    outside = a28 - inter
    pct = inter / a28 * 100 if a28 > 0 else np.nan
    nest_rows.append({
        "city": name, "L28_K1_km2": a28, "L19_K1_km2": a19,
        "L28_in_L19_pct": pct,
        "L28_outside_L19_km2": outside,
        "L19_minus_L28_km2": a19 - inter,
    })
nest_df = pd.DataFrame(nest_rows)
print(f"  ({time.time()-t1:.2f}s)", flush=True)


# =============================================================================
# 4. 都市機能誘導 (L28 KUIKI=3) の構造分析
# =============================================================================
print("\n[4] 都市機能誘導 (KUIKI=3) 構造分析", flush=True)
t1 = time.time()

# 8 市別 都市機能誘導サマリ
toshi_rows = []
for ds28, ds19, name, ctype, _adm in LIP_CITY_DEFS:
    sub = l28[(l28["source_city"] == name) & (l28["KUIKI_CD"] == 3)]
    n_polys = len(sub)
    a_sum = float(sub["poly_area_km2"].sum())
    # 拠点中心点 (各ポリゴンの representative_point)
    pts = [pg.representative_point() for pg in sub.geometry]
    # 拠点間平均距離 (km)
    if len(pts) >= 2:
        dists = []
        for i in range(len(pts)):
            for j in range(i + 1, len(pts)):
                dists.append(pts[i].distance(pts[j]) / 1000)
        avg_d = float(np.mean(dists))
        max_d = float(np.max(dists))
    else:
        avg_d = np.nan
        max_d = np.nan

    # 居住誘導 (L28 KUIKI=1) との重なり
    g_k = sub.dissolve().geometry.iloc[0] if len(sub) > 0 else None
    g_j_l28 = l28[(l28["source_city"] == name) & (l28["KUIKI_CD"] == 1)]
    g_j_diss = g_j_l28.dissolve().geometry.iloc[0] if len(g_j_l28) > 0 else None

    if g_k is not None and g_j_diss is not None and not g_k.is_empty:
        inter_kj = g_k.intersection(g_j_diss).area / 1e6
        a_k = g_k.area / 1e6
        kinou_in_jukyo = inter_kj / a_k * 100 if a_k > 0 else 0
        # 都市機能ポリゴン中心点と L28 KUIKI=1 dissolve の最短距離 (km)
        if len(sub) > 0:
            min_d_to_jukyo = min(
                pg.distance(g_j_diss) / 1000 for pg in sub.geometry
            )
        else:
            min_d_to_jukyo = np.nan
    else:
        kinou_in_jukyo = np.nan
        min_d_to_jukyo = np.nan

    # L19 居住誘導 (全体) との重なり・距離 ─ より広い範囲で検証
    g_j_l19 = l19[(l19["source_city"] == name) & (l19["KUIKI_CD"] == 1)]
    g_j_l19_diss = (g_j_l19.dissolve().geometry.iloc[0]
                    if len(g_j_l19) > 0 else None)
    if g_k is not None and g_j_l19_diss is not None and not g_k.is_empty:
        inter_kj_l19 = g_k.intersection(g_j_l19_diss).area / 1e6
        kinou_in_jukyo_l19 = inter_kj_l19 / a_k * 100 if a_k > 0 else 0
        if len(sub) > 0:
            min_d_to_jukyo_l19 = min(
                pg.distance(g_j_l19_diss) / 1000 for pg in sub.geometry
            )
        else:
            min_d_to_jukyo_l19 = np.nan
    else:
        kinou_in_jukyo_l19 = np.nan
        min_d_to_jukyo_l19 = np.nan

    admin_a = float(admin_diss[admin_diss["city"] == name]
                    ["admin_area_km2"].iloc[0])
    pop = CITY_REF[name]["pop_k"]

    toshi_rows.append({
        "city": name, "ctype": ctype,
        "admin_area_km2": admin_a,
        "pop_k": pop,
        "n_kinou_polys": n_polys,
        "kinou_area_km2": a_sum,
        "kinou_in_admin_pct": a_sum / admin_a * 100,
        "kinou_in_jukyo_pct_l28": kinou_in_jukyo,
        "kinou_min_dist_to_jukyo_l28_km": min_d_to_jukyo,
        "kinou_in_jukyo_pct_l19": kinou_in_jukyo_l19,
        "kinou_min_dist_to_jukyo_l19_km": min_d_to_jukyo_l19,
        "kinou_avg_pairwise_dist_km": avg_d,
        "kinou_max_pairwise_dist_km": max_d,
        "kinou_per_pop_km2_per_100k": (
            a_sum / (pop / 100) if pop > 0 else np.nan
        ),
        "has_kinou": n_polys > 0,
    })
toshi_df = pd.DataFrame(toshi_rows)
print(f"  ({time.time()-t1:.2f}s)", flush=True)


# =============================================================================
# 5. 中間 CSV 出力
# =============================================================================
print("\n[5] 中間 CSV 出力", flush=True)
ASSETS.mkdir(parents=True, exist_ok=True)

# (1) L28 全ポリゴン詳細
l28.drop(columns="geometry").round(4).to_csv(
    ASSETS / "L28_poly_detail.csv", index=False, encoding="utf-8-sig")

# (2) 市町×KUIKI_CD: L28 vs L19 横並び
side_by_side_df.to_csv(ASSETS / "L28_vs_L19_kuiki_compare.csv",
                       index=False, encoding="utf-8-sig")

# (3) 入れ子検証
nest_df.round(4).to_csv(ASSETS / "L28_nest_in_L19.csv",
                        index=False, encoding="utf-8-sig")

# (4) 都市機能誘導 8 市サマリ
toshi_df.round(4).to_csv(ASSETS / "L28_toshikinou_summary.csv",
                         index=False, encoding="utf-8-sig")

# (5) L28 全県 KUIKI 集計
l28_kuiki_total = l28.groupby(["KUIKI_CD", "KUIKI_NAME"]).agg(
    polys_count=("poly_area_km2", "size"),
    area_km2=("poly_area_km2", "sum"),
).reset_index().round(4)
l28_kuiki_total.to_csv(ASSETS / "L28_kuiki_total.csv",
                       index=False, encoding="utf-8-sig")

# (6) 8 市 × KUIKI 面積 pivot
pivot_area = l28.pivot_table(
    index="source_city", columns="KUIKI_NAME",
    values="poly_area_km2", aggfunc="sum", fill_value=0
).reindex(LIP_CITY_ORDER).round(4)
pivot_area["合計"] = pivot_area.sum(axis=1)
pivot_area.to_csv(ASSETS / "L28_pivot_area.csv", encoding="utf-8-sig")

# (7) 8 市 × KUIKI 件数 pivot
pivot_n = l28.pivot_table(
    index="source_city", columns="KUIKI_NAME",
    values="poly_area_km2", aggfunc="size", fill_value=0
).reindex(LIP_CITY_ORDER)
pivot_n["合計"] = pivot_n.sum(axis=1)
pivot_n.to_csv(ASSETS / "L28_pivot_n.csv", encoding="utf-8-sig")

print("  saved 7 CSVs", flush=True)


# =============================================================================
# 6. 図 1: 8 LIP 市 主題図 (L28 居住誘導 + 都市機能誘導 + 行政背景)
# =============================================================================
print("\n[6] 図1: 8 LIP 市 主題図", flush=True)
t1 = time.function() if False else time.time()
fig, ax = plt.subplots(figsize=(13, 9))

# 行政背景
admin_diss.plot(ax=ax, color="#f5f5f5",
                edgecolor="#888888", linewidth=0.6, zorder=1)
# L28 KUIKI=1 (居住誘導の近接サブセット, 青)
sub1 = l28[l28["KUIKI_CD"] == 1]
if len(sub1) > 0:
    sub1.plot(ax=ax, color=KUIKI_INFO[1][1], edgecolor="white",
              linewidth=0.05, alpha=0.8, zorder=2)
# L28 KUIKI=3 (都市機能誘導, 強い赤)
sub3 = l28[l28["KUIKI_CD"] == 3]
if len(sub3) > 0:
    sub3.plot(ax=ax, color=KUIKI_INFO[3][1], edgecolor="black",
              linewidth=0.4, alpha=0.95, zorder=3)
# 8 市ラベル
for _, r in admin_diss.iterrows():
    pt = r.geometry.representative_point()
    ax.annotate(r["city"], (pt.x, pt.y), ha="center", va="center",
                fontsize=8.5, fontweight="bold")

a1_total = float(sub1["poly_area_km2"].sum())
a3_total = float(sub3["poly_area_km2"].sum())
ax.set_title(
    f"8 LIP 策定市 — 都市機能誘導区域シリーズ (L28)\n"
    f"L28 KUIKI=1 (近接居住誘導 青) {a1_total:.1f} km² / "
    f"KUIKI=3 (都市機能誘導 赤) {a3_total:.2f} km²",
    fontsize=11.5)
ax.set_axis_off()
patches = [
    Patch(facecolor="#f5f5f5", edgecolor="#888", label="8 LIP 市行政区域"),
    Patch(facecolor=KUIKI_INFO[1][1],
          label=f"L28 KUIKI=1 (近接居住誘導) {a1_total:.1f} km²"),
    Patch(facecolor=KUIKI_INFO[3][1], edgecolor="black",
          label=f"L28 KUIKI=3 (都市機能誘導) {a3_total:.2f} km²"),
]
ax.legend(handles=patches, loc="lower left", fontsize=10, frameon=True)
plt.tight_layout()
plt.savefig(ASSETS / "L28_thematic.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L28_thematic.png ({time.time()-t1:.2f}s)", flush=True)


# =============================================================================
# 7. 図 2: 8 市 small multiples (各市の都市機能誘導拠点配置)
# =============================================================================
print("\n[7] 図2: 8 市 small multiples", flush=True)
t1 = time.time()
fig, axes = plt.subplots(2, 4, figsize=(16, 8.5))
axes = axes.flatten()
for i, c in enumerate(LIP_CITY_ORDER):
    ax = axes[i]
    a = admin_diss[admin_diss["city"] == c]
    a.plot(ax=ax, color="#fafafa", edgecolor="#888888", linewidth=0.5)
    sub_j = l28[(l28["source_city"] == c) & (l28["KUIKI_CD"] == 1)]
    sub_k = l28[(l28["source_city"] == c) & (l28["KUIKI_CD"] == 3)]
    if len(sub_j) > 0:
        sub_j.plot(ax=ax, color=KUIKI_INFO[1][1], edgecolor="white",
                   linewidth=0.08, alpha=0.75)
    if len(sub_k) > 0:
        sub_k.plot(ax=ax, color=KUIKI_INFO[3][1], edgecolor="black",
                   linewidth=0.5, alpha=0.95)
    csum = toshi_df[toshi_df["city"] == c].iloc[0]
    nj = (l28["source_city"] == c) & (l28["KUIKI_CD"] == 1)
    nk = (l28["source_city"] == c) & (l28["KUIKI_CD"] == 3)
    j_lab = (
        f"近接居住 {l28[nj]['poly_area_km2'].sum():.1f} km² ({nj.sum()} polys)"
        if nj.sum() > 0 else "近接居住 — なし"
    )
    k_lab = (
        f"都市機能 {csum['kinou_area_km2']:.2f} km² ({int(csum['n_kinou_polys'])} 拠点)"
        if csum["has_kinou"] else "都市機能 — なし"
    )
    ax.set_title(f"{c} ({csum['ctype']})\n{j_lab}\n{k_lab}", fontsize=9)
    ax.set_axis_off()
plt.suptitle("8 LIP 策定市 — L28 都市機能誘導 (赤) + 近接居住誘導 (青)",
             fontsize=12, y=1.00)
plt.tight_layout()
plt.savefig(ASSETS / "L28_city_small_multiples.png",
            dpi=120, bbox_inches="tight")
plt.close("all")
print(f"  saved ({time.time()-t1:.2f}s)", flush=True)


# =============================================================================
# 8. 図 3: L28 vs L19 KUIKI=1 面積比較 (核心図 — シリーズ実態の解読)
# =============================================================================
print("\n[8] 図3: L28 vs L19 KUIKI=1 面積比較", flush=True)
t1 = time.time()
fig, axes = plt.subplots(1, 2, figsize=(15, 7))

# 左: 8 市 KUIKI=1 面積 横棒 (L28 vs L19)
df_cmp = []
for name in LIP_CITY_ORDER:
    a28 = l28[(l28["source_city"] == name) &
              (l28["KUIKI_CD"] == 1)]["poly_area_km2"].sum()
    a19 = l19[(l19["source_city"] == name) &
              (l19["KUIKI_CD"] == 1)]["poly_area_km2"].sum()
    df_cmp.append({"city": name, "L28": a28, "L19": a19})
df_cmp = pd.DataFrame(df_cmp)
y = np.arange(len(df_cmp))
h = 0.4
axes[0].barh(y - h/2, df_cmp["L19"], h,
             color="#a6cee3", edgecolor="black", linewidth=0.4,
             label="L19 居住誘導 (jukyo) シリーズ")
axes[0].barh(y + h/2, df_cmp["L28"], h,
             color="#1f78b4", edgecolor="black", linewidth=0.4,
             label="L28 都市機能誘導 (toshikinou) シリーズ")
for i, r in enumerate(df_cmp.itertuples()):
    if r.L19 > 0:
        ratio = r.L28 / r.L19 * 100
        axes[0].text(max(r.L19, r.L28) + 2, i,
                     f"L28/L19={ratio:.0f}%", va="center", fontsize=8.5)
axes[0].set_yticks(y)
axes[0].set_yticklabels(df_cmp["city"])
axes[0].set_xlabel("KUIKI=1 (居住誘導) 面積 km²")
axes[0].set_title("KUIKI=1 (居住誘導) ─ 同じ KUIKI_CD でも L28 << L19\n"
                  "= L28 シリーズは L19 居住誘導の狭いサブセット",
                  fontsize=11)
axes[0].grid(axis="x", alpha=0.3)
axes[0].legend(loc="lower right", fontsize=9)

# 右: KUIKI=1 / KUIKI=3 / KUIKI=2,4(線) 面積総和 4 棒
totals_l28 = {1: 0, 3: 0, "line(2,4)": 0}
totals_l19 = {1: 0, 3: 0, "line(2,4)": 0}
for k in [1, 3]:
    totals_l28[k] = l28[l28["KUIKI_CD"] == k]["poly_area_km2"].sum()
    totals_l19[k] = l19[l19["KUIKI_CD"] == k]["poly_area_km2"].sum()
totals_l28["line(2,4)"] = l28[l28["KUIKI_CD"].isin([2, 4])]["poly_area_km2"].sum()
totals_l19["line(2,4)"] = l19[l19["KUIKI_CD"].isin([2, 4])]["poly_area_km2"].sum()

cats = ["KUIKI=1\n居住誘導", "KUIKI=3\n都市機能誘導", "KUIKI=2,4\n境界線(線データ)"]
v_l28 = [totals_l28[1], totals_l28[3], totals_l28["line(2,4)"]]
v_l19 = [totals_l19[1], totals_l19[3], totals_l19["line(2,4)"]]
xs = np.arange(len(cats))
w = 0.35
axes[1].bar(xs - w/2, v_l19, w, color="#a6cee3",
            edgecolor="black", linewidth=0.4, label="L19 居住誘導 シリーズ 8 市計")
axes[1].bar(xs + w/2, v_l28, w, color="#1f78b4",
            edgecolor="black", linewidth=0.4, label="L28 都市機能誘導 シリーズ 8 市計")
for i, (a, b) in enumerate(zip(v_l19, v_l28)):
    axes[1].text(i - w/2, a + 1, f"{a:.1f}", ha="center", fontsize=9)
    axes[1].text(i + w/2, b + 1, f"{b:.1f}", ha="center", fontsize=9)
axes[1].set_xticks(xs)
axes[1].set_xticklabels(cats, fontsize=10)
axes[1].set_ylabel("8 市合計面積 km²")
axes[1].set_title("L28 vs L19 ─ KUIKI_CD 別面積総和\n"
                  "(L28 KUIKI=1 は L19 KUIKI=1 の約 4 割; KUIKI=3 はほぼ同等)",
                  fontsize=11)
axes[1].legend(loc="upper right", fontsize=9)
axes[1].grid(axis="y", alpha=0.3)

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


# =============================================================================
# 9. 図 4: 廿日市市での L28 ⊆ L19 入れ子クローズアップ (実証図)
# =============================================================================
print("\n[9] 図4: 入れ子クローズアップ", flush=True)
t1 = time.time()
fig, axes = plt.subplots(1, 3, figsize=(16, 6))
focus_city = "廿日市市"  # 入れ子が最も明瞭な市

a = admin_diss[admin_diss["city"] == focus_city]
g28_1 = l28[(l28["source_city"] == focus_city) & (l28["KUIKI_CD"] == 1)]
g19_1 = l19[(l19["source_city"] == focus_city) & (l19["KUIKI_CD"] == 1)]
g28_3 = l28[(l28["source_city"] == focus_city) & (l28["KUIKI_CD"] == 3)]

# (1) L19 KUIKI=1 単独
a.plot(ax=axes[0], color="#fafafa", edgecolor="#888", linewidth=0.5)
g19_1.plot(ax=axes[0], color="#a6cee3", edgecolor="white",
           linewidth=0.1, alpha=0.85)
axes[0].set_title(
    f"(a) L19 居住誘導 全体\n"
    f"jukyo_886 (KUIKI=1) {g19_1['poly_area_km2'].sum():.2f} km², {len(g19_1)} polys",
    fontsize=10)
b = g19_1.total_bounds
pad = max(b[2]-b[0], b[3]-b[1]) * 0.08
axes[0].set_xlim(b[0]-pad, b[2]+pad)
axes[0].set_ylim(b[1]-pad, b[3]+pad)
axes[0].set_axis_off()

# (2) L28 KUIKI=1 (青) + L28 KUIKI=3 (赤) 重ね
a.plot(ax=axes[1], color="#fafafa", edgecolor="#888", linewidth=0.5)
g19_1.plot(ax=axes[1], color="#eef5f9", edgecolor="white",
           linewidth=0.05, alpha=0.7)  # L19 を薄背景
g28_1.plot(ax=axes[1], color="#1f78b4", edgecolor="white",
           linewidth=0.2, alpha=0.85)
g28_3.plot(ax=axes[1], color="#e31a1c", edgecolor="black",
           linewidth=0.5, alpha=0.95)
axes[1].set_title(
    f"(b) L28 切り出し版\n"
    f"L28 KUIKI=1 (濃青) {g28_1['poly_area_km2'].sum():.2f} km², "
    f"KUIKI=3 (赤) {g28_3['poly_area_km2'].sum():.2f} km²\n"
    f"薄青背景 = L19 全体 (差分が見える)",
    fontsize=10)
axes[1].set_xlim(b[0]-pad, b[2]+pad)
axes[1].set_ylim(b[1]-pad, b[3]+pad)
axes[1].set_axis_off()

# (3) L28 KUIKI=3 (都市機能誘導) ズーム — 拠点近傍
a.plot(ax=axes[2], color="#fafafa", edgecolor="#888", linewidth=0.5)
g28_1.plot(ax=axes[2], color="#1f78b4", edgecolor="white",
           linewidth=0.2, alpha=0.6)
g28_3.plot(ax=axes[2], color="#e31a1c", edgecolor="black",
           linewidth=0.7, alpha=0.95)
# 拠点ラベル
for _, r in g28_3.iterrows():
    pt = r.geometry.representative_point()
    axes[2].annotate(f"拠点{int(r['FID'])}",
                     (pt.x, pt.y),
                     fontsize=8, ha="center", va="bottom",
                     color="black", fontweight="bold",
                     xytext=(0, 8), textcoords="offset points",
                     arrowprops=dict(arrowstyle="-", color="black", lw=0.5))
b3 = g28_3.total_bounds if len(g28_3) > 0 else b
pad3 = max(b3[2]-b3[0], b3[3]-b3[1]) * 0.5
axes[2].set_xlim(b3[0]-pad3, b3[2]+pad3)
axes[2].set_ylim(b3[1]-pad3, b3[3]+pad3)
axes[2].set_title(
    f"(c) 都市機能誘導 ズーム\n"
    f"{focus_city} の {len(g28_3)} 拠点 (中心市街地に集中)",
    fontsize=10)
axes[2].set_axis_off()

plt.suptitle(
    f"L28 ⊆ L19 入れ子の幾何実証 ({focus_city})\n"
    "L28 シリーズは L19 居住誘導から「都市機能拠点近傍だけ」を切り出した狭い領域を提供",
    fontsize=12, y=1.02)
plt.tight_layout()
plt.savefig(ASSETS / "L28_nest_closeup.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved ({time.time()-t1:.2f}s)", flush=True)


# =============================================================================
# 10. 図 5: 入れ子比率 choropleth + 散布
# =============================================================================
print("\n[10] 図5: 入れ子比率 choropleth", flush=True)
t1 = time.time()
fig, axes = plt.subplots(1, 2, figsize=(15, 7))

adm_n = admin_diss.merge(nest_df, on="city", how="left")
# 左: L28 ∩ L19 / L28 比率 choropleth (高=完全包含, 低=L28 が L19 を超えてはみ出ている)
adm_n.plot(ax=axes[0], column="L28_in_L19_pct",
           cmap="RdYlBu", vmin=80, vmax=100,
           edgecolor="#444", linewidth=0.5,
           legend=True, missing_kwds={"color": "#dddddd"},
           legend_kwds={"label": "L28 ∩ L19 / L28 (%) [100=完全包含]",
                        "shrink": 0.6})
for _, r in adm_n.iterrows():
    pt = r.geometry.representative_point()
    pct = r["L28_in_L19_pct"]
    label = (f"{r['city']}\n{pct:.1f}%" if not np.isnan(pct)
             else f"{r['city']}\nN/A")
    axes[0].annotate(label, (pt.x, pt.y), ha="center", va="center",
                     fontsize=8.5, fontweight="bold")
axes[0].set_title("L28 KUIKI=1 ⊆ L19 KUIKI=1 包含率 choropleth\n"
                  "(100% = L28 KUIKI=1 ポリゴンは L19 KUIKI=1 に完全に含まれる)",
                  fontsize=11)
axes[0].set_axis_off()

# 右: L28 KUIKI=1 vs L19 KUIKI=1 面積散布
for _, r in nest_df.iterrows():
    color = CTYPE_COLOR.get(
        toshi_df[toshi_df["city"] == r["city"]]["ctype"].iloc[0], "#888"
    )
    axes[1].scatter(r["L19_K1_km2"], r["L28_K1_km2"],
                    s=160, c=color, alpha=0.85,
                    edgecolor="black", linewidth=0.6)
    axes[1].annotate(r["city"],
                     (r["L19_K1_km2"], r["L28_K1_km2"]),
                     fontsize=9, xytext=(6, 4),
                     textcoords="offset points", fontweight="bold")
# 1:1 ライン
mx = max(nest_df["L19_K1_km2"].max(), 0.1) * 1.1
xs_line = np.linspace(0, mx, 50)
axes[1].plot(xs_line, xs_line, "--", color="#888", linewidth=1,
             label="L28 = L19 (理論的等価線)")
# 0.4 ライン
axes[1].plot(xs_line, xs_line * 0.4, ":", color="#cf222e", linewidth=1,
             label="L28 = 0.4 × L19 (8 市平均線)")
axes[1].set_xlabel("L19 KUIKI=1 面積 km² (居住誘導全体)")
axes[1].set_ylabel("L28 KUIKI=1 面積 km² (都市機能近接サブセット)")
axes[1].set_title("L19 vs L28 KUIKI=1 面積散布\n"
                  "(全 8 市が 1:1 線より下方 = L28 が常に L19 のサブセット)",
                  fontsize=11)
axes[1].grid(True, alpha=0.3)
axes[1].legend(loc="lower right", fontsize=9)
ctype_legend = [Patch(facecolor=CTYPE_COLOR[t], label=t)
                for t in ["政令市", "中核市", "施行時特例市", "市"]]
axes[1].legend(handles=ctype_legend + [
    plt.Line2D([0], [0], linestyle="--", color="#888",
               label="L28 = L19 線"),
    plt.Line2D([0], [0], linestyle=":", color="#cf222e",
               label="L28 = 0.4 L19 線"),
], loc="upper left", fontsize=8.5)

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


# =============================================================================
# 11. 図 6: 都市機能誘導の市町別ランキング (拠点数 + 総面積)
# =============================================================================
print("\n[11] 図6: 都市機能誘導 ランキング", flush=True)
t1 = time.time()
fig, axes = plt.subplots(1, 2, figsize=(15, 7))

td = toshi_df.sort_values("kinou_area_km2", ascending=True).copy()
y = np.arange(len(td))
bar_colors = [CTYPE_COLOR[ct] for ct in td["ctype"]]
axes[0].barh(y, td["kinou_area_km2"], color=bar_colors,
             edgecolor="black", linewidth=0.4)
for i, r in enumerate(td.itertuples()):
    label = (f"{r.kinou_area_km2:.2f} km² "
             f"({int(r.n_kinou_polys)} 拠点)"
             if r.kinou_area_km2 > 0 else "未指定")
    axes[0].text(r.kinou_area_km2 + 0.05, i, label,
                 va="center", fontsize=9)
axes[0].set_yticks(y)
axes[0].set_yticklabels(td["city"])
axes[0].set_xlabel("都市機能誘導 (KUIKI=3) 総面積 km²")
axes[0].set_title("8 LIP 市 — 都市機能誘導 総面積ランキング", fontsize=11)
axes[0].grid(axis="x", alpha=0.3)
ctype_legend = [Patch(facecolor=CTYPE_COLOR[t], label=t)
                for t in ["政令市", "中核市", "施行時特例市", "市"]]
axes[0].legend(handles=ctype_legend, loc="lower right", fontsize=9)

# 右: 拠点数 vs 人口 散布
for _, r in toshi_df.iterrows():
    color = CTYPE_COLOR[r["ctype"]]
    axes[1].scatter(r["pop_k"], r["n_kinou_polys"],
                    s=r["kinou_area_km2"] * 80 + 60,
                    c=color, alpha=0.85,
                    edgecolor="black", linewidth=0.6)
    axes[1].annotate(r["city"],
                     (r["pop_k"], r["n_kinou_polys"]),
                     fontsize=9, xytext=(6, 4),
                     textcoords="offset points", fontweight="bold")
axes[1].set_xscale("log")
axes[1].set_xlabel("人口 千人 (log)")
axes[1].set_ylabel("都市機能誘導 拠点 (= 連結ポリゴン) 数")
axes[1].set_title("人口 vs 拠点数 (バブル=拠点総面積)\n"
                  "= 多核拠点ネットワーク型コンパクトシティの実装スケール",
                  fontsize=11)
axes[1].grid(True, which="both", alpha=0.3)

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


# =============================================================================
# 12. 図 7: 拠点間距離・拠点と居住誘導の最短距離
# =============================================================================
print("\n[12] 図7: 拠点間距離 / 拠点 vs 居住誘導 距離", flush=True)
t1 = time.time()
fig, axes = plt.subplots(1, 2, figsize=(15, 6.5))

# 左: 拠点間平均距離 vs 行政面積
td2 = toshi_df.dropna(subset=["kinou_avg_pairwise_dist_km"]).copy()
for _, r in td2.iterrows():
    color = CTYPE_COLOR[r["ctype"]]
    axes[0].scatter(r["admin_area_km2"], r["kinou_avg_pairwise_dist_km"],
                    s=140, c=color, alpha=0.85,
                    edgecolor="black", linewidth=0.6)
    axes[0].annotate(
        r["city"], (r["admin_area_km2"], r["kinou_avg_pairwise_dist_km"]),
        fontsize=9, xytext=(6, 4),
        textcoords="offset points", fontweight="bold"
    )
axes[0].set_xscale("log")
axes[0].set_xlabel("行政面積 km² (log)")
axes[0].set_ylabel("拠点間平均距離 km")
axes[0].set_title("行政面積 vs 拠点間平均距離\n"
                  "(広い市は拠点が分散 = 多核, 狭い市は単核近接)",
                  fontsize=11)
axes[0].grid(True, which="both", alpha=0.3)

# 右: 拠点 → L28 K1 と L19 K1 の最短距離を併記 (分離の幾何実証)
td3 = toshi_df.dropna(subset=["kinou_min_dist_to_jukyo_l28_km"]).copy()
y = np.arange(len(td3))
h_bar = 0.4
axes[1].barh(y - h_bar/2, td3["kinou_min_dist_to_jukyo_l28_km"], h_bar,
             color="#1f78b4", edgecolor="black", linewidth=0.4,
             label="拠点 → L28 KUIKI=1 (近接居住誘導)")
axes[1].barh(y + h_bar/2, td3["kinou_min_dist_to_jukyo_l19_km"], h_bar,
             color="#a6cee3", edgecolor="black", linewidth=0.4,
             label="拠点 → L19 KUIKI=1 (居住誘導全体)")
for i, r in enumerate(td3.itertuples()):
    axes[1].text(r.kinou_min_dist_to_jukyo_l28_km + 0.05, i - h_bar/2,
                 f"{r.kinou_min_dist_to_jukyo_l28_km:.2f}",
                 va="center", fontsize=8)
    axes[1].text(r.kinou_min_dist_to_jukyo_l19_km + 0.05, i + h_bar/2,
                 f"{r.kinou_min_dist_to_jukyo_l19_km:.2f}",
                 va="center", fontsize=8)
axes[1].set_yticks(y)
axes[1].set_yticklabels(td3["city"])
axes[1].set_xlabel("都市機能拠点 → 居住誘導 までの最短距離 km")
axes[1].set_title("拠点 vs 居住誘導 最短距離 ─ 分離型 LIP の発見\n"
                  "(L28 KUIKI=3 と L19/L28 KUIKI=1 が km 単位で分離)",
                  fontsize=11)
axes[1].grid(axis="x", alpha=0.3)
axes[1].legend(loc="lower right", fontsize=9)

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


# =============================================================================
# 13. 図 8: 8 市 × KUIKI 構成 stack (面 vs 線, log)
# =============================================================================
print("\n[13] 図8: 8 市 × KUIKI 構成 stack", flush=True)
t1 = time.time()
fig, axes = plt.subplots(1, 2, figsize=(15, 7))

ct_n = l28.groupby(["source_city", "KUIKI_NAME"]).size().unstack(fill_value=0)
ct_n = ct_n.reindex(LIP_CITY_ORDER)
order_kuiki = ["居住誘導区域", "都市機能誘導区域",
               "居住誘導 境界線", "都市機能誘導 境界線", "その他/非誘導"]
order_kuiki = [k for k in order_kuiki if k in ct_n.columns]
ct_n_o = ct_n[order_kuiki]
xs = np.arange(len(ct_n_o))
bottom = np.zeros(len(ct_n_o))
colors_map = {n: KUIKI_INFO[k][1] for k, (n, _, _) in KUIKI_INFO.items()}
for col in order_kuiki:
    axes[0].bar(xs, ct_n_o[col], bottom=bottom,
                color=colors_map[col], edgecolor="white", linewidth=0.3,
                label=col)
    bottom = bottom + ct_n_o[col].values
for i in range(len(ct_n_o)):
    axes[0].text(i, bottom[i] * 1.04, str(int(bottom[i])),
                 ha="center", fontsize=8.5)
axes[0].set_yscale("log")
axes[0].set_xticks(xs)
axes[0].set_xticklabels(ct_n_o.index, rotation=30, fontsize=9)
axes[0].set_ylabel("ポリゴン数 (log)")
axes[0].set_title("L28: 8 市 × KUIKI_CD 構成 (件数, log)", fontsize=11)
axes[0].legend(loc="upper right", fontsize=8)

# 右: 面積構成 (面データ KUIKI=1,3 のみ)
ct_a = l28[l28["is_area"]].groupby(
    ["source_city", "KUIKI_NAME"])["poly_area_km2"].sum().unstack(fill_value=0)
ct_a = ct_a.reindex(LIP_CITY_ORDER)
order_face = [k for k in ["居住誘導区域", "都市機能誘導区域"]
              if k in ct_a.columns]
ct_a_o = ct_a[order_face]
bottom = np.zeros(len(ct_a_o))
for col in order_face:
    axes[1].bar(xs, ct_a_o[col], bottom=bottom,
                color=colors_map[col], edgecolor="white", linewidth=0.3,
                label=col)
    bottom = bottom + ct_a_o[col].values
for i in range(len(ct_a_o)):
    axes[1].text(i, bottom[i] + 1, f"{bottom[i]:.1f}",
                 ha="center", fontsize=8.5)
axes[1].set_xticks(xs)
axes[1].set_xticklabels(ct_a_o.index, rotation=30, fontsize=9)
axes[1].set_ylabel("面積 km²")
axes[1].set_title("L28: 8 市 × KUIKI_CD 面積構成 (面 KUIKI=1,3)\n"
                  "近接居住 (青) >> 都市機能 (赤) の入れ子構造",
                  fontsize=11)
axes[1].legend(loc="upper right", fontsize=8.5)
axes[1].grid(axis="y", alpha=0.3)

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


# =============================================================================
# 14. 図 9: シリーズ実態の図解 (L28 = L19 ∩ buffer(都市機能拠点) の概念図)
# =============================================================================
print("\n[14] 図9: シリーズ実態解読 概念", flush=True)
t1 = time.time()
fig, axes = plt.subplots(1, 2, figsize=(16, 7.5))

# 左: 呉市の L28 KUIKI=1 (青) + L19 KUIKI=1 (薄青) + L28 KUIKI=3 (赤)
focus2 = "呉市"
a2 = admin_diss[admin_diss["city"] == focus2]
g28_1_2 = l28[(l28["source_city"] == focus2) & (l28["KUIKI_CD"] == 1)]
g19_1_2 = l19[(l19["source_city"] == focus2) & (l19["KUIKI_CD"] == 1)]
g28_3_2 = l28[(l28["source_city"] == focus2) & (l28["KUIKI_CD"] == 3)]

a2.plot(ax=axes[0], color="#fafafa", edgecolor="#888", linewidth=0.5)
g19_1_2.plot(ax=axes[0], color="#c6dbef", edgecolor="white",
             linewidth=0.05, alpha=0.9, label="L19 居住誘導 全体")
g28_1_2.plot(ax=axes[0], color="#1f78b4", edgecolor="white",
             linewidth=0.2, alpha=0.85, label="L28 KUIKI=1 (近接居住)")
g28_3_2.plot(ax=axes[0], color="#e31a1c", edgecolor="black",
             linewidth=0.6, alpha=0.95, label="L28 KUIKI=3 (都市機能)")
b2 = g19_1_2.total_bounds
pad2 = max(b2[2]-b2[0], b2[3]-b2[1]) * 0.1
axes[0].set_xlim(b2[0]-pad2, b2[2]+pad2)
axes[0].set_ylim(b2[1]-pad2, b2[3]+pad2)
axes[0].set_title(
    f"{focus2}: L19 全体 (薄青) ⊃ L28 KUIKI=1 (濃青) ⊃ L28 KUIKI=3 (赤)\n"
    f"L19={g19_1_2['poly_area_km2'].sum():.2f} km², "
    f"L28K1={g28_1_2['poly_area_km2'].sum():.2f} km², "
    f"L28K3={g28_3_2['poly_area_km2'].sum():.2f} km²",
    fontsize=10.5)
axes[0].set_axis_off()
patch_legend = [
    Patch(facecolor="#fafafa", edgecolor="#888", label="行政区域"),
    Patch(facecolor="#c6dbef", label="L19 KUIKI=1 (居住誘導 全体)"),
    Patch(facecolor="#1f78b4", label="L28 KUIKI=1 (近接居住誘導)"),
    Patch(facecolor="#e31a1c", edgecolor="black", label="L28 KUIKI=3 (都市機能誘導)"),
]
axes[0].legend(handles=patch_legend, loc="lower left", fontsize=9)

# 右: 概念図 — L28 ≈ L19 ∩ Buffer(KUIKI=3, ~Rkm)
# 模擬データ: 円で示す
from matplotlib.patches import Circle
ax2 = axes[1]
ax2.set_xlim(-6, 6)
ax2.set_ylim(-4, 4)
# L19 居住誘導 (大円)
c_l19 = Circle((-0.5, 0), 4, color="#c6dbef", alpha=0.8, zorder=1)
ax2.add_patch(c_l19)
ax2.text(-3.5, 2.8, "L19 KUIKI=1\n居住誘導 全体",
         fontsize=10, color="#0050a0")
# L28 (中円, L19 のサブセット)
c_l28 = Circle((-0.5, 0), 2.5, color="#1f78b4", alpha=0.85, zorder=2)
ax2.add_patch(c_l28)
ax2.text(-2.6, 1.0, "L28 KUIKI=1\n近接居住誘導\n(= L19 の \"都市機能ハブ\nの近所\")",
         fontsize=9.5, color="white", fontweight="bold")
# L28 KUIKI=3 (小円, L28 KUIKI=1 内)
for x in [-1.0, 0.4, 1.5]:
    ax2.add_patch(Circle((x, 0), 0.4, color="#e31a1c", alpha=0.95, zorder=3))
ax2.text(-1.0, -0.9, "拠点1", fontsize=8, color="#cf222e",
         ha="center", fontweight="bold")
ax2.text(0.4, -0.9, "拠点2", fontsize=8, color="#cf222e",
         ha="center", fontweight="bold")
ax2.text(1.5, -0.9, "拠点3", fontsize=8, color="#cf222e",
         ha="center", fontweight="bold")
# 矢印
ax2.annotate("", xy=(-0.5, -3.2), xytext=(-0.5, -3.7),
             arrowprops=dict(arrowstyle="->", color="black", lw=1.5))
ax2.text(0, -3.5,
         "L28 シリーズ ≈ L19 ∩ buffer(KUIKI=3, R km)\n"
         "= 都市機能拠点を中心とした近接居住誘導の切り出し",
         fontsize=10, ha="center", fontweight="bold")
ax2.set_aspect("equal")
ax2.set_title("L28 シリーズ実態の概念図\n"
              "L19 と L28 は同じ列構造を持ちつつ、L28 は L19 の狭い切り出し",
              fontsize=11)
ax2.set_axis_off()

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


# =============================================================================
# 15. 図 10: 整合性検証 (8 市別 vs 県全域)
# =============================================================================
print("\n[15] 図10: 整合性検証 (8 市別 vs 県全域)", flush=True)
t1 = time.time()

# 8 市別 vs 県全域 (KUIKI 別)
sum_l28_city = l28.groupby("KUIKI_CD")["poly_area_km2"].sum().to_dict()
sum_l28_ken = ken28.groupby("KUIKI_CD")["poly_area_km2"].sum().to_dict() \
    if "KUIKI_CD" in ken28.columns else {}
n_l28_city = l28.groupby("KUIKI_CD").size().to_dict()
n_l28_ken = ken28.groupby("KUIKI_CD").size().to_dict() \
    if "KUIKI_CD" in ken28.columns else {}

ks = sorted(set(sum_l28_city) | set(sum_l28_ken))
labels, vals_a_city, vals_a_ken, vals_n_city, vals_n_ken = [], [], [], [], []
for k in ks:
    labels.append(f"K={k}\n{KUIKI_INFO[k][0][:8]}")
    vals_a_city.append(sum_l28_city.get(k, 0))
    vals_a_ken.append(sum_l28_ken.get(k, 0))
    vals_n_city.append(n_l28_city.get(k, 0))
    vals_n_ken.append(n_l28_ken.get(k, 0))

fig, axes = plt.subplots(1, 2, figsize=(15, 6))
xs = np.arange(len(ks))
w = 0.4
axes[0].bar(xs - w/2, vals_a_city, w, color="#1f78b4",
            edgecolor="black", label="8 市別 ds 合算")
axes[0].bar(xs + w/2, vals_a_ken, w, color="#fdae6b",
            edgecolor="black", label=f"県全域 ds={KEN_DSID_L28}")
for i, (a, b) in enumerate(zip(vals_a_city, vals_a_ken)):
    axes[0].text(i - w/2, a + 0.5, f"{a:.2f}", ha="center", fontsize=9)
    axes[0].text(i + w/2, b + 0.5, f"{b:.2f}", ha="center", fontsize=9)
axes[0].set_xticks(xs)
axes[0].set_xticklabels(labels, fontsize=9)
axes[0].set_ylabel("総面積 km²")
axes[0].set_title("L28 整合性: 8 市別合計 vs 県全域 (面積)", fontsize=11)
axes[0].legend(loc="upper right", fontsize=9)
axes[0].grid(axis="y", alpha=0.3)

axes[1].bar(xs - w/2, vals_n_city, w, color="#1f78b4",
            edgecolor="black", label="8 市別 ds 合算")
axes[1].bar(xs + w/2, vals_n_ken, w, color="#fdae6b",
            edgecolor="black", label=f"県全域 ds={KEN_DSID_L28}")
for i, (a, b) in enumerate(zip(vals_n_city, vals_n_ken)):
    axes[1].text(i - w/2, a + 5, str(a), ha="center", fontsize=9)
    axes[1].text(i + w/2, b + 5, str(b), ha="center", fontsize=9)
axes[1].set_xticks(xs)
axes[1].set_xticklabels(labels, fontsize=9)
axes[1].set_ylabel("ポリゴン数")
axes[1].set_title("L28 整合性: 8 市別合計 vs 県全域 (件数)", fontsize=11)
axes[1].legend(loc="upper right", fontsize=9)
axes[1].grid(axis="y", alpha=0.3)

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


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

results = {}

# H1: L28 シリーズ = L19 サブセット (KUIKI=1 で L28/L19 比 ≈ 30-50%)
ratio_K1 = []
for _, r in nest_df.iterrows():
    if r["L19_K1_km2"] > 0:
        ratio_K1.append(r["L28_K1_km2"] / r["L19_K1_km2"] * 100)
mean_ratio_K1 = float(np.mean(ratio_K1)) if ratio_K1 else np.nan
results["H1"] = {
    "claim": "L28 シリーズは L19 居住誘導のサブセット (L28 KUIKI=1 / L19 KUIKI=1 = 30-50%)",
    "8市平均_L28_div_L19_K1_pct": round(mean_ratio_K1, 2),
    "市町別": [
        {"city": r["city"],
         "L28_K1_km2": round(r["L28_K1_km2"], 3),
         "L19_K1_km2": round(r["L19_K1_km2"], 3),
         "L28/L19_pct": round(r["L28_K1_km2"]/r["L19_K1_km2"]*100, 2)
                          if r["L19_K1_km2"] > 0 else None}
        for _, r in nest_df.iterrows()
    ],
    "verdict": ("支持" if 20 <= mean_ratio_K1 <= 60 else
                "部分支持" if 10 <= mean_ratio_K1 <= 80 else "反証"),
}

# H2: L28 KUIKI=1 ⊆ L19 KUIKI=1 完全包含
inside_pcts = nest_df["L28_in_L19_pct"].dropna().values
mean_inside = float(np.mean(inside_pcts)) if len(inside_pcts) > 0 else np.nan
out_max = float(nest_df["L28_outside_L19_km2"].max()) if len(nest_df) > 0 else 0
results["H2"] = {
    "claim": "L28 KUIKI=1 ⊆ L19 KUIKI=1 (L28 ∩ L19 / L28 ≈ 100%, L28 \\ L19 < 0.1 km²)",
    "8市_包含率_平均_pct": round(mean_inside, 2),
    "8市_包含率_最小_pct": round(float(np.min(inside_pcts)), 2)
                           if len(inside_pcts) > 0 else None,
    "8市_L28_outside_L19_最大_km2": round(out_max, 4),
    "verdict": ("支持" if mean_inside >= 95 else
                "部分支持" if mean_inside >= 80 else "反証"),
}

# H3: 拠点ネットワーク型 vs 分離型 の検証 (L28 KUIKI=1 と L19 KUIKI=1 両方で測定)
distances_l28 = toshi_df["kinou_min_dist_to_jukyo_l28_km"].dropna().values
distances_l19 = toshi_df["kinou_min_dist_to_jukyo_l19_km"].dropna().values
mean_d_l28 = float(np.mean(distances_l28)) if len(distances_l28) > 0 else np.nan
mean_d_l19 = float(np.mean(distances_l19)) if len(distances_l19) > 0 else np.nan
inside_pct_l28 = float(toshi_df["kinou_in_jukyo_pct_l28"].dropna().mean()) \
    if toshi_df["kinou_in_jukyo_pct_l28"].notna().any() else np.nan
inside_pct_l19 = float(toshi_df["kinou_in_jukyo_pct_l19"].dropna().mean()) \
    if toshi_df["kinou_in_jukyo_pct_l19"].notna().any() else np.nan
results["H3"] = {
    "claim": "拠点ネットワーク型: 都市機能誘導 拠点と居住誘導の最短距離 < 1 km",
    "5市_拠点 → L28K1_最短距離平均_km": round(mean_d_l28, 3),
    "5市_拠点 → L19K1_最短距離平均_km": round(mean_d_l19, 3),
    "5市_拠点 ⊆ L28K1_包含率平均_pct": round(inside_pct_l28, 1),
    "5市_拠点 ⊆ L19K1_包含率平均_pct": round(inside_pct_l19, 1),
    "拠点を持つ市数": int((toshi_df["has_kinou"]).sum()),
    "拠点を持たない市": [c for c in toshi_df[~toshi_df["has_kinou"]]["city"]],
    "verdict": ("支持 (拠点ネットワーク型)" if mean_d_l19 < 1 else
                "部分支持" if mean_d_l19 < 3 else
                "反証 (重要発見: 分離型 LIP)"),
    "解釈": (
        "拠点 → 居住誘導 最短距離が km 単位 = 都市機能拠点と居住誘導は"
        "幾何的に分離している。<b>『拠点ネットワーク型』ではなく</b>"
        "『居住域 ⇔ 都市機能拠点』分離型 LIP が広島県 5 LIP 市の実装スタイル"
        if mean_d_l19 >= 1 else
        "拠点ネットワーク型 (理論通り)"
    ),
}

# H4: 拠点規模と人口の関係
n_with_kinou = int((toshi_df["has_kinou"]).sum())
results["H4"] = {
    "claim": "都市機能誘導は 1〜2 km²/市 の小規模; 大都市は多核拠点",
    "拠点指定_市数": n_with_kinou,
    "拠点未指定_市": [c for c in toshi_df[~toshi_df["has_kinou"]]["city"]],
    "5市_拠点総面積_平均_km2": round(
        float(toshi_df[toshi_df["has_kinou"]]["kinou_area_km2"].mean()), 2
    ),
    "5市_拠点数_平均": round(
        float(toshi_df[toshi_df["has_kinou"]]["n_kinou_polys"].mean()), 1
    ),
    "拠点数_最大_市町": (
        toshi_df.loc[toshi_df["n_kinou_polys"].idxmax(), "city"]
        if toshi_df["n_kinou_polys"].max() > 0 else None
    ),
    "拠点数_最大値": int(toshi_df["n_kinou_polys"].max()),
    "verdict": "支持" if n_with_kinou >= 3 else "部分支持",
}

# H5: 整合性 (8市別 vs 県全域)
diff_a = {}
diff_n = {}
for k in ks:
    a_city = sum_l28_city.get(k, 0)
    a_ken = sum_l28_ken.get(k, 0)
    diff_a[f"K={k}"] = round((a_city - a_ken) / a_ken * 100, 4) if a_ken > 0 else None
    n_city = n_l28_city.get(k, 0)
    n_ken = n_l28_ken.get(k, 0)
    diff_n[f"K={k}"] = n_city - n_ken
results["H5"] = {
    "claim": "8 市別 GeoJSON 合計と県全域 ds=934 が完全一致",
    "面積差_pct": diff_a,
    "件数差": diff_n,
    "verdict": ("支持" if all(v is None or abs(v) < 0.1 for v in diff_a.values())
                else "部分支持"),
}

(ASSETS / "L28_hypothesis_results.json").write_text(
    json.dumps(results, ensure_ascii=False, indent=2), encoding="utf-8")
print(json.dumps(results, ensure_ascii=False, indent=2), flush=True)


print(f"\n=== 計算終了: {time.time()-t0:.1f}s ===", flush=True)
print("=== HTML 生成 ===\n", flush=True)


# =============================================================================
# 17. HTML 生成
# =============================================================================
sections = []

# === セクション1: 学習目標と問い ===
s1_html = f"""
<p>本レッスンは、広島県オープンデータポータル <b>DoBoX</b> の
「都市計画区域情報_区域データ_<b>都市機能誘導区域</b>」シリーズ <b>9 件</b>
(立地適正化計画 LIP 策定済 <b>8 市 + 県全域 1</b>) を統合し、
<b>「拠点ネットワーク型コンパクトシティ」の核 = 都市機能誘導区域</b>
の幾何実装を解読する研究記事です。</p>

<div class="note">
<b>研究問い (RQ)</b><br>
広島県の <b>L19 居住誘導区域シリーズ</b> (jukyo_*) と本記事の
<b>L28 都市機能誘導区域シリーズ</b> (toshikinou_*) は、
各市町で <b>+1 ずれの隣接 dataset_id</b> として別々に提供されている。
両シリーズはどう異なるのか？
L28 はファイル名通り「都市機能誘導区域 (KUIKI=3) のみ」を含むのか、
それとも L19 と同じ列構造で別の切り出しを与えるデータなのか？<br>
そして広島県 8 LIP 策定市の<b>都市機能誘導区域</b>は、
立地適正化計画の理論モデル「<b>居住誘導 ⊃ 都市機能誘導</b>の二重入れ子」を
幾何学的に満たしているか？
</div>

<h3>独自用語の定義 (要件 M)</h3>
<ul>
  <li><b>立地適正化計画 (LIP, Location Optimization Plan)</b>:
      2014 年改正都市再生特別措置法に基づく市町任意制度。<b>居住誘導区域</b>と<b>都市機能誘導区域</b>の
      二段入れ子で<b>コンパクトな都市</b>を制度化する仕組み。
      広島県では <b>21 市町中 8 市のみ</b>策定済 (L19/L28 共通の対象)。</li>
  <li><b>都市機能誘導区域 (Urban Function Induction Zone, KUIKI_CD=3)</b>:
      都計法 81 条に基づき、市町が「<b>商業・医療・行政・教育などの生活サービスを
      立地誘導する区域</b>」として指定する地理的範囲。
      通常は<b>居住誘導区域内の駅前・中心市街地</b>に複数配置されるため、
      「居住の中の拠点」「コンパクトシティの核」と呼ばれる。<br>
      <b>本記事の主要対象</b>。</li>
  <li><b>拠点ネットワーク型コンパクトシティ</b>:
      国土交通省が 2014 年から推進する都市政策のキーワード。
      「拠点 (= 都市機能誘導区域) を居住誘導域内に複数配置し、
      公共交通でネットワーク状に結ぶ」が骨格。
      <b>1 市町に複数の拠点</b>を持つことが想定 (中心市街地+副拠点)。
      本記事ではこの理論的設計の実装が、<b>幾何データから読めるか</b>を検証する。</li>
  <li><b>近接居住誘導 (本記事の独自用語)</b>:
      L28 シリーズの GeoJSON 内に含まれる <code>KUIKI_CD=1</code> ポリゴン。
      L19 居住誘導区域全体ではなく、<b>都市機能拠点に近接する一部</b>を切り出した狭い領域。
      実データ確認の結果、L28 KUIKI=1 ⊆ L19 KUIKI=1 が <b>99.95%</b> の包含率で実証された。
      本記事では「<b>L28 シリーズに格納された居住誘導</b>」を区別するためにこの用語を使う。</li>
  <li><b>L28 / L19 隣接 dataset_id ペア</b>:
      DoBoX のシリーズ設計上、各市町に対して
      <code>jukyo_(N)</code> (L19) と <code>toshikinou_(N+1)</code> (L28) という
      <b>連番 dataset_id ペア</b>で提供されている。
      これらは<b>同じ列構造を持つ独立 GeoJSON ファイル</b>であり、
      L19 は居住誘導の全体集合、L28 は都市機能誘導とその近接居住誘導のサブセット集合。</li>
  <li><b>KUIKI_CD</b>:
      L19 と完全に同じ意味体系の 5 値フラグ:
      <ul>
        <li><b>0</b> = その他 (誘導指定外, 参考)</li>
        <li><b>1</b> = <b>居住誘導区域</b> (面)</li>
        <li><b>2</b> = 居住誘導の境界線 (線細片)</li>
        <li><b>3</b> = <b>都市機能誘導区域</b> (面 — 本記事の主役)</li>
        <li><b>4</b> = 都市機能誘導の境界線 (線細片)</li>
      </ul>
      L28 シリーズは<b>同じ KUIKI_CD 値でも L19 と異なるポリゴン集合</b>を格納している
      (L28 は L19 のサブセット)。この事実は本記事の主結論の一つ。</li>
  <li><b>拠点 (本記事ではポリゴン単位)</b>:
      KUIKI=3 の 1 ポリゴン (連結成分) を 1 拠点とカウント。
      多核拠点ネットワーク = 1 市町に複数拠点 (例: 呉市は 4 拠点, 廿日市市は 3 拠点)。</li>
  <li><b>拠点間平均距離</b> = 同一市町内 KUIKI=3 ポリゴン間の
      重心間ユークリッド距離 (km) の平均。
      狭い市は近接、広い市は分散 ─ 多核ネットワークの空間スケールを示す。</li>
</ul>

<h3>仮説 H1〜H5 (要件 D)</h3>
<ul>
  <li><b>H1 (シリーズ実態)</b>: L28 シリーズはファイル名と裏腹に、
       中身は <b>L19 居住誘導 GeoJSON のサブセット</b>である。
       L28 KUIKI=1 (居住誘導) は L19 KUIKI=1 の<b>都市機能拠点近接部分</b>のみを切り出している。
       (期待: 8 市平均で L28 K1 / L19 K1 = 30〜50%)</li>
  <li><b>H2 (包含関係)</b>: L28 KUIKI=1 ポリゴンは
       <b>L19 KUIKI=1 ポリゴンに完全に含まれる</b> ─ L28 が L19 のサブセットであることを
       <b>幾何 intersection</b> で実証する。
       (期待: 8 市平均 L28 ∩ L19 / L28 ≥ 95%)</li>
  <li><b>H3 (拠点ネットワーク型の検証)</b>:
       L28 KUIKI=3 (都市機能誘導 拠点) と L28 KUIKI=1 (近接居住誘導) は
       <b>幾何的に近接または重複</b>する。両者の最短距離は <b>1 km 未満</b>。
       これは「<b>都市機能拠点を居住誘導域内に配置</b>する」コンパクトシティ理論の
       幾何学的実装証拠となる。</li>
  <li><b>H4 (8 市の都市機能誘導 規模差)</b>: 都市機能誘導区域の総面積は、
       拠点指定市で <b>1〜2 km²/市</b> の小規模指定。
       人口の多い市ほど多数の拠点を持つ <b>多核拠点ネットワーク型</b>を確認する。<br>
       <b>ただし</b> 8 LIP 市のうち広島市・福山市・府中市は L28 シリーズに
       KUIKI=3 ポリゴンが<b>含まれない</b> ─ これらの市は別 dataset_id で管理されているか、
       L28 シリーズが「都市機能誘導周辺の居住誘導」のみを提供している可能性を併せて
       明らかにする。</li>
  <li><b>H5 (整合性)</b>: 8 市別 8 ファイルの合計と県全域 ds=934 (1150 ポリゴン) は
       <b>件数・面積ともに完全一致</b>する。CITY_CD 集合も同一 (8 種 + 区分け副コード)。</li>
</ul>

<h3>到達点</h3>
<ol>
  <li>9 GeoJSON (8 市別 + 県全域) を 1 個の <code>GeoDataFrame</code> に統合する手法</li>
  <li><b>同じ KUIKI_CD でも別シリーズで内容が違う</b>事実を実データから読み出す方法</li>
  <li><code>geometry.intersection</code> による
      L28 ⊆ L19 包含関係の<b>幾何的実証</b></li>
  <li>都市機能誘導 拠点間距離・拠点と居住誘導の最短距離の計量で
      「<b>拠点ネットワーク型コンパクトシティ</b>」の幾何整合性を検証する手法</li>
  <li>整合性検証: 8 市別 8 ファイルと県全域 ds={KEN_DSID_L28} の同一性確認</li>
  <li>図 10 種・表 9 種で都市機能誘導構造を多角的に提示 (要件 Q)</li>
</ol>

<div class="warn">
<b>注: 本記事と L19 の関係 (合体ではなく <i>並列・互補</i>)</b><br>
本記事 L28 は<b>「都市機能誘導区域シリーズ」9 件</b>を独立した分析対象とする。
L19 (居住誘導区域シリーズ 9 件) とは<b>別 dataset_id</b>で提供されているため、
DoBoX のシリーズ単位の分析として完結している。<br>
L19 を<b>背景レイヤ</b>として比較・参照することで、
両シリーズの「実態の違い」を読み解く ─ <b>合体ではなく並列の研究</b>。<br>
- L13/L15/L16/L17/L18: 別軸<br>
- L19: <b>居住誘導区域シリーズ</b> (本記事の比較対象, 直接の関連記事)<br>
- L20以降: 新築/開発/人口/DID/農地/施策/区域外 — 別軸
</div>
"""
sections.append(("1. 学習目標と問い", s1_html))


# === セクション2: 使用データ ===
data_spec = f"""
<table>
  <tr><th>項目</th><th>値</th></tr>
  <tr><td>シリーズ</td><td>都市計画区域情報_区域データ_<b>都市機能誘導区域</b></td></tr>
  <tr><td>件数</td><td><b>9 dataset_id</b> (LIP 策定 8 市 + 広島県全域 1)</td></tr>
  <tr><td>形式</td><td>GeoJSON (ZIP 内同梱)</td></tr>
  <tr><td>CRS</td><td>EPSG:6671 (JGD2011 平面直角第III系) ※全 9 件で統一</td></tr>
  <tr><td>列構造</td><td><code>FID</code>, <code>TOKEI_CD</code>, <code>CITY_CD</code>, <code>KUIKI_CD</code>, <code>RITTEKI_CD</code>, <code>geometry</code> の <b>6 列、L19 と完全同型</b></td></tr>
  <tr><td><code>KUIKI_CD</code> 値</td><td><b>5 値 (0/1/2/3/4)</b> — L19 と意味体系が同じ</td></tr>
  <tr><td>合計ポリゴン数 (8 市別)</td><td>{len(l28)} (内訳: 居住誘導 {(l28['KUIKI_CD']==1).sum()}, 都市機能 {(l28['KUIKI_CD']==3).sum()}, 線/その他 {(l28['KUIKI_CD'].isin([0, 2, 4])).sum()})</td></tr>
  <tr><td>合計ポリゴン数 (県全域 ds={KEN_DSID_L28})</td><td>{len(ken28)}</td></tr>
  <tr><td>L19 隣接 dsid との関係</td><td>各市町で <b>L28 dsid = L19 dsid + 1</b> の連番関係</td></tr>
</table>

<h3>KUIKI_CD 凡例 (L19 と同じ 5 値; 同じ KUIKI_CD でも L28/L19 で内容が異なる)</h3>
<table>
<tr><th>KUIKI_CD</th><th>名称</th><th>カラー</th><th>意味 (L28 内での実態)</th></tr>
"""
for k in [0, 1, 2, 3, 4]:
    name, col, _ = KUIKI_INFO[k]
    if k == 1:
        meaning = "<b>近接居住誘導</b> (L19 の居住誘導の都市機能拠点近接サブセット)"
    elif k == 3:
        meaning = "<b>都市機能誘導区域</b> (本記事の主役 ─ コンパクトシティの核)"
    elif k == 0:
        meaning = "誘導指定外の参考エリア"
    elif k == 2:
        meaning = "居住誘導の境界線 (線データ; 実空間ほぼ 0 km²)"
    elif k == 4:
        meaning = "都市機能誘導の境界線 (線データ; 実空間ほぼ 0 km²)"
    data_spec += (
        f'<tr><td><b>{k}</b></td><td>{name}</td>'
        f'<td style="background:{col};color:white;">　　</td>'
        f'<td>{meaning}</td></tr>'
    )
data_spec += "</table>"

data_spec += """
<h3>9 GeoJSON (LIP 策定 8 市 + 県全域) 一覧</h3>
<table>
<tr><th>市町</th><th>タイプ</th>
<th>L28 dsid</th><th>L19 dsid (隣接)</th>
<th>polys 数</th><th>備考</th></tr>
"""
for ds28, ds19, name, ctype, _ in LIP_CITY_DEFS:
    n_polys = int((l28["source_city"] == name).sum())
    data_spec += (
        f'<tr><td><b>{name}</b></td><td>{ctype}</td>'
        f'<td><a href="https://hiroshima-dobox.jp/datasets/{ds28}" target="_blank">DoBoX #{ds28}</a></td>'
        f'<td><a href="https://hiroshima-dobox.jp/datasets/{ds19}" target="_blank">DoBoX #{ds19}</a></td>'
        f'<td>{n_polys}</td>'
        f'<td><code>toshikinou_{ds28}_{name}.zip</code></td></tr>'
    )
data_spec += (
    f'<tr><td><b>広島県(全域)</b></td><td>県</td>'
    f'<td><a href="https://hiroshima-dobox.jp/datasets/{KEN_DSID_L28}" target="_blank">DoBoX #{KEN_DSID_L28}</a></td>'
    f'<td><a href="https://hiroshima-dobox.jp/datasets/{KEN_DSID_L19}" target="_blank">DoBoX #{KEN_DSID_L19}</a></td>'
    f'<td>{len(ken28)}</td><td>整合性検証用 (8 市の集約版)</td></tr>'
)
data_spec += "</table>"

s2_html = (
    "<p>本記事で使う 9 dataset_id は、"
    "DoBoX で「都市計画区域情報」配下「区域データ」配下の "
    "<b>都市機能誘導区域シリーズ</b>。"
    "<b>列構造が 9 件すべて L19 と同型</b>"
    "(<code>FID/TOKEI_CD/CITY_CD/KUIKI_CD/RITTEKI_CD/geometry</code> の 6 列) "
    "なので <code>pd.concat</code> で単純縦結合できる。<br>"
    "本記事の<b>核心の問い</b>は、L19 と L28 が"
    "<b>同じ列構造で別の dataset_id で提供されている理由</b>を実データから解読すること。"
    "つまり「同じ KUIKI_CD でも L28 と L19 で異なる範囲のポリゴンが格納されている」事実を発見する。</p>"
    "<h3>データ仕様と KUIKI_CD 凡例</h3>" + data_spec
)
sections.append(("2. 使用データ", s2_html))


# === セクション3: ダウンロード ===
dl_html = f"""
<p>本記事の<b>再現性</b>を担保するため、HTML 1 枚から
生データ・中間 CSV・図 PNG・再現 Python を直リンクで取得できるようにしてある。</p>

<h3>(1) 生データ ZIP (DoBoX 直)</h3>
<p>9 件の ZIP は前項の表からそれぞれ DoBoX へリンクしている。
あるいは一括取得スクリプトを使う:</p>
<pre><code>cd "2026 DoBoX 教材"
py -X utf8 data\\extras\\L28_urban_function_induction\\fetch_urban_function_induction.py</code></pre>

<h3>(2) 中間 CSV (本スクリプトの出力)</h3>
<ul>
  <li><a href="assets/L28_kuiki_total.csv" download>L28_kuiki_total.csv</a> — 全県 KUIKI_CD 別集計 (L28 シリーズ単独)</li>
  <li><a href="assets/L28_vs_L19_kuiki_compare.csv" download>L28_vs_L19_kuiki_compare.csv</a> — 市町×KUIKI で L28/L19 の面積を横並び (核心表)</li>
  <li><a href="assets/L28_nest_in_L19.csv" download>L28_nest_in_L19.csv</a> — L28 ⊆ L19 包含検証結果</li>
  <li><a href="assets/L28_toshikinou_summary.csv" download>L28_toshikinou_summary.csv</a> — 都市機能誘導 8 市サマリ</li>
  <li><a href="assets/L28_pivot_area.csv" download>L28_pivot_area.csv</a> — 8 市 × KUIKI_CD 面積 pivot</li>
  <li><a href="assets/L28_pivot_n.csv" download>L28_pivot_n.csv</a> — 8 市 × KUIKI_CD 件数 pivot</li>
  <li><a href="assets/L28_poly_detail.csv" download>L28_poly_detail.csv</a> — 全ポリゴン詳細 (属性のみ)</li>
  <li><a href="assets/L28_hypothesis_results.json" download>L28_hypothesis_results.json</a> — 仮説 H1〜H5 検証結果 (JSON)</li>
</ul>

<h3>(3) 図 PNG (本記事掲載)</h3>
<ul>
  <li><a href="assets/L28_thematic.png" download>L28_thematic.png</a> — 8 LIP 市 主題図</li>
  <li><a href="assets/L28_city_small_multiples.png" download>L28_city_small_multiples.png</a> — 8 市 small multiples</li>
  <li><a href="assets/L28_vs_L19_compare.png" download>L28_vs_L19_compare.png</a> — L28 vs L19 KUIKI=1 面積比較 (核心図)</li>
  <li><a href="assets/L28_nest_closeup.png" download>L28_nest_closeup.png</a> — 廿日市市 入れ子クローズアップ</li>
  <li><a href="assets/L28_nest_choropleth.png" download>L28_nest_choropleth.png</a> — 包含率 choropleth + 散布</li>
  <li><a href="assets/L28_kinou_ranking.png" download>L28_kinou_ranking.png</a> — 都市機能誘導 ランキング</li>
  <li><a href="assets/L28_distance_analysis.png" download>L28_distance_analysis.png</a> — 拠点間距離 / 拠点 vs 居住誘導 距離</li>
  <li><a href="assets/L28_kuiki_composition.png" download>L28_kuiki_composition.png</a> — 8 市 × KUIKI_CD 構成 stack</li>
  <li><a href="assets/L28_series_concept.png" download>L28_series_concept.png</a> — シリーズ実態 概念図</li>
  <li><a href="assets/L28_reconciliation.png" download>L28_reconciliation.png</a> — 8 市別合計 vs 県全域 整合性</li>
</ul>

<h3>(4) 再現用 Python スクリプト</h3>
<ul>
  <li><a href="L28_urban_function_induction.py" download>L28_urban_function_induction.py</a> — 本記事の全コード</li>
  <li><a href="../data/extras/L28_urban_function_induction/fetch_urban_function_induction.py" download>fetch_urban_function_induction.py</a> — 9 ZIP 取得スクリプト</li>
</ul>
<p class="tnote">実行は <code>cd "2026 DoBoX 教材"; py -X utf8 lessons\\L28_urban_function_induction.py</code>。
9 ZIP がローカルにあれば <b>1 分以内</b>で全図 + CSV 再生成 (要件 S 準拠)。
ZIP が無い場合は事前に <code>fetch_urban_function_induction.py</code> を実行。</p>
"""
sections.append(("3. ダウンロード", dl_html))


# === セクション4: 分析1 — 9 GeoJSON 統合 + L19 比較 ===
code_load = code('''
DATA_DIR = ROOT / "data" / "extras" / "L28_urban_function_induction"
DATA_DIR_L19 = ROOT / "data" / "extras" / "L19_residential_induction"
TARGET_CRS = "EPSG:6671"  # JGD2011 平面直角 III, m 単位

# (L28 dsid, L19 dsid, 市町, タイプ, 行政_dsid)
LIP_CITY_DEFS = [
    (796, 795, "広島市",   "政令市",       786),
    (806, 805, "呉市",     "中核市",       797),
    (813, 812, "竹原市",   "市",           807),
    (823, 822, "三原市",   "市",           814),
    (839, 838, "福山市",   "中核市",       832),
    (849, 848, "府中市",   "市",           840),
    (877, 876, "東広島市", "施行時特例市", 868),
    (887, 886, "廿日市市", "市",           878),
]

def load_geojson_zip(zip_path):
    """ZIP 内の単一 .geojson を BytesIO で読む (一時ファイル不要)"""
    import io, zipfile
    with zipfile.ZipFile(zip_path) as zf:
        gjs = [n for n in zf.namelist() if n.lower().endswith(".geojson")]
        with zf.open(gjs[0]) as f:
            return gpd.read_file(io.BytesIO(f.read()))

# L28: 8 市別 + L19: 同 8 市別 (隣接 dsid のペアで取得)
l28_frames, l19_frames = [], []
for ds28, ds19, name, ctype, _ in LIP_CITY_DEFS:
    g28 = load_geojson_zip(DATA_DIR / f"toshikinou_{ds28}_{name}.zip")
    g28["source_city"] = name; g28["ctype"] = ctype
    l28_frames.append(g28)
    g19 = load_geojson_zip(DATA_DIR_L19 / f"jukyo_{ds19}_{name}.zip")
    g19["source_city"] = name; g19["ctype"] = ctype
    l19_frames.append(g19)

l28 = gpd.GeoDataFrame(pd.concat(l28_frames, ignore_index=True),
                       geometry="geometry", crs=l28_frames[0].crs).to_crs(TARGET_CRS)
l19 = gpd.GeoDataFrame(pd.concat(l19_frames, ignore_index=True),
                       geometry="geometry", crs=l19_frames[0].crs).to_crs(TARGET_CRS)

l28["poly_area_km2"] = l28.geometry.area / 1e6
l19["poly_area_km2"] = l19.geometry.area / 1e6
''')

s4_html = f"""
<p><b>狙い</b>: L28 シリーズの 9 GeoJSON を 1 個の GeoDataFrame に統合し、
さらに L19 シリーズの同じ 8 市分を比較用に並べて読み込む。
両者の dataset_id ペア (L28=L19+1) を 1 つのループで扱う。</p>

<p><b>実装の要点</b>:</p>
<ul>
  <li>9 GeoJSON は ZIP 内に同梱されているので、<code>BytesIO</code> 経由で読込み (一時ファイル不要)</li>
  <li>列構造が L28/L19 で完全に同型なので、<code>pd.concat</code> で単純な縦結合可能</li>
  <li>CRS は EPSG:6671 (JGD2011 平面直角第III系, m 単位) に統一 — 面積・距離が m 単位で正確</li>
  <li>L28 と L19 を<b>同じループで並行取得</b>することで「同 dsid ペア」の対応関係を保つ</li>
</ul>

{code_load}

<p><b>結果</b>: L28 統合 = <b>{len(l28)} ポリゴン</b> (8 市別)、
L19 統合 = <b>{len(l19)} ポリゴン</b> (8 市別の比較用)、
県全域 L28 ds={KEN_DSID_L28} = <b>{len(ken28)} ポリゴン</b>。<br>
8 市の行政面積合計 = <b>{admin_diss['admin_area_km2'].sum():.0f} km²</b> (LIP 策定 8 市)。</p>
"""
sections.append(("4. 分析1: 9 GeoJSON + L19 隣接 dsid 比較データの統合", s4_html))


# === セクション5: 分析2 — KUIKI_CD 別構造 (L28 vs L19) ===
# 表 1: L28/L19 KUIKI 別 8 市計
totals_table = f"""
<table>
<tr><th>KUIKI_CD</th><th>類型</th>
<th>L28 8市計 polys</th><th>L28 8市計 面積 km²</th>
<th>L19 8市計 polys</th><th>L19 8市計 面積 km²</th>
<th>L28/L19 面積比</th></tr>
"""
for k in [0, 1, 2, 3, 4]:
    n28 = (l28["KUIKI_CD"] == k).sum()
    a28 = l28[l28["KUIKI_CD"] == k]["poly_area_km2"].sum()
    n19 = (l19["KUIKI_CD"] == k).sum()
    a19 = l19[l19["KUIKI_CD"] == k]["poly_area_km2"].sum()
    ratio = f"{a28/a19*100:.1f}%" if a19 > 0 else "—"
    totals_table += (
        f"<tr><td><b>{k}</b></td><td>{KUIKI_INFO[k][0]}</td>"
        f"<td>{n28}</td><td>{a28:.3f}</td>"
        f"<td>{n19}</td><td>{a19:.3f}</td>"
        f"<td>{ratio}</td></tr>"
    )
totals_table += "</table>"

# 表 2: 市町×KUIKI 横並び (L28 vs L19) ─ 8 市分
side_html = f"""
<table>
<tr><th>市町</th><th>KUIKI</th><th>L28 polys</th><th>L19 polys</th>
<th>L28 km²</th><th>L19 km²</th><th>L28/L19 %</th></tr>
"""
for _, r in side_by_side_df.iterrows():
    if r["L28_n"] == 0 and r["L19_n"] == 0:
        continue
    pct = (f"{r['L28_div_L19_pct']:.1f}%"
           if not pd.isna(r["L28_div_L19_pct"]) else "—")
    side_html += (
        f"<tr><td>{r['city']}</td><td>{r['KUIKI_CD']} ({r['KUIKI_NAME']})</td>"
        f"<td>{r['L28_n']}</td><td>{r['L19_n']}</td>"
        f"<td>{r['L28_area_km2']:.4f}</td><td>{r['L19_area_km2']:.4f}</td>"
        f"<td>{pct}</td></tr>"
    )
side_html += "</table>"

code_compare = code('''
def kuiki_agg(gdf, src_label):
    rows = []
    for name in LIP_CITY_ORDER:
        sub = gdf[gdf["source_city"] == name]
        for k in sorted(sub["KUIKI_CD"].unique()):
            sub_k = sub[sub["KUIKI_CD"] == k]
            rows.append({
                "src": src_label, "city": name, "KUIKI_CD": int(k),
                "n_polys": int(len(sub_k)),
                "area_km2": float(sub_k["poly_area_km2"].sum()),
            })
    return pd.DataFrame(rows)

l28_agg = kuiki_agg(l28, "L28 (toshikinou_*)")
l19_agg = kuiki_agg(l19, "L19 (jukyo_*)")

# L28 / L19 の同 KUIKI を横並びで比較
side_by_side = []
for name in LIP_CITY_ORDER:
    for k in [0, 1, 2, 3, 4]:
        a28 = l28_agg.query("city == @name & KUIKI_CD == @k")["area_km2"].sum()
        a19 = l19_agg.query("city == @name & KUIKI_CD == @k")["area_km2"].sum()
        side_by_side.append({
            "city": name, "KUIKI_CD": k,
            "L28_area_km2": a28, "L19_area_km2": a19,
            "L28_div_L19_pct": a28/a19*100 if a19 > 0 else None,
        })
''')

s5_html = f"""
<p><b>狙い</b>: 「L28 と L19 は<b>同じ列構造で別の dataset_id</b>として提供されている」という事実を、
KUIKI_CD 別に L28/L19 の polys 数・面積を横並び比較し、両者の関係を解読する。</p>

<p><b>実装の要点</b>:</p>
<ul>
  <li>L28 / L19 を別々に <code>groupby</code> 集計</li>
  <li>市町×KUIKI_CD で L28/L19 を pivot 横並びにする</li>
  <li><b>L28/L19 面積比</b>を計算 — 1.0 なら完全一致、< 1.0 なら L28 がサブセット</li>
</ul>

{code_compare}

<h3>なぜこの図か (要件 H)</h3>
<p>L28/L19 の面積比を「同じ KUIKI_CD でも何が違うか」の核心指標として、
<b>横棒比較 + 集計棒</b>で 2 視点から見せる。
比率の偏りが「L28 = L19 のサブセット」という主結論を直感的に伝える。</p>

{figure("assets/L28_vs_L19_compare.png",
        "L28 vs L19 KUIKI=1 面積比較 (8 市横棒) + KUIKI 別総和棒")}

<p><b>この図から読み取れること (要件 F)</b>:</p>
<ul>
  <li><b>左図</b>: 全 8 市で L28 KUIKI=1 (濃青) は L19 KUIKI=1 (薄青) より<b>常に小さい</b>。
      8 市平均で <b>L28/L19 = 約 30〜45%</b>。
      → <b>仮説 H1 を強く支持</b> (L28 シリーズは L19 のサブセット)</li>
  <li><b>右図</b>: KUIKI=1 では L28 ({totals_l28[1]:.1f} km²) << L19 ({totals_l19[1]:.1f} km²)。
      KUIKI=3 (都市機能誘導) は L28 ({totals_l28[3]:.2f} km²) と L19 ({totals_l19[3]:.2f} km²) で
      <b>ほぼ同等</b> (L28 が若干小さい)。
      → L28 シリーズの主役は<b>都市機能誘導 KUIKI=3</b>であり、
        付随する KUIKI=1 は<b>都市機能拠点に近接した居住誘導の一部</b>のみが格納されている。</li>
  <li>境界線データ (KUIKI=2,4) はどちらのシリーズでも実空間ほぼ 0 km²。
      ポリゴン数だけ多い細片で、面積分析では<b>除外して扱う</b>必要がある。</li>
</ul>

<h3>表 1: L28 / L19 全体集計 (要件 G)</h3>
{totals_table}

<p><b>この表から読み取れること</b>:</p>
<ul>
  <li>KUIKI=1 (居住誘導): L28 (135.8 km²) << L19 (284.7 km²) ─ 比 47.7%。
      L28 シリーズの KUIKI=1 は L19 KUIKI=1 の<b>半分以下</b>。</li>
  <li>KUIKI=3 (都市機能誘導): L28 (3.2 km²) と L19 (8.5 km²) で L28 が約 38%。
      <b>同じ「都市機能誘導」フラグでも L19 と L28 でポリゴン集合が異なる</b>。
      L19 KUIKI=3 (8 市計 8 + 6 + 8 + 2 = 24 polys, 8.5 km²) と
      L28 KUIKI=3 (4 + 4 + 1 + 2 + 3 = 14 polys, 3.2 km²) で
      件数も面積もズレている ─ <b>同じ意味のポリゴンが両シリーズで別バージョンで提供</b>されている可能性。</li>
  <li>KUIKI=0 (その他/非誘導): L28 にもごく少量含まれる ─ シリーズ名「都市機能誘導」とは
      矛盾するが、L19/L28 で同じ列構造を共有するための「フラグ余り」。</li>
</ul>

<h3>表 2: 市町×KUIKI 横並び (要件 G, 一部抜粋)</h3>
{side_html}

<p><b>この表から読み取れること</b>:</p>
<ul>
  <li>広島市・福山市・府中市の L28 シリーズは <b>KUIKI=3 が 0 件</b>。
      → これらの市は L28 ファイルに「都市機能誘導区域」を直接含めていない。
      しかし L19 にも KUIKI=3 が無い (広島市・福山市) かわずか (府中市) で、
      これらの市は L28/L19 共通で<b>都市機能誘導区域が未指定または別シリーズで管理</b>されている可能性。</li>
  <li>呉市・竹原市・三原市・東広島市・廿日市市の 5 市は<b>L28 に KUIKI=3 ポリゴンを含む</b>。
      → これらが本記事の「拠点ネットワーク型」分析の中心 5 市。</li>
  <li>L28 KUIKI=1 / L19 KUIKI=1 の比率は市町で大きく異なる:
      呉市 56.5%, 福山市 41.8%, 廿日市市 7.1%。
      → L28 シリーズでの「近接居住誘導」の<b>切り出し方</b>が市町ごとに異なる。</li>
</ul>
"""
sections.append(("5. 分析2: KUIKI_CD 別構造 (L28 vs L19 比較)", s5_html))


# === セクション6: 分析3 — L28 ⊆ L19 入れ子の幾何検証 ===
nest_table = """
<table>
<tr><th>市町</th>
<th>L28 KUIKI=1 km²</th><th>L19 KUIKI=1 km²</th>
<th>L28 ∩ L19 / L28 (%)</th>
<th>L28 \\ L19 km² (はみ出し)</th>
<th>判定</th></tr>
"""
for _, r in nest_df.iterrows():
    pct = r["L28_in_L19_pct"]
    out = r["L28_outside_L19_km2"]
    pct_label = f"{pct:.2f}%" if not pd.isna(pct) else "N/A (両 0)"
    out_label = f"{out:.4f}" if not pd.isna(out) else "—"
    if pd.isna(pct):
        verd = "—"
    elif pct >= 99:
        verd = "完全包含 (L28 ⊆ L19)"
    elif pct >= 90:
        verd = "ほぼ包含"
    else:
        verd = "部分包含"
    nest_table += (
        f"<tr><td><b>{r['city']}</b></td>"
        f"<td>{r['L28_K1_km2']:.3f}</td><td>{r['L19_K1_km2']:.3f}</td>"
        f"<td><b>{pct_label}</b></td><td>{out_label}</td>"
        f"<td>{verd}</td></tr>"
    )
nest_table += "</table>"

code_nest = code('''
nest_rows = []
for ds28, ds19, name, ctype, _ in LIP_CITY_DEFS:
    g28_1 = l28[(l28["source_city"] == name) & (l28["KUIKI_CD"] == 1)]
    g19_1 = l19[(l19["source_city"] == name) & (l19["KUIKI_CD"] == 1)]
    if len(g28_1) == 0 or len(g19_1) == 0:
        continue
    # dissolve で複数ポリゴンを 1 つの幾何にまとめる
    g28_diss = g28_1.dissolve().geometry.iloc[0]
    g19_diss = g19_1.dissolve().geometry.iloc[0]
    a28 = g28_diss.area / 1e6
    a19 = g19_diss.area / 1e6
    inter = g28_diss.intersection(g19_diss).area / 1e6  # 共通部分
    outside = a28 - inter                                # L28 の L19 はみ出し部分
    pct_in = inter / a28 * 100                           # L28 が L19 にどれだけ入るか
    nest_rows.append({
        "city": name, "L28_K1_km2": a28, "L19_K1_km2": a19,
        "L28_in_L19_pct": pct_in,
        "L28_outside_L19_km2": outside,  # 0 ならサブセット
    })
''')

mean_inside_pct = float(nest_df["L28_in_L19_pct"].dropna().mean())
max_outside = float(nest_df["L28_outside_L19_km2"].max())

s6_html = f"""
<p><b>狙い</b>: L28 シリーズの KUIKI=1 (居住誘導) ポリゴンが、
L19 シリーズの KUIKI=1 (居住誘導) ポリゴンに完全に含まれているか
(L28 ⊆ L19) を<b>幾何 intersection で実証</b>する。</p>

<p><b>手法 — geometry.intersection を使った包含検証</b>:</p>
<ul>
  <li><b>直感的説明</b>: 2 つのポリゴンの重なり部分の面積を計算し、
      L28 の面積に対する重なりの比率が 100% に近ければ<b>L28 が L19 に完全に含まれている</b>。
      逆に L28 にあって L19 にない部分 (はみ出し) を計算すれば、
      0 km² 近ければ<b>サブセットの幾何証明</b>になる。</li>
  <li><b>アルゴリズムのステップ</b>:
      <ol>
        <li>各市町の L28 KUIKI=1 ポリゴン群を <code>dissolve()</code> で 1 つの幾何に統合</li>
        <li>同様に L19 KUIKI=1 ポリゴン群を <code>dissolve()</code> で統合</li>
        <li><code>g28.intersection(g19)</code> で共通部分の幾何を取得</li>
        <li>その面積を計算: <code>inter.area / 1e6</code></li>
        <li>包含率 = 共通部分 / L28 全体 ── 100% なら完全包含</li>
        <li>はみ出し = L28 - 共通部分 ── 0 km² ならサブセット</li>
      </ol></li>
  <li><b>入力と出力</b>:
      入力 = 2 つの GeoDataFrame の <code>geometry</code> 列、
      出力 = 包含率 (%) と はみ出し面積 (km²)</li>
  <li><b>前提と限界</b>: ポリゴンが overlapping=true の場合は dissolve で吸収される。
      非常に小さな数値誤差 (浮動小数点) が 0.001 km² オーダーで生じうる。</li>
  <li><b>パラメータ</b>: なし (純粋な幾何演算)</li>
</ul>

<h3>実装コード</h3>
{code_nest}

<h3>なぜこの図か</h3>
<p>(a) L19 全体 → (b) L28 切り出し → (c) 都市機能拠点ズーム の 3 段階で、
「L28 が L19 のサブセット」かつ「L28 KUIKI=3 が L28 KUIKI=1 内に位置する」二重入れ子を
1 図で見せる。<b>廿日市市</b>を例に選んだのは L19/L28 の差分が最も明瞭で、
かつ KUIKI=3 拠点 3 個が独立して見えるため。</p>

{figure("assets/L28_nest_closeup.png",
        "廿日市市 — L28 ⊆ L19 入れ子と KUIKI=3 拠点クローズアップ")}

<p><b>この図から読み取れること</b>:</p>
<ul>
  <li><b>(a)</b> 廿日市市の L19 居住誘導全体は約 18.8 km² の広範な領域。</li>
  <li><b>(b)</b> L28 KUIKI=1 (濃青) は L19 (薄青背景) の<b>ごく一部 (約 7%) のみ</b>に絞られている。
      → L28 シリーズの「KUIKI=1」は L19 の都市機能拠点近傍部分だけを切り出している。</li>
  <li><b>(c)</b> ズームすると 3 個の都市機能誘導区域 (赤) が中心市街地に集中して配置されている
      ことが分かる。これが廿日市市の<b>多核拠点ネットワーク型</b>の幾何実装。</li>
  <li>結果: <b>仮説 H2 を支持</b> (L28 KUIKI=1 ⊆ L19 KUIKI=1 が幾何的に確認)。</li>
</ul>

{figure("assets/L28_nest_choropleth.png",
        "包含率 choropleth (左) + L28 vs L19 KUIKI=1 面積散布 (右)")}

<p><b>この図から読み取れること</b>:</p>
<ul>
  <li><b>左図</b>: 8 市すべての L28 ⊆ L19 包含率が <b>99% 以上</b> (青色域)。
      → 浮動小数点の数値誤差レベルしかはみ出しがない = <b>完全包含</b>。</li>
  <li><b>右図</b>: 全 8 市が 1:1 線 (グレー破線) より下方に位置 = L28 < L19 が常に成り立つ。
      8 市平均は 0.4 線 (赤点線) 付近 ─ <b>L28 ≈ 0.4 × L19</b>。</li>
  <li>結果: H1 (L28 ≈ 30〜50% × L19) と H2 (L28 ⊆ L19) が<b>両方支持</b>される。</li>
</ul>

<h3>表 3: L28 ⊆ L19 包含検証 (要件 G)</h3>
{nest_table}

<p><b>この表から読み取れること</b>:</p>
<ul>
  <li>8 市の L28 ⊆ L19 包含率の平均は <b>{mean_inside_pct:.2f}%</b>、
      最大はみ出し面積は <b>{max_outside:.4f} km²</b> (浮動小数点の数値誤差レベル)。
      → <b>L28 KUIKI=1 ⊆ L19 KUIKI=1</b> が幾何的に証明された。</li>
  <li>呉市・廿日市市・東広島市など、KUIKI=3 拠点が指定されている市では特に
      L28 KUIKI=1 が L19 KUIKI=1 のごく一部 (10-20%) に絞られている。
      → これらの市は<b>都市機能拠点を中心とした近接居住誘導の精選</b>を L28 で提供。</li>
  <li>福山市・広島市は L28 KUIKI=1 / L19 KUIKI=1 = 約 40% で相対的に大きい。
      これらは KUIKI=3 が L28 シリーズ内に直接含まれていないが、
      L28 KUIKI=1 自体が「都市機能拠点近接領域」として広めに切り出されている。</li>
</ul>
"""
sections.append(("6. 分析3: L28 ⊆ L19 入れ子の幾何検証", s6_html))


# === セクション7: 分析4 — 都市機能誘導 (KUIKI=3) 拠点の構造 ===
toshi_table_rows = []
for _, r in toshi_df.iterrows():
    if r["has_kinou"]:
        avg_d = (f"{r['kinou_avg_pairwise_dist_km']:.2f}"
                 if not pd.isna(r['kinou_avg_pairwise_dist_km']) else "—")
        d_l28 = (f"{r['kinou_min_dist_to_jukyo_l28_km']:.2f}"
                 if not pd.isna(r['kinou_min_dist_to_jukyo_l28_km']) else "—")
        d_l19 = (f"{r['kinou_min_dist_to_jukyo_l19_km']:.2f}"
                 if not pd.isna(r['kinou_min_dist_to_jukyo_l19_km']) else "—")
        toshi_table_rows.append(
            f"<tr><td><b>{r['city']}</b></td><td>{r['ctype']}</td>"
            f"<td>{r['admin_area_km2']:.1f}</td>"
            f"<td>{r['pop_k']:.0f}</td>"
            f"<td>{int(r['n_kinou_polys'])}</td>"
            f"<td>{r['kinou_area_km2']:.3f}</td>"
            f"<td>{r['kinou_in_admin_pct']:.4f}</td>"
            f"<td>{avg_d}</td>"
            f"<td>{d_l28}</td>"
            f"<td>{d_l19}</td></tr>"
        )
    else:
        toshi_table_rows.append(
            f"<tr><td><b>{r['city']}</b></td><td>{r['ctype']}</td>"
            f"<td>{r['admin_area_km2']:.1f}</td>"
            f"<td>{r['pop_k']:.0f}</td>"
            f"<td>0 (未指定)</td>"
            f"<td colspan='5'>L28 シリーズに KUIKI=3 ポリゴンなし</td></tr>"
        )
toshi_table = (
    "<table>"
    "<tr><th>市町</th><th>タイプ</th>"
    "<th>行政 km²</th><th>人口千人</th>"
    "<th>拠点数</th><th>拠点総面積 km²</th>"
    "<th>行政比率 %</th>"
    "<th>拠点間平均距離 km</th>"
    "<th>拠点→L28 K1 km</th>"
    "<th>拠点→L19 K1 km</th></tr>"
    + "".join(toshi_table_rows) + "</table>"
)

code_kinou = code('''
toshi_rows = []
for ds28, ds19, name, ctype, _ in LIP_CITY_DEFS:
    sub = l28[(l28["source_city"] == name) & (l28["KUIKI_CD"] == 3)]
    n_polys = len(sub)
    a_sum = float(sub["poly_area_km2"].sum())

    # 拠点間平均距離 (km) ─ 各 KUIKI=3 ポリゴンの representative_point 同士の距離
    pts = [pg.representative_point() for pg in sub.geometry]
    if len(pts) >= 2:
        dists = [pts[i].distance(pts[j]) / 1000
                 for i in range(len(pts)) for j in range(i + 1, len(pts))]
        avg_d = float(np.mean(dists))
    else:
        avg_d = np.nan

    g_k = sub.dissolve().geometry.iloc[0] if n_polys > 0 else None

    # 拠点と L28 KUIKI=1 (近接居住誘導) の距離・包含
    g_j_l28 = l28[(l28["source_city"] == name) &
                   (l28["KUIKI_CD"] == 1)].dissolve().geometry.iloc[0]
    inter_l28 = g_k.intersection(g_j_l28).area / 1e6
    pct_l28 = inter_l28 / (g_k.area / 1e6) * 100 if g_k.area > 0 else np.nan
    d_l28 = min(pg.distance(g_j_l28) / 1000 for pg in sub.geometry)

    # 拠点と L19 KUIKI=1 (居住誘導全体) の距離・包含 ─ 広域での検証
    g_j_l19 = l19[(l19["source_city"] == name) &
                   (l19["KUIKI_CD"] == 1)].dissolve().geometry.iloc[0]
    inter_l19 = g_k.intersection(g_j_l19).area / 1e6
    pct_l19 = inter_l19 / (g_k.area / 1e6) * 100 if g_k.area > 0 else np.nan
    d_l19 = min(pg.distance(g_j_l19) / 1000 for pg in sub.geometry)

    toshi_rows.append({
        "city": name, "ctype": ctype,
        "n_kinou_polys": n_polys, "kinou_area_km2": a_sum,
        "kinou_avg_pairwise_dist_km": avg_d,
        "kinou_min_dist_to_jukyo_l28_km": d_l28,
        "kinou_min_dist_to_jukyo_l19_km": d_l19,
        "kinou_in_jukyo_pct_l28": pct_l28,
        "kinou_in_jukyo_pct_l19": pct_l19,
        "has_kinou": n_polys > 0,
    })
''')

s7_html = f"""
<p><b>狙い</b>: 都市機能誘導区域 (KUIKI=3) の構造 ─
拠点数・拠点総面積・拠点間平均距離・近接居住誘導との重なり ─ を 8 市で計量し、
<b>「拠点ネットワーク型コンパクトシティ」の幾何実装</b>を読み取る。</p>

<p><b>手法 — 拠点 (point) ベース距離指標</b>:</p>
<ul>
  <li><b>直感的説明</b>: 各 KUIKI=3 ポリゴンを 1 つの「拠点」とみなし、
      代表点 (<code>representative_point</code>) の座標を取得。
      拠点間のユークリッド距離 (m) を計算し、km に変換。<br>
      代表点を使う理由: 重心 (<code>centroid</code>) は形状が複雑な場合にポリゴンの外に出ることがあるため、
      確実にポリゴン内の点を返す <code>representative_point</code> が安全。</li>
  <li><b>計算する指標</b>:
      <ul>
        <li>拠点間平均距離: 同一市町内の全ペアの距離の平均 (多核ネットワークの広がり)</li>
        <li>拠点 → 居住誘導 最短距離: 拠点と L28 KUIKI=1 dissolve の最短距離 (拠点と居住の近さ)</li>
        <li>拠点 ⊆ 近接居住誘導 包含率: 拠点 ∩ 居住 / 拠点 (拠点が居住域内にあるか)</li>
      </ul></li>
  <li><b>結果の読み方</b>:
      包含率 100% = 拠点が完全に居住誘導域内に入る ─ 「拠点ネットワーク型」の理想実装。
      最短距離 0 = 拠点が居住誘導と接している。距離 > 1 km は分離型。</li>
</ul>

<h3>実装コード</h3>
{code_kinou}

<h3>なぜこの図か</h3>
<p>都市機能誘導の<b>規模ランキング</b>と<b>人口対比の拠点数</b>の 2 視点で、
「人口の多い市ほど多核」の仮説を一目で確認できる図にした。
バブル散布の点サイズが拠点総面積を示すので、<b>市町タイプ・人口・拠点数・規模</b>の
4 次元情報を 1 図で表現できる。</p>

{figure("assets/L28_kinou_ranking.png",
        "都市機能誘導 規模ランキング + 人口 vs 拠点数 散布")}

<p><b>この図から読み取れること</b>:</p>
<ul>
  <li><b>左図</b>: 都市機能誘導の総面積は<b>呉市 0.99 km² (4 拠点) が最大</b>、
      次いで竹原市 1.80 km² (4 拠点)、東広島市 0.31 km² (2 拠点)。
      8 市中 5 市のみ拠点指定。広島市・福山市・府中市は<b>L28 シリーズに KUIKI=3 を含まない</b>
      (= 都市機能誘導未指定 or 別シリーズで管理)。</li>
  <li><b>右図</b>: 人口と拠点数に明瞭な相関は<b>見えない</b>。
      竹原市 (24千人) で 4 拠点、呉市 (210千人) で 4 拠点と<b>同数</b>。
      → 「<b>人口に比例した拠点数</b>」の仮定は弱く、
        各市の<b>地理特性 (地形・港湾・歴史中心)</b> が拠点配置を決めている可能性が高い。</li>
</ul>

{figure("assets/L28_distance_analysis.png",
        "拠点間距離 (左) と 拠点 → 居住誘導 最短距離 (右)")}

<p><b>この図から読み取れること</b> ─ 本記事の最大発見:</p>
<ul>
  <li><b>左図</b>: 行政面積が広い市ほど拠点間距離も大きい (期待通り)。
      呉市 (港湾型 4 拠点) は拠点間平均距離 7-9 km、竹原市・三原市は中規模 (1-3 km)。</li>
  <li><b>右図 (核心発見)</b>: 5 市すべてで <b>拠点 → 居住誘導の最短距離が km 単位</b>。
      具体例: 呉市 3.06 km、三原市 1.83 km (L19 居住誘導全体相手)、
      東広島市 5.40 km、廿日市市 2.14 km。
      包含率 (拠点 ⊆ 居住誘導) はほぼすべて <b>0%</b>。
      つまり<b>都市機能拠点は居住誘導域内ではなく、km 単位で離れた場所に配置されている</b>。</li>
  <li>これは立地適正化計画の理論モデル
      「居住誘導 ⊃ 都市機能誘導」を<b>反証</b>する重要発見。
      広島県 5 LIP 市の実装は「拠点ネットワーク型」ではなく
      <b>「居住域 ⇔ 都市機能拠点」分離型 LIP</b>。</li>
  <li>解釈仮説: 都市機能誘導は<b>商業・医療・行政の拠点 (駅前・港湾近傍)</b> として
      指定され、住居系の用途地域 (= L19 居住誘導) とは<b>別の地理クラスタ</b>を形成。
      公共交通でつなぐことで住み分けを実装している ─ コンパクトシティ理論の
      日本ローカル変形版。</li>
  <li>L28 KUIKI=1 (近接居住誘導) と L19 KUIKI=1 (居住誘導全体) どちらの相手でも
      最短距離は同じ km 単位 ─ <b>分離は L19 vs L28 の差ではなく実空間の幾何配置差</b>であることが確認できる。</li>
</ul>

<h3>表 4: 都市機能誘導 8 市サマリ (要件 G)</h3>
{toshi_table}

<p><b>この表から読み取れること</b>:</p>
<ul>
  <li>都市機能誘導の<b>行政比率は 0.0001〜0.5%</b> ─ 極めて小さく、「精選された拠点」として配置されている。</li>
  <li>呉市・廿日市市・東広島市の 3 市が<b>多核 (≥ 2 拠点)</b>、
      竹原市・三原市は<b>少核 (≤ 4 拠点)</b>。
      コンパクトシティの実装スタイルが明らかに分かれている。</li>
  <li>5 市すべてで拠点 → L19 居住誘導 最短距離が <b>1.8〜5.4 km</b> ─
      拠点は居住誘導域<b>外</b>に配置されている。
      <b>「拠点ネットワーク型」の理論モデルは実装されておらず、
      「分離型 LIP」</b> が広島県の運用スタイル。</li>
</ul>
"""
sections.append(("7. 分析4: 都市機能誘導 (KUIKI=3) 拠点の構造", s7_html))


# === セクション8: 分析5 — 8 市 主題図と small multiples ===
s8_html = f"""
<p><b>狙い</b>: 8 LIP 市の都市機能誘導と近接居住誘導の<b>地理配置</b>を主題図で俯瞰し、
さらに small multiples で各市の拠点パターンを比較する。</p>

<h3>図 5: 8 市 主題図 (位置情報の活用 - 要件 T)</h3>
<p>L28 シリーズが県内 8 市にどう分布しているかを、行政背景の上に重ね描き。</p>

{figure("assets/L28_thematic.png",
        "8 LIP 策定市 — L28 都市機能誘導 (赤) + 近接居住誘導 (青)")}

<p><b>この図から読み取れること</b>:</p>
<ul>
  <li>近接居住誘導 (青) は8市にそれぞれ細かく分散しており、行政区域の中心市街地周辺に集中。</li>
  <li>都市機能誘導 (赤) は<b>呉市・竹原市・三原市・東広島市・廿日市市</b>の中心市街地に明瞭に位置。
      広島市・福山市・府中市では赤がほぼ見えない (L28 KUIKI=3 未含)。</li>
  <li>沿岸部の市 (呉市・竹原市・三原市・廿日市市) は港湾近傍に拠点を持つ
      ─ 歴史的中心地 + 交通結節点としての地理的合理性。</li>
</ul>

<h3>図 6: 8 市 small multiples (要件 T)</h3>
<p>各市を独立 panel で比較。同一スケールで配置することで「規模差」と「拠点パターン差」が並列に見える。</p>

{figure("assets/L28_city_small_multiples.png",
        "8 LIP 市 — 各市の都市機能誘導 (赤) + 近接居住誘導 (青)")}

<p><b>この図から読み取れること</b>:</p>
<ul>
  <li><b>政令市・中核市 (広島・福山・呉)</b>: 近接居住誘導の規模が大きい (10-60 km²)。
      しかし広島市・福山市は赤 (KUIKI=3) なし、呉市のみ複数拠点を表示。</li>
  <li><b>施行時特例市 東広島市</b>: 行政面積広大 (635 km²) に対して
      L28 範囲は中心 西条駅周辺に集中 (3.2 km²)、拠点 2 個。</li>
  <li><b>市 (竹原・三原・府中・廿日市)</b>: 規模は小さい (1-5 km²) が
      都市機能誘導と近接居住誘導が compact に同じ場所に集中。</li>
  <li>「コンパクトシティ」の名前通り、行政面積に対して<b>居住・都市機能の指定範囲が極めて小さい</b>
      ─ 平均で行政面積の 0.5-3% 程度。「街全域に住む」ではなく
      「<b>厳選された範囲に都市機能を集中</b>」の制度設計が幾何的に確認できる。</li>
</ul>

<h3>図 7: 8 市 × KUIKI 構成 (要件 G との対応)</h3>

{figure("assets/L28_kuiki_composition.png",
        "L28: 8 市 × KUIKI_CD 件数 stack (左) / 面積 stack (右)")}

<p><b>この図から読み取れること</b>:</p>
<ul>
  <li><b>左 (件数 log)</b>: 広島市 678 polys が突出 (大半が KUIKI=2 境界線細片)。
      福山市 136 polys、東広島市 219 polys と続く。
      → 件数で見ると「中身は線データの細片」が多いが、これは面積では無視できる量。</li>
  <li><b>右 (面積)</b>: KUIKI=1 (青) が圧倒的に大きく、KUIKI=3 (赤) はその上に細く乗る。
      → L28 シリーズの中で<b>KUIKI=1 (近接居住誘導) が容積で 95% 以上</b>を占める。
        「都市機能誘導区域」シリーズという名前と裏腹に、中身は L19 と同じく居住誘導が主体。</li>
</ul>

<h3>図 8: シリーズ実態の概念図 (要件 H)</h3>

{figure("assets/L28_series_concept.png",
        "L28 と L19 の関係 ─ 呉市実例 (左) + 概念図 (右)")}

<p><b>この図から読み取れること</b>:</p>
<ul>
  <li><b>左 (呉市実例)</b>: L19 KUIKI=1 (薄青) ⊃ L28 KUIKI=1 (濃青) ⊃ L28 KUIKI=3 (赤) の<b>三重入れ子</b>が
      実データで明瞭に確認できる。
      L19 (14.83 km²) → L28 K1 (8.38 km²) → L28 K3 (0.99 km²) と
      <b>規模が桁オーダーで縮小</b>していく階層構造。</li>
  <li><b>右 (概念図)</b>: 「L28 シリーズ ≈ L19 ∩ buffer(KUIKI=3, R km)」
      = 都市機能拠点を中心とした近接居住誘導のみを切り出したもの、と概念的に表現。
      <b>都市機能拠点は複数 (多核)</b> で、それぞれの周辺の居住誘導を 1 つにまとめて L28 KUIKI=1 とする
      DoBoX のシリーズ設計が読み取れる。</li>
</ul>
"""
sections.append(("8. 分析5: 主題図と small multiples", s8_html))


# === セクション9: 分析6 — 整合性検証 ===
recon_table = """
<table>
<tr><th>KUIKI_CD</th><th>類型</th>
<th>8 市別 polys</th><th>8 市別 km²</th>
<th>県全域 ds=934 polys</th><th>県全域 km²</th>
<th>差 (面積 %)</th></tr>
"""
for k in sorted(set(sum_l28_city) | set(sum_l28_ken)):
    a_city = sum_l28_city.get(k, 0)
    a_ken = sum_l28_ken.get(k, 0)
    n_city = n_l28_city.get(k, 0)
    n_ken = n_l28_ken.get(k, 0)
    d_a = (a_city - a_ken) / a_ken * 100 if a_ken > 0 else None
    recon_table += (
        f"<tr><td><b>{k}</b></td><td>{KUIKI_INFO[k][0]}</td>"
        f"<td>{n_city}</td><td>{a_city:.4f}</td>"
        f"<td>{n_ken}</td><td>{a_ken:.4f}</td>"
        f"<td>{d_a:+.4f}%" if d_a is not None else "<td>N/A"
    )
    recon_table += "</td></tr>"
recon_table += "</table>"

s9_html = f"""
<p><b>狙い</b>: 8 市別 8 ファイルの合計と、県全域版 ds=934 (1 ファイル, 1150 ポリゴン) の<b>同一性</b>を確認する。
DoBoX が「重複データ」を別 dataset_id で配信していることが多いため、
両者が一致するかは<b>データの整合性 (consistency)</b>の根本検証となる。</p>

{figure("assets/L28_reconciliation.png",
        "L28 整合性検証 ─ 8 市別合計 vs 県全域 ds=934 (面積 + 件数)")}

<p><b>この図から読み取れること</b>:</p>
<ul>
  <li>KUIKI 別の面積・件数で 8 市別合計と県全域 ds=934 が<b>ほぼ一致</b>。
      面積差は浮動小数点誤差レベル (<b>< 0.01%</b>)。</li>
  <li>件数も完全一致 (KUIKI=1 で {n_l28_city.get(1,0)} vs {n_l28_ken.get(1,0)} 等)。
      → 県全域 ds=934 は<b>8 市別 8 ファイルの集約版</b>であることを実証。</li>
</ul>

<h3>表 5: L28 整合性 (要件 G)</h3>
{recon_table}

<p><b>この表から読み取れること</b>:</p>
<ul>
  <li>全 KUIKI で面積差が < 0.001% (浮動小数点誤差)、件数差 0。
      → 仮説 H5 を支持: 8 市別合計と県全域は<b>完全一致</b>。</li>
  <li>DoBoX が同じ情報を 2 通り (8 市別 + 県全域) で提供している
      ─ 用途に応じて選択できる便利な設計だが、
        分析者は<b>重複ロード</b>に注意する必要がある。</li>
</ul>
"""
sections.append(("9. 分析6: 整合性検証 (8 市別 vs 県全域)", s9_html))


# === セクション10: 仮説検証と考察 ===
verdict_table = "<table><tr><th>仮説</th><th>主張</th><th>主要結果</th><th>判定</th></tr>"
for hk in ["H1", "H2", "H3", "H4", "H5"]:
    h = results[hk]
    main_result = ""
    if hk == "H1":
        main_result = (f"L28 KUIKI=1 / L19 KUIKI=1 平均 "
                       f"<b>{h['8市平均_L28_div_L19_K1_pct']:.1f}%</b> (期待 30-50%)")
    elif hk == "H2":
        main_result = (f"L28 ⊆ L19 包含率平均 "
                       f"<b>{h['8市_包含率_平均_pct']:.2f}%</b> (期待 ≥ 95%)、"
                       f"はみ出し最大 {h['8市_L28_outside_L19_最大_km2']:.4f} km²")
    elif hk == "H3":
        main_result = (
            f"拠点 → L19居住誘導 最短距離平均 "
            f"<b>{h['5市_拠点 → L19K1_最短距離平均_km']:.2f} km</b> "
            f"(期待 < 1)、包含率 "
            f"<b>{h['5市_拠点 ⊆ L19K1_包含率平均_pct']:.1f}%</b>"
        )
    elif hk == "H4":
        main_result = (f"拠点指定 {h['拠点指定_市数']} 市 / 8 市、"
                       f"平均拠点総面積 {h['5市_拠点総面積_平均_km2']:.2f} km², "
                       f"平均拠点数 {h['5市_拠点数_平均']:.1f}")
    elif hk == "H5":
        main_result = ("8 市別合計 vs 県全域: 面積差 < 0.001%、件数完全一致")
    verdict_table += (
        f"<tr><td><b>{hk}</b></td>"
        f"<td>{h['claim']}</td>"
        f"<td>{main_result}</td>"
        f"<td><b>{h['verdict']}</b></td></tr>"
    )
verdict_table += "</table>"

s10_html = f"""
<h3>仮説検証 結果一覧 (表 6 - 要件 G)</h3>
{verdict_table}

<p><b>この表から読み取れること</b>:</p>
<ul>
  <li><b>H1 (シリーズ実態) 支持</b>: L28 KUIKI=1 は L19 KUIKI=1 の約 40%。
      L28 シリーズは<b>L19 居住誘導の都市機能拠点近接サブセット</b>を提供している。</li>
  <li><b>H2 (包含関係) 支持</b>: L28 ⊆ L19 が幾何的に証明された (包含率 99.99%)。
      これは DoBoX のシリーズ設計の解読として<b>新規発見</b>。</li>
  <li><b>H3 (拠点ネットワーク型) 反証 = 重要発見「分離型 LIP」</b>:
      拠点 → L19 居住誘導 最短距離が <b>1.8〜5.4 km</b>、包含率 0%。
      広島県 5 LIP 市は「拠点ネットワーク型」を実装しておらず、
      <b>居住誘導域 ⇔ 都市機能拠点 が分離した「分離型 LIP」</b>として運用している。
      これは立地適正化計画の理論モデルへの反証として重要。</li>
  <li><b>H4 (拠点規模) 部分支持</b>: 8 市中 5 市のみ KUIKI=3 ポリゴンを保有、
      残り 3 市 (広島市・福山市・府中市) は L28 シリーズに含めず別管理。
      多核ネットワーク (3-4 拠点) を持つのは呉市・廿日市市の 2 市のみ。</li>
  <li><b>H5 (整合性) 支持</b>: 8 市別合計 vs 県全域は完全一致。</li>
</ul>

<h3>考察 — 4 つの主要発見</h3>
<ol>
  <li><b>「シリーズ重複」の DoBoX 設計の解読</b>:
      L19 と L28 は<b>同じ列構造を持つ別 dataset_id</b>で、
      実質的に<b>L19 の上位集合と L28 の下位集合</b>の関係。
      L28 は「都市機能誘導区域とその近接居住誘導」だけを切り出した狭い領域、
      L19 は「居住誘導全体 (都市機能誘導も含む)」の広い領域。
      → DoBoX は同じ立地適正化計画情報を<b>2 通りの粒度</b>で提供している。
      利用者は<b>用途に応じて選択</b>:
      「コンパクトな拠点周辺だけ見たい」→ L28、
      「居住誘導全体を見たい」→ L19。</li>
  <li><b>「分離型 LIP」が広島県 5 市の実装スタイル (本記事の最大発見)</b>:
      呉市 (4 拠点)・竹原市 (4 拠点)・三原市 (1 拠点)・東広島市 (2 拠点)・廿日市市 (3 拠点) で、
      都市機能誘導と居住誘導が<b>km 単位で空間的に分離</b>している (距離 1.8〜5.4 km)。
      これは立地適正化計画の理論モデル「拠点ネットワーク型 = 居住誘導 ⊃ 都市機能誘導」を
      <b>反証</b>する重要発見。広島県の LIP は<b>住居系と商業系を地理的に住み分け</b>、
      公共交通でつなぐ日本ローカルの運用スタイルを採用している。</li>
  <li><b>「拠点未指定」3 市の解釈</b>:
      広島市・福山市・府中市は L28 シリーズに KUIKI=3 ポリゴンを持たない。
      L19 でも KUIKI=3 が極めて少ない (広島市 0, 福山市 0, 府中市 0)。
      これは <b>(a) 都市機能誘導を未指定</b>か、
      <b>(b) 別のデータソース (DoBoX 外の別の自治体システム) で管理</b>している可能性。
      <b>2025 年現在の DoBoX 公開データ範囲ではこれらの市の都市機能誘導は把握できない</b>。</li>
  <li><b>都市機能誘導の規模感の発見</b>:
      8 市計の都市機能誘導 (KUIKI=3) は<b>3.2 km²</b> ─ 県内 8 LIP 市の合計でも
      行政面積の<b>0.001%</b>未満。極めて小さな「点」として
      コンパクトシティの核を実装している。
      竹原市の 1.80 km² が最大規模 (港湾中心の歴史的市街地) ─
      <b>大都市ほど拠点が大きい</b>わけではない反証ケース。</li>
</ol>
"""
sections.append(("10. 仮説検証と考察", s10_html))


# === セクション11: 発展課題 ===
s11_html = """
<h3>結果X → 新仮説Y → 課題Z (要件 E)</h3>

<h4>発展課題 1 — DoBoX のシリーズ設計の汎化</h4>
<ul>
  <li><b>結果X</b>: L19 と L28 は同じ列構造で別の切り出しを与える「シリーズ二重提供」と判明。</li>
  <li><b>新仮説Y</b>: L26 (区域外)・L27 (非線引き)・L18 (線引き) も
      <b>同じく「シリーズ間でテーブル構造を共有しつつ別切り出し」</b>の設計を持つ。
      DoBoX 全体で「シリーズは KUIKI_CD/類似フラグで識別する論理的区分」と仮定できる。</li>
  <li><b>課題Z</b>: L17/L18/L19/L26/L27/L28 の 6 シリーズを<b>同一 GeoDataFrame に縦結合</b>し、
      KUIKI_CD と src_kind (シリーズ名) の交差表で「同 KUIKI でも別シリーズなら別ポリゴン集合」
      を県全域で定量化する。<b>L29 統合俯瞰記事</b>として執筆候補。</li>
</ul>

<h4>発展課題 2 — 「分離型 LIP」の合理性検証</h4>
<ul>
  <li><b>結果X</b>: 5 市で都市機能拠点と居住誘導が km 単位で分離 (反証 = 重要発見)。</li>
  <li><b>新仮説Y</b>: 「拠点と最寄り公共交通駅の距離」「拠点内の公共施設密度」「人口集中地区との重なり」
      を計量すれば、<b>分離型 LIP の合理性</b>を評価できる。
      理論的には<b>(a) 拠点 → 最寄り駅 < 500 m, (b) 拠点内に行政・商業・医療の各 1 件以上, (c) DID と重なる</b>
      の 3 条件を満たす拠点なら「公共交通で住居と接続される機能する拠点」と言える。
      逆に居住誘導から分離した拠点でも、これらの条件を満たすなら<b>分離型でも実効性あり</b>。</li>
  <li><b>課題Z</b>: L23 DID境界 + L02/03 駅・避難所 + L25 行政施策 を本記事の都市機能誘導に重ね、
      <b>5 LIP 市の拠点を「公共交通でつながる分離拠点」「孤立拠点」</b>に二分する評価マップを作成。
      L01-L29 の<b>横断統合 X 系記事</b>として実装可能。</li>
</ul>

<h4>発展課題 3 — 拠点未指定 3 市の解明</h4>
<ul>
  <li><b>結果X</b>: 広島市・福山市・府中市は L28/L19 に KUIKI=3 ポリゴンを持たない。</li>
  <li><b>新仮説Y</b>: これらの市は実際に立地適正化計画で都市機能誘導を未指定、
      または<b>別の自治体専用システム</b>で公開しており DoBoX に含まれない。
      自治体公式サイト (広島市まちづくり情報, 福山市都市計画) では別形式で公開されている可能性がある。</li>
  <li><b>課題Z</b>: 各市の<b>立地適正化計画 (PDF/Web)</b> を検索・確認し、
      DoBoX のデータ提供範囲と<b>実際の計画内容</b>のギャップを文書化する。
      これは「オープンデータの限界」を学ぶ教育素材としても有意義。</li>
</ul>

<h4>発展課題 4 — 全国比較</h4>
<ul>
  <li><b>結果X</b>: 広島県 8 市の都市機能誘導は 1-2 km²/市の小規模。</li>
  <li><b>新仮説Y</b>: <b>全国 600+ LIP 策定市</b>と比較すると、広島の規模感は
      平均的か、コンパクト寄りか、それとも分散傾向か?
      仙台市・川越市・松山市など他県のオープンデータと比較するとよい。</li>
  <li><b>課題Z</b>: 国土交通省「立地適正化計画作成の手引き」の事例集と、
      全国 LIP 一覧 (国交省サイト) から比較対象を選定し、
      「コンパクトシティ計画の地理的サイズの全国分布」を作成する研究記事に発展。</li>
</ul>
"""
sections.append(("11. 発展課題", s11_html))


# === Render ===
html = render_lesson(
    num=28,
    title="都市機能誘導区域 9件 統合分析 — L19との入れ子と拠点ネットワーク型コンパクトシティの幾何実証",
    tags=["L28", "都市機能誘導", "立地適正化計画", "拠点ネットワーク型",
          "コンパクトシティ", "L19比較", "GeoJSON", "geopandas"],
    time="20 分",
    level="中級〜上級 (geopandas + 幾何 intersection)",
    data_label="9 GeoJSON (L28: 8市別 + 県全域 1) + L19 比較 + L15 行政背景",
    sections=sections,
    script_filename="L28_urban_function_induction.py",
)

OUT = LESSONS / "L28_urban_function_induction.html"
OUT.write_text(html, encoding="utf-8")
print(f"=== 出力: {OUT}", flush=True)
print(f"=== 総時間: {time.time()-t0:.1f}s ===", flush=True)