# -*- coding: utf-8 -*-
"""L34 臨港交通施設 2 件統合分析
       — 広島県 24 港湾 + 7 漁港の 296 道路・橋・駐車場が描く港湾アクセス第 3 層

カバー宣言:
  本記事は DoBoX のシリーズ <b>「臨港交通施設_港湾/漁港」 2 件</b>
  (dataset_id = 1252, 1256) を統合し、広島県内における
  <b>港湾アクセス交通施設</b> (港湾道路・港湾橋りょう・港湾駐車場) の
  <b>陸海接続構造</b>を分析する研究記事である。

  L32 (外郭施設、dsid 1250+1254) と L33 (係留施設、dsid 1251+1255) の続編で、
  港湾施設の<b>3 層階層 (外郭→係留→臨港交通)</b> の最後のレイヤを扱う。
  すなわち、外郭で守られ、岸壁・物揚場で船を泊めた後、その<b>陸側で貨物・人を
  内陸へ運ぶ「最終 1 km インフラ」</b>。

研究の問い (RQ):
  広島県の臨港交通施設 (合計 296 件) は、L32/L33 の 41 港集合のうち
  <b>どれだけの港にしか整備されていない</b>のか? 整備された港では道路・橋・駐車場が
  <b>どんな構成</b>で陸海接続を支え、<b>外郭・係留と空間的にどう重なる</b>のか?
    (1) L32/L33 で同定された 27 港湾 + 14 漁港のうち、臨港交通が整備された港は何港か?
        <b>整備格差</b> (=港湾交通の有無による港の階層化) はどれくらい大きいか?
    (2) 施設種類 (道路-車道 / 駐車場 / 橋りょう) はカテゴリ別にどう分布しているか?
        港湾は道路+駐車場+橋の3形式、漁港は道路+橋のみ?
    (3) 大規模商業港 (広島港・福山港・尾道糸崎港) は<b>道路面積</b>でも他を引き離すか?
    (4) GIS 情報の<b>欠損パターン</b>: なぜ橋りょう (漁港 6 件) は全て GIS なしか?
        これは「データ品質欠落」か「未測量」か「占用許可ベースで地番管理」か?
    (5) 臨港道路 (POLYGON 表現) は L32 外郭・L33 係留と空間的にどう関係するか?
        係留からの 100 m 圏内に道路がある = 「岸壁直結アクセス」が読めるか?
    (6) 駐車場は港湾のみに整備 = フェリー・客船利用者の自家用車対応 = 旅客港の指標か?

仮説 H1〜H7:
  H1 (整備格差・極端): 臨港交通整備済み港は 24 港湾 + 7 漁港 = 31 港。
       L32/L33 の 41 港のうち <b>10 港 (≈24%) は臨港交通施設なし</b>。
       特に漁港は <b>14 中 7 のみ (50%)</b>整備で、半分は車道・橋すらデータ上ない。
  H2 (件数規模差・最大): 港湾 278 (94%) vs 漁港 18 (6%)、比 <b>15.4:1</b>。
       L32 外郭 (1.33:1)、L33 係留 (1.9:1) と比べ、<b>圧倒的な港湾偏重</b>。
       理由: 漁港は伝統的に既存集落道路を併用するためデータ上の独立交通施設が極少。
  H3 (施設種類の分化): 港湾は 3 種 (道路 172 / 駐車場 93 / 橋 13)、
       漁港は 2 種 (道路 12 / 橋 6) で<b>駐車場が漁港にゼロ</b>。
       これは「漁港利用者は徒歩・近隣駐車・漁協車両のみ = 不特定多数の旅客なし」
       という機能的差異を反映。駐車場 = 旅客港マーカー。
  H4 (大規模港集中・最強): 広島港 79 件 (28.4%) + 尾道糸崎 39 + 福山 21 で
       港湾上位 3 港で <b>50%</b>。L32 (44%)・L33 (46%) よりさらに集中。
       道路網は商業ハブ港でしか面的に整備されない。
  H5 (GIS 欠損の構造的偏り): GIS 情報があるのは港湾 62/278 (22%)、漁港 3/18 (17%)。
       <b>「道路-車道」と「橋りょう」で欠損率が異なる</b>。漁港の橋 6 件は全て GIS なし。
       これは<b>データ整備の優先度</b>の現れ (主要港の道路から測量、橋は後回し)。
  H6 (3 層整合): 臨港交通の道路 POLYGON は、L33 係留施設から <b>100 m 以内</b>に
       <b>80% 以上</b>位置する。すなわち岸壁の真裏に道路が走る「岸壁直結」が原則。
       残りはバックヤード駐車場・連絡路で、外郭からも内側 (=港内陸地) に位置する。
  H7 (3 層比例関係): 1 港あたりの 3 層件数比は概ね
       <b>外郭 : 係留 : 交通 ≈ 5 : 7 : 2</b> で、交通が最も少ない。
       これは「交通施設は港の規模ではなく、貨物量の関数」であることを示唆 (発展)。

要件 S 準拠: 1 分以内完走。GIS 欠損が多いので空間処理は軽量 (POLYGON の centroid と sjoin)。
要件 T 準拠: 県全域マップ + 上位港ズーム + 3 層重ね合わせマップ + 整備格差ヒートマップ。
要件 Q 準拠: 図 9 種以上、表 9 種以上 (296 施設 × 3 種 × 31 港 × 3 層)。

データ仕様 (2 件は同一スキーマ):
  A. 港湾臨港交通 1252: CSV, 278 行 × 14 列, GIS情報 = POLYGON (面ジオメトリ),
     24 港湾, 3 施設種類 (道路-車道 / 駐車場 / 橋りょう)
  B. 漁港臨港交通 1256: CSV, 18 行 × 14 列, 同一スキーマ, 7 漁港,
     2 施設種類 (道路-車道 / 橋りょう)
  共通列: 事業, 所管, 施設分類, 施設種類, 港湾名, 事務所, 市区町村１, 市区町村２,
          施設番号, 施設名称, 管理者名等, GIS情報, 開始位置緯度, 開始位置経度

GIS タイプの違い (重要):
  - L32 外郭 / L33 係留: LINESTRING (線) → 延長 m を計算
  - L34 臨港交通: POLYGON (面) → 面積 m² を計算 (道路や駐車場は「面」で記録)

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

sys.path.insert(0, str(Path(__file__).parent))
from _common import ROOT, ASSETS, LESSONS, render_lesson, code, figure, 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
from matplotlib.lines import Line2D
import geopandas as gpd
from shapely.geometry import Point, Polygon, MultiPolygon
from shapely.wkt import loads as wkt_loads

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

t0 = time.time()
print("=== L34 臨港交通施設 2 件統合分析 ===", flush=True)

# =============================================================================
# 0. 定数: 2 dataset_id と 3 施設種類の凡例
# =============================================================================
DATA_DIR = ROOT / "data" / "extras" / "L34_port_traffic"
DATA_DIR.mkdir(parents=True, exist_ok=True)
ADMIN_DIR = ROOT / "data" / "extras" / "L15_admin_zones"
L32_CSV = LESSONS / "assets" / "L32_all_facilities.csv"
L33_CSV = LESSONS / "assets" / "L33_all_facilities.csv"
TARGET_CRS = "EPSG:6671"  # JGD2011 平面直角 III, m 単位

SERIES = [
    (1252, "臨港交通施設（港湾）", "港湾", "harbor_traffic_facility.csv",       32496),
    (1256, "臨港交通施設（漁港）", "漁港", "fishing_port_traffic_facility.csv", 32500),
]

# 施設種類 (港湾3+漁港2 の合計 union = 3 種)
KIND_ORDER = ["道路-車道", "駐車場", "橋りょう"]
# 物理機能でグルーピング
KIND_FUNCTION = {
    "道路-車道":  "港内車道 (車両 + 貨物移動)",
    "駐車場":    "旅客・貨物の一時駐車",
    "橋りょう":  "水路・運河を跨ぐ短橋",
}
# 施設種類別の色
KIND_COLOR = {
    "道路-車道":  "#1f77b4",  # 青 (主役 = 道路網)
    "駐車場":    "#e7ba52",  # 砂色 (旅客対応)
    "橋りょう":  "#cf222e",  # 赤 (構造物)
}
# カテゴリ色 (港湾 / 漁港) — L32/L33 と統一
CAT_COLOR = {"港湾": "#0969da", "漁港": "#1a7f37"}


# =============================================================================
# 1. 2 dataset の存在確認とロード
# =============================================================================
print("\n[1] 2 dataset の存在確認とロード", flush=True)
t1 = time.time()

dfs = []
for dsid, label, cat, fname, rid in SERIES:
    p = DATA_DIR / fname
    try:
        ensure_dataset(p, dataset_id=dsid, resource_id=rid, label=f"L34/{dsid} {label}")
    except Exception as e:
        print(f"  WARN dsid={dsid}: {e}", flush=True)
    df = pd.read_csv(p, encoding="utf-8-sig")
    df["port_category"] = cat
    df["dsid"] = dsid
    dfs.append(df)
    print(f"  {label} (dsid={dsid}): {df.shape}", flush=True)

# 連結
ALL = pd.concat(dfs, ignore_index=True)
print(f"  合計: {len(ALL)} 行 × {ALL.shape[1]} 列", flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 2. WKT パース → GeoDataFrame (POLYGON ベース)
# =============================================================================
print("\n[2] GIS情報 WKT を Shapely geometry にパース", flush=True)
t1 = time.time()

def parse_wkt(s):
    if pd.isna(s):
        return None
    try:
        return wkt_loads(str(s))
    except Exception:
        return None

ALL["geometry"] = ALL["GIS情報"].apply(parse_wkt)
gdf = ALL.dropna(subset=["geometry"]).copy()
gdf = gpd.GeoDataFrame(gdf, geometry="geometry", crs="EPSG:4326").to_crs(TARGET_CRS)
gdf["area_m2"] = gdf.geometry.area
gdf["perimeter_m"] = gdf.geometry.length  # POLYGON の周長
print(f"  geom 有効: {len(gdf)} / {len(ALL)} ({len(gdf)/len(ALL)*100:.1f}%)", flush=True)
print(f"  geom 欠損: 港湾={int((ALL[ALL['port_category']=='港湾']['geometry'].isna()).sum())}, "
      f"漁港={int((ALL[ALL['port_category']=='漁港']['geometry'].isna()).sum())}", flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 3. カテゴリ × 施設種類のクロス集計
# =============================================================================
print("\n[3] カテゴリ × 施設種類 クロス", flush=True)
t1 = time.time()

# 全件 (geom 欠損も含む) の集計
pv_kind = ALL.groupby(["port_category", "施設種類"]).size().unstack(fill_value=0)\
    .reindex(columns=KIND_ORDER, fill_value=0)
pv_kind["合計"] = pv_kind.sum(axis=1)
pv_kind.loc["合計"] = pv_kind.sum(axis=0)
print(pv_kind.to_string())

# 面積 m² ベース (POLYGON のみ) → ha (10000 m²) に変換
pv_area = pd.pivot_table(
    gdf, index="port_category", columns="施設種類",
    values="area_m2", aggfunc="sum", fill_value=0,
).reindex(columns=KIND_ORDER, fill_value=0) / 10000.0
pv_area = pv_area.round(3)
pv_area["合計ha"] = pv_area.sum(axis=1).round(3)
pv_area.loc["合計"] = pv_area.sum(axis=0).round(3)
print()
print("面積 (ha):")
print(pv_area.to_string())
print(f"  ({time.time()-t1:.1f}s)", flush=True)

pv_kind.to_csv(ASSETS / "L34_kind_count_cross.csv", encoding="utf-8-sig")
pv_area.to_csv(ASSETS / "L34_kind_area_cross.csv", encoding="utf-8-sig")


# =============================================================================
# 4. GIS 欠損パターン分析 (本記事固有の重要分析)
# =============================================================================
print("\n[4] GIS 欠損パターン分析", flush=True)
t1 = time.time()

gis_pattern = ALL.groupby(["port_category", "施設種類"]).agg(
    n_total=("GIS情報", "size"),
    n_with_gis=("GIS情報", lambda x: x.notna().sum()),
).reset_index()
gis_pattern["欠損率(%)"] = ((1 - gis_pattern["n_with_gis"] / gis_pattern["n_total"]) * 100).round(1)
gis_pattern.to_csv(ASSETS / "L34_gis_missing_pattern.csv", index=False, encoding="utf-8-sig")
print(gis_pattern.to_string(index=False))
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 5. 港別集計
# =============================================================================
print("\n[5] 港別の施設数・面積集計", flush=True)
t1 = time.time()

port_agg = ALL.groupby(["port_category", "港湾名"]).size().reset_index(name="n_facilities")
# 面積 (geom ありのみ)
area_agg = gdf.groupby(["port_category", "港湾名"]).agg(
    total_area_m2=("area_m2", "sum"),
    n_with_geom=("area_m2", "size"),
).reset_index()
port_agg = port_agg.merge(area_agg, on=["port_category", "港湾名"], how="left")
port_agg["total_area_m2"] = port_agg["total_area_m2"].fillna(0)
port_agg["n_with_geom"] = port_agg["n_with_geom"].fillna(0).astype(int)
port_agg["total_area_ha"] = (port_agg["total_area_m2"] / 10000).round(3)
port_agg = port_agg.sort_values(["port_category", "n_facilities"], ascending=[True, False])
port_agg.to_csv(ASSETS / "L34_port_aggregate.csv", index=False, encoding="utf-8-sig")
print(port_agg.to_string(index=False))

n_ports_h_l34 = port_agg[port_agg["port_category"] == "港湾"]["港湾名"].nunique()
n_ports_f_l34 = port_agg[port_agg["port_category"] == "漁港"]["港湾名"].nunique()
print(f"  L34 整備済み: 港湾 {n_ports_h_l34} 港 / 漁港 {n_ports_f_l34} 港", flush=True)

# 上位
top5_h = port_agg[port_agg["port_category"] == "港湾"].head(5)
top5_f = port_agg[port_agg["port_category"] == "漁港"].head(5)
print()
print("港湾 上位 5:")
print(top5_h.to_string(index=False))
print("漁港 上位 5:")
print(top5_f.to_string(index=False))
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 6. 整備格差: L32/L33 の 41 港集合と比較
# =============================================================================
print("\n[6] 整備格差: L32/L33 41 港集合との比較", flush=True)
t1 = time.time()

# L32/L33 の港集合を取得 (L33_all_facilities.csv は L32 と同集合のはず)
ports_l32_l33_h = set()
ports_l32_l33_f = set()
GAP_OK = False
if L33_CSV.exists():
    l33 = pd.read_csv(L33_CSV, encoding="utf-8-sig", usecols=["port_category", "港湾名"])
    ports_l32_l33_h = set(l33[l33["port_category"] == "港湾"]["港湾名"].unique())
    ports_l32_l33_f = set(l33[l33["port_category"] == "漁港"]["港湾名"].unique())
    GAP_OK = True
elif L32_CSV.exists():
    l32 = pd.read_csv(L32_CSV, encoding="utf-8-sig", usecols=["port_category", "港湾名"])
    ports_l32_l33_h = set(l32[l32["port_category"] == "港湾"]["港湾名"].unique())
    ports_l32_l33_f = set(l32[l32["port_category"] == "漁港"]["港湾名"].unique())
    GAP_OK = True
ports_l34_h = set(ALL[ALL["port_category"] == "港湾"]["港湾名"].unique())
ports_l34_f = set(ALL[ALL["port_category"] == "漁港"]["港湾名"].unique())

if GAP_OK:
    only_l32_l33_h = ports_l32_l33_h - ports_l34_h
    only_l32_l33_f = ports_l32_l33_f - ports_l34_f
    both_h = ports_l32_l33_h & ports_l34_h
    both_f = ports_l32_l33_f & ports_l34_f
    print(f"  港湾 41 港のうち: L34 整備済み {len(both_h)} 港、L32/L33 のみ {len(only_l32_l33_h)} 港", flush=True)
    print(f"    L34 未整備の港湾: {sorted(only_l32_l33_h)}", flush=True)
    print(f"  漁港 14 漁港のうち: L34 整備済み {len(both_f)}、L32/L33 のみ {len(only_l32_l33_f)}", flush=True)
    print(f"    L34 未整備の漁港: {sorted(only_l32_l33_f)}", flush=True)
    n_total_pref = len(ports_l32_l33_h) + len(ports_l32_l33_f)
    n_l34_covered = len(both_h) + len(both_f)
    coverage_rate = n_l34_covered / max(1, n_total_pref) * 100
    print(f"  L34 のカバー率: {n_l34_covered}/{n_total_pref} = {coverage_rate:.1f}%", flush=True)
    # 表として保存
    gap_rows = []
    for p in sorted(ports_l32_l33_h | ports_l34_h):
        in_l32 = p in ports_l32_l33_h
        in_l34 = p in ports_l34_h
        gap_rows.append({
            "port_category": "港湾", "港湾名": p,
            "L32_L33": "○" if in_l32 else "−",
            "L34": "○" if in_l34 else "−",
            "状態": "全層整備" if (in_l32 and in_l34) else ("交通のみ" if in_l34 else "交通なし"),
        })
    for p in sorted(ports_l32_l33_f | ports_l34_f):
        in_l32 = p in ports_l32_l33_f
        in_l34 = p in ports_l34_f
        gap_rows.append({
            "port_category": "漁港", "港湾名": p,
            "L32_L33": "○" if in_l32 else "−",
            "L34": "○" if in_l34 else "−",
            "状態": "全層整備" if (in_l32 and in_l34) else ("交通のみ" if in_l34 else "交通なし"),
        })
    gap_df = pd.DataFrame(gap_rows)
    gap_df.to_csv(ASSETS / "L34_layer_coverage_gap.csv", index=False, encoding="utf-8-sig")
else:
    only_l32_l33_h = set()
    only_l32_l33_f = set()
    both_h = ports_l34_h
    both_f = ports_l34_f
    n_total_pref = len(ports_l34_h) + len(ports_l34_f)
    n_l34_covered = n_total_pref
    coverage_rate = 100.0
    print(f"  WARN: L32/L33 CSV not found, gap analysis skipped", flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 7. 県境 polygon の load
# =============================================================================
print("\n[7] 県境の load", flush=True)
t1 = time.time()

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)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 8. 3 層空間関係: 臨港交通 × L33 係留 × L32 外郭
# =============================================================================
print("\n[8] 3 層 (外郭/係留/交通) の空間関係", flush=True)
t1 = time.time()

# L32 外郭と L33 係留の WKT を再ジオ化
L32_OK = L32_CSV.exists()
L33_OK = L33_CSV.exists()
PROXIMITY_M = 100.0  # 交通 → 係留 の近接距離 (岸壁直結を見るため 100m)

l33_gdf = None
l32_gdf = None
n_near_moor_h = n_near_moor_f = 0
n_total_geo_h = int((gdf["port_category"] == "港湾").sum())
n_total_geo_f = int((gdf["port_category"] == "漁港").sum())

if L33_OK:
    l33 = pd.read_csv(L33_CSV, encoding="utf-8-sig")
    if "GIS情報" in l33.columns:
        l33["geometry"] = l33["GIS情報"].apply(parse_wkt)
        l33_gdf = gpd.GeoDataFrame(
            l33.dropna(subset=["geometry"]),
            geometry="geometry", crs="EPSG:4326",
        ).to_crs(TARGET_CRS)
        # 係留全体の MultiLineString
        moor_diss = l33_gdf.geometry.union_all()
        # 各臨港交通施設 (POLYGON) の境界から係留への最近距離
        gdf["dist_to_moor_m"] = gdf.geometry.distance(moor_diss)
        gdf["near_moor"] = gdf["dist_to_moor_m"] <= PROXIMITY_M
        n_near_moor_h = int(((gdf["port_category"]=="港湾") & gdf["near_moor"]).sum())
        n_near_moor_f = int(((gdf["port_category"]=="漁港") & gdf["near_moor"]).sum())
        print(f"  港湾 交通施設が係留から {PROXIMITY_M:.0f}m 以内: "
              f"{n_near_moor_h}/{n_total_geo_h} ({n_near_moor_h/max(1,n_total_geo_h)*100:.1f}%)", flush=True)
        print(f"  漁港 交通施設が係留から {PROXIMITY_M:.0f}m 以内: "
              f"{n_near_moor_f}/{n_total_geo_f} ({n_near_moor_f/max(1,n_total_geo_f)*100:.1f}%)", flush=True)
        print(f"  交通 全体の係留距離 中央値: {gdf['dist_to_moor_m'].median():.1f} m", flush=True)
    else:
        gdf["dist_to_moor_m"] = np.nan
        gdf["near_moor"] = False
        L33_OK = False
else:
    gdf["dist_to_moor_m"] = np.nan
    gdf["near_moor"] = False

# L32 外郭との距離も
n_near_outer_h = n_near_outer_f = 0
PROXIMITY_OUTER_M = 300.0
if L32_OK:
    l32 = pd.read_csv(L32_CSV, encoding="utf-8-sig")
    if "GIS情報" in l32.columns:
        l32["geometry"] = l32["GIS情報"].apply(parse_wkt)
        l32_gdf = gpd.GeoDataFrame(
            l32.dropna(subset=["geometry"]),
            geometry="geometry", crs="EPSG:4326",
        ).to_crs(TARGET_CRS)
        outer_diss = l32_gdf.geometry.union_all()
        gdf["dist_to_outer_m"] = gdf.geometry.distance(outer_diss)
        gdf["near_outer"] = gdf["dist_to_outer_m"] <= PROXIMITY_OUTER_M
        n_near_outer_h = int(((gdf["port_category"]=="港湾") & gdf["near_outer"]).sum())
        n_near_outer_f = int(((gdf["port_category"]=="漁港") & gdf["near_outer"]).sum())
        print(f"  港湾 交通施設が外郭から {PROXIMITY_OUTER_M:.0f}m 以内: "
              f"{n_near_outer_h}/{n_total_geo_h} ({n_near_outer_h/max(1,n_total_geo_h)*100:.1f}%)", flush=True)
    else:
        gdf["dist_to_outer_m"] = np.nan
        gdf["near_outer"] = False
        L32_OK = False
else:
    gdf["dist_to_outer_m"] = np.nan
    gdf["near_outer"] = False
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 9. 3 層件数比 (港別 x 層) - 整備済み港のみ
# =============================================================================
print("\n[9] 3 層件数比", flush=True)
t1 = time.time()

three_layer = []
if GAP_OK:
    if L32_OK and l32_gdf is not None:
        # l32 ALL (geom 含めず) も読む
        l32_all = pd.read_csv(L32_CSV, encoding="utf-8-sig", usecols=["port_category", "港湾名"])
        l33_all = pd.read_csv(L33_CSV, encoding="utf-8-sig", usecols=["port_category", "港湾名"])
        for cat, ports_set in [("港湾", ports_l32_l33_h), ("漁港", ports_l32_l33_f)]:
            for p in sorted(ports_set):
                n_outer = int(((l32_all["port_category"]==cat) & (l32_all["港湾名"]==p)).sum())
                n_moor = int(((l33_all["port_category"]==cat) & (l33_all["港湾名"]==p)).sum())
                n_traffic = int(((ALL["port_category"]==cat) & (ALL["港湾名"]==p)).sum())
                three_layer.append({
                    "port_category": cat, "港湾名": p,
                    "外郭(L32)": n_outer, "係留(L33)": n_moor, "交通(L34)": n_traffic,
                    "合計": n_outer + n_moor + n_traffic,
                })
three_layer_df = pd.DataFrame(three_layer)
if len(three_layer_df) > 0:
    three_layer_df = three_layer_df.sort_values(["port_category", "合計"], ascending=[True, False])
    three_layer_df.to_csv(ASSETS / "L34_three_layer_by_port.csv", index=False, encoding="utf-8-sig")
    print(three_layer_df.head(15).to_string(index=False))
    # 全体合計比
    total_outer = three_layer_df["外郭(L32)"].sum()
    total_moor = three_layer_df["係留(L33)"].sum()
    total_traffic = three_layer_df["交通(L34)"].sum()
    print(f"\n  総計 比: 外郭 {total_outer} : 係留 {total_moor} : 交通 {total_traffic}", flush=True)
    if total_traffic > 0:
        print(f"  正規化 (交通=1): "
              f"{total_outer/total_traffic:.1f} : {total_moor/total_traffic:.1f} : 1", flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 10. 中間 CSV: 全件統合
# =============================================================================
print("\n[10] 中間 CSV 出力", flush=True)
t1 = time.time()

out = gdf.drop(columns=["geometry"]).copy()
out.to_csv(ASSETS / "L34_all_facilities.csv", index=False, encoding="utf-8-sig")

series_df = pd.DataFrame([
    {"dsid": d, "label": lab, "category": cat, "file": f, "n_rows": int((ALL["dsid"]==d).sum())}
    for (d, lab, cat, f, _) in SERIES
])
series_df.to_csv(ASSETS / "L34_series.csv", index=False, encoding="utf-8-sig")
print(f"  saved csv ({time.time()-t1:.1f}s)", flush=True)


# 値の事前整理 (HTML から参照)
n_total_all = len(ALL)
n_total_geo = len(gdf)
n_h_total = int((ALL["port_category"] == "港湾").sum())
n_f_total = int((ALL["port_category"] == "漁港").sum())
total_ha_h = pv_area.loc["港湾", "合計ha"] if "港湾" in pv_area.index else 0
total_ha_f = pv_area.loc["漁港", "合計ha"] if "漁港" in pv_area.index else 0
fish_pct = n_f_total / max(1, n_h_total + n_f_total) * 100
count_ratio = n_h_total / max(1, n_f_total)
total_road = int(pv_kind.loc["合計", "道路-車道"]) if "道路-車道" in pv_kind.columns else 0
total_park = int(pv_kind.loc["合計", "駐車場"]) if "駐車場" in pv_kind.columns else 0
total_bridge = int(pv_kind.loc["合計", "橋りょう"]) if "橋りょう" in pv_kind.columns else 0


# =============================================================================
# 11. 図 1: 2 dataset 構造マトリクス + カード
# =============================================================================
print("\n[11] 図 1: dataset 構造", flush=True)
t1 = time.time()

n_h_kinds_str = " / ".join(
    f"{k} {int(pv_kind.loc['港湾',k])}" for k in KIND_ORDER if pv_kind.loc['港湾',k] > 0
)
n_f_kinds_str = " / ".join(
    f"{k} {int(pv_kind.loc['漁港',k])}" for k in KIND_ORDER if pv_kind.loc['漁港',k] > 0
)

fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))

# 左: dataset カード
ax = axes[0]
ax.set_xlim(0, 10); ax.set_ylim(0, 6)
ax.set_axis_off()
ax.text(5, 5.5, "DoBoX 臨港交通施設シリーズ 2 件", ha="center",
        fontsize=15, fontweight="bold")
# 港湾カード
ax.add_patch(plt.Rectangle((0.3, 0.7), 4.4, 4.0, fill=True,
            facecolor=CAT_COLOR["港湾"], alpha=0.13, edgecolor=CAT_COLOR["港湾"], linewidth=2.5))
ax.text(2.5, 4.2, "dsid 1252\n臨港交通施設（港湾）", ha="center", fontsize=12,
        fontweight="bold", color=CAT_COLOR["港湾"])
ax.text(2.5, 3.05, f"{n_h_total} 件 / {n_ports_h_l34} 港湾", ha="center", fontsize=14)
ax.text(2.5, 2.05, n_h_kinds_str, ha="center", fontsize=9.5)
ax.text(2.5, 1.0, "L33 27 港湾のうち 24 港のみ整備\n(3 港 = 三高港・横田港等は交通なし)",
        ha="center", fontsize=8.5, color="#666")

# 漁港カード
ax.add_patch(plt.Rectangle((5.3, 0.7), 4.4, 4.0, fill=True,
            facecolor=CAT_COLOR["漁港"], alpha=0.13, edgecolor=CAT_COLOR["漁港"], linewidth=2.5))
ax.text(7.5, 4.2, "dsid 1256\n臨港交通施設（漁港）", ha="center", fontsize=12,
        fontweight="bold", color=CAT_COLOR["漁港"])
ax.text(7.5, 3.05, f"{n_f_total} 件 / {n_ports_f_l34} 漁港", ha="center", fontsize=14)
ax.text(7.5, 2.05, n_f_kinds_str, ha="center", fontsize=9.5)
ax.text(7.5, 1.0, "L33 14 漁港のうち 7 漁港のみ整備\n(半数は交通施設データなし)",
        ha="center", fontsize=8.5, color="#666")

ax.text(5, 0.25,
        f"合計 {n_total_all} 件 / {n_ports_h_l34 + n_ports_f_l34} 港 — "
        f"L32/L33 の 41 港のうち {n_total_pref - n_l34_covered} 港は L34 未整備",
        ha="center", fontsize=11, fontweight="bold")

# 右: 施設種類比較バー
ax = axes[1]
labels = KIND_ORDER
x = np.arange(len(labels))
w = 0.4
n_h = [int(pv_kind.loc["港湾", k]) if k in pv_kind.columns else 0 for k in labels]
n_f = [int(pv_kind.loc["漁港", k]) if k in pv_kind.columns else 0 for k in labels]
ax.bar(x - w/2, n_h, w, color=CAT_COLOR["港湾"], label="港湾", edgecolor="#222")
ax.bar(x + w/2, n_f, w, color=CAT_COLOR["漁港"], label="漁港", edgecolor="#222")
for i, (a, b) in enumerate(zip(n_h, n_f)):
    if a > 0: ax.text(i - w/2, a + 3, str(a), ha="center", fontsize=10)
    if b > 0: ax.text(i + w/2, b + 3, str(b), ha="center", fontsize=10)
ax.set_xticks(x)
ax.set_xticklabels(labels, fontsize=11)
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 / "L34_fig1_dataset_overview.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L34_fig1_dataset_overview.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 12. 図 2: 整備格差 (3 層カバー比較) ヒートマップ
# =============================================================================
print("\n[12] 図 2: 3 層カバー", flush=True)
t1 = time.time()

if GAP_OK and len(three_layer_df) > 0:
    # 上位 20 港 (合計件数で) のヒートマップ
    top_for_heatmap = three_layer_df.sort_values("合計", ascending=False).head(25).copy()
    top_for_heatmap["label"] = top_for_heatmap["港湾名"] + " (" + top_for_heatmap["port_category"] + ")"
    mat = top_for_heatmap[["外郭(L32)", "係留(L33)", "交通(L34)"]].values
    mat_norm = mat / mat.max(axis=0, keepdims=True).clip(min=1)  # 列ごと正規化

    fig, axes = plt.subplots(1, 2, figsize=(15, 8))

    # 左: 整備格差ヒートマップ (絶対件数)
    ax = axes[0]
    im = ax.imshow(mat, aspect="auto", cmap="YlOrRd")
    ax.set_xticks(range(3))
    ax.set_xticklabels(["外郭 L32", "係留 L33", "交通 L34"], fontsize=11)
    ax.set_yticks(range(len(top_for_heatmap)))
    ax.set_yticklabels(top_for_heatmap["label"], fontsize=9)
    for i in range(len(top_for_heatmap)):
        for j in range(3):
            v = mat[i, j]
            txt_color = "white" if mat_norm[i, j] > 0.5 else "black"
            ax.text(j, i, str(int(v)), ha="center", va="center",
                    color=txt_color, fontsize=9)
    ax.set_title(f"上位 25 港の 3 層件数 — ヒートマップ\n(数値=件数、色=列内相対)", fontsize=12)
    plt.colorbar(im, ax=ax, label="件数", shrink=0.7)

    # 右: 3 層比のスタックバー (上位 15 港)
    ax = axes[1]
    top15 = three_layer_df.sort_values("合計", ascending=False).head(15).copy()
    top15["label"] = top15["港湾名"] + "(" + top15["port_category"].str[0] + ")"
    y = np.arange(len(top15))
    layers = [("外郭(L32)", "#9467bd"), ("係留(L33)", "#1f77b4"), ("交通(L34)", "#cf222e")]
    left = np.zeros(len(top15))
    for col, color in layers:
        vals = top15[col].values
        ax.barh(y, vals, left=left, color=color, label=col, edgecolor="#222")
        left = left + vals
    for i, total_v in enumerate(top15["合計"]):
        ax.text(total_v + 3, i, f" {int(total_v)}", va="center", fontsize=9)
    ax.set_yticks(y)
    ax.set_yticklabels(top15["label"], fontsize=9.5)
    ax.invert_yaxis()
    ax.set_xlabel("施設件数 (3 層合計)")
    ax.set_title("上位 15 港の 3 層構成 (外郭 + 係留 + 交通)", fontsize=12)
    ax.legend(loc="lower right", fontsize=9)
    ax.grid(axis="x", alpha=0.3)

    plt.tight_layout()
    fig.savefig(ASSETS / "L34_fig2_three_layer_heatmap.png", dpi=130, bbox_inches="tight")
    plt.close("all")
    print(f"  saved L34_fig2_three_layer_heatmap.png ({time.time()-t1:.1f}s)", flush=True)
else:
    print(f"  SKIP: 3 層分析不可", flush=True)


# =============================================================================
# 13. 図 3: 港別ランキング (件数 + 面積)
# =============================================================================
print("\n[13] 図 3: 港別ランキング", flush=True)
t1 = time.time()

top_all = port_agg.sort_values("n_facilities", ascending=False).head(15).copy()
top_all["label_with_cat"] = top_all["港湾名"] + " (" + top_all["port_category"] + ")"

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

# 左: 件数バー
ax = axes[0]
y = np.arange(len(top_all))
colors = [CAT_COLOR[c] for c in top_all["port_category"]]
ax.barh(y, top_all["n_facilities"], color=colors, edgecolor="#222")
for i, n in enumerate(top_all["n_facilities"]):
    ax.text(n + 0.5, i, f" {int(n)}", va="center", fontsize=10)
ax.set_yticks(y)
ax.set_yticklabels(top_all["label_with_cat"], fontsize=10)
ax.invert_yaxis()
ax.set_xlabel("施設数")
ax.set_title(f"臨港交通施設 件数 上位 15 港 (整備済 {n_ports_h_l34+n_ports_f_l34} 港から)")
legend_handles = [Patch(facecolor=CAT_COLOR["港湾"], label="港湾"),
                  Patch(facecolor=CAT_COLOR["漁港"], label="漁港")]
ax.legend(handles=legend_handles, loc="lower right")
ax.grid(axis="x", alpha=0.3)

# 右: 面積バー (geom ありの分のみ)
ax = axes[1]
ax.barh(y, top_all["total_area_ha"], color=colors, edgecolor="#222")
for i, v in enumerate(top_all["total_area_ha"]):
    ax.text(v + 0.05, i, f" {v:.2f}", va="center", fontsize=10)
ax.set_yticks(y)
ax.set_yticklabels(top_all["label_with_cat"], fontsize=10)
ax.invert_yaxis()
ax.set_xlabel("臨港交通 面積 (ha, GIS有のみ)")
ax.set_title("臨港交通 整備面積 上位 15 港")
ax.legend(handles=legend_handles, loc="lower right")
ax.grid(axis="x", alpha=0.3)

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


# =============================================================================
# 14. 図 4: 県全域 臨港交通マップ (POLYGON 描画 + カテゴリ色)
# =============================================================================
print("\n[14] 図 4: 県全域マップ", flush=True)
t1 = time.time()

fig, ax = plt.subplots(figsize=(14, 8))
g_admin_pref.boundary.plot(ax=ax, color="#aaa", linewidth=0.4)
pref_diss.boundary.plot(ax=ax, color="#222", linewidth=1.0)

# POLYGON は塗りつぶしで描画 (面なので)
for cat, color in CAT_COLOR.items():
    sub = gdf[gdf["port_category"] == cat]
    if len(sub) > 0:
        sub.plot(ax=ax, color=color, edgecolor=color, alpha=0.7, linewidth=0.6)

# 全体の bbox は L33 (係留) と同じくらいの瀬戸内海岸線中心
fbox = gdf.total_bounds
pad = 5000
ax.set_xlim(fbox[0]-pad, fbox[2]+pad)
ax.set_ylim(fbox[1]-pad, fbox[3]+pad)

# 上位 8 港の名前ラベル
top_label = port_agg[port_agg["n_with_geom"] > 0].sort_values("n_facilities", ascending=False).head(8)
for _, r in top_label.iterrows():
    sub_p = gdf[(gdf["港湾名"] == r["港湾名"]) & (gdf["port_category"] == r["port_category"])]
    if len(sub_p) == 0: continue
    cx = sub_p.geometry.union_all().centroid
    ax.annotate(r["港湾名"], xy=(cx.x, cx.y), fontsize=10, color="#000",
                ha="center", va="center",
                bbox=dict(boxstyle="round,pad=0.2", fc="white", ec="#444", alpha=0.9))

legend_handles = [
    Patch(facecolor=CAT_COLOR["港湾"], edgecolor=CAT_COLOR["港湾"], alpha=0.7,
          label=f"港湾 ({n_h_total} 件 / {n_ports_h_l34} 港)"),
    Patch(facecolor=CAT_COLOR["漁港"], edgecolor=CAT_COLOR["漁港"], alpha=0.7,
          label=f"漁港 ({n_f_total} 件 / {n_ports_f_l34} 漁港)"),
]
ax.legend(handles=legend_handles, loc="lower left", fontsize=11, title="行政カテゴリ")
ax.set_title(
    f"広島県 臨港交通施設マップ — GIS 有 {len(gdf)} 件 ({n_total_all - len(gdf)} 件は GIS 欠損)",
    fontsize=14)
ax.set_xticks([]); ax.set_yticks([])
ax.set_aspect("equal")
plt.tight_layout()
fig.savefig(ASSETS / "L34_fig4_pref_overview.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L34_fig4_pref_overview.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 15. 図 5: 施設種類 small multiples マップ (3 種)
# =============================================================================
print("\n[15] 図 5: 施設種類 small multiples", flush=True)
t1 = time.time()

# 3 形式 + 全種重ね = 4 panel
fig, axes = plt.subplots(1, 4, figsize=(20, 6))
fbox = gdf.total_bounds
pad = 5000

for i, k in enumerate(KIND_ORDER):
    ax = axes[i]
    g_admin_pref.boundary.plot(ax=ax, color="#bbb", linewidth=0.3)
    pref_diss.boundary.plot(ax=ax, color="#222", linewidth=0.7)
    sub = gdf[gdf["施設種類"] == k]
    n_total_kind = int((ALL["施設種類"] == k).sum())
    if len(sub) > 0:
        sub.plot(ax=ax, color=KIND_COLOR[k], edgecolor=KIND_COLOR[k], linewidth=0.6, alpha=0.85)
    ax.set_xlim(fbox[0]-pad, fbox[2]+pad); ax.set_ylim(fbox[1]-pad, fbox[3]+pad)
    ax.set_title(f"{k}\n(全 {n_total_kind} 件 / GIS有 {len(sub)})",
                 fontsize=12, color=KIND_COLOR[k], fontweight="bold")
    ax.set_xticks([]); ax.set_yticks([]); ax.set_aspect("equal")

# 4 番目: 全種重ね
ax = axes[3]
g_admin_pref.boundary.plot(ax=ax, color="#bbb", linewidth=0.3)
pref_diss.boundary.plot(ax=ax, color="#222", linewidth=0.7)
for k in KIND_ORDER:
    sub = gdf[gdf["施設種類"] == k]
    if len(sub) > 0:
        sub.plot(ax=ax, color=KIND_COLOR[k], edgecolor=KIND_COLOR[k], linewidth=0.5, alpha=0.85)
ax.set_xlim(fbox[0]-pad, fbox[2]+pad); ax.set_ylim(fbox[1]-pad, fbox[3]+pad)
ax.set_title(f"3 種重ね\n(GIS有 {len(gdf)} 件)", fontsize=12, fontweight="bold")
ax.set_xticks([]); ax.set_yticks([]); ax.set_aspect("equal")

fig.suptitle("臨港交通施設 種類別 small multiples マップ", fontsize=14, y=1.02)
plt.tight_layout()
fig.savefig(ASSETS / "L34_fig5_kind_smallmult.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L34_fig5_kind_smallmult.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 16. 図 6: 上位 4 港の詳細マップ (3 層重ね)
# =============================================================================
print("\n[16] 図 6: 上位 4 港 3 層重ね", flush=True)
t1 = time.time()

# 港湾上位 3 + 漁港上位 1
top_h_ports = port_agg[(port_agg["port_category"]=="港湾") & (port_agg["n_with_geom"] > 0)].head(3)["港湾名"].tolist()
top_f_ports = port_agg[(port_agg["port_category"]=="漁港") & (port_agg["n_with_geom"] > 0)].head(1)["港湾名"].tolist()
detail_ports = [(p, "港湾") for p in top_h_ports] + [(p, "漁港") for p in top_f_ports]

fig, axes = plt.subplots(2, 2, figsize=(14, 12))
axes_flat = axes.flatten()

for i, (port, cat) in enumerate(detail_ports):
    ax = axes_flat[i]
    sub = gdf[(gdf["港湾名"] == port) & (gdf["port_category"] == cat)]
    if len(sub) == 0:
        ax.text(0.5, 0.5, f"{port}: GIS 無", transform=ax.transAxes, ha="center", fontsize=14)
        continue
    bbox = sub.total_bounds
    pad_p = max((bbox[2]-bbox[0])*0.15, (bbox[3]-bbox[1])*0.15, 200)
    g_admin_pref.boundary.plot(ax=ax, color="#aaa", linewidth=0.3)
    legend_handles_port = []

    # 1. L32 外郭 (背景: 紫線)
    if l32_gdf is not None:
        l32_local = l32_gdf.cx[bbox[0]-pad_p:bbox[2]+pad_p, bbox[1]-pad_p:bbox[3]+pad_p]
        if len(l32_local):
            l32_local.plot(ax=ax, color="#9467bd", linewidth=2.0, alpha=0.85)
            legend_handles_port.append(
                Line2D([0], [0], color="#9467bd", lw=3, label=f"外郭(L32) {len(l32_local)}")
            )
    # 2. L33 係留 (中景: 青線)
    if l33_gdf is not None:
        l33_local = l33_gdf.cx[bbox[0]-pad_p:bbox[2]+pad_p, bbox[1]-pad_p:bbox[3]+pad_p]
        if len(l33_local):
            l33_local.plot(ax=ax, color="#0969da", linewidth=2.0, alpha=0.95)
            legend_handles_port.append(
                Line2D([0], [0], color="#0969da", lw=3, label=f"係留(L33) {len(l33_local)}")
            )
    # 3. L34 交通 (前景: 種類色で塗り)
    for k in KIND_ORDER:
        sk = sub[sub["施設種類"] == k]
        if len(sk):
            sk.plot(ax=ax, color=KIND_COLOR[k], edgecolor=KIND_COLOR[k], linewidth=0.6, alpha=0.85)
            legend_handles_port.append(
                Patch(facecolor=KIND_COLOR[k], edgecolor=KIND_COLOR[k],
                      alpha=0.85, label=f"{k}({len(sk)})")
            )

    ax.set_xlim(bbox[0]-pad_p, bbox[2]+pad_p)
    ax.set_ylim(bbox[1]-pad_p, bbox[3]+pad_p)
    n_total_port = int(((ALL["港湾名"]==port) & (ALL["port_category"]==cat)).sum())
    total_area = sub["area_m2"].sum() / 10000.0
    ax.set_title(f"{port} ({cat}) — 交通全{n_total_port}件 / GIS有{len(sub)} / {total_area:.2f}ha",
                 fontsize=11, color=CAT_COLOR[cat], fontweight="bold")
    ax.legend(handles=legend_handles_port, loc="lower right", fontsize=8)
    ax.set_xticks([]); ax.set_yticks([])
    ax.set_aspect("equal")

fig.suptitle("上位 4 港の 3 層重ね詳細マップ — 紫:外郭 青:係留 色:交通", fontsize=14, y=1.00)
plt.tight_layout()
fig.savefig(ASSETS / "L34_fig6_top_ports_3layer.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L34_fig6_top_ports_3layer.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 17. 図 7: 面積分布 (boxplot + ヒスト + Lorenz)
# =============================================================================
print("\n[17] 図 7: 面積分布", flush=True)
t1 = time.time()

fig, axes = plt.subplots(2, 2, figsize=(13, 9))

# 左上: 種類別の面積 boxplot
ax = axes[0, 0]
data_box = []
labels_box = []
colors_box = []
for k in KIND_ORDER:
    arr = gdf[gdf["施設種類"]==k]["area_m2"].values
    if len(arr) >= 1:
        data_box.append(arr)
        labels_box.append(f"{k}\n(n={len(arr)})")
        colors_box.append(KIND_COLOR[k])
if data_box:
    bp = ax.boxplot(data_box, tick_labels=labels_box, patch_artist=True, showfliers=True)
    for patch, c in zip(bp["boxes"], colors_box):
        patch.set_facecolor(c); patch.set_alpha(0.7)
ax.set_ylabel("面積 [m²]")
ax.set_yscale("log")
ax.set_title("施設種類別 面積分布 (log スケール)")
ax.grid(axis="y", alpha=0.3)
ax.tick_params(axis="x", labelsize=10)

# 右上: 道路 vs 駐車場 のヒスト (港湾 vs 漁港 でダブル)
ax = axes[0, 1]
road_h = gdf[(gdf["port_category"]=="港湾") & (gdf["施設種類"]=="道路-車道")]["area_m2"]
road_f = gdf[(gdf["port_category"]=="漁港") & (gdf["施設種類"]=="道路-車道")]["area_m2"]
park_h = gdf[(gdf["port_category"]=="港湾") & (gdf["施設種類"]=="駐車場")]["area_m2"]
# 全データの最大を見て bin を決める
all_areas = gdf["area_m2"].values
if len(all_areas):
    bins = np.logspace(np.log10(max(1, all_areas.min())), np.log10(all_areas.max()), 25)
else:
    bins = 25
if len(road_h):
    ax.hist(road_h, bins=bins, alpha=0.55, color=CAT_COLOR["港湾"],
            label=f"港湾 道路 (n={len(road_h)})", edgecolor="#222")
if len(road_f):
    ax.hist(road_f, bins=bins, alpha=0.7, color=CAT_COLOR["漁港"],
            label=f"漁港 道路 (n={len(road_f)})", edgecolor="#222")
if len(park_h):
    ax.hist(park_h, bins=bins, alpha=0.55, color="#e7ba52",
            label=f"港湾 駐車場 (n={len(park_h)})", edgecolor="#222")
ax.set_xscale("log")
ax.set_xlabel("面積 [m²] (log)")
ax.set_ylabel("施設数")
ax.set_title("道路 vs 駐車場 の面積分布")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)

# 左下: 港別 累積面積 (Lorenz)
ax = axes[1, 0]
sorted_ports = port_agg[port_agg["total_area_ha"] > 0].sort_values("total_area_ha", ascending=False).reset_index(drop=True)
if len(sorted_ports):
    sorted_ports["cum_pct"] = sorted_ports["total_area_ha"].cumsum() / sorted_ports["total_area_ha"].sum() * 100
    sorted_ports["rank_pct"] = (np.arange(len(sorted_ports)) + 1) / len(sorted_ports) * 100
    ax.plot(sorted_ports["rank_pct"], sorted_ports["cum_pct"], marker="o", color="#cf222e", lw=2)
    ax.plot([0, 100], [0, 100], "--", color="#888", label="完全均等")
    ax.fill_between(sorted_ports["rank_pct"], sorted_ports["cum_pct"], alpha=0.2, color="#cf222e")
    s10 = sorted_ports[sorted_ports["rank_pct"] <= 10]["cum_pct"].iloc[-1] if (sorted_ports["rank_pct"]<=10).any() else None
    s25 = sorted_ports[sorted_ports["rank_pct"] <= 25]["cum_pct"].iloc[-1] if (sorted_ports["rank_pct"]<=25).any() else None
    if s10 is not None:
        ax.axvline(10, color="#444", linestyle=":", alpha=0.5)
        ax.text(11, s10-3, f"上位10%\n{s10:.0f}%集中", fontsize=9, color="#444")
    if s25 is not None:
        ax.axvline(25, color="#444", linestyle=":", alpha=0.5)
        ax.text(26, s25-3, f"上位25%\n{s25:.0f}%", fontsize=9, color="#444")
ax.set_xlabel("港のランク (%)")
ax.set_ylabel("累積面積 (%)")
ax.set_title("臨港交通面積の港間集中 (Lorenz)")
ax.legend()
ax.grid(alpha=0.3)
ax.set_xlim(0, 100); ax.set_ylim(0, 105)

# 右下: 駐車場 1 件あたり面積 散布図 (港 × 件数 × 平均面積)
ax = axes[1, 1]
park_per_port = gdf[gdf["施設種類"]=="駐車場"].groupby(["port_category", "港湾名"]).agg(
    n=("area_m2", "size"),
    mean_m2=("area_m2", "mean"),
    total_m2=("area_m2", "sum"),
).reset_index()
park_per_port["total_ha"] = park_per_port["total_m2"] / 10000
if len(park_per_port):
    for cat, color in CAT_COLOR.items():
        sub_p = park_per_port[park_per_port["port_category"] == cat]
        if len(sub_p):
            ax.scatter(sub_p["n"], sub_p["mean_m2"], s=sub_p["total_ha"]*200+30, alpha=0.65,
                       color=color, edgecolor="#222", label=f"{cat} (n={len(sub_p)})")
    for _, r in park_per_port.sort_values("total_m2", ascending=False).head(5).iterrows():
        ax.annotate(r["港湾名"], xy=(r["n"], r["mean_m2"]), fontsize=9,
                    xytext=(4, 3), textcoords="offset points")
ax.set_xlabel("駐車場の本数")
ax.set_ylabel("駐車場 1 件あたり平均面積 [m²]")
ax.set_yscale("log")
ax.set_title("港ごとの駐車場プロファイル\n(円サイズ = 総面積 ha)")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)

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


# =============================================================================
# 18. 図 8: GIS 欠損パターン
# =============================================================================
print("\n[18] 図 8: GIS 欠損", flush=True)
t1 = time.time()

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

# 左: 種類×カテゴリの欠損率ヒート
ax = axes[0]
miss_pct_pivot = ALL.groupby(["port_category", "施設種類"]).apply(
    lambda x: (x["GIS情報"].isna().sum() / max(1, len(x)) * 100), include_groups=False
).unstack(fill_value=np.nan).reindex(columns=KIND_ORDER)
n_pivot = ALL.groupby(["port_category", "施設種類"]).size().unstack(fill_value=0).reindex(columns=KIND_ORDER, fill_value=0)
mat = miss_pct_pivot.values
im = ax.imshow(mat, cmap="Reds", aspect="auto", vmin=0, vmax=100)
for i in range(mat.shape[0]):
    for j in range(mat.shape[1]):
        v = mat[i, j]
        n = n_pivot.iloc[i, j]
        if not np.isnan(v) and n > 0:
            color = "white" if v > 50 else "black"
            ax.text(j, i, f"{v:.0f}%\n(n={int(n)})", ha="center", va="center",
                    fontsize=11, color=color)
        elif n == 0:
            ax.text(j, i, "(なし)", ha="center", va="center", fontsize=10, color="#999")
ax.set_xticks(range(len(KIND_ORDER)))
ax.set_xticklabels(KIND_ORDER, fontsize=11)
ax.set_yticks(range(len(miss_pct_pivot.index)))
ax.set_yticklabels(miss_pct_pivot.index, fontsize=11)
ax.set_title("GIS 欠損率 — カテゴリ × 施設種類")
plt.colorbar(im, ax=ax, label="欠損率 (%)", shrink=0.7)

# 右: 港別の欠損率バー (件数 5+ 港のみ)
ax = axes[1]
gis_per_port = ALL.groupby(["port_category", "港湾名"]).apply(
    lambda x: pd.Series({"n_total": len(x), "n_missing": x["GIS情報"].isna().sum()}),
    include_groups=False
).reset_index()
gis_per_port["miss_pct"] = (gis_per_port["n_missing"] / gis_per_port["n_total"] * 100).round(1)
gis_per_port = gis_per_port[gis_per_port["n_total"] >= 3].sort_values("miss_pct", ascending=False).head(15)
y = np.arange(len(gis_per_port))
colors_bar = [CAT_COLOR[c] for c in gis_per_port["port_category"]]
ax.barh(y, gis_per_port["miss_pct"], color=colors_bar, edgecolor="#222")
for i, (mp, n_t, n_m) in enumerate(zip(gis_per_port["miss_pct"], gis_per_port["n_total"], gis_per_port["n_missing"])):
    ax.text(mp + 1, i, f" {mp:.0f}% ({int(n_m)}/{int(n_t)})", va="center", fontsize=9)
ax.set_yticks(y)
ax.set_yticklabels([f"{r.港湾名}({r.port_category[0]})" for r in gis_per_port.itertuples()], fontsize=10)
ax.invert_yaxis()
ax.set_xlabel("GIS 欠損率 (%)")
ax.set_xlim(0, 110)
ax.set_title("港別 GIS 欠損率 (件数 3+ の港、上位 15)")
ax.legend(handles=[Patch(facecolor=CAT_COLOR["港湾"], label="港湾"),
                   Patch(facecolor=CAT_COLOR["漁港"], label="漁港")], loc="lower right")
ax.grid(axis="x", alpha=0.3)

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


# =============================================================================
# 19. 図 9: 3 層空間関係 (広島港 + 福山港 + 尾道糸崎港)
# =============================================================================
print("\n[19] 図 9: 3 層空間関係", flush=True)
t1 = time.time()

if L32_OK and L33_OK and l33_gdf is not None and l32_gdf is not None:
    fig, axes = plt.subplots(1, 2, figsize=(16, 7))
    fbox = gdf.total_bounds
    pad = 5000

    # 左: 県全域 — 3 層重ね
    ax = axes[0]
    g_admin_pref.boundary.plot(ax=ax, color="#aaa", linewidth=0.3)
    pref_diss.boundary.plot(ax=ax, color="#222", linewidth=0.8)
    l32_gdf.plot(ax=ax, color="#9467bd", linewidth=1.0, alpha=0.55)  # 外郭=紫
    l33_gdf.plot(ax=ax, color="#0969da", linewidth=1.5, alpha=0.7)   # 係留=青
    gdf.plot(ax=ax, color="#cf222e", edgecolor="#cf222e", linewidth=0.6, alpha=0.85)  # 交通=赤
    ax.set_xlim(fbox[0]-pad, fbox[2]+pad); ax.set_ylim(fbox[1]-pad, fbox[3]+pad)
    near_pct_total = (n_near_moor_h + n_near_moor_f) / max(1, n_total_geo_h + n_total_geo_f) * 100
    ax.set_title(f"県全域 3 層重ね (紫=外郭 / 青=係留 / 赤=交通)\n"
                 f"交通の {PROXIMITY_M:.0f}m 圏内に係留: {near_pct_total:.0f}%",
                 fontsize=11)
    ax.set_xticks([]); ax.set_yticks([]); ax.set_aspect("equal")

    # 右: 広島港 ズーム
    ax = axes[1]
    target = "広島港"
    sub_t = gdf[gdf["港湾名"] == target]
    if len(sub_t) > 0:
        bbox = sub_t.total_bounds
        zoom_pad = max((bbox[2]-bbox[0])*0.15, (bbox[3]-bbox[1])*0.15, 200)
        l32_local = l32_gdf.cx[bbox[0]-zoom_pad:bbox[2]+zoom_pad, bbox[1]-zoom_pad:bbox[3]+zoom_pad]
        l33_local = l33_gdf.cx[bbox[0]-zoom_pad:bbox[2]+zoom_pad, bbox[1]-zoom_pad:bbox[3]+zoom_pad]
        g_admin_pref.boundary.plot(ax=ax, color="#aaa", linewidth=0.3)
        l32_local.plot(ax=ax, color="#9467bd", linewidth=2.0, alpha=0.8)
        l33_local.plot(ax=ax, color="#0969da", linewidth=2.0, alpha=0.95)
        legend_handles_zoom = [
            Line2D([0], [0], color="#9467bd", lw=3, label=f"外郭(L32) {len(l32_local)}"),
            Line2D([0], [0], color="#0969da", lw=3, label=f"係留(L33) {len(l33_local)}"),
        ]
        for k in KIND_ORDER:
            sk = sub_t[sub_t["施設種類"] == k]
            if len(sk):
                sk.plot(ax=ax, color=KIND_COLOR[k], edgecolor=KIND_COLOR[k],
                        linewidth=0.6, alpha=0.85)
                legend_handles_zoom.append(
                    Patch(facecolor=KIND_COLOR[k], alpha=0.85, label=f"{k}({len(sk)})")
                )
        ax.set_xlim(bbox[0]-zoom_pad, bbox[2]+zoom_pad)
        ax.set_ylim(bbox[1]-zoom_pad, bbox[3]+zoom_pad)
        ax.set_title(f"{target} ズーム — 紫:外郭 青:係留 色:交通\n"
                     f"交通 {len(sub_t)} / 係留 {len(l33_local)} / 外郭 {len(l32_local)}",
                     fontsize=11)
        ax.legend(handles=legend_handles_zoom, loc="lower right", fontsize=8)
        ax.set_xticks([]); ax.set_yticks([]); ax.set_aspect("equal")

    plt.tight_layout()
    fig.savefig(ASSETS / "L34_fig9_three_layer_overlay.png", dpi=130, bbox_inches="tight")
    plt.close("all")
    print(f"  saved L34_fig9_three_layer_overlay.png ({time.time()-t1:.1f}s)", flush=True)
else:
    print(f"  SKIP: L32 or L33 not available", flush=True)


# =============================================================================
# 20. 図 10: 事務所別の管理件数
# =============================================================================
print("\n[20] 図 10: 事務所別管理件数", flush=True)
t1 = time.time()

office_pv = ALL.groupby(["事務所", "port_category"]).size().unstack(fill_value=0)
office_pv["合計"] = office_pv.sum(axis=1)
office_pv = office_pv.sort_values("合計", ascending=False)
office_pv.to_csv(ASSETS / "L34_office_breakdown.csv", encoding="utf-8-sig")

fig, ax = plt.subplots(figsize=(10, 5.5))
y = np.arange(len(office_pv))
left = np.zeros(len(office_pv))
for cat in ["港湾", "漁港"]:
    if cat in office_pv.columns:
        vals = office_pv[cat].values
        ax.barh(y, vals, left=left, color=CAT_COLOR[cat], label=cat, edgecolor="#222")
        left = left + vals
for i, total_v in enumerate(office_pv["合計"]):
    ax.text(total_v + 1, i, f" {int(total_v)}", va="center", fontsize=10)
ax.set_yticks(y)
ax.set_yticklabels(office_pv.index, fontsize=11)
ax.invert_yaxis()
ax.set_xlabel("管理件数")
ax.set_title("管理事務所別の臨港交通施設数 (港湾 + 漁港)")
ax.legend(loc="lower right")
ax.grid(axis="x", alpha=0.3)
plt.tight_layout()
fig.savefig(ASSETS / "L34_fig10_office_breakdown.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L34_fig10_office_breakdown.png ({time.time()-t1:.1f}s)", flush=True)


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

hypos = []

# H1: 整備格差
n_gap_h = len(only_l32_l33_h) if GAP_OK else 0
n_gap_f = len(only_l32_l33_f) if GAP_OK else 0
gap_pct = (n_gap_h + n_gap_f) / max(1, n_total_pref) * 100 if GAP_OK else 0
h1 = GAP_OK and (n_gap_h + n_gap_f >= 5)
hypos.append({
    "H": "H1",
    "claim": "L32/L33 41 港のうち 10 港 (≈24%) は L34 未整備、漁港の半数は未整備",
    "result": (f"未整備: 港湾 {n_gap_h}/{len(ports_l32_l33_h) if GAP_OK else 0}, "
               f"漁港 {n_gap_f}/{len(ports_l32_l33_f) if GAP_OK else 0}, "
               f"全体 {gap_pct:.0f}%") if GAP_OK else "L32/L33 不在のため未検証",
    "verdict": "支持" if h1 else ("部分支持" if (n_gap_h + n_gap_f >= 3) else "反証/未検証"),
})

# H2: 件数規模差 (港湾偏重)
h2 = (count_ratio >= 10)
hypos.append({
    "H": "H2",
    "claim": "件数比 港湾:漁港 ≥ 10:1 (L32 1.33:1、L33 1.9:1 より圧倒的偏重)",
    "result": f"港湾 {n_h_total} ({100-fish_pct:.1f}%) / 漁港 {n_f_total} ({fish_pct:.1f}%), 比 {count_ratio:.1f}:1",
    "verdict": "支持" if h2 else ("部分支持" if count_ratio >= 5 else "反証"),
})

# H3: 施設種類の分化 (駐車場が漁港にゼロ)
park_in_h = int(pv_kind.loc["港湾", "駐車場"]) if "駐車場" in pv_kind.columns else 0
park_in_f = int(pv_kind.loc["漁港", "駐車場"]) if "駐車場" in pv_kind.columns else 0
h3 = (park_in_h > 0) and (park_in_f == 0)
hypos.append({
    "H": "H3",
    "claim": "駐車場は港湾のみ整備、漁港には 0 件 (旅客港マーカー)",
    "result": f"駐車場 港湾 {park_in_h} / 漁港 {park_in_f}",
    "verdict": "支持" if h3 else "部分支持",
})

# H4: 上位 3 港集中
top3_h_cnt = int(top5_h.head(3)["n_facilities"].sum())
top3_h_pct = top3_h_cnt / max(1, n_h_total) * 100
h4 = top3_h_pct >= 50
hypos.append({
    "H": "H4",
    "claim": "港湾上位 3 港で 50%+ (L32/L33 より集中)",
    "result": f"上位 3 港 ({', '.join(top5_h.head(3)['港湾名'].tolist())}): "
              f"{top3_h_cnt}/{n_h_total} = {top3_h_pct:.1f}%",
    "verdict": "支持" if h4 else "部分支持",
})

# H5: GIS 欠損の構造的偏り (橋りょう特に多い)
miss_road = (((ALL["施設種類"]=="道路-車道") & ALL["GIS情報"].isna()).sum() / max(1, (ALL["施設種類"]=="道路-車道").sum()) * 100)
miss_park = (((ALL["施設種類"]=="駐車場") & ALL["GIS情報"].isna()).sum() / max(1, (ALL["施設種類"]=="駐車場").sum()) * 100)
miss_bridge = (((ALL["施設種類"]=="橋りょう") & ALL["GIS情報"].isna()).sum() / max(1, (ALL["施設種類"]=="橋りょう").sum()) * 100)
miss_bridge_f = (((ALL["施設種類"]=="橋りょう") & (ALL["port_category"]=="漁港") & ALL["GIS情報"].isna()).sum() /
                  max(1, ((ALL["施設種類"]=="橋りょう") & (ALL["port_category"]=="漁港")).sum()) * 100)
h5 = (miss_bridge > miss_road) and (miss_bridge_f >= 80)
hypos.append({
    "H": "H5",
    "claim": "GIS 欠損は橋りょうで顕著、漁港の橋は 100% 欠損",
    "result": f"道路欠損 {miss_road:.0f}% / 駐車場 {miss_park:.0f}% / 橋りょう {miss_bridge:.0f}% (漁港橋りょう {miss_bridge_f:.0f}%)",
    "verdict": "支持" if h5 else "部分支持",
})

# H6: 交通 → 係留 100m 圏内 (岸壁直結)
near_moor_pct = (n_near_moor_h + n_near_moor_f) / max(1, n_total_geo_h + n_total_geo_f) * 100 if L33_OK else 0
h6 = L33_OK and (near_moor_pct >= 80)
hypos.append({
    "H": "H6",
    "claim": f"交通施設の 80%+ が L33 係留から {PROXIMITY_M:.0f}m 圏内 (岸壁直結)",
    "result": (f"全体 {near_moor_pct:.1f}% (港湾 {n_near_moor_h/max(1,n_total_geo_h)*100:.1f}%, "
               f"漁港 {n_near_moor_f/max(1,n_total_geo_f)*100:.1f}%)" if L33_OK
               else "L33 不在のため未検証"),
    "verdict": "支持" if h6 else ("部分支持" if near_moor_pct >= 60 else ("反証" if L33_OK else "未検証")),
})

# H7: 3 層比例 (外郭:係留:交通)
if GAP_OK and len(three_layer_df) > 0:
    total_outer = int(three_layer_df["外郭(L32)"].sum())
    total_moor = int(three_layer_df["係留(L33)"].sum())
    total_traffic = int(three_layer_df["交通(L34)"].sum())
    ratio_str = f"外郭:係留:交通 = {total_outer}:{total_moor}:{total_traffic}"
    norm = f"正規化 (交通=1): {total_outer/max(1,total_traffic):.1f} : {total_moor/max(1,total_traffic):.1f} : 1"
    # H7: 比は ~5:7:2 → 交通最少
    h7 = (total_traffic < total_outer) and (total_traffic < total_moor)
else:
    ratio_str = "未検証"
    norm = ""
    h7 = False
hypos.append({
    "H": "H7",
    "claim": "3 層比 外郭:係留:交通 ≈ 5:7:2 で交通が最少",
    "result": f"{ratio_str} ({norm})",
    "verdict": "支持" if h7 else ("反証" if GAP_OK else "未検証"),
})

hypos_df = pd.DataFrame(hypos)
hypos_df.to_csv(ASSETS / "L34_hypothesis_results.csv", index=False, encoding="utf-8-sig")
with (ASSETS / "L34_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))
print(f"  ({time.time()-t1:.1f}s)", flush=True)


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

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

sections = []

# 値の整理
n_gap_total = len(only_l32_l33_h) + len(only_l32_l33_f) if GAP_OK else 0
gap_h_list = ", ".join(sorted(only_l32_l33_h)) if GAP_OK else ""
gap_f_list = ", ".join(sorted(only_l32_l33_f)) if GAP_OK else ""
miss_total_pct = ALL["GIS情報"].isna().mean() * 100
near_moor_pct = (n_near_moor_h + n_near_moor_f) / max(1, n_total_geo_h + n_total_geo_f) * 100 if L33_OK else 0


# === Section 1: 学習目標と問い ===
sec1 = f"""
<div class="note">
<b>カバー宣言</b>: 本記事は DoBoX のシリーズ <b>「臨港交通施設_港湾/漁港」 2 件</b>
(dataset_id = 1252, 1256) を統合し、広島県内における<b>港湾アクセス交通インフラ</b>
(港湾道路・港湾橋・港湾駐車場) の<b>陸海接続構造</b>を分析する研究記事です。
L32 (外郭施設) と L33 (係留施設) の続編として、
港湾施設の<b>3 層階層 (外郭→係留→臨港交通)</b> の最後の層を扱います。
</div>

<h3>シリーズ構造判定: (a) 行政カテゴリ分割型 + (b) 整備格差型のハイブリッド</h3>
<p>このシリーズの 2 件は、L32/L33 と同様<b>「港湾」「漁港」の 2 行政カテゴリ</b>で
分割される<b>相補型</b>構造ですが、L32/L33 とは決定的に違う点があります。</p>
<ul>
<li><b>L32 外郭</b>: 27 港湾 + 14 漁港 = 41 港 ({842} 件) <b>すべての港に整備</b></li>
<li><b>L33 係留</b>: 同 41 港 (1,224 件) <b>すべての港に整備</b></li>
<li><b>L34 臨港交通</b>: <b>{n_ports_h_l34} 港湾 + {n_ports_f_l34} 漁港 = {n_ports_h_l34 + n_ports_f_l34} 港のみ</b>
({n_total_all} 件)。<b>残り {n_gap_total} 港は交通施設データなし</b></li>
</ul>
<p>すなわち L34 は<b>「整備された港のみが入っているサブセット型」</b>であり、
これ自体が研究的に重要な発見です。<b>「外郭+係留はあるが臨港交通はない港」</b>とは何か?
これらの港は<b>独立した港湾道路を持たず、既存の集落道路・市道に依存している</b>と推察されます。</p>

<h3>独自用語の定義 (本記事内で固定)</h3>
<table>
<tr><th>用語</th><th>定義</th></tr>
<tr><td>臨港交通施設</td><td>港湾区域・漁港区域に整備された<b>陸上交通インフラの総称</b>。
本データでは<b>「道路-車道」「駐車場」「橋りょう」</b>の 3 種を含む。
<b>港湾の海側 (外郭) と内側 (係留) ではなく、係留の陸側にある</b>第 3 層インフラ。</td></tr>
<tr><td>道路-車道</td><td>港湾区域内に管理者が整備した<b>港湾道路</b>。
一般市道とは別途、港湾管理者 (広島県) が占用許可を出して整備・維持する。
本データに <b>{total_road} 件</b> (港湾 {int(pv_kind.loc['港湾','道路-車道'])} + 漁港 {int(pv_kind.loc['漁港','道路-車道'])})。
<b>POLYGON (面)</b> として記録され、道路の<b>幅 × 長さ</b>の 2 次元範囲を保持。</td></tr>
<tr><td>駐車場</td><td>港湾区域内の<b>船客・荷主・職員用駐車場</b>。
フェリー・客船利用者の自家用車置き場、貨物運送の待機場、港湾労働者の通勤駐車場など。
本データに <b>{total_park} 件</b>、<b>港湾のみ ({int(pv_kind.loc['港湾','駐車場'])} 件)</b> で漁港にはゼロ。
これは「漁港利用者は徒歩・近隣駐車・漁協車両のみ」という機能差を反映。</td></tr>
<tr><td>橋りょう</td><td>港湾区域内に整備された<b>短い橋</b> (陸地と桟橋を繋ぐ橋、
水路を跨ぐ連絡橋など)。一般道路橋とは別途、港湾管理者が整備。
本データに <b>{total_bridge} 件</b>。<b>GIS 欠損率が極めて高い</b>
(全 {total_bridge} 件中 {int(((ALL['施設種類']=='橋りょう') & ALL['GIS情報'].isna()).sum())} 件が GIS 無)。</td></tr>
<tr><td>整備格差</td><td>本記事独自概念。「外郭+係留はあるが、臨港交通はない」港の存在。
広島県では<b>{n_gap_total} 港 ({gap_pct:.0f}%) が L34 未整備</b>。
これらは<b>独立した港湾道路を持たず、既存の市道・町道・集落道路をそのまま利用</b>している港。</td></tr>
<tr><td>3 層階層</td><td>L32 (外郭) → L33 (係留) → L34 (臨港交通) の港湾施設の階層。
<b>外洋から陸地に向かう設計順序</b>。波を外郭で遮り、係留で船を泊め、交通で陸へ繋ぐ。
本記事はこの 3 層が空間的に同心円状に重なることを実証する。</td></tr>
<tr><td>岸壁直結</td><td>本記事独自の判定基準。臨港交通施設 (POLYGON) の境界から
L33 係留施設 (LineString) への最近距離が <b>{PROXIMITY_M:.0f} m 以内</b>のもの。
これは<b>岸壁に船が泊まり、すぐ後ろの道路で貨物が運ばれる</b>という陸海直結の物理的指標。</td></tr>
<tr><td>POLYGON 面積</td><td>L32/L33 が LineString (線) で延長を計算するのに対し、
L34 は POLYGON (面) で<b>面積 (m² / ha)</b>を計算する。これは
道路や駐車場が「線ではなく幅広な面で整備される」という現実を反映した記録方式。</td></tr>
</table>

<h3>研究の問い (RQ)</h3>
<p>広島県の臨港交通施設 ({n_total_all} 件) は、L32/L33 で同定された 41 港集合のうち
<b>どれだけの港にしか整備されていない</b>か? 整備された {n_ports_h_l34 + n_ports_f_l34} 港では、
道路・橋・駐車場が<b>どんな構成</b>で陸海接続を支え、
<b>L32 外郭・L33 係留と空間的にどう重なる</b>のか?</p>

<ol>
<li>整備済み港は何港か? <b>整備格差</b> (L32/L33 にあって L34 にない港) はどれくらいあるか?</li>
<li>施設種類 (道路 / 駐車場 / 橋) はカテゴリ別にどう分布し、駐車場は港湾のみか?</li>
<li>件数規模差は? L32 (1.33:1)、L33 (1.9:1) と比べて港湾偏重はどれくらい強いか?</li>
<li>大規模商業港 (広島港・福山港・尾道糸崎港) は道路面積でも他を引き離すか?</li>
<li>GIS 情報の欠損パターン (なぜ橋は欠損が多いか?)</li>
<li>臨港交通の道路は L33 係留から {PROXIMITY_M:.0f}m 圏内に何 % か (岸壁直結)?</li>
<li>3 層比 (外郭 : 係留 : 交通) はどれくらいの数比か?</li>
</ol>

<h3>仮説 H1〜H7</h3>
<ul>
<li><b>H1 (整備格差・極端)</b>: L32/L33 41 港のうち <b>{n_gap_total} 港 (≈{gap_pct:.0f}%) は L34 未整備</b>。
特に漁港は 14 中 7 のみ整備で、<b>半分は車道・橋すらデータ上ない</b>。
これは「漁港は集落道路を併用」という現実を反映。</li>
<li><b>H2 (件数規模差・最大)</b>: 港湾 {n_h_total} ({100-fish_pct:.0f}%) vs 漁港 {n_f_total} ({fish_pct:.0f}%)、
比 <b>{count_ratio:.1f}:1</b>。L32 (1.33:1)、L33 (1.9:1) と比べ<b>圧倒的な港湾偏重</b>。
理由: 漁港は集落道路併用でデータ上の独立交通施設が極少。</li>
<li><b>H3 (施設種類の分化)</b>: 港湾は 3 種、漁港は 2 種。<b>駐車場は漁港にゼロ</b>。
駐車場 = フェリー・客船客の自家用車対応 = <b>旅客港マーカー</b>。
漁港は不特定多数旅客なし → 駐車場不要。</li>
<li><b>H4 (大規模港集中・最強)</b>: 港湾上位 3 港 (広島・尾道糸崎・福山) で
全体の <b>50%+</b>。L32 (44%)・L33 (46%) よりさらに集中。
道路網は商業ハブ港でしか面的に整備されない。</li>
<li><b>H5 (GIS 欠損の構造的偏り)</b>: GIS 情報あり 港湾 22% / 漁港 17% のみ。
<b>「道路-車道」 vs 「橋りょう」で欠損率が異なる</b>。
漁港の橋 6 件は<b>全て GIS なし</b>。データ整備の優先度の現れ。</li>
<li><b>H6 (3 層整合)</b>: 臨港交通 POLYGON は L33 係留から <b>{PROXIMITY_M:.0f}m 以内</b>に
<b>80%+</b> 位置。岸壁の真裏に道路 = 「岸壁直結」が原則。</li>
<li><b>H7 (3 層比例関係)</b>: 1 港の 3 層件数比 ≈ <b>外郭 : 係留 : 交通 = 5 : 7 : 2</b>。
交通が最少。これは「交通施設は港の規模ではなく、貨物量の関数」を示唆。</li>
</ul>

<h3>到達点</h3>
<p>2 dataset_id を「港湾 + 漁港」の相補的な 2 カテゴリとして読み解き、
{n_total_all} 件の臨港交通施設の地理分布・施設種類構成・面積分布・整備格差・
<b>L32 外郭・L33 係留との 3 層階層関係</b>を統合分析する。これにより、
広島県の港湾交通整備が
<b>「外郭で守り、係留で泊め、交通で陸へ繋ぐ」3 層階層設計</b>であり、
かつ<b>整備された港湾の中でもさらに大規模商業港に強く偏在している</b>ことを実データで裏付ける。
さらに、「臨港交通なし」の 10 港の存在自体が<b>港湾整備の階層性</b>を逆照射する。</p>

<div class="warn"><b>本記事のスコープ外</b>:
本記事は臨港交通の<b>地理構造</b>に集中する。
道路の<b>舗装種別・幅員・通行制限</b>、駐車場の<b>収容台数・料金</b>、
橋の<b>桁長・架設年・耐荷重</b>、貨物量・通行台数などは
本記事では扱わず発展課題とする。</div>
"""
sections.append(("学習目標と問い", sec1))


# === Section 2: 使用データ ===
sec2_table = "<table><tr><th>dsid</th><th>カテゴリ</th><th>件数</th><th>形式</th><th>列数</th><th>DoBoX</th></tr>"
for r in series_df.itertuples():
    n_rows = int(ALL[ALL["dsid"] == r.dsid].shape[0])
    sec2_table += (
        f"<tr><td><b>{r.dsid}</b></td><td>{r.category}</td>"
        f"<td>{n_rows}</td>"
        f"<td>CSV (UTF-8 BOM)</td><td>14</td>"
        f"<td><a href='https://hiroshima-dobox.jp/datasets/{r.dsid}' target='_blank'>#{r.dsid}</a></td></tr>"
    )
sec2_table += "</table>"

# 列スキーマ表
schema = [
    ("事業",       "港湾 / 漁港"),
    ("所管",       "港湾 / 漁港 (=事業と同じ)"),
    ("施設分類",   "臨港交通施設 (定数)"),
    ("施設種類",   "道路-車道 / 駐車場 / 橋りょう (港湾 3 種、漁港 2 種)"),
    ("港湾名",     f"港または漁港の名前 ({n_ports_h_l34 + n_ports_f_l34} 種、L33 41 港の部分集合)"),
    ("事務所",     "管理事務所 (例: 三原支所、広島港湾振興事務所、呉支所等)"),
    ("市区町村１/２",  "市町名 (NaN がほとんど)"),
    ("施設番号",   "施設番号 (D-1-3 等の体系。D は道路系 = Road の D ではなく、港湾整備事業区分 D)"),
    ("施設名称",   "施設の固有名 (例: 宇和部北臨港道路、小用桟橋駐車場)"),
    ("管理者名等", "広島県 / NaN 等"),
    ("GIS情報",    "WKT 形式の POLYGON (経度・緯度)、欠損が多い"),
    ("開始位置緯度/経度", "POLYGON の起点の緯度経度"),
]
schema_html = "<table><tr><th>列名</th><th>意味</th></tr>" + "".join(
    f"<tr><td><code>{c}</code></td><td>{m}</td></tr>" for c, m in schema
) + "</table>"

sec2 = f"""
<p>本記事が使用する 2 dataset_id の一覧。L32/L33 と完全に同一スキーマ (14 列) で、
カテゴリ列 <code>事業 = 港湾 / 漁港</code> が二分の唯一の指標です。</p>

{sec2_table}

<h3>列スキーマ (両 dataset 共通)</h3>
{schema_html}

<h3>L32/L33 との重要な相違</h3>
<table>
<tr><th>軸</th><th>L32 外郭</th><th>L33 係留</th><th>L34 臨港交通</th></tr>
<tr><td>件数 (港湾+漁港)</td><td>480 + 362 = 842</td><td>802 + 422 = 1,224</td>
<td>{n_h_total} + {n_f_total} = {n_total_all}</td></tr>
<tr><td>カバー港数</td><td>27 + 14 = 41</td><td>27 + 14 = 41</td>
<td>{n_ports_h_l34} + {n_ports_f_l34} = {n_ports_h_l34 + n_ports_f_l34} <b>(差 {n_gap_total} 港)</b></td></tr>
<tr><td>件数比 港:漁</td><td>1.33:1</td><td>1.9:1</td>
<td><b>{count_ratio:.1f}:1</b> (最大偏重)</td></tr>
<tr><td>GIS 形式</td><td>LINESTRING (線)</td><td>LINESTRING (線)</td>
<td><b>POLYGON (面)</b></td></tr>
<tr><td>主要指標</td><td>延長 m</td><td>延長 m</td><td><b>面積 m²</b></td></tr>
<tr><td>GIS 有効率</td><td>~99% (ほぼ完備)</td><td>~98% (ほぼ完備)</td>
<td><b>{(1-miss_total_pct/100)*100:.0f}% (大量欠損)</b></td></tr>
</table>

<h3>データ品質メモ</h3>
<ul>
<li><b>geom 欠損</b>: 全 {n_total_all} 件のうち <b>{int(ALL['GIS情報'].isna().sum())} 件 ({miss_total_pct:.1f}%)</b> が GIS 欠損。
<b>L32/L33 は 1-2% しか欠損がなかった</b>のと対照的。
これは「臨港交通施設は地番ベースで管理され、空間範囲を測量して登録するインセンティブが弱い」可能性を示唆 (発展課題 Z3)。</li>
<li><b>施設番号体系</b>: <code>D-1-3</code> のような番号。<b>D</b> = 道路系 (港湾整備事業区分の Code D)。
L32 外郭は B 系統、L33 係留は C 系統。これでアルファベット A〜D の 4 体系のうち本記事は D を扱う。</li>
<li><b>geom タイプ</b>: 全て <code>POLYGON</code> または <code>MULTIPOLYGON</code> (面)。
道路・駐車場の幅と長さを面で記録するため。L33 (LineString) と<b>計算する指標が違う</b>:
延長 m → 面積 m² へ。</li>
<li><b>整備格差データ</b>: L32/L33 で扱われた 41 港のうち、L34 に登場しない港は
<b>港湾 {len(only_l32_l33_h)} 港</b> ({gap_h_list[:80] + ('...' if len(gap_h_list) > 80 else '')})、
<b>漁港 {len(only_l32_l33_f)} 港</b> ({gap_f_list[:80] + ('...' if len(gap_f_list) > 80 else '')})。
これは欠損ではなく<b>「臨港交通施設が物理的に存在しない (=既存道路で代替されている)」</b>港。</li>
</ul>

<h3>本記事の主要分析テーブル</h3>
<p>2 dataset を縦結合した <code>L34_all_facilities.csv</code> ({n_total_geo} 行 × geom 有効分のみ) を主軸に、
各分析セクションでクロス集計と GIS 操作 (主に L33 との空間結合) を重ねる。</p>
"""
sections.append(("使用データ", sec2))


# === Section 3: ダウンロード ===
dl_buttons = []
for r in series_df.itertuples():
    rid = [s[4] for s in SERIES if s[0]==r.dsid][0]
    dl_buttons.append(
        f'<li>{r.label} (dsid={r.dsid}, {r.category}): '
        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("L34_series.csv")} — 2 dataset 一覧</li>
<li>{dl("L34_all_facilities.csv")} — 統合臨港交通施設 ({n_total_geo} 行) + area_m2 + dist_to_moor_m + near_moor</li>
<li>{dl("L34_kind_count_cross.csv")} — カテゴリ × 施設種類 (件数)</li>
<li>{dl("L34_kind_area_cross.csv")} — カテゴリ × 施設種類 (面積 ha)</li>
<li>{dl("L34_port_aggregate.csv")} — 港別集計 (件数 + 総面積)</li>
<li>{dl("L34_gis_missing_pattern.csv")} — GIS 欠損率 (種類×カテゴリ)</li>
<li>{dl("L34_layer_coverage_gap.csv")} — L32/L33 vs L34 整備格差表 (41 港の状態)</li>
<li>{dl("L34_three_layer_by_port.csv")} — 港別 3 層件数 (外郭・係留・交通)</li>
<li>{dl("L34_office_breakdown.csv")} — 事務所別管理件数</li>
<li>{dl("L34_hypothesis_results.csv")} — H1〜H7 検証結果</li>
</ul>

<h3>図 (本記事生成 PNG)</h3>
<ul>
<li>{dl("L34_fig1_dataset_overview.png")} / {dl("L34_fig2_three_layer_heatmap.png")} / {dl("L34_fig3_port_ranking.png")}</li>
<li>{dl("L34_fig4_pref_overview.png")} / {dl("L34_fig5_kind_smallmult.png")} / {dl("L34_fig6_top_ports_3layer.png")}</li>
<li>{dl("L34_fig7_area_distribution.png")} / {dl("L34_fig8_gis_missing.png")} / {dl("L34_fig9_three_layer_overlay.png")}</li>
<li>{dl("L34_fig10_office_breakdown.png")}</li>
</ul>

<h3>再現用 Python スクリプト</h3>
<p><a href="L34_port_traffic.py" download>L34_port_traffic.py</a> を取得して
プロジェクトルートで <code>py -X utf8 lessons/L34_port_traffic.py</code> を実行。
データが無ければ自動取得します (L32/L33 の中間 CSV が <code>lessons/assets/</code> に
あれば 3 層関係も自動分析、無ければそれらのセクションのみスキップ)。</p>
"""
sections.append(("ダウンロード", sec3))


# === Section 4: 分析 1 — dataset 構造 ===
sec4 = f"""
<h3>狙い</h3>
<p>「港湾」「漁港」の 2 dataset が、件数・カバー範囲・施設種類構成でどう分化しているかを
1 枚の絵で示す。L32/L33 と比較し、<b>臨港交通固有の極端な偏重</b>を識別する。</p>

<h3>手法 (簡潔に)</h3>
<p>L32/L33 と全く同じ手順: 2 dataset を縦結合し、<code>port_category</code> 列 (港湾 / 漁港) で
件数・港数・施設種類分布を集計。施設種類別件数を二系列バーで比較する。</p>

<h3>実装</h3>
{code('''
SERIES = [
    (1252, "臨港交通施設（港湾）", "港湾", "harbor_traffic_facility.csv",       32496),
    (1256, "臨港交通施設（漁港）", "漁港", "fishing_port_traffic_facility.csv", 32500),
]
dfs = []
for dsid, label, cat, fname, rid in SERIES:
    p = DATA_DIR / fname
    ensure_dataset(p, dataset_id=dsid, resource_id=rid)
    df = pd.read_csv(p, encoding="utf-8-sig")
    df["port_category"] = cat
    df["dsid"] = dsid
    dfs.append(df)

ALL = pd.concat(dfs, ignore_index=True)  # 全 296 行
pv_kind = pd.pivot_table(ALL, index="port_category", columns="施設種類",
                         values="施設番号", aggfunc="count", fill_value=0)
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 2 dataset の規模を文字 + バーで同時に伝える。
左 (カード) は「カテゴリ別の件数+港数+構造形式の内訳」をテキストで、
右 (バー) は「3 種の絶対件数」を視覚化。両方を 1 枚に置くことで
「カテゴリの極端な規模差 (H2)」と「施設種類の分化 (H3)」が同時に読める。</p>
{figure("assets/L34_fig1_dataset_overview.png", "2 dataset の構造概観 — カードビュー (左) と施設種類比較 (右)")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>港湾は {n_h_total} 件 / {n_ports_h_l34} 港、漁港は {n_f_total} 件 / {n_ports_f_l34} 漁港</b>。
件数比 <b>港湾:漁港 = {count_ratio:.1f}:1</b>。
これは L32 外郭 (1.33:1) ・ L33 係留 (1.9:1) と比べて<b>圧倒的な偏重</b>で、H2 を強支持。</li>
<li><b>主役は道路-車道</b>: 港湾 {int(pv_kind.loc['港湾','道路-車道'])} + 漁港 {int(pv_kind.loc['漁港','道路-車道'])}
= {total_road} 件で全体の <b>{total_road/n_total_all*100:.0f}%</b>。
車両通行が港湾アクセスの主機能であることを反映。</li>
<li><b>駐車場は港湾のみ</b>: {int(pv_kind.loc['港湾','駐車場'])} 件すべて港湾、漁港は<b>ゼロ</b>。
これは H3 を支持し、駐車場 = <b>旅客港マーカー</b> という解釈が成り立つ。</li>
<li><b>橋りょうは少数</b>: 全体で {total_bridge} 件 ({total_bridge/n_total_all*100:.1f}%) のみ。
水路・運河を跨ぐ短橋に限定。本格的な道路橋は別途国道・県道として管理されるため。</li>
<li><b>L33 41 港集合との比較</b>: L34 整備済みは {n_ports_h_l34 + n_ports_f_l34} 港のみで、
<b>{n_gap_total} 港 ({gap_pct:.0f}%) は L34 未整備</b>。
これらの未整備港は<b>独立港湾道路を持たず既存道路に依存</b>している (H1 支持)。</li>
</ul>

<h3>表と読み取り (件数クロス)</h3>
{pv_kind.to_html(classes='', border=0)}
<p><b>読み取り</b>: 道路-車道は港湾・漁港両方に多数あり、橋りょうは両方に少数あるが、
<b>駐車場は港湾のみ</b>。漁港の施設構成は「道路 + 橋」の 2 機能に絞られ、
旅客機能 (駐車場) を持たない。これは<b>漁港 = 漁業者専用 + 集落道路で代替</b>という運用思想の現れ。</p>

<h3>表と読み取り (面積クロス)</h3>
{pv_area.to_html(classes='', border=0)}
<p><b>読み取り</b>:</p>
<ul>
<li>港湾の道路-車道は <b>{pv_area.loc['港湾','道路-車道']:.2f} ha</b>。
広島県の港湾内に整備された車道の総面積で、サッカーコート約 {pv_area.loc['港湾','道路-車道']*1.4:.0f} 個分。</li>
<li>港湾の駐車場は <b>{pv_area.loc['港湾','駐車場']:.2f} ha</b>。
1 件あたり平均 {pv_area.loc['港湾','駐車場']*10000/max(1,int(pv_kind.loc['港湾','駐車場'])):.0f} m² で、
標準的なフェリー駐車場のサイズ感。</li>
<li>漁港の合計面積は <b>{total_ha_f:.2f} ha</b> で港湾 ({total_ha_h:.2f} ha) の {total_ha_f/max(0.01,total_ha_h)*100:.1f}%。
件数比 ({count_ratio:.1f}:1) を上回る面積比 {total_ha_h/max(0.01,total_ha_f):.1f}:1 となり、
<b>漁港の交通施設は数だけでなく面積も極小</b>。</li>
</ul>
"""
sections.append(("分析 1: 2 dataset の構造を可視化", sec4))


# === Section 5: 分析 2 — 整備格差 (3 層カバー) ===
gap_html = ""
if GAP_OK:
    gap_table_h = pd.DataFrame({
        "L32/L33 港湾名 (整備済 41 港中)": [", ".join(sorted(only_l32_l33_h)) if only_l32_l33_h else "(なし)"],
    }).T
    # 整備格差表をシンプルに
    gap_summary = f"""
<table>
<tr><th>カテゴリ</th><th>L32/L33 港数</th><th>L34 整備済み</th><th>L34 未整備</th><th>未整備の港名</th></tr>
<tr><td>港湾</td><td>{len(ports_l32_l33_h)}</td><td>{len(both_h)}</td>
<td><b>{len(only_l32_l33_h)} 港</b></td><td>{', '.join(sorted(only_l32_l33_h)) if only_l32_l33_h else '(なし)'}</td></tr>
<tr><td>漁港</td><td>{len(ports_l32_l33_f)}</td><td>{len(both_f)}</td>
<td><b>{len(only_l32_l33_f)} 漁港</b></td><td>{', '.join(sorted(only_l32_l33_f)) if only_l32_l33_f else '(なし)'}</td></tr>
<tr><td>合計</td><td>{n_total_pref}</td><td>{n_l34_covered}</td>
<td><b>{n_gap_total} 港 ({gap_pct:.1f}%)</b></td><td>—</td></tr>
</table>
"""
else:
    gap_summary = "<p>L32/L33 中間 CSV 不在のためスキップ。先に <code>L32_port_breakwaters.py</code> と <code>L33_port_mooring.py</code> を実行してください。</p>"

sec5 = f"""
<h3>狙い</h3>
<p>L34 のシリーズ構造は<b>L32/L33 と異なる</b> — 同じ 41 港集合のうち
<b>整備済の港のサブセット</b>のみが入っている。
どの港が整備済み / 未整備かを特定し、<b>「整備格差」の地理</b>を読む。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p><b>集合演算 (差集合)</b>: 単純な 2 集合の差を取る:</p>
<ul>
<li><b>(L32/L33 港集合) − (L34 港集合) = L34 未整備港</b></li>
<li><b>入力</b>: L33_all_facilities.csv の港名一覧、L34 自身の港名一覧</li>
<li><b>出力</b>: 整備状態 (全層整備 / L32-L33 のみ / 交通のみ) を持つ港リスト</li>
<li><b>限界</b>: 「データに登録されていない」 ≠ 「物理的に存在しない」。
DoBoX の登録漏れの可能性も否定できないが、シリーズの公式定義からは「整備格差」と解釈するのが自然。</li>
</ul>

<h3>実装</h3>
{code('''
l33 = pd.read_csv("lessons/assets/L33_all_facilities.csv", encoding="utf-8-sig")
ports_l33_h = set(l33[l33["port_category"]=="港湾"]["港湾名"].unique())
ports_l34_h = set(ALL[ALL["port_category"]=="港湾"]["港湾名"].unique())

only_l33_h = ports_l33_h - ports_l34_h  # L34 未整備の港湾
both_h = ports_l33_h & ports_l34_h      # 全層整備の港湾
''')}

<h3>表と読み取り</h3>
{gap_summary}
<p><b>読み取り</b>:</p>
<ul>
<li><b>港湾 {len(only_l32_l33_h)} 港</b>が L34 未整備。
これらは「外郭 + 係留」はあるのに「臨港交通」がない港で、<b>独立した港湾道路を持たない</b>。
<b>{', '.join(sorted(only_l32_l33_h)) if only_l32_l33_h else 'なし'}</b> など、規模の小さい島嶼港・地方港湾が中心。</li>
<li><b>漁港 {len(only_l32_l33_f)} 漁港</b>が L34 未整備、これは <b>14 漁港中 {len(only_l32_l33_f)/14*100:.0f}%</b>。
漁港は半数以上が独立交通施設を持たない。
<b>{', '.join(sorted(only_l32_l33_f)) if only_l32_l33_f else 'なし'}</b> などの小規模漁港。
集落道路で代替され、漁港専用道として登記されていない。</li>
<li><b>整備格差は明確に「規模に従う」</b>: 大規模港 (広島・福山・尾道糸崎) は完全整備、
小規模港・島嶼漁港は道路を持たない。
これは<b>「貨物量に応じた段階的整備」</b>という政策ロジックの実証。</li>
</ul>

<h3>図と読み取り — 3 層件数ヒートマップ</h3>
<p><b>なぜこの図か</b>: 港ごとの 3 層 (外郭/係留/交通) 件数を 1 表で並べ、
「どの港が全層豊富か」「交通だけ薄い港はどこか」を視覚的に発見できる。
ヒートマップ (左) で全体像、スタックバー (右) で 3 層比率を見る。</p>
{figure("assets/L34_fig2_three_layer_heatmap.png", "上位 25 港の 3 層件数ヒートマップ + 上位 15 港の 3 層スタックバー")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>広島港</b>: 3 層全てで首位 (外郭○ + 係留○ + 交通 79 件)。<b>多層豊富</b>な総合港湾。</li>
<li><b>尾道糸崎港</b>: 外郭・係留も多いが、交通も <b>39 件</b>。多島港のため複数の島面に道路。</li>
<li><b>福山港</b>: 工業港湾。外郭・係留に対して交通 21 件は<b>相対的に少ない</b> = 道路網が集約されている。</li>
<li><b>倉橋 (漁港)</b>: 係留が圧倒的に多い (倉橋島の漁港集積) が、<b>交通は 4 件</b>と極少。
漁港の係留は集落道路で繋がるため、独立した「漁港道路」を持たない。</li>
<li>3 層比 <b>外郭 : 係留 : 交通 = {int(three_layer_df['外郭(L32)'].sum()) if len(three_layer_df) else 0} : {int(three_layer_df['係留(L33)'].sum()) if len(three_layer_df) else 0} : {int(three_layer_df['交通(L34)'].sum()) if len(three_layer_df) else 0}</b>。
正規化すると <b>≈ {int(three_layer_df['外郭(L32)'].sum())/max(1,int(three_layer_df['交通(L34)'].sum())):.1f} : {int(three_layer_df['係留(L33)'].sum())/max(1,int(three_layer_df['交通(L34)'].sum())):.1f} : 1</b> で、
交通が最少 (H7 を {hypos[6]['verdict']})。</li>
</ul>
"""
sections.append(("分析 2: 整備格差と 3 層構造 — どの港が 3 層揃うか", sec5))


# === Section 6: 分析 3 — 港別ランキング ===
top10_html = port_agg.sort_values("n_facilities", ascending=False).head(10).to_html(index=False, classes='', border=0)
sec6 = f"""
<h3>狙い</h3>
<p>整備済 {n_ports_h_l34 + n_ports_f_l34} 港の中で「どこが交通施設を最も豊富に持つか」を
件数・面積の 2 軸で順位付け。H4 (上位港集中) を検証。</p>

<h3>手法</h3>
<p>港湾名 + カテゴリでグループ集計。施設数・総面積 (POLYGON 面積、ha 単位) を取り、
件数で降順ソート。上位 15 港を 2 軸バーで比較。</p>

<h3>実装</h3>
{code('''
port_agg = ALL.groupby(["port_category", "港湾名"]).size().reset_index(name="n_facilities")
area_agg = gdf.groupby(["port_category", "港湾名"]).agg(
    total_area_m2=("area_m2", "sum"),
    n_with_geom=("area_m2", "size"),
).reset_index()
port_agg = port_agg.merge(area_agg, on=["port_category", "港湾名"], how="left")
port_agg["total_area_ha"] = port_agg["total_area_m2"].fillna(0) / 10000
port_agg = port_agg.sort_values("n_facilities", ascending=False)
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 件数だけでは「小さい施設多数の港」が上位、面積だけでは
「大型 1 件の港」が上位になる。両方を並べることで両ベクトルの整合・不整合が読める。
特に L34 は GIS 欠損が多いので、面積はサンプル分のみ。</p>
{figure("assets/L34_fig3_port_ranking.png", "臨港交通施設 上位 15 港 — 件数 (左) と総面積 (右、GIS有のみ)")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>{top_all.iloc[0]['港湾名']} ({top_all.iloc[0]['port_category']})</b> が件数で首位 —
<b>{int(top_all.iloc[0]['n_facilities'])} 件</b> ({int(top_all.iloc[0]['n_facilities'])/n_total_all*100:.1f}%)。
これだけで全体の {int(top_all.iloc[0]['n_facilities'])/n_total_all*100:.0f}% を 1 港で占有 —
極端な集中。</li>
<li><b>上位 3 港 ({', '.join(top_all.head(3)['港湾名'].tolist())})</b> で全体の <b>{top_all.head(3)['n_facilities'].sum()/n_total_all*100:.1f}%</b>。
H4 (上位 3 港で 50%+) を <b>{hypos[3]['verdict']}</b>。</li>
<li><b>カテゴリ混在ランキング</b>: 上位 15 港のうち漁港は<b>{(top_all['port_category']=='漁港').sum()} 港のみ</b>。
H2 (港湾偏重) の視覚的確認。</li>
<li><b>L33 ランキングとの比較</b>: L33 上位は広島港・倉橋・福山・尾道糸崎で、
L34 上位は {", ".join(top_all.head(4)['港湾名'].tolist())} で<b>順序が変わる</b>。
特に倉橋 (L33 で 248 件の係留豊富) が L34 では大きく順位を落とすのが特徴的。
これは漁港が「係留はあるが交通はほぼ無い」という構造の証拠。</li>
</ul>

<h3>表と読み取り (上位 10 港)</h3>
{top10_html}

<h3>カテゴリ別の上位集中 (H4 検証)</h3>
<table>
<tr><th>カテゴリ</th><th>上位 N 港</th><th>件数 / カテゴリ計</th><th>シェア (%)</th></tr>
<tr><td>港湾</td><td>上位 3 ({', '.join(top5_h.head(3)['港湾名'].tolist())})</td>
<td>{top5_h.head(3)['n_facilities'].sum()} / {n_h_total}</td>
<td><b>{top5_h.head(3)['n_facilities'].sum()/max(1,n_h_total)*100:.1f}%</b></td></tr>
<tr><td>港湾</td><td>上位 5</td>
<td>{top5_h['n_facilities'].sum()} / {n_h_total}</td>
<td>{top5_h['n_facilities'].sum()/max(1,n_h_total)*100:.1f}%</td></tr>
<tr><td>漁港</td><td>上位 3</td>
<td>{top5_f.head(3)['n_facilities'].sum()} / {n_f_total}</td>
<td>{top5_f.head(3)['n_facilities'].sum()/max(1,n_f_total)*100:.1f}%</td></tr>
</table>
<p><b>読み取り</b>: 港湾上位 3 で全体の <b>{top5_h.head(3)['n_facilities'].sum()/max(1,n_h_total)*100:.0f}%</b>。
これは L32 (44%)・L33 (46%) より<b>強い集中</b>。
道路網は<b>商業ハブ港でしか面的に整備されない</b>という法則の数値的裏付け。</p>
"""
sections.append(("分析 3: 港別ランキングと上位集中", sec6))


# === Section 7: 分析 4 — 県全域マップ ===
sec7 = f"""
<h3>狙い</h3>
<p>件数や面積だけでは見えない<b>地理偏在</b>を、県全域の POLYGON マップで一望する。
瀬戸内海岸 + 主要島嶼のどこに臨港交通が密集し、どこが手薄かを可視化。</p>

<h3>手法</h3>
<p>WKT パース済 GeoDataFrame を EPSG:6671 (平面直角 III 系) に投影し、
<code>geopandas.plot()</code> で<b>面塗り</b>。
L33 (LineString = 線で描画) と違い、L34 は POLYGON で<b>面で描画</b>するため、
広域マップでは小さな道路 POLYGON が点状に見えることに注意。
県全域 polygon (L15 admin_922) を背景に、カテゴリ・施設種類で色分け。</p>

<h3>実装</h3>
{code('''
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=(14, 8))
pref_diss.boundary.plot(ax=ax, color="#222", linewidth=1.0)
for cat, color in CAT_COLOR.items():
    sub = gdf[gdf["port_category"] == cat]
    sub.plot(ax=ax, color=color, edgecolor=color, alpha=0.7, linewidth=0.6)
''')}

<h3>図と読み取り — 県全域</h3>
<p><b>なぜこの図か</b>: 県全域に対する臨港交通 POLYGON の分布密度を、青 (港湾) と緑 (漁港) で
直接見ることで「どこに何カテゴリが集中しているか」を一目で把握できる。
個別の港名ラベルで上位 8 港を識別。</p>
{figure("assets/L34_fig4_pref_overview.png", f"広島県 臨港交通施設マップ — GIS 有 {len(gdf)} 件をカテゴリで色分け")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>広島湾岸 (広島市〜廿日市)</b> に港湾 (青) のクラスタ。広島港の道路網と駐車場群。</li>
<li><b>福山〜尾道〜三原の東部</b>: 港湾 (青) が連続。尾道糸崎港・福山港の連鎖。
ただし<b>L33 と比べて密度がはるかに低い</b> (件数規模が 1/4 程度)。</li>
<li><b>呉湾岸 + 倉橋島・江田島</b>: <b>漁港 (緑) はほとんど見えない</b>。
L33 では緑 (漁港係留) が密集していたエリアだが、L34 では沈黙。
<b>「漁港は集落道路で代替」</b>の視覚的確認。</li>
<li><b>瀬戸内海島嶼 (大崎上島・生口島・因島)</b>: 中田港・小用港・蒲刈港など中規模港の道路が点在。</li>
<li>L32/L33 と<b>同じ瀬戸内海岸線</b>に分布するが、<b>密度が大幅に薄い</b>のが特徴。</li>
</ul>

<h3>図と読み取り — 施設種類 small multiples</h3>
<p><b>なぜこの図か</b>: 県全域マップでは施設種類の偏在が見えない。
3 種を別パネルにすることで「道路は全域、駐車場は港湾のみ、橋は極稀」のような
形式ごとの地理パターンが分離して見える。</p>
{figure("assets/L34_fig5_kind_smallmult.png", "施設種類 3 種別 small multiples — それぞれの偏在パターン")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>道路-車道 ({total_road} 件)</b>: 県全域に広く分布。整備済 31 港すべてに見られる主要施設。</li>
<li><b>駐車場 ({total_park} 件)</b>: <b>主要商業港 (広島・尾道糸崎・福山) に集中</b>。
漁港エリアにはほぼ無い (H3 視覚的支持)。</li>
<li><b>橋りょう ({total_bridge} 件)</b>: GIS 有が極少 (4 件のみ) で、ほぼ点でしか見えない。
データ品質欠落 (H5) を視覚的に確認。</li>
<li>3 種重ね合わせのパネルでも、<b>L33 のような「面的な瀬戸内海岸の係留密集」</b>は見えない。
臨港交通は<b>「特定の主要港にしか面整備されない」</b>のが視覚化されている。</li>
</ul>
"""
sections.append(("分析 4: 県全域マップで臨港交通分布を観る", sec7))


# === Section 8: 分析 5 — 上位 4 港の 3 層詳細 ===
detail_text = "<ul>\n"
for port, cat in detail_ports:
    sub_p = gdf[(gdf["港湾名"]==port) & (gdf["port_category"]==cat)]
    n_total_p = int(((ALL["港湾名"]==port) & (ALL["port_category"]==cat)).sum())
    n_geo_p = len(sub_p)
    area_p = sub_p["area_m2"].sum() / 10000.0
    kinds_p = ALL[(ALL["港湾名"]==port) & (ALL["port_category"]==cat)]["施設種類"].value_counts().to_dict()
    kinds_str = ", ".join(f"{k} {v}" for k, v in kinds_p.items())
    detail_text += (
        f"<li><b>{port} ({cat}, 全{n_total_p}件 / GIS有{n_geo_p}件 / {area_p:.2f}ha)</b>: {kinds_str}。"
    )
    if port == "広島港":
        detail_text += "県内最大の総合港湾。<b>道路+駐車場+橋の 3 種共存</b>。フェリー乗り場を中心に同心円状に駐車場群が配置され、その外側に港湾道路が走る典型構造。3 層 (外郭=北防波堤・宇品防波堤、係留=岸壁・物揚場群、交通=道路+駐車場) が同一エリアに共存し、<b>港湾の標準形</b>を形成。</li>\n"
    elif port == "尾道糸崎港":
        detail_text += "多島港ゆえ<b>道路が島ごとに分散</b>。橋りょうもこの港に集中 (島と島を繋ぐ短橋)。フェリー航路用の駐車場も整備済み。L32/L33 と組み合わせると、各島に小さな「外郭+係留+交通」セットが点在する<b>群島型港湾</b>。</li>\n"
    elif port == "福山港":
        detail_text += "工業港湾。<b>大型岸壁背後の貨物搬出道路</b>が中心。駐車場は労働者用と少数の旅客フェリー対応。鉄鋼・機械貨物の物流動線として道路面積が大きい。</li>\n"
    elif port == "小用港":
        n_park_konoura = int(((ALL['港湾名']=='小用港')&(ALL['施設種類']=='駐車場')).sum())
        n_road_konoura = int(((ALL['港湾名']=='小用港')&(ALL['施設種類']=='道路-車道')).sum())
        detail_text += (f"江田島の主要港。広島港と江田島を結ぶフェリー航路の起点で、"
                        f"<b>駐車場 {n_park_konoura} 件 + 道路 {n_road_konoura} 件</b>"
                        f"が施設の中心。フェリー旅客の自家用車預け場として機能。</li>\n")
    elif port == "蒲刈港":
        detail_text += "島嶼中規模港。下蒲刈・上蒲刈島の連絡。道路が分散配置で、橋も小規模あり。</li>\n"
    elif port == "倉橋":
        detail_text += "倉橋島の漁港集積地 (L33 では係留 248 件で県最大級だった)。L34 では<b>道路 4 件のみ</b>と極少。これは漁港特有で、係留岸壁の背後は集落道路がそのまま機能するため独立な漁港道路がほぼ不要。</li>\n"
    else:
        detail_text += "中規模港の典型構成。</li>\n"
detail_text += "</ul>\n"

sec8 = f"""
<h3>狙い</h3>
<p>上位 4 港 ({", ".join([p for p, c in detail_ports])}) の<b>3 層重ね詳細マップ</b>を並べ、
<b>各港で外郭 (紫線) + 係留 (青線) + 交通 (種類色塗り)</b>がどう配置されているかを精読する。
本記事の主役分析。</p>

<h3>手法</h3>
<p>各港について bbox を取り、L32 外郭 + L33 係留 + L34 交通 の 3 レイヤを<b>同一座標系で重ね描画</b>。
これにより<b>「外郭で守られた港内に係留が並び、その背後に道路と駐車場」</b>という階層が
1 枚の絵で読める。</p>

<h3>実装</h3>
{code('''
# 3 層重ね描画
l32_local = l32_gdf.cx[bbox[0]-pad:bbox[2]+pad, bbox[1]-pad:bbox[3]+pad]
l33_local = l33_gdf.cx[bbox[0]-pad:bbox[2]+pad, bbox[1]-pad:bbox[3]+pad]

l32_local.plot(ax=ax, color="#9467bd", linewidth=2.0)  # 紫: 外郭 (線)
l33_local.plot(ax=ax, color="#0969da", linewidth=2.0)  # 青: 係留 (線)
for k in KIND_ORDER:
    sk = sub[sub["施設種類"] == k]
    sk.plot(ax=ax, color=KIND_COLOR[k], alpha=0.85)  # 色: 交通 (面)
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 3 層を別々に見るのではなく<b>同じ地図上に重ねる</b>ことで、
階層関係 (外郭 → 係留 → 交通) が空間的にどう同心円配置されているかが直感的に読める。
これが本記事の最大の発見ポイント。</p>
{figure("assets/L34_fig6_top_ports_3layer.png", "上位 4 港の 3 層重ね詳細マップ — 紫:外郭 青:係留 色:交通")}
<p><b>読み取り — 1 港ずつ</b>:</p>
{detail_text}
<div class="note"><b>3 層配置の階層パターン (発見)</b>: 上位 4 港の比較から、
3 層の空間配置パターンは少なくとも以下 3 つに分類できる:
<ol>
<li><b>同心円型 (広島港型)</b>: 外側 (外洋) → 外郭 → 係留 → 交通 → 内側 (陸地) と
<b>同心円状に階層配置</b>。総合港湾の標準形。</li>
<li><b>群島分散型 (尾道糸崎型)</b>: 各島に小さな 3 層セットが分散配置。
1 つ 1 つは小さいが、群として大規模港湾と同等の機能を持つ。</li>
<li><b>係留偏重型 (倉橋型)</b>: 係留が圧倒的多数だが、交通はほぼゼロ。
<b>陸への接続は集落道路に依存</b>する漁港特有のパターン。
これが H1 (整備格差) の物理的実態。</li>
</ol>
これらは港の<b>機能</b> (商業 / 漁業) と<b>地形</b> (本土 / 群島) の積で決まる。</div>
"""
sections.append(("分析 5: 上位 4 港の 3 層重ね詳細マップ — 階層配置の発見", sec8))


# === Section 9: 分析 6 — 面積分布 ===
sec9 = f"""
<h3>狙い</h3>
<p>臨港交通施設の<b>面積分布</b>を多角的に観察し、施設種類ごとの
標準サイズ・港間集中度・駐車場プロファイルを統計的に検証する。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p>4 つの可視化手法を統合する:</p>
<ul>
<li><b>(1) Box plot (log 軸)</b>: 3 種の面積分布を 1 枚に集約。
log スケールで小さな道路片 (10 m²) と大型駐車場 (5,000 m²) を同時可視化。</li>
<li><b>(2) ヒストグラム</b>: 道路と駐車場の面積分布を港湾・漁港で重ねる。
log 軸で広い面積レンジを表示。</li>
<li><b>(3) Lorenz 曲線</b>: x 軸 = 港のランク (%)、y 軸 = 累積面積 (%)。
完全均等線 (対角線) からの偏差が集中度。</li>
<li><b>(4) 散布図 (バブル)</b>: 各港の駐車場プロファイル。
1 港 = 1 バブルで「件数 × 平均面積」を散布。</li>
</ul>

<h3>実装</h3>
{code('''
# POLYGON 面積を計算 (EPSG:6671 で m² 単位)
gdf["area_m2"] = gdf.geometry.area

# 駐車場プロファイル
park_per_port = gdf[gdf["施設種類"]=="駐車場"].groupby(["port_category", "港湾名"]).agg(
    n=("area_m2", "size"),
    mean_m2=("area_m2", "mean"),
    total_m2=("area_m2", "sum"),
).reset_index()

# Lorenz 曲線
sorted_ports = port_agg.sort_values("total_area_ha", ascending=False)
sorted_ports["cum_pct"] = sorted_ports["total_area_ha"].cumsum() / sorted_ports["total_area_ha"].sum() * 100
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 4 視点を 2x2 で並べることで、分布形状・分散・集中・港単位プロファイルを
1 枚に統合できる。</p>
{figure("assets/L34_fig7_area_distribution.png", "面積分布 4 視点 — boxplot、ヒスト、Lorenz、駐車場散布図")}
<p><b>読み取り (左上 — 種類別 boxplot)</b>:</p>
<ul>
<li><b>道路-車道</b>: 中央値が最も大きく、面積分布も広い。port_road の長尺な POLYGON。</li>
<li><b>駐車場</b>: 中央値はやや小さいが、ばらつきが大きい (フェリー大型駐車場 vs 港湾事務所駐車場)。</li>
<li><b>橋りょう</b>: サンプル数が少なく統計が不安定だが、中央値は中規模。</li>
</ul>
<p><b>読み取り (右上 — 道路 vs 駐車場 ヒスト)</b>:</p>
<ul>
<li>港湾道路と港湾駐車場は<b>面積レンジが重なる</b> (どちらも数百〜数千 m²)。</li>
<li>漁港道路は港湾道路より<b>小さめにシフト</b>。漁港の道路は短く狭い。</li>
</ul>
<p><b>読み取り (左下 — Lorenz 曲線)</b>:</p>
<ul>
<li>曲線は対角線から<b>大きく上に膨らむ</b> = 強い上位集中。</li>
<li>上位 10% の港 (整備済 31 港のうち 3 港) で面積の<b>大半</b>を占有。
L33 (係留) より集中度が高い。</li>
<li>これは<b>「臨港交通は特定大型港でしか整備されない」</b>べき分布構造。</li>
</ul>
<p><b>読み取り (右下 — 駐車場プロファイル散布図)</b>:</p>
<ul>
<li>右上 (件数多 + 平均大) に<b>広島港・尾道糸崎港</b>: 多数の大型駐車場を抱えるフェリー旅客港。</li>
<li>左下 (件数少 + 平均小) は<b>中規模港</b>: 駐車場 1〜3 件を最小限に保有。</li>
<li>漁港は<b>1 件もプロットなし</b> (駐車場ゼロ、H3 視覚確認)。</li>
</ul>
"""
sections.append(("分析 6: 面積分布の 4 視点解析", sec9))


# === Section 10: 分析 7 — GIS 欠損パターン ===
sec10 = f"""
<h3>狙い</h3>
<p>L34 の<b>GIS 欠損率 {miss_total_pct:.0f}%</b> は L32 (~1%)・L33 (~2%) と比べて異常に高い。
この欠損が<b>「ランダム」</b>か<b>「構造的偏り」</b>かを判定する。
構造的偏りなら、データ整備の優先度・歴史的経緯が読めるはず。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p><b>2 軸クロス集計 + ヒートマップ</b>:</p>
<ol>
<li>カテゴリ × 施設種類で欠損率を計算。</li>
<li>ヒートマップで「どこが欠損集中か」を視覚化。</li>
<li>港別の欠損率もバーで並べ、整備差が「種類差」か「港差」かを切り分け。</li>
</ol>
<ul>
<li><b>入力</b>: 全 296 件の施設、各々の GIS 情報の有無</li>
<li><b>出力</b>: 種類×カテゴリの欠損率 + 港別欠損率</li>
<li><b>仮説</b>: 「データ整備優先度は<b>主要港 → 地方港、道路 → 駐車場 → 橋</b>の順」</li>
</ul>

<h3>実装</h3>
{code('''
gis_pattern = ALL.groupby(["port_category", "施設種類"]).agg(
    n_total=("GIS情報", "size"),
    n_with_gis=("GIS情報", lambda x: x.notna().sum()),
).reset_index()
gis_pattern["欠損率(%)"] = (1 - gis_pattern["n_with_gis"] / gis_pattern["n_total"]) * 100

# 港別の欠損率
gis_per_port = ALL.groupby(["port_category", "港湾名"]).apply(
    lambda x: pd.Series({{"n_total": len(x), "n_missing": x["GIS情報"].isna().sum()}})
).reset_index()
gis_per_port["miss_pct"] = gis_per_port["n_missing"] / gis_per_port["n_total"] * 100
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: ヒートマップ (左) で種類×カテゴリの<b>2 次元偏り</b>を、
バー (右) で港別の<b>1 次元偏り</b>を、両方提示することで欠損の構造的パターンが分離できる。</p>
{figure("assets/L34_fig8_gis_missing.png", "GIS 欠損パターン — 種類×カテゴリ ヒートマップ + 港別バー")}
<p><b>読み取り (左ヒートマップ)</b>:</p>
<ul>
<li><b>道路-車道</b>: 欠損率 港湾 {((ALL[(ALL['港湾名']!='dummy')&(ALL['施設種類']=='道路-車道')&(ALL['port_category']=='港湾')]['GIS情報'].isna().mean())*100):.0f}% / 漁港 {((ALL[(ALL['施設種類']=='道路-車道')&(ALL['port_category']=='漁港')]['GIS情報'].isna().mean())*100):.0f}%。
道路が比較的整備度が高い (主要対象)。</li>
<li><b>駐車場</b>: 欠損率 {((ALL[(ALL['施設種類']=='駐車場')]['GIS情報'].isna().mean())*100):.0f}% (港湾のみ)。
道路と同程度に欠損が多い。</li>
<li><b>橋りょう</b>: <b>欠損率が顕著に高い</b>。漁港の橋 6 件は<b>100% GIS 欠損</b>。
港湾の橋 13 件中 {(ALL[(ALL['施設種類']=='橋りょう')&(ALL['port_category']=='港湾')]['GIS情報'].isna().sum())} 件が欠損。
これは橋が<b>桁・橋台ベースの「点的」管理</b>で、面積測量されていないためと推察。</li>
<li>H5 (橋りょう欠損が顕著) を<b>{hypos[4]['verdict']}</b>。</li>
</ul>
<p><b>読み取り (右バー)</b>:</p>
<ul>
<li>港別では欠損率 0% (完全整備) から 100% (全件欠損) まで広く分布。
特定の港に整備が集中する構造的偏りが見える。</li>
<li>大規模商業港は欠損率が低く (測量済)、小規模港・漁港は欠損率が高い。
これは<b>「重要港から優先測量」の整備順序</b>の現れ。</li>
<li>整備格差 (どの港か) と GIS 整備度 (測量されたか) は<b>強く相関</b>する可能性が高い。</li>
</ul>

<div class="warn"><b>解釈の注意</b>: GIS 欠損 = 「データに緯度経度がない」だけで、
<b>施設自体は実存する</b>。施設番号・名称はあるが、
shapefile での空間範囲が登録されていないだけ。
これは「DoBoX に整備優先度がある」現実の反映であり、
重要港 (広島港・福山港) の道路から先に空間整備が進む。</div>
"""
sections.append(("分析 7: GIS 欠損パターン分析 — データ整備優先度の発見", sec10))


# === Section 11: 分析 8 — 3 層空間関係 (本記事の主結果) ===
if L33_OK:
    near_pct_total = (n_near_moor_h + n_near_moor_f) / max(1, n_total_geo_h + n_total_geo_f) * 100
    near_pct_h = n_near_moor_h / max(1, n_total_geo_h) * 100
    near_pct_f = n_near_moor_f / max(1, n_total_geo_f) * 100
    median_dist_moor = gdf["dist_to_moor_m"].median()
    sec11 = f"""
<h3>狙い</h3>
<p>本記事の<b>主結果</b>: L32 (外郭) と L33 (係留) と L34 (交通) の<b>3 層空間関係</b>を分析する。
仮説 H6: 臨港交通の道路 POLYGON は L33 係留 LineString から <b>{PROXIMITY_M:.0f}m 以内</b>に
80%+ 位置するはず (=岸壁直結アクセス)。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p><b>距離計算 + 閾値判定 (POLYGON ↔ LineString)</b>:</p>
<ol>
<li><b>L33 係留 LineString を 1 つの巨大 MultiLineString に統合</b> (<code>union_all()</code>)。
これにより距離計算が高速化される (各交通施設 vs 1 つの統合線)。</li>
<li>各交通施設の<b>POLYGON 境界</b>から係留 LineString への最近距離を <code>geometry.distance()</code> で計算。
POLYGON が LineString と重なる場合は距離 = 0 m。</li>
<li><b>{PROXIMITY_M:.0f} m 以内</b>を「岸壁直結」と判定。100 m は岸壁の背後 1 ブロック程度の距離で、
「岸壁から運び出された貨物が短距離で道路に乗る」物理的範囲。</li>
</ol>

<ul>
<li><b>入力</b>: 交通 POLYGON {n_total_geo} 件 + 係留 LineString (L33 由来) </li>
<li><b>出力</b>: 各交通施設の <code>dist_to_moor_m</code> + <code>near_moor</code> ブール</li>
<li><b>限界</b>: 「近い」 ≠ 「物理的に直結」。実際は<b>岸壁のヘリ部 (船から降ろす場所) と道路の入口 (車に積む場所)</b>
が連続する必要がある。本分析は<b>位置の必要条件</b>を見ている。</li>
</ul>

<h3>実装</h3>
{code('''
l33 = pd.read_csv("lessons/assets/L33_all_facilities.csv", encoding="utf-8-sig")
l33["geometry"] = l33["GIS情報"].apply(parse_wkt)
l33_gdf = gpd.GeoDataFrame(l33.dropna(subset=["geometry"]),
                           geometry="geometry", crs="EPSG:4326").to_crs("EPSG:6671")

# 係留 LineString を 1 つの巨大 MultiLineString に統合 (高速化)
moor_diss = l33_gdf.geometry.union_all()

# 各交通施設 (POLYGON) → 係留 への最近距離
gdf["dist_to_moor_m"] = gdf.geometry.distance(moor_diss)
gdf["near_moor"] = gdf["dist_to_moor_m"] <= 100.0  # 100m 圏 = 岸壁直結
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 県全域 (左) で 3 層全体の重なりを、広島港ズーム (右) で個別の港内構造を精査する。
紫 (L32 外郭線) → 青 (L33 係留線) → 色塗り (L34 交通面) の 3 色階層が読める。</p>
{figure("assets/L34_fig9_three_layer_overlay.png", "3 層空間関係 — 県全域 (左) + 広島港ズーム (右)")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>交通施設の {near_pct_total:.1f}%</b> が L33 係留から {PROXIMITY_M:.0f}m 圏内 (港湾 {near_pct_h:.1f}% / 漁港 {near_pct_f:.1f}%)。
H6 (80%+) を <b>{hypos[5]['verdict']}</b>。</li>
<li>交通 → 係留の<b>最近距離 中央値: {median_dist_moor:.0f} m</b>。
ほとんどの交通施設は係留のすぐ近く (= 岸壁直後の道路) に位置する。</li>
<li>広島港ズームでは<b>「外郭の内側に係留、係留の陸側に道路+駐車場」</b>という
3 層階層が一目で見える。これは港湾工学の教科書設計の実データ確認。</li>
<li><b>3 層階層の主結論</b>: 41 港の港湾施設は、
<b>外洋から陸へ向かって 外郭 → 係留 → 交通</b>の 3 層が同心円状に重なる。
ただし<b>交通層は半数以上の港で欠落</b> (整備格差 H1)、
存在する港でも<b>主要 3 港に集中</b> (H4)。</li>
<li><b>「外郭→係留→交通」は時系列順序でもある</b>: 整備の歴史としても
先に防波堤 (外郭) を築き、次に岸壁 (係留)、最後に道路 (交通) が整備される。
施設番号 (B → C → D) の<b>アルファベット順</b>がこれを反映 (発展課題 Z2)。</li>
</ul>

<div class="warn"><b>解釈の注意</b>: 「100m 圏内」 = 「物理的に近い」だけで、
実際の貨物動線は<b>岸壁のヘリ位置 × 道路の出入口位置 × 港内動線</b>に依存する。
本分析は階層配置の<b>位置検証</b>であり、動線評価ではない。
本格的な動線評価は<b>ヤード設計 + 貨物量シミュレーション</b>で別途行う必要がある (発展課題 Z3)。</div>
"""
else:
    sec11 = f"""
<h3>狙い</h3>
<p>L33 (係留施設) と本記事の臨港交通の空間関係を分析する予定でしたが、
L33 中間 CSV が見つからないためスキップ。</p>
<p>L33 を先に実行してから本記事を再実行することで、3 層関係分析が有効化されます。</p>
"""
sections.append(("分析 8: 3 層空間関係 (外郭 → 係留 → 交通)", sec11))


# === Section 12: 仮説検証と考察 ===
sec12_table = hypos_df.to_html(index=False, classes='', border=0)
sec12 = f"""
<p>H1〜H7 の検証結果を 1 表で示す。</p>
{sec12_table}

<h3>事務所別の管理件数</h3>
{figure("assets/L34_fig10_office_breakdown.png", "管理事務所別の臨港交通施設数 (港湾 + 漁港)")}
<p><b>読み取り</b>:</p>
<ul>
<li>事務所別では <b>{office_pv.index[0]}</b> が首位 ({int(office_pv.iloc[0]['合計'])} 件)。
これは {office_pv.index[0]} の管轄エリアに大規模商業港が含まれるため。</li>
<li>事務所が<b>港湾 + 漁港の両方</b>を管理する地域もあれば、片方しか持たない地域もある。
管轄ごとの設計思想・整備順序を読む手がかりになる。</li>
</ul>

<h3>総括: 広島県の臨港交通整備思想</h3>
<p>2 dataset から再構成した臨港交通施設の構造分析により、以下の<b>4 つの設計思想</b>が読み取れる。</p>
<ul>
<li><b>(1) 整備格差・極端</b>: L32/L33 41 港のうち <b>{n_gap_total} 港 ({gap_pct:.0f}%) が L34 未整備</b>。
特に漁港は半数以上が独立交通施設なし。これは<b>「漁港は集落道路で代替」</b>という運用思想の物理的実体。</li>
<li><b>(2) 機能的二相設計</b>: 港湾は<b>道路 + 駐車場 + 橋の 3 機能</b>共存、
漁港は<b>道路 + 橋の 2 機能</b>のみ (駐車場ゼロ)。
これは<b>「旅客対応か漁業者専用か」</b>の機能差を反映。
L32 (外郭の二相)、L33 (係留の二相) と整合する一貫した設計。</li>
<li><b>(3) 商業ハブ集中</b>: 整備済 {n_ports_h_l34 + n_ports_f_l34} 港のうち<b>上位 3 港</b>
({", ".join(top_all.head(3)['港湾名'].tolist())}) で
全交通施設の <b>{top_all.head(3)['n_facilities'].sum()/n_total_all*100:.0f}%</b> を占有。
L32 (44%)・L33 (46%) より集中度が高く、<b>道路網は商業ハブでしか面整備されない</b>法則。</li>
<li><b>(4) 3 層階層整備</b>: 交通施設の <b>{near_pct_total if L33_OK else 0:.0f}%</b>
が L33 係留から {PROXIMITY_M:.0f}m 圏内に位置。
<b>「外郭で守り、係留で泊め、交通で陸へ繋ぐ」</b>同心円状階層が実データで確認できた。</li>
</ul>

<p>本記事は<b>「臨港交通施設は港湾施設階層の最終 1 km インフラ」</b>という視点を実データで裏付けた。
{n_total_all} 件の交通施設が、<b>2 つの法的カテゴリ</b>で {n_ports_h_l34 + n_ports_f_l34} 港の<b>陸海接続容量</b>を形成し、
<b>L32 で示した外郭施設の防護網と L33 で示した係留施設の接岸網の陸側</b>に階層配置されている。
この網の幾何構造を理解することは、<b>港湾の物流動線・防災避難・観光客対応の総合分析</b>の出発点である。</p>

<h3>L32 / L33 / L34 の 3 層比較表 (本記事の到達点)</h3>
<table>
<tr><th>軸</th><th>L32 外郭</th><th>L33 係留</th><th>L34 臨港交通</th></tr>
<tr><td>カバー港数</td><td>27 + 14 = 41</td><td>27 + 14 = 41</td>
<td><b>{n_ports_h_l34} + {n_ports_f_l34} = {n_ports_h_l34 + n_ports_f_l34}</b> (10 港未整備)</td></tr>
<tr><td>件数 (港湾+漁港)</td><td>480 + 362 = 842</td><td>802 + 422 = 1,224</td>
<td><b>{n_h_total} + {n_f_total} = {n_total_all}</b></td></tr>
<tr><td>件数比 港:漁</td><td>1.33:1</td><td>1.9:1</td>
<td><b>{count_ratio:.1f}:1</b> (最大偏重)</td></tr>
<tr><td>主役構造</td><td>防波堤 67%</td><td>物揚場 56%</td>
<td><b>道路 {total_road/n_total_all*100:.0f}% + 駐車場 {total_park/n_total_all*100:.0f}%</b></td></tr>
<tr><td>主役 (港湾)</td><td>防波堤</td><td>物揚場 + 浮さん橋 + 岸壁</td>
<td>道路 + 駐車場</td></tr>
<tr><td>主役 (漁港)</td><td>防波堤 + 護岸</td><td>物揚場 (圧倒的)</td>
<td>道路のみ</td></tr>
<tr><td>GIS 形式</td><td>LINESTRING</td><td>LINESTRING</td>
<td><b>POLYGON</b></td></tr>
<tr><td>主要指標</td><td>延長 km</td><td>延長 km</td><td><b>面積 ha</b></td></tr>
<tr><td>GIS 完備率</td><td>~99%</td><td>~98%</td>
<td><b>{(1-miss_total_pct/100)*100:.0f}%</b> (大量欠損)</td></tr>
<tr><td>上位 3 港集中率 (港湾)</td><td>~44%</td><td>~46%</td>
<td><b>{top5_h.head(3)['n_facilities'].sum()/max(1,n_h_total)*100:.0f}%</b> (最強)</td></tr>
<tr><td>下層との空間関係</td><td>—</td><td>外郭から 500m 圏内 89%</td>
<td><b>係留から 100m 圏内 {near_pct_total if L33_OK else 0:.0f}%</b></td></tr>
<tr><td>施設番号 prefix</td><td>B</td><td>C</td><td><b>D</b></td></tr>
</table>
"""
sections.append(("仮説検証と考察", sec12))


# === Section 13: 発展課題 ===
sec13 = f"""
<h3>結果 X1 → 新仮説 Y1 → 課題 Z1: 整備格差の歴史的経緯と政策評価</h3>
<ul>
<li><b>結果 X1</b>: L32/L33 41 港のうち {n_gap_total} 港 ({gap_pct:.0f}%) が L34 未整備。
特に漁港は半数以上が独立交通施設なし。</li>
<li><b>新仮説 Y1</b>: これは<b>整備の優先度</b>の現れであり、
「貨物量・旅客量が一定基準を超えた港のみが交通施設整備の対象になる」運用基準が存在するはず。
未整備港は<b>「交通量がしきい値未満」</b>として独立道路整備が見送られている可能性が高い。</li>
<li><b>課題 Z1</b>: 国交省<b>港湾統計年報</b>で各港の年間貨物量・旅客数を取得し、
本記事の整備済 / 未整備フラグと散布。<b>「貨物量しきい値」</b>を ROC 曲線で同定。
これは港湾整備政策の<b>定量評価</b>の入口になる。
さらに、過去 30 年の整備変遷 (整備計画書) と照合すれば<b>整備順序の歴史</b>が読める。</li>
</ul>

<h3>結果 X2 → 新仮説 Y2 → 課題 Z2: 施設番号 (B/C/D) の付番順と整備時系列</h3>
<ul>
<li><b>結果 X2</b>: 施設番号は L32 が B 系統、L33 が C 系統、L34 が D 系統。
これは港湾整備のアルファベット順と推察される。</li>
<li><b>新仮説 Y2</b>: 港湾整備の歴史的順序は<b>「先に防波堤 (B) を築き、次に岸壁 (C)、最後に道路 (D)」</b>。
施設番号の数字部分 (<code>D-1-3</code> の「1-3」) は整備順序を示すかもしれない。
広島港の場合、明治期に北防波堤 (B-1)、大正期に宇品の岸壁 (C-1)、戦後に港湾道路 (D-1) という順。</li>
<li><b>課題 Z2</b>: 港湾管理者の<b>港湾整備史</b> (記念誌・年報) と本データの施設番号を照合し、
時系列マップで「外郭 → 係留 → 交通」の整備順序を可視化。
3 層の整備時期差 (外郭の方が古いか、同時か) を実証。
これは<b>港湾発達の地理史</b>研究の核となる。</li>
</ul>

<h3>結果 X3 → 新仮説 Y3 → 課題 Z3: GIS 欠損 {miss_total_pct:.0f}% の構造的理由</h3>
<ul>
<li><b>結果 X3</b>: L34 の GIS 欠損率は {miss_total_pct:.0f}% で L32/L33 の 1-2% と比べ異常に高い。
特に漁港の橋りょう 6 件は<b>全て欠損</b>。</li>
<li><b>新仮説 Y3</b>: 道路・駐車場・橋は<b>「占用許可ベースで地番管理」</b>される性質があり、
shapefile での面積測量が後回しになる。一方、防波堤・岸壁は<b>港湾施設構造物として直接測量</b>される。
これが GIS 整備度の差を生む。</li>
<li><b>課題 Z3</b>: 港湾管理者 (広島県) の管理システム文書を調査し、
施設データの<b>取得経路</b>を明らかにする。占用許可台帳 vs 構造物台帳の差を確認し、
GIS 欠損が<b>「データソースの違い」</b>であることを実証する。
これは公共データ整備の<b>メタ研究</b>の入口。</li>
</ul>

<h3>結果 X4 → 新仮説 Y4 → 課題 Z4: 駐車場面積と旅客数の対応</h3>
<ul>
<li><b>結果 X4</b>: 駐車場 {total_park} 件 ({pv_area.loc['港湾','駐車場']:.2f} ha) は<b>港湾のみ</b>に整備され、
特に広島港・尾道糸崎港・小用港 (江田島) などフェリー航路港に集中。</li>
<li><b>新仮説 Y4</b>: 駐車場面積は<b>港の年間フェリー旅客数</b>と強く相関する。
具体的には<b>「ピーク時 1 時間旅客数 × 1 台あたり 25 m²」</b>でほぼ説明できるはず。</li>
<li><b>課題 Z4</b>: 国交省<b>旅客船統計</b>または各フェリー会社の航路時刻表 + 1 便あたり乗用車収容数から
ピーク旅客数を推定し、本記事の駐車場面積と散布。
<b>「ピーク旅客 × 25 m²」モデル</b>の適合度を測る。
過大な駐車場 = 計画時旅客数より少ない実数 = 計画見直し対象 を識別できる。</li>
</ul>

<h3>結果 X5 → 新仮説 Y5 → 課題 Z5: 道路の管理者責任分界</h3>
<ul>
<li><b>結果 X5</b>: 港湾道路は道路法でなく<b>港湾法</b>で整備され、管理者は広島県 (港湾管理者)。
一般市道とは別物。</li>
<li><b>新仮説 Y5</b>: 港湾道路と一般市道の<b>境界線 (=管理責任分界)</b> は実地でどう引かれているか?
災害時の通行止め判断・補修費用負担などで、両者の境界が問題になることがあるはず。</li>
<li><b>課題 Z5</b>: 国交省 GeoBase または OSM の市道データと本記事の港湾道路 POLYGON を空間結合し、
<b>境界エリア</b> (オーバーラップ・空隙) を抽出。
広島港では国道 487 号と港湾道路の境界がどこかを特定。
これは<b>災害時通行管理</b>の実務に直結する課題で、
災害復旧計画 (L11 ハザード重ね) との統合分析の入口にもなる。</li>
</ul>
"""
sections.append(("発展課題", sec13))


# レンダリング
html = render_lesson(
    num=34,
    title=f"臨港交通施設 2 件統合分析 — 広島県 {n_ports_h_l34} 港湾 + {n_ports_f_l34} 漁港の {n_total_all} 道路・橋・駐車場が描く港湾アクセス第 3 層",
    tags=["臨港交通", "GIS", "港湾", "漁港", "geopandas", "POLYGON", "道路", "駐車場", "3 層階層"],
    time="40 分",
    level="リテラシ",
    data_label="DoBoX 臨港交通施設 2 dataset (1252 港湾, 1256 漁港)",
    sections=sections,
    script_filename="L34_port_traffic.py",
)

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


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