# -*- coding: utf-8 -*-
"""L33 係留施設 2 件統合分析
       — 広島県 27 港湾 + 14 漁港の 1,224 係留施設に見る船舶受入構造

カバー宣言:
  本記事は DoBoX のシリーズ <b>「係留施設_港湾/漁港」 2 件</b>
  (dataset_id = 1251, 1255) を統合し、
  広島県内における<b>港湾係留施設 (=岸壁・物揚場・浮さん橋・船揚場・係船くい・係船浮標・さん橋)</b>
  の<b>船舶受入構造</b>を分析する研究記事である。

  L32 (外郭施設 dsid 1250+1254) と同じく、2 件は <b>「港湾」(=県管理一般港湾、27 港) と
  「漁港」(=県管理漁港、14 漁港)</b> という <b>2 つの行政カテゴリ</b>で同じ 41 港集合を
  二分する<b>相補型</b>構造。本記事は<b>外郭 (L32) と係留 (L33) の 2 層関係</b>も主題とし、
  「外郭で守られた港内のどこに、どんな係留施設が配置されるか」を空間結合で読み解く。

研究の問い (RQ):
  広島県内に整備された係留施設 (合計 1,224 件) は、
  「港湾 27 港」と「漁港 14 漁港」の <b>2 行政カテゴリ</b>でどう分布し、
  <b>どの構造形式</b> (岸壁 / 物揚場 / 浮さん橋 / 船揚場 / さん橋 / 係船くい / 係船浮標)
  で<b>船舶接岸・荷役・漁船係留</b>を支えているか? L32 の外郭施設が形成する
  <b>防護網の内側</b>にこれら係留がどう配置されているか?
    (1) 港湾 vs 漁港 で<b>施設数・延長・構造形式の構成比</b>はどう違うか?
    (2) <b>主役構造</b>は岸壁 (大型船受入) か、物揚場 (小型漁船・荷役) か?
        L32 では「防波堤」が主役 (67%) だったが、係留では?
    (3) 大型岸壁 (>200m) は<b>どの港にどれだけ</b>集中しているか?
        商業港 (広島・福山・尾道糸崎) で長大岸壁の優位は確認できるか?
    (4) <b>浮さん橋</b> (フロート式桟橋) は瀬戸内の潮位差 (3 m+) に対応する構造で、
        港湾と漁港のどちらに多いか? 設計選択の論理は何か?
    (5) <b>外郭 (L32) と係留 (L33) の空間関係</b>:
        係留施設は外郭施設の<b>陸側 (=守られた港内)</b>にどれくらい配置されているか?
        外郭で守られていない「裸の係留」(沖合や外海面) はあるか?
    (6) 港の<b>係留密度</b> (1 港あたり係留延長 / 海岸線延長) は港湾と漁港でどう違うか?

仮説 H1〜H6:
  H1 (カテゴリ規模差・反転): 港湾 802 件 (66%) vs 漁港 422 件 (34%) で
       港湾優位が <b>外郭 (1.33:1) より強い</b> (1.9:1)。
       これは<b>商業港湾は 1 港あたり多数のバースが必要</b>な一方、
       漁港は規模が小さく係留施設も少なくて済むため。延長比は 2 倍以上。
  H2 (構造形式の意外な主役): <b>岸壁は 80 件 (6.5%) と少数</b>、
       <b>主役は物揚場 (689 件 / 56.3%)</b>。これは港湾工学の常識 (大型岸壁が主役)
       とは異なり、<b>「広島県の係留 = 小型船受入の物揚場が圧倒的多数」</b>という構造。
       岸壁は大型船専用の高規格構造で、相対的に少数。
  H3 (浮さん橋の港湾偏在): 浮さん橋 370 件は<b>港湾 (253) > 漁港 (117)</b>。
       潮汐変動 (瀬戸内で 3-4 m) に対応する<b>フロート式桟橋</b>は、
       旅客船・小型船バースとして港湾に多い。漁港は伝統的に物揚場で陸揚げ。
  H4 (主要港集中): 港湾 27 港のうち上位 5 港 (広島・福山・尾道糸崎・呉地区・厳島)
       で全 802 件の <b>50% 以上</b>。漁港 14 漁港のうち上位 3 で 50%+。
       L32 と同じ偏在構造。係留も外郭と同じ「ごく少数の大規模港 + 多数の小規模港」べき分布。
  H5 (大型岸壁の都市港集中): 岸壁 80 件のうち<b>長さ 200 m 超の長大岸壁</b>は、
       広島港・福山港・尾道糸崎港の 3 大商業港に集中。漁港の岸壁 12 件は
       全て 100 m 未満の中型岸壁。
  H6 (外郭の内側に係留): 係留施設の<b>80% 以上</b>が L32 外郭施設から
       <b>500 m 以内</b>に位置する = 外郭で守られた港内に係留が配置されている。
       裸の係留 (外郭 1 km 以遠) は 5% 未満。これは<b>外郭→係留の 2 層整備</b>
       という港湾の階層設計を実データで裏付ける。

要件 S 準拠: 1 分以内完走。1,224 LineString + 外郭 730 LineString の空間結合は
            空間インデックス (gdf.sindex) で対応。
要件 T 準拠: 県全域マップ + 港別小マルチプル + 外郭重ね + 大型岸壁マップ。
要件 Q 準拠: 図 11 種、表 9 種以上 (1,224 施設 × 7 構造形式 × 41 港)。

データ仕様 (2 件は同一スキーマ):
  A. 港湾係留 1251: CSV, 802 行 × 14 列, GIS情報 (LINESTRING / MULTILINESTRING),
     27 港湾, 7 構造形式 (物揚場/浮さん橋/岸壁/船揚場/係船くい/さん橋/係船浮標)
  B. 漁港係留 1255: CSV, 422 行 × 14 列, 同一スキーマ, 14 漁港, 5 構造形式
     (物揚場/浮さん橋/船揚場/岸壁/さん橋)
  共通列: 事業, 所管, 施設分類, 施設種類, 港湾名, 事務所, 市区町村１, 市区町村２,
          施設番号, 施設名称, 管理者名等, GIS情報, 開始位置緯度, 開始位置経度

メモリ対策: 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, LineString, MultiLineString
from shapely.wkt import loads as wkt_loads

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

t0 = time.time()
print("=== L33 係留施設 2 件統合分析 ===", flush=True)

# =============================================================================
# 0. 定数: 2 dataset_id と 7 構造形式の凡例
# =============================================================================
DATA_DIR = ROOT / "data" / "extras" / "L33_port_mooring"
DATA_DIR.mkdir(parents=True, exist_ok=True)
ADMIN_DIR = ROOT / "data" / "extras" / "L15_admin_zones"
L32_DIR = ROOT / "data" / "extras" / "L32_port_breakwaters"
TARGET_CRS = "EPSG:6671"  # JGD2011 平面直角 III, m 単位

SERIES = [
    (1251, "係留施設（港湾）",  "港湾", "harbor_mooring_facility.csv",       32495),
    (1255, "係留施設（漁港）",  "漁港", "fishing_port_mooring_facility.csv", 32499),
]

# 構造形式 (港湾7+漁港5の合計 union = 7 種)
KIND_ORDER = ["岸壁", "物揚場", "浮さん橋", "船揚場", "さん橋", "係船くい", "係船浮標"]
# 物理機能でグルーピング
KIND_FUNCTION = {
    "岸壁":     "大型船バース (商業)",
    "物揚場":   "小型船陸揚げ・荷役",
    "浮さん橋": "フロート桟橋 (潮汐対応)",
    "船揚場":   "小型船陸上引き揚げ",
    "さん橋":   "固定桟橋 (旅客・小型)",
    "係船くい": "杭打ち係留点",
    "係船浮標": "沖合ブイ係留",
}
# 構造形式別の色 (機能ごとに色相を変える)
KIND_COLOR = {
    "岸壁":     "#cf222e",  # 赤 (大型バース、最重要)
    "物揚場":   "#1f77b4",  # 青 (主役 = 数が多い)
    "浮さん橋": "#2ca02c",  # 緑 (フロート構造)
    "船揚場":   "#e7ba52",  # 砂色 (陸揚げ、漁港の伝統)
    "さん橋":   "#9467bd",  # 紫
    "係船くい": "#8c6d31",  # 茶
    "係船浮標": "#17becf",  # 水色 (沖合)
}
# カテゴリ色 (港湾 / 漁港) — L32 と統一
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"L33/{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
# =============================================================================
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["length_m"] = gdf.geometry.length
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())

# 延長 km ベース
pv_len = pd.pivot_table(
    gdf, index="port_category", columns="施設種類",
    values="length_m", aggfunc="sum", fill_value=0,
).reindex(columns=KIND_ORDER, fill_value=0) / 1000.0
pv_len = pv_len.round(2)
pv_len["合計km"] = pv_len.sum(axis=1).round(2)
pv_len.loc["合計"] = pv_len.sum(axis=0).round(2)
print()
print("延長 (km):")
print(pv_len.to_string())
print(f"  ({time.time()-t1:.1f}s)", flush=True)

pv_kind.to_csv(ASSETS / "L33_kind_count_cross.csv", encoding="utf-8-sig")
pv_len.to_csv(ASSETS / "L33_kind_length_cross.csv", encoding="utf-8-sig")


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

port_agg = gdf.groupby(["port_category", "港湾名"]).agg(
    n_facilities=("length_m", "size"),
    total_length_m=("length_m", "sum"),
    n_kinds=("施設種類", "nunique"),
).reset_index()
port_agg["total_length_km"] = (port_agg["total_length_m"] / 1000).round(3)
port_agg = port_agg.sort_values(["port_category", "total_length_km"], ascending=[True, False])
port_agg.to_csv(ASSETS / "L33_port_aggregate.csv", index=False, encoding="utf-8-sig")
print(port_agg.to_string(index=False))

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


# =============================================================================
# 5. 構造形式別の延長統計 (1 施設あたり)
# =============================================================================
print("\n[5] 構造形式別の延長統計", flush=True)
t1 = time.time()

len_stats = gdf.groupby(["port_category", "施設種類"]).agg(
    n=("length_m", "size"),
    median_m=("length_m", "median"),
    mean_m=("length_m", "mean"),
    p25_m=("length_m", lambda x: np.percentile(x, 25)),
    p75_m=("length_m", lambda x: np.percentile(x, 75)),
    max_m=("length_m", "max"),
).round(1).reset_index()
len_stats.to_csv(ASSETS / "L33_length_stats.csv", index=False, encoding="utf-8-sig")
print(len_stats.to_string(index=False))
print(f"  ({time.time()-t1:.1f}s)", flush=True)

# 岸壁 vs 物揚場 の延長中央値
qw_h_med = gdf[(gdf["port_category"]=="港湾") & (gdf["施設種類"]=="岸壁")]["length_m"].median()
qw_f_med = gdf[(gdf["port_category"]=="漁港") & (gdf["施設種類"]=="岸壁")]["length_m"].median()
mp_h_med = gdf[(gdf["port_category"]=="港湾") & (gdf["施設種類"]=="物揚場")]["length_m"].median()
mp_f_med = gdf[(gdf["port_category"]=="漁港") & (gdf["施設種類"]=="物揚場")]["length_m"].median()
print(f"  岸壁 延長中央値: 港湾 {qw_h_med:.1f} m vs 漁港 {qw_f_med:.1f} m", flush=True)
print(f"  物揚場 延長中央値: 港湾 {mp_h_med:.1f} m vs 漁港 {mp_f_med:.1f} m", flush=True)


# =============================================================================
# 6. 大型岸壁 (>=200 m) の検出
# =============================================================================
print("\n[6] 大型岸壁 (>=200 m) の検出", flush=True)
t1 = time.time()

LARGE_QUAY_THRESH = 200.0  # m
quays = gdf[gdf["施設種類"] == "岸壁"].copy()
quays_large = quays[quays["length_m"] >= LARGE_QUAY_THRESH].copy()
quays_large_pv = quays_large.groupby(["port_category", "港湾名"]).agg(
    n=("length_m", "size"),
    total_m=("length_m", "sum"),
    max_m=("length_m", "max"),
).reset_index().sort_values("total_m", ascending=False)
quays_large_pv.to_csv(ASSETS / "L33_large_quays_by_port.csv", index=False, encoding="utf-8-sig")
print(f"  全岸壁 {len(quays)} 件のうち {len(quays_large)} 件 ({len(quays_large)/max(1,len(quays))*100:.1f}%) が {LARGE_QUAY_THRESH:.0f} m 以上", flush=True)
print(quays_large_pv.to_string(index=False))
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. 外郭施設 (L32) との空間関係解析
# =============================================================================
print("\n[8] 外郭施設 (L32) との空間関係", flush=True)
t1 = time.time()

# L32 で生成した CSV を再利用 (外郭施設の WKT を含む)
l32_csv = LESSONS / "assets" / "L32_all_facilities.csv"
PROXIMITY_M = 500.0  # 外郭から係留までの近接距離

if l32_csv.exists():
    l32 = pd.read_csv(l32_csv, encoding="utf-8-sig")
    print(f"  L32 外郭 CSV: {len(l32)} 行", flush=True)
    # WKT が CSV に含まれている前提で再ジオ化
    if "GIS情報" in l32.columns:
        l32["geometry"] = l32["GIS情報"].apply(parse_wkt)
        l32_gdf_raw = gpd.GeoDataFrame(
            l32.dropna(subset=["geometry"]),
            geometry="geometry", crs="EPSG:4326"
        ).to_crs(TARGET_CRS)
        # 外郭の dissolve で巨大 MultiLine 1 つにし、距離計算を高速化
        outer_diss = l32_gdf_raw.geometry.union_all()
        # 各係留施設の中心点から外郭への最近距離 (m)
        # geom がポリゴンでなく LineString なので distance() 使用
        gdf["dist_to_outer_m"] = gdf.geometry.distance(outer_diss)
        gdf["near_outer"] = gdf["dist_to_outer_m"] <= PROXIMITY_M
        n_near_h = int(gdf[(gdf["port_category"]=="港湾") & gdf["near_outer"]].shape[0])
        n_near_f = int(gdf[(gdf["port_category"]=="漁港") & gdf["near_outer"]].shape[0])
        n_total_h = int((gdf["port_category"]=="港湾").sum())
        n_total_f = int((gdf["port_category"]=="漁港").sum())
        print(f"  港湾 係留が外郭から {PROXIMITY_M:.0f} m 以内: "
              f"{n_near_h}/{n_total_h} ({n_near_h/max(1,n_total_h)*100:.1f}%)", flush=True)
        print(f"  漁港 係留が外郭から {PROXIMITY_M:.0f} m 以内: "
              f"{n_near_f}/{n_total_f} ({n_near_f/max(1,n_total_f)*100:.1f}%)", flush=True)
        print(f"  係留 全体の距離 中央値: {gdf['dist_to_outer_m'].median():.1f} m", flush=True)
        L32_OK = True
    else:
        print("  WARN: L32 CSV に GIS情報 列なし、L32 関連分析スキップ", flush=True)
        gdf["dist_to_outer_m"] = np.nan
        gdf["near_outer"] = False
        L32_OK = False
        n_near_h = n_near_f = 0
        n_total_h = int((gdf["port_category"]=="港湾").sum())
        n_total_f = int((gdf["port_category"]=="漁港").sum())
else:
    print(f"  WARN: L32 CSV not found at {l32_csv}, skipping outer-mooring relation", flush=True)
    gdf["dist_to_outer_m"] = np.nan
    gdf["near_outer"] = False
    L32_OK = False
    n_near_h = n_near_f = 0
    n_total_h = int((gdf["port_category"]=="港湾").sum())
    n_total_f = int((gdf["port_category"]=="漁港").sum())
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 9. 係留密度: 港別 1 港あたり係留延長
# =============================================================================
print("\n[9] 係留密度の港別計算", flush=True)
t1 = time.time()

# 各港の bounding box の対角線長 (m) で港の物理サイズを近似
def port_diag(group_geom):
    bbox = group_geom.total_bounds  # [minx, miny, maxx, maxy]
    return ((bbox[2]-bbox[0])**2 + (bbox[3]-bbox[1])**2)**0.5

port_size = []
for (cat, name), grp in gdf.groupby(["port_category", "港湾名"]):
    diag = port_diag(grp)
    total_len = grp["length_m"].sum()
    port_size.append({
        "port_category": cat,
        "港湾名": name,
        "n_facilities": len(grp),
        "total_length_km": total_len / 1000.0,
        "diag_km": diag / 1000.0,
        "density_m_per_km": (total_len / max(1.0, diag/1000.0)),  # m / km diagonal
    })
port_size_df = pd.DataFrame(port_size).sort_values("density_m_per_km", ascending=False)
port_size_df["density_m_per_km"] = port_size_df["density_m_per_km"].round(0)
port_size_df["total_length_km"] = port_size_df["total_length_km"].round(3)
port_size_df["diag_km"] = port_size_df["diag_km"].round(2)
port_size_df.to_csv(ASSETS / "L33_port_density.csv", index=False, encoding="utf-8-sig")
print(port_size_df.head(15).to_string(index=False))
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 / "L33_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 / "L33_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_ports_h = ALL[ALL["port_category"] == "港湾"]["港湾名"].nunique()
n_ports_f = ALL[ALL["port_category"] == "漁港"]["港湾名"].nunique()
n_h_total = int((ALL["port_category"] == "港湾").sum())
n_f_total = int((ALL["port_category"] == "漁港").sum())
total_km_h = pv_len.loc["港湾", "合計km"] if "港湾" in pv_len.index else 0
total_km_f = pv_len.loc["漁港", "合計km"] if "漁港" in pv_len.index else 0
fish_pct = n_f_total / (n_h_total + n_f_total) * 100
length_ratio = total_km_h / max(0.01, total_km_f)


# =============================================================================
# 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 1251\n係留施設（港湾）", ha="center", fontsize=12,
        fontweight="bold", color=CAT_COLOR["港湾"])
ax.text(2.5, 3.05, f"{n_h_total} 件 / {n_ports_h} 港湾", ha="center", fontsize=14)
ax.text(2.5, 2.05, n_h_kinds_str, ha="center", fontsize=8.5)
ax.text(2.5, 1.0, "事務所: 広島港湾・三原・呉・東広島・東部・廿日市", 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 1255\n係留施設（漁港）", ha="center", fontsize=12,
        fontweight="bold", color=CAT_COLOR["漁港"])
ax.text(7.5, 3.05, f"{n_f_total} 件 / {n_ports_f} 漁港", ha="center", fontsize=14)
ax.text(7.5, 2.05, n_f_kinds_str, ha="center", fontsize=8.5)
ax.text(7.5, 1.0, "事務所: 呉・東部・廿日市・東広島・三原・広島港湾", ha="center",
        fontsize=8.5, color="#666")

ax.text(5, 0.25,
        f"合計 {n_total_all} 件 / 41 港 (= {n_ports_h} 港湾 + {n_ports_f} 漁港) — "
        f"L32 外郭と同一の港集合",
        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 + 6, str(a), ha="center", fontsize=9)
    if b > 0: ax.text(i + w/2, b + 6, str(b), ha="center", fontsize=9)
ax.set_xticks(x)
ax.set_xticklabels(labels, fontsize=10, rotation=20)
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 / "L33_fig1_dataset_overview.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L33_fig1_dataset_overview.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 12. 図 2: 構造形式 × カテゴリ 構成比 (100% スタック)
# =============================================================================
print("\n[12] 図 2: 構造形式構成比", flush=True)
t1 = time.time()

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

# 左: 件数 100% スタック
ax = axes[0]
n_h_sum = sum(n_h); n_f_sum = sum(n_f)
ratios_h = [c / max(1, n_h_sum) * 100 for c in n_h]
ratios_f = [c / max(1, n_f_sum) * 100 for c in n_f]
left_h = 0; left_f = 0
for i, k in enumerate(labels):
    ax.barh(["港湾", "漁港"], [ratios_h[i], ratios_f[i]],
            left=[left_h, left_f], color=KIND_COLOR[k], label=k, edgecolor="#222", linewidth=0.5)
    if ratios_h[i] > 4:
        ax.text(left_h + ratios_h[i]/2, 0, f"{ratios_h[i]:.0f}%",
                ha="center", va="center", fontsize=10, color="white", fontweight="bold")
    if ratios_f[i] > 4:
        ax.text(left_f + ratios_f[i]/2, 1, f"{ratios_f[i]:.0f}%",
                ha="center", va="center", fontsize=10, color="white", fontweight="bold")
    left_h += ratios_h[i]; left_f += ratios_f[i]
ax.set_xlabel("構成比 (%)")
ax.set_title(f"構造形式 構成比 (件数ベース)\n港湾 n={n_h_sum} / 漁港 n={n_f_sum}")
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.13), ncol=4, fontsize=9)
ax.set_xlim(0, 100)

# 右: 延長 km 100% スタック
ax = axes[1]
len_h = [pv_len.loc["港湾", k] if k in pv_len.columns and "港湾" in pv_len.index else 0 for k in labels]
len_f = [pv_len.loc["漁港", k] if k in pv_len.columns and "漁港" in pv_len.index else 0 for k in labels]
total_h = sum(len_h); total_f = sum(len_f)
ratios_h2 = [v / max(0.01, total_h) * 100 for v in len_h]
ratios_f2 = [v / max(0.01, total_f) * 100 for v in len_f]
left_h = 0; left_f = 0
for i, k in enumerate(labels):
    ax.barh(["港湾", "漁港"], [ratios_h2[i], ratios_f2[i]],
            left=[left_h, left_f], color=KIND_COLOR[k], edgecolor="#222", linewidth=0.5)
    if ratios_h2[i] > 4:
        ax.text(left_h + ratios_h2[i]/2, 0, f"{ratios_h2[i]:.0f}%",
                ha="center", va="center", fontsize=10, color="white", fontweight="bold")
    if ratios_f2[i] > 4:
        ax.text(left_f + ratios_f2[i]/2, 1, f"{ratios_f2[i]:.0f}%",
                ha="center", va="center", fontsize=10, color="white", fontweight="bold")
    left_h += ratios_h2[i]; left_f += ratios_f2[i]
ax.set_xlabel("延長構成比 (%)")
ax.set_title(f"構造形式 構成比 (延長 km ベース)\n港湾 {total_h:.1f} km / 漁港 {total_f:.1f} km")
ax.set_xlim(0, 100)

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


# =============================================================================
# 13. 図 3: 港別ランキング (上位 15 港、件数 + 延長)
# =============================================================================
print("\n[13] 図 3: 港別ランキング", flush=True)
t1 = time.time()

top_all = port_agg.sort_values("total_length_km", 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 + 1, 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("係留施設 件数 上位 15 港 (カテゴリ色)")
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)

# 右: 延長バー
ax = axes[1]
ax.barh(y, top_all["total_length_km"], color=colors, edgecolor="#222")
for i, v in enumerate(top_all["total_length_km"]):
    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("係留線 総延長 (km)")
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 / "L33_fig3_port_ranking.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L33_fig3_port_ranking.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 14. 図 4: 県全域 係留施設マップ (LineString 描画 + カテゴリ色)
# =============================================================================
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)

for cat, color in CAT_COLOR.items():
    sub = gdf[gdf["port_category"] == cat]
    sub.plot(ax=ax, color=color, linewidth=2.5, alpha=0.95)

fbox = gdf.total_bounds
pad = 5000
ax.set_xlim(fbox[0]-pad, fbox[2]+pad)
ax.set_ylim(fbox[1]-pad, fbox[3]+pad)

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

legend_handles = [
    Line2D([0], [0], color=CAT_COLOR["港湾"], lw=3, label=f"港湾 ({n_total_h} 件 / 27 港)"),
    Line2D([0], [0], color=CAT_COLOR["漁港"], lw=3, label=f"漁港 ({n_total_f} 件 / 14 漁港)"),
]
ax.legend(handles=legend_handles, loc="lower left", fontsize=11, title="行政カテゴリ")
ax.set_title(
    f"広島県 係留施設マップ — 全 {len(gdf)} 件 ({n_total_h} 港湾 + {n_total_f} 漁港)",
    fontsize=14)
ax.set_xticks([]); ax.set_yticks([])
ax.set_aspect("equal")
plt.tight_layout()
fig.savefig(ASSETS / "L33_fig4_pref_overview.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L33_fig4_pref_overview.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 15. 図 5: 構造形式 small multiples マップ (7 種)
# =============================================================================
print("\n[15] 図 5: 構造形式 small multiples", flush=True)
t1 = time.time()

# 7 形式 + 全種重ね = 8 panel
fig, axes = plt.subplots(2, 4, figsize=(18, 9))
axes_flat = axes.flatten()
fbox = gdf.total_bounds
pad = 5000

for i, k in enumerate(KIND_ORDER):
    ax = axes_flat[i]
    g_admin_pref.boundary.plot(ax=ax, color="#bbb", linewidth=0.3)
    pref_diss.boundary.plot(ax=ax, color="#222", linewidth=0.7)
    sub = gdf[gdf["施設種類"] == k]
    if len(sub) > 0:
        # 数が極端に少ない (<10) ものは線太く
        lw = 4.0 if len(sub) < 10 else 2.6
        sub.plot(ax=ax, color=KIND_COLOR[k], linewidth=lw, alpha=0.95)
    ax.set_xlim(fbox[0]-pad, fbox[2]+pad); ax.set_ylim(fbox[1]-pad, fbox[3]+pad)
    ax.set_title(f"{k} ({len(sub)} 件)", fontsize=12, color=KIND_COLOR[k], fontweight="bold")
    ax.set_xticks([]); ax.set_yticks([]); ax.set_aspect("equal")

# 8 番目: 全種重ね
ax = axes_flat[7]
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], linewidth=2.0, 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"7 種重ね ({len(gdf)} 件)", fontsize=12, fontweight="bold")
ax.set_xticks([]); ax.set_yticks([]); ax.set_aspect("equal")

fig.suptitle("係留施設 構造形式別 small multiples マップ — 7 種の地理偏在", fontsize=14, y=1.00)
plt.tight_layout()
fig.savefig(ASSETS / "L33_fig5_kind_smallmult.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L33_fig5_kind_smallmult.png ({time.time()-t1:.1f}s)", flush=True)


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

# 港湾上位 2 + 漁港上位 2 を選定 (実測でランクを取り直す)
top_h_ports = port_agg[port_agg["port_category"]=="港湾"].head(2)["港湾名"].tolist()
top_f_ports = port_agg[port_agg["port_category"]=="漁港"].head(2)["港湾名"].tolist()
detail_ports = [(p, "港湾") for p in top_h_ports] + [(p, "漁港") for p in top_f_ports]

fig, axes = plt.subplots(2, 2, figsize=(13, 11))
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}: 0 件", transform=ax.transAxes, ha="center")
        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 = []
    for k in KIND_ORDER:
        sk = sub[sub["施設種類"] == k]
        if len(sk):
            sk.plot(ax=ax, color=KIND_COLOR[k], linewidth=2.4, alpha=0.95)
            legend_handles_port.append(
                Line2D([0], [0], color=KIND_COLOR[k], lw=3, 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)
    total_len_km = sub["length_m"].sum() / 1000.0
    ax.set_title(f"{port} ({cat}) — {len(sub)} 件 / {total_len_km:.2f} km",
                 fontsize=12, 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 港の係留施設 詳細マップ", fontsize=14, y=1.00)
plt.tight_layout()
fig.savefig(ASSETS / "L33_fig6_top_ports_detail.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L33_fig6_top_ports_detail.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))

# 左上: 岸壁の延長分布 (港湾 vs 漁港 ヒストグラム)
ax = axes[0, 0]
qw_h = gdf[(gdf["port_category"]=="港湾") & (gdf["施設種類"]=="岸壁")]["length_m"]
qw_f = gdf[(gdf["port_category"]=="漁港") & (gdf["施設種類"]=="岸壁")]["length_m"]
mp_h = gdf[(gdf["port_category"]=="港湾") & (gdf["施設種類"]=="物揚場")]["length_m"]
mp_f = gdf[(gdf["port_category"]=="漁港") & (gdf["施設種類"]=="物揚場")]["length_m"]
bins = np.linspace(0, 500, 30)
ax.hist(qw_h, bins=bins, alpha=0.7, color=CAT_COLOR["港湾"], label=f"港湾 岸壁 (n={len(qw_h)})", edgecolor="#222")
if len(qw_f) > 0:
    ax.hist(qw_f, bins=bins, alpha=0.7, color=CAT_COLOR["漁港"], label=f"漁港 岸壁 (n={len(qw_f)})", edgecolor="#222")
if len(qw_h) > 0:
    ax.axvline(qw_h.median(), color=CAT_COLOR["港湾"], linestyle="--",
               label=f"港湾岸壁 中央値 {qw_h.median():.0f}m")
if len(qw_f) > 0:
    ax.axvline(qw_f.median(), color=CAT_COLOR["漁港"], linestyle="--",
               label=f"漁港岸壁 中央値 {qw_f.median():.0f}m")
ax.axvline(LARGE_QUAY_THRESH, color="#cf222e", linestyle=":",
           label=f"大型岸壁閾値 {LARGE_QUAY_THRESH:.0f}m")
ax.set_xlabel("岸壁 1 施設あたり延長 [m]")
ax.set_ylabel("施設数")
ax.set_title("岸壁の延長分布 — 港湾 vs 漁港")
ax.legend(fontsize=8.5)
ax.grid(alpha=0.3)

# 右上: 全構造形式の boxplot
ax = axes[0, 1]
data_box = []
labels_box = []
colors_box = []
for k in KIND_ORDER:
    arr = gdf[gdf["施設種類"]==k]["length_m"].values
    if len(arr) >= 1:
        data_box.append(arr)
        labels_box.append(f"{k}\n(n={len(arr)})")
        colors_box.append(KIND_COLOR[k])
bp_plot = ax.boxplot(data_box, tick_labels=labels_box, patch_artist=True, showfliers=True)
for patch, c in zip(bp_plot["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", rotation=15)

# 左下: 港別 累積延長 (Lorenz 曲線)
ax = axes[1, 0]
all_sorted = port_agg.sort_values("total_length_km", ascending=False).reset_index(drop=True)
all_sorted["cum_pct"] = all_sorted["total_length_km"].cumsum() / all_sorted["total_length_km"].sum() * 100
all_sorted["rank_pct"] = (np.arange(len(all_sorted)) + 1) / len(all_sorted) * 100
ax.plot(all_sorted["rank_pct"], all_sorted["cum_pct"], marker="o", color="#cf222e", lw=2)
ax.plot([0, 100], [0, 100], "--", color="#888", label="完全均等")
ax.fill_between(all_sorted["rank_pct"], all_sorted["cum_pct"], alpha=0.18, color="#cf222e")
share_at_10 = all_sorted[all_sorted["rank_pct"] <= 10]["cum_pct"].iloc[-1] if (all_sorted["rank_pct"]<=10).any() else None
ax.axvline(10, color="#444", linestyle=":", alpha=0.5)
share_at_25 = all_sorted[all_sorted["rank_pct"] <= 25]["cum_pct"].iloc[-1] if (all_sorted["rank_pct"]<=25).any() else None
ax.axvline(25, color="#444", linestyle=":", alpha=0.5)
if share_at_10:
    ax.text(11, share_at_10-3, f"上位 10%\n{share_at_10:.0f}% 集中",
            fontsize=10, color="#444")
if share_at_25:
    ax.text(26, share_at_25-3, f"上位 25%\n{share_at_25:.0f}%",
            fontsize=10, 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]
mp_per_port = gdf[gdf["施設種類"]=="物揚場"].groupby(["port_category", "港湾名"]).agg(
    n=("length_m", "size"),
    mean_m=("length_m", "mean"),
    total_m=("length_m", "sum"),
).reset_index()
mp_per_port["total_km"] = mp_per_port["total_m"] / 1000
for cat, color in CAT_COLOR.items():
    sub_p = mp_per_port[mp_per_port["port_category"] == cat]
    ax.scatter(sub_p["n"], sub_p["mean_m"], s=sub_p["total_km"]*100+30, alpha=0.65,
               color=color, edgecolor="#222", label=f"{cat} (n={len(sub_p)})")
for _, r in mp_per_port.sort_values("total_m", ascending=False).head(5).iterrows():
    ax.annotate(r["港湾名"], xy=(r["n"], r["mean_m"]), fontsize=9,
                xytext=(4, 3), textcoords="offset points")
ax.set_xlabel("物揚場の本数")
ax.set_ylabel("物揚場 1 本あたり平均延長 [m]")
ax.set_title("港ごとの物揚場プロファイル\n(円サイズ = 総延長 km)")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)

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


# =============================================================================
# 18. 図 8: 大型岸壁マップ
# =============================================================================
print("\n[18] 図 8: 大型岸壁マップ", flush=True)
t1 = time.time()

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

# 全係留 (背景の薄灰)
gdf.plot(ax=ax, color="#bbb", linewidth=0.6, alpha=0.6)

# 岸壁 (中型) 100<=L<200
quays_mid = quays[(quays["length_m"] >= 100) & (quays["length_m"] < LARGE_QUAY_THRESH)]
quays_mid.plot(ax=ax, color="#e7ba52", linewidth=2.5, alpha=0.95, label=f"中型岸壁 100-200m (n={len(quays_mid)})")

# 大型岸壁 (>=200)
quays_large.plot(ax=ax, color="#cf222e", linewidth=4.0, alpha=0.95, label=f"大型岸壁 ≥{LARGE_QUAY_THRESH:.0f}m (n={len(quays_large)})")

# 大型岸壁にラベル (上位 10)
for _, r in quays_large.sort_values("length_m", ascending=False).head(10).iterrows():
    cx = r.geometry.centroid
    ax.annotate(f"{r['港湾名']} {r['length_m']:.0f}m",
                xy=(cx.x, cx.y), fontsize=8,
                xytext=(5, 5), textcoords="offset points",
                bbox=dict(boxstyle="round,pad=0.15", fc="white", ec="#444", alpha=0.85))

fbox = gdf.total_bounds
pad = 5000
ax.set_xlim(fbox[0]-pad, fbox[2]+pad)
ax.set_ylim(fbox[1]-pad, fbox[3]+pad)
ax.set_title(
    f"大型岸壁 (≥{LARGE_QUAY_THRESH:.0f} m) の地理分布 — "
    f"全岸壁 {len(quays)} 件のうち {len(quays_large)} 件 ({len(quays_large)/max(1,len(quays))*100:.1f}%)",
    fontsize=13)
ax.legend(loc="lower left", fontsize=10)
ax.set_xticks([]); ax.set_yticks([]); ax.set_aspect("equal")
plt.tight_layout()
fig.savefig(ASSETS / "L33_fig8_large_quays.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L33_fig8_large_quays.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 19. 図 9: 外郭 (L32) と係留 (L33) の 2 層関係マップ
# =============================================================================
print("\n[19] 図 9: 外郭と係留の 2 層関係", flush=True)
t1 = time.time()

if L32_OK:
    fig, axes = plt.subplots(1, 2, figsize=(16, 7))
    fbox = gdf.total_bounds
    pad = 5000

    # 左: 県全域 — 外郭 (灰) + 係留 (色)
    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_raw.plot(ax=ax, color="#9467bd", linewidth=1.6, alpha=0.7)  # 外郭=紫
    near = gdf[gdf["near_outer"]]
    far = gdf[~gdf["near_outer"]]
    near.plot(ax=ax, color="#1f77b4", linewidth=2.0, alpha=0.92)  # 守られた係留=青
    far.plot(ax=ax, color="#cf222e", linewidth=2.5, alpha=0.95)  # 裸の係留=赤
    ax.set_xlim(fbox[0]-pad, fbox[2]+pad); ax.set_ylim(fbox[1]-pad, fbox[3]+pad)
    ax.set_title(f"県全域 — 外郭 (紫) と係留 (青=守られた / 赤=外郭から{PROXIMITY_M:.0f}m超)\n"
                 f"守られた係留: {n_near_h+n_near_f}/{n_total_h+n_total_f} "
                 f"({(n_near_h+n_near_f)/max(1,n_total_h+n_total_f)*100:.1f}%)",
                 fontsize=11)
    ax.set_xticks([]); ax.set_yticks([]); ax.set_aspect("equal")

    # 右: 広島港 ズーム
    ax = axes[1]
    target = "広島港"
    sub_m = gdf[gdf["港湾名"] == target]
    if len(sub_m) > 0:
        bbox = sub_m.total_bounds
        zoom_pad = max((bbox[2]-bbox[0])*0.1, (bbox[3]-bbox[1])*0.1, 200)
        # L32 のうち同じ bbox 圏 (cx で空間インデックス活用、高速)
        l32_local = l32_gdf_raw.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.5, alpha=0.85)
        # 構造形式ごとに色分け (Patch handle で凡例を作る — geopandas Legend 警告回避)
        legend_handles_zoom = [Line2D([0], [0], color="#9467bd", lw=3,
                                      label=f"外郭 (L32, {len(l32_local)})")]
        for k in KIND_ORDER:
            sk = sub_m[sub_m["施設種類"] == k]
            if len(sk):
                sk.plot(ax=ax, color=KIND_COLOR[k], linewidth=2.5, alpha=0.95)
                legend_handles_zoom.append(
                    Line2D([0], [0], color=KIND_COLOR[k], lw=3, 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_m)} 件、外郭 {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 / "L33_fig9_outer_mooring_layers.png", dpi=130, bbox_inches="tight")
    plt.close("all")
    print(f"  saved L33_fig9_outer_mooring_layers.png ({time.time()-t1:.1f}s)", flush=True)
else:
    print(f"  SKIP: L32 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 / "L33_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 + 5, 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 / "L33_fig10_office_breakdown.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L33_fig10_office_breakdown.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 21. 図 11: 機能グルーピング マップ + 港湾 vs 漁港 比較
# =============================================================================
print("\n[21] 図 11: 機能グルーピング + 比較", flush=True)
t1 = time.time()

# 機能 4 グループ
FUNC_GROUPS = {
    "大型船バース": ["岸壁"],
    "小型船・荷役": ["物揚場", "船揚場"],
    "桟橋系": ["浮さん橋", "さん橋"],
    "杭・浮標": ["係船くい", "係船浮標"],
}
FUNC_COLOR = {
    "大型船バース":  "#cf222e",
    "小型船・荷役":  "#1f77b4",
    "桟橋系":       "#2ca02c",
    "杭・浮標":      "#9467bd",
}

fig, axes = plt.subplots(1, 2, figsize=(16, 7.5))

# 左: 機能別マップ
ax = axes[0]
g_admin_pref.boundary.plot(ax=ax, color="#aaa", linewidth=0.4)
pref_diss.boundary.plot(ax=ax, color="#222", linewidth=1.0)
for func, kinds in FUNC_GROUPS.items():
    sub = gdf[gdf["施設種類"].isin(kinds)]
    if len(sub):
        sub.plot(ax=ax, color=FUNC_COLOR[func], linewidth=2.6, alpha=0.95)
fbox = gdf.total_bounds
pad = 5000
ax.set_xlim(fbox[0]-pad, fbox[2]+pad)
ax.set_ylim(fbox[1]-pad, fbox[3]+pad)
legend_handles = [
    Line2D([0], [0], color=FUNC_COLOR[f], lw=3,
           label=f"{f} (n={int((gdf['施設種類'].isin(ks)).sum())})")
    for f, ks in FUNC_GROUPS.items()
]
ax.legend(handles=legend_handles, loc="lower left", fontsize=10, title="機能グループ")
ax.set_title(f"係留施設の機能別マップ — {len(gdf)} 件を 4 機能で色分け", fontsize=13)
ax.set_xticks([]); ax.set_yticks([]); ax.set_aspect("equal")

# 右: 港湾 vs 漁港 機能別構成比
ax = axes[1]
func_pv = pd.DataFrame(0.0, index=["港湾", "漁港"], columns=list(FUNC_GROUPS.keys()))
for func, kinds in FUNC_GROUPS.items():
    for cat in ["港湾", "漁港"]:
        c = ((ALL["port_category"]==cat) & (ALL["施設種類"].isin(kinds))).sum()
        func_pv.loc[cat, func] = c
func_pv_pct = func_pv.div(func_pv.sum(axis=1), axis=0) * 100
left_h = 0; left_f = 0
for func in FUNC_GROUPS:
    rh = func_pv_pct.loc["港湾", func]
    rf = func_pv_pct.loc["漁港", func]
    ax.barh(["港湾", "漁港"], [rh, rf], left=[left_h, left_f],
            color=FUNC_COLOR[func], label=func, edgecolor="#222", linewidth=0.5)
    if rh > 4: ax.text(left_h + rh/2, 0, f"{rh:.0f}%", ha="center", va="center", color="white", fontweight="bold", fontsize=11)
    if rf > 4: ax.text(left_f + rf/2, 1, f"{rf:.0f}%", ha="center", va="center", color="white", fontweight="bold", fontsize=11)
    left_h += rh; left_f += rf
ax.set_xlim(0, 100)
ax.set_xlabel("構成比 (%)")
ax.set_title("機能グループ構成比 — 港湾 vs 漁港")
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.13), ncol=4, fontsize=10)

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


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

hypos = []

# H1: カテゴリ規模差 (港湾 ≥ 漁港、漁港 30%+、延長 1.8x+)
h1 = (n_h_total >= n_f_total) and (length_ratio >= 1.8)
hypos.append({
    "H": "H1",
    "claim": "港湾 ≥ 漁港、外郭 (1.33:1) より港湾優位が強い (1.8x+)、延長は 2 倍以上",
    "result": (f"港湾={n_h_total} ({100-fish_pct:.1f}%), 漁港={n_f_total} ({fish_pct:.1f}%), "
               f"延長 港湾={total_km_h} km / 漁港={total_km_f} km, 比 {length_ratio:.2f}x"),
    "verdict": "支持" if h1 else ("部分支持" if length_ratio >= 1.5 else "反証"),
})

# H2: 構造形式の意外な主役 (岸壁 < 10%, 物揚場 50%+)
total_mp = pv_kind.loc["合計", "物揚場"] if "物揚場" in pv_kind.columns else 0
total_qw = pv_kind.loc["合計", "岸壁"] if "岸壁" in pv_kind.columns else 0
n_all = n_h_total + n_f_total
mp_pct = total_mp / max(1, n_all) * 100
qw_pct = total_qw / max(1, n_all) * 100
h2 = (qw_pct < 10) and (mp_pct >= 50)
hypos.append({
    "H": "H2",
    "claim": "岸壁 < 10%、主役は物揚場 (50%+)",
    "result": f"岸壁 {qw_pct:.1f}% ({total_qw}件) / 物揚場 {mp_pct:.1f}% ({total_mp}件)",
    "verdict": "支持" if h2 else "部分支持",
})

# H3: 浮さん橋 港湾 > 漁港 (港湾の方が件数多い)
fp_h = int(pv_kind.loc["港湾", "浮さん橋"]) if "浮さん橋" in pv_kind.columns else 0
fp_f = int(pv_kind.loc["漁港", "浮さん橋"]) if "浮さん橋" in pv_kind.columns else 0
fp_h_pct = fp_h / max(1, n_h_total) * 100
fp_f_pct = fp_f / max(1, n_f_total) * 100
h3 = fp_h > fp_f
hypos.append({
    "H": "H3",
    "claim": "浮さん橋: 港湾 > 漁港",
    "result": f"港湾 {fp_h} 件 ({fp_h_pct:.1f}%) vs 漁港 {fp_f} 件 ({fp_f_pct:.1f}%)",
    "verdict": "支持" if h3 else "反証",
})

# H4: 上位港集中
top5_h = port_agg[port_agg["port_category"]=="港湾"].head(5)
top5_h_sum = int(top5_h["n_facilities"].sum())
top5_h_pct = top5_h_sum / max(1, port_agg[port_agg["port_category"]=="港湾"]["n_facilities"].sum()) * 100
top3_f = port_agg[port_agg["port_category"]=="漁港"].head(3)
top3_f_sum = int(top3_f["n_facilities"].sum())
top3_f_pct = top3_f_sum / max(1, port_agg[port_agg["port_category"]=="漁港"]["n_facilities"].sum()) * 100
h4 = (top5_h_pct >= 50) and (top3_f_pct >= 50)
hypos.append({
    "H": "H4",
    "claim": "港湾上位 5 港で 50%+、漁港上位 3 漁港で 50%+",
    "result": f"港湾上位 5: {top5_h_pct:.1f}% / 漁港上位 3: {top3_f_pct:.1f}%",
    "verdict": "支持" if h4 else "部分支持",
})

# H5: 大型岸壁 (>=200m) は 3 大商業港に集中
big_quay_ports = quays_large_pv["港湾名"].tolist() if len(quays_large_pv) > 0 else []
expected_big_ports = ["広島港", "福山港", "尾道糸崎港"]
matched = [p for p in big_quay_ports if p in expected_big_ports]
h5 = len(matched) >= 2 and len(quays_large) > 0
hypos.append({
    "H": "H5",
    "claim": "200m 超岸壁は広島港・福山港・尾道糸崎港の 3 大商業港に集中",
    "result": (f"大型岸壁 {len(quays_large)} 件 / 該当港 {big_quay_ports} / "
               f"3 大商業港 ({len(matched)} 港マッチ)"),
    "verdict": "支持" if h5 else ("部分支持" if len(matched) >= 1 else "反証"),
})

# H6: 外郭の内側に係留 (80%+ が外郭 500m 圏)
near_pct_total = (n_near_h + n_near_f) / max(1, n_total_h + n_total_f) * 100
h6 = near_pct_total >= 80 if L32_OK else False
hypos.append({
    "H": "H6",
    "claim": f"係留 80%+ が L32 外郭から {PROXIMITY_M:.0f}m 圏内",
    "result": (f"全体 {near_pct_total:.1f}% (港湾 {n_near_h/max(1,n_total_h)*100:.1f}%, "
               f"漁港 {n_near_f/max(1,n_total_f)*100:.1f}%)" if L32_OK
               else "L32 データ不在のためスキップ"),
    "verdict": "支持" if h6 else ("部分支持" if near_pct_total >= 60 else ("反証" if L32_OK else "未検証")),
})

hypos_df = pd.DataFrame(hypos)
hypos_df.to_csv(ASSETS / "L33_hypothesis_results.csv", index=False, encoding="utf-8-sig")
with (ASSETS / "L33_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)


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

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

sections = []

# === Section 1: 学習目標と問い ===
sec1 = f"""
<div class="note">
<b>カバー宣言</b>: 本記事は DoBoX のシリーズ <b>「係留施設_港湾/漁港」 2 件</b>
(dataset_id = 1251, 1255) を統合し、広島県内における<b>港湾係留施設</b>
の<b>船舶受入構造</b>を分析する研究記事です。L32 (外郭施設) の続編として、
<b>「外郭で守られた港内のどこに、どんな係留が配置されるか」</b>を空間結合で読み解きます。
</div>

<h3>シリーズ構造判定: (a) 行政カテゴリ分割型 (L32 と同型)</h3>
<p>このシリーズの 2 件は、<b>「2 港分の同一データ」ではなく</b>、
<b>「港湾」(=県管理一般港湾、{n_ports_h} 港 {n_h_total} 件) と「漁港」(=県管理漁港、{n_ports_f} 漁港 {n_f_total} 件)</b>
という <b>2 つの行政カテゴリ</b>で全県の係留施設を二分する<b>相補型</b>構造。
合計 <b>{n_total_all} 件 / 41 港</b>、対象港集合は <b>L32 (外郭施設) と完全一致</b>。
両者は<b>完全に同一スキーマ</b> (14 列、GIS情報は LineString WKT) で、
本記事ではこれを縦結合して横断分析します。</p>

<h3>独自用語の定義 (本記事内で固定)</h3>
<table>
<tr><th>用語</th><th>定義</th></tr>
<tr><td>係留施設</td><td>船舶が<b>陸 (または陸近傍構造物) に接岸して係留</b>するための土木構造物の総称。
本記事では岸壁・物揚場・浮さん橋・船揚場・さん橋・係船くい・係船浮標 の 7 種を含む。
<b>外郭施設 (=防波堤等、L32 で扱った) や臨港交通施設 (=道路、L34) とは別カテゴリ</b>。</td></tr>
<tr><td>岸壁 (quay wall)</td><td>大型船 (貨物船・客船・フェリー) が直接接岸する高規格構造。
水深 4 m 以上の前面水域を持ち、垂直壁面を備える。本データに <b>{total_qw} 件のみ</b>
(港湾 {int(pv_kind.loc['港湾','岸壁'])} + 漁港 {int(pv_kind.loc['漁港','岸壁'])}) で<b>少数</b>。</td></tr>
<tr><td>物揚場 (mono-age-ba)</td><td>小型船・漁船が接岸して<b>水揚げ・荷役</b>する構造。
水深は浅く (1〜3 m)、岸壁より低規格。<b>本データの主役 ({total_mp} 件 / {mp_pct:.0f}%)</b>。
「もの (=漁獲・荷物) を揚げる場」が語源。</td></tr>
<tr><td>浮さん橋 (floating pier)</td><td><b>フロート (浮き) の上に乗る桟橋</b>。
潮位の上下動 (瀬戸内で 3 m 級) に追従するため、ガイド杭で水平位置を固定し垂直方向に動く。
本データに {fp_h+fp_f} 件 (港湾 {fp_h} + 漁港 {fp_f})。漢字表記の「浮さん橋」は古い土木用語。</td></tr>
<tr><td>船揚場</td><td>小型漁船・伝馬船を<b>陸に引き揚げる</b>傾斜面 (スリップウェイ)。
波打ち際の傾斜板で、エンジン故障時の引き揚げや干物乾燥に利用。
本データに {int(pv_kind.loc['合計','船揚場'])} 件。</td></tr>
<tr><td>さん橋 (pier)</td><td><b>固定杭の上に板床を架けた桟橋</b>。
浮さん橋と違い水位変化に追従せず、潮位差の少ない場所で旅客船・小型船バースとして使う。
本データに {int(pv_kind.loc['合計','さん橋'])} 件と少数。</td></tr>
<tr><td>係船くい</td><td>港内に立てた<b>1 本の杭 (or 数本の杭群)</b>に船が直接ロープを取って係留する。
桟橋を作るほどの予算が無いか、一時係留用。本データに {int(pv_kind.loc['合計','係船くい'])} 件。</td></tr>
<tr><td>係船浮標 (mooring buoy)</td><td>沖合に浮かべた<b>大型ブイ</b>に船をつなぐ。
錨地 (anchorage) 代わりに使う簡易係留。本データに {int(pv_kind.loc['合計','係船浮標'])} 件 (極少)。</td></tr>
<tr><td>係留線 総延長</td><td>1 港の全係留施設の幾何長 (LineString length) の合計。
港の<b>船舶受入容量</b>を表す指標。バース数の代理指標としても解釈可能。</td></tr>
<tr><td>外郭から {PROXIMITY_M:.0f} m 以内</td><td>本記事独自の<b>「守られた係留」</b>判定基準。
L32 で計算した外郭施設 (防波堤・突堤等) の任意点から <b>{PROXIMITY_M:.0f} m 以内</b>に位置する係留を
「外郭で守られた」とみなす。{PROXIMITY_M:.0f}m は港湾の典型的サイズ (港口幅 100-200m + 港内奥行き 200-300m) の概算上限。</td></tr>
</table>

<h3>研究の問い (RQ)</h3>
<p>広島県内に整備された係留施設 ({n_total_all} 件) は、
<b>「港湾 {n_ports_h} 港」</b>と<b>「漁港 {n_ports_f} 漁港」</b>の 2 行政カテゴリでどう分布し、
<b>どの構造形式</b>で<b>船舶接岸・荷役・漁船係留</b>を支えているか?
L32 で示した<b>外郭施設の防護網の内側</b>に、係留はどう配置されているか?</p>

<ol>
<li>港湾 vs 漁港 で<b>施設数・延長・構造形式の構成比</b>はどう違うか? L32 (外郭) との比較で何が見えるか?</li>
<li><b>主役構造</b>は岸壁か物揚場か? L32 では「防波堤」が主役 (67%) だったが、係留では?</li>
<li>大型岸壁 (>200m) は<b>どの港にどれだけ</b>集中しているか?</li>
<li><b>浮さん橋</b>は港湾と漁港のどちらに多いか? 設計選択の論理は?</li>
<li><b>外郭 (L32) と係留 (L33) の空間関係</b>: 係留は外郭から {PROXIMITY_M:.0f} m 圏内に何 % あるか?</li>
<li>港の<b>係留密度</b> (1 港あたり延長 / 港の物理サイズ) は港湾と漁港でどう違うか?</li>
</ol>

<h3>仮説 H1〜H6</h3>
<ul>
<li><b>H1 (カテゴリ規模差・反転)</b>: 港湾 {n_h_total} ({100-fish_pct:.0f}%) vs 漁港 {n_f_total} ({fish_pct:.0f}%)。
港湾優位が <b>L32 外郭 (1.33:1) より強く、約 1.9:1</b>。これは<b>商業港湾は 1 港あたり多数のバースが必要</b>な一方、
漁港は規模が小さく係留施設も少なくて済むため。延長比は <b>2 倍以上</b>。</li>
<li><b>H2 (構造形式の意外な主役)</b>: <b>岸壁は少数 ({qw_pct:.1f}%)</b>、
<b>主役は物揚場 ({mp_pct:.1f}%)</b>。「港湾係留 = 大型岸壁が主役」という素朴な予想に反し、
広島県の係留は<b>小型船受入の物揚場が圧倒的</b>。岸壁は大型船専用の高規格構造で相対的に少数。</li>
<li><b>H3 (浮さん橋の港湾偏在)</b>: 浮さん橋 {fp_h+fp_f} 件は<b>港湾 ({fp_h}) > 漁港 ({fp_f})</b>。
潮汐変動 (瀬戸内で 3 m 級) に対応する<b>フロート式桟橋</b>は、
旅客船・小型船バースとして港湾に多い。漁港は伝統的に物揚場で陸揚げ。</li>
<li><b>H4 (主要港集中)</b>: 港湾上位 5 で全体の <b>50%+</b>、漁港上位 3 で <b>50%+</b>。
L32 と同じ偏在構造。係留も外郭と同じ「ごく少数の大規模港 + 多数の小規模港」べき分布。</li>
<li><b>H5 (大型岸壁の都市港集中)</b>: 200 m 超の長大岸壁は<b>広島港・福山港・尾道糸崎港の 3 大商業港</b>に集中。
漁港の岸壁 12 件は全て 100 m 未満の中型岸壁。</li>
<li><b>H6 (外郭の内側に係留)</b>: 係留施設の<b>80% 以上</b>が L32 外郭施設から
<b>{PROXIMITY_M:.0f} m 以内</b>に位置する = 外郭で守られた港内に係留が配置されている。
裸の係留は <b>5% 未満</b>。これは<b>「外郭→係留」の 2 層整備</b>という階層設計の実証。</li>
</ul>

<h3>到達点</h3>
<p>2 dataset_id を「港湾 + 漁港」の相補的な 2 カテゴリとして読み解き、
{n_total_all} 件の係留施設の地理分布・構造形式構成・延長分布・大型岸壁集中・<b>外郭との 2 層関係</b>を
統合分析する。これにより、広島県の係留整備が
<b>「港湾は岸壁+物揚場+浮さん橋の混合」「漁港は物揚場主体の小規模分散」</b>という
<b>カテゴリで分化した二相設計</b>であり、かつ<b>外郭施設の防護網の内側に階層的に配置</b>されていることを実データで裏付ける。</p>

<div class="warn"><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 = [
    ("事業",       "港湾 / 漁港"),
    ("所管",       "港湾 / 漁港 (=事業と同じ)"),
    ("施設分類",   "係留施設 (定数)"),
    ("施設種類",   "岸壁 / 物揚場 / 浮さん橋 / 船揚場 / さん橋 / 係船くい / 係船浮標 (港湾 7 種、漁港 5 種)"),
    ("港湾名",     "港または漁港の名前 (41 種、L32 と同集合)"),
    ("事務所",     "管理事務所 (6 種): 三原・広島港湾・呉・東広島・東部・廿日市"),
    ("市区町村１/２",  "市町名 (NaN がほとんど)"),
    ("施設番号",   "施設番号 (C-1-01 等の体系。C は係留系 = mooring/quay 略)"),
    ("施設名称",   "施設の固有名 (例: 第1物揚場、北防波堤岸壁)"),
    ("管理者名等", "港湾管理者 / 漁港管理者 / 広島県 等"),
    ("GIS情報",    "WKT 形式の LineString または MultiLineString (経度・緯度)"),
    ("開始位置緯度/経度", "WKT の起点の緯度経度"),
]
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 の一覧。完全に同一スキーマ (14 列) で、
カテゴリ列 <code>事業 = 港湾 / 漁港</code> が二分の唯一の指標です。</p>

{sec2_table}

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

<h3>データ品質メモ</h3>
<ul>
<li><b>geom 欠損</b>: 港湾 {int(ALL[(ALL['port_category']=='港湾')&(ALL['GIS情報'].isna())].shape[0])} 件
+ 漁港 {int(ALL[(ALL['port_category']=='漁港')&(ALL['GIS情報'].isna())].shape[0])} 件
= <b>{int(ALL['GIS情報'].isna().sum())} 件 ({ALL['GIS情報'].isna().mean()*100:.1f}%)</b>。
これは古い施設・廃止施設・未測量施設の可能性。緯度経度列はあるが LineString が無いので延長計算には使えない。</li>
<li><b>施設番号体系</b>: <code>C-1-01</code> のような番号。<b>C</b> = 係留系 (Mooring の M ではなく、
港湾整備事業区分の Code C)。L32 外郭は <b>B</b> 系統。</li>
<li><b>geom タイプ</b>: 大半は <code>LINESTRING</code> (短い 2-4 点のもの多数)、
複合線が <code>MULTILINESTRING</code>。点 (POINT) は無く、線状で記録。</li>
<li><b>注意</b>: 「横田港」は港湾と漁港の両方に存在する。これは同名の異なる港 (港湾の横田港=広島港の一部の島しょ部、
漁港の横田=安芸郡蒲刈の漁港など) で、<code>(港湾名, port_category)</code> の 2 列で識別する。</li>
</ul>

<h3>本記事の主要分析テーブル</h3>
<p>2 dataset を縦結合した <code>L33_all_facilities.csv</code> ({n_total_geo} 行 × geom 有効分のみ) を主軸に、
各分析セクションでクロス集計と GIS 操作 (主に L32 との空間結合) を重ねる。</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("L33_series.csv")} — 2 dataset 一覧</li>
<li>{dl("L33_all_facilities.csv")} — 統合係留施設 ({n_total_geo} 行) + length_m + dist_to_outer_m + near_outer</li>
<li>{dl("L33_kind_count_cross.csv")} — カテゴリ × 構造形式 (件数)</li>
<li>{dl("L33_kind_length_cross.csv")} — カテゴリ × 構造形式 (延長 km)</li>
<li>{dl("L33_port_aggregate.csv")} — 港別集計 (件数 + 総延長 + 種別数)</li>
<li>{dl("L33_length_stats.csv")} — 構造形式別の延長統計</li>
<li>{dl("L33_large_quays_by_port.csv")} — 大型岸壁 (≥{LARGE_QUAY_THRESH:.0f}m) の港別集計</li>
<li>{dl("L33_port_density.csv")} — 港別 係留密度</li>
<li>{dl("L33_office_breakdown.csv")} — 事務所別管理件数</li>
<li>{dl("L33_hypothesis_results.csv")} — H1〜H6 検証結果</li>
</ul>

<h3>図 (本記事生成 PNG)</h3>
<ul>
<li>{dl("L33_fig1_dataset_overview.png")} / {dl("L33_fig2_kind_share.png")} / {dl("L33_fig3_port_ranking.png")}</li>
<li>{dl("L33_fig4_pref_overview.png")} / {dl("L33_fig5_kind_smallmult.png")} / {dl("L33_fig6_top_ports_detail.png")}</li>
<li>{dl("L33_fig7_length_distribution.png")} / {dl("L33_fig8_large_quays.png")} / {dl("L33_fig9_outer_mooring_layers.png")}</li>
<li>{dl("L33_fig10_office_breakdown.png")} / {dl("L33_fig11_function_compare.png")}</li>
</ul>

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


# === Section 4: 分析 1 — dataset 構造 ===
sec4 = f"""
<h3>狙い</h3>
<p>「港湾」と「漁港」の 2 dataset が、件数・カバー範囲・構造形式構成で
どう分化しているかを 1 枚の絵で示す。L32 (外郭) と比較し、係留固有の構造を識別する。</p>

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

<h3>実装</h3>
{code('''
SERIES = [
    (1251, "係留施設（港湾）", "港湾", "harbor_mooring_facility.csv",       32495),
    (1255, "係留施設（漁港）", "漁港", "fishing_port_mooring_facility.csv", 32499),
]
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)  # 全 1,224 行
pv_kind = pd.pivot_table(ALL, index="port_category", columns="施設種類",
                          values="施設番号", aggfunc="count", fill_value=0)
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 2 dataset の規模を文字 + バーで同時に伝える。
左 (カード) は「カテゴリ別の件数+港数+構造形式の内訳」をテキストで、
右 (バー) は「7 構造形式の絶対件数」を視覚化する。両方を 1 枚に置くことで
「カテゴリの規模差 (H1)」と「構造形式の意外な主役 (H2)」が同時に読める。</p>
{figure("assets/L33_fig1_dataset_overview.png", "2 dataset の構造概観 — カードビュー (左) と構造形式比較 (右)")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>港湾は {n_h_total} 件 / {n_ports_h} 港、漁港は {n_f_total} 件 / {n_ports_f} 漁港</b>。
件数比は <b>港湾:漁港 ≈ {length_ratio:.2f}:1</b>。
これは L32 外郭 (1.33:1) より港湾優位が強く、<b>商業港湾は 1 港あたり係留が多く必要</b>な反映。</li>
<li><b>主役は物揚場</b>: 港湾 {int(pv_kind.loc['港湾','物揚場'])} + 漁港 {int(pv_kind.loc['漁港','物揚場'])}
= {total_mp} 件で全体の <b>{mp_pct:.0f}%</b>。これは予想 (大型岸壁が主役) を覆す結果。
広島県の係留は<b>小型船・荷役の物揚場が圧倒的多数</b>。</li>
<li><b>岸壁は少数</b>: 全体で {total_qw} 件 ({qw_pct:.0f}%) のみ。
うち港湾 {int(pv_kind.loc['港湾','岸壁'])} 件、漁港 {int(pv_kind.loc['漁港','岸壁'])} 件。
これは「岸壁は大型船バース専用の高規格構造で、絶対数は少ない」という工学設計の現実を反映。</li>
<li><b>浮さん橋</b>: 港湾 {fp_h} ({fp_h_pct:.0f}%) > 漁港 {fp_f} ({fp_f_pct:.0f}%)。
H3 を支持。フロート構造は<b>潮汐対応が必要な客船・小型船バース</b>として港湾に多い。</li>
<li>港湾 7 種の構造形式 (係船くい・さん橋・係船浮標を含む) vs 漁港 5 種。
港湾の方が<b>構造形式の多様性</b>が高い = 多様な船種を受け入れる必要がある。</li>
</ul>

<h3>表と読み取り</h3>
{pv_kind.to_html(classes='', border=0)}
<p><b>読み取り</b>: 構造形式 7 種のうち、<b>港湾と漁港の両方に多数あるのは「物揚場」「浮さん橋」「岸壁」「船揚場」</b>。
「係船くい」「さん橋」「係船浮標」は港湾固有 (漁港にはほぼ無い)。
これは漁港側の構造選択が<b>「物揚場主体 + 補助的な浮さん橋・船揚場」</b>に絞られていることを示す。</p>
"""
sections.append(("分析 1: 2 dataset の構造を可視化", sec4))


# === Section 5: 分析 2 — 構造形式の構成比と延長 ===
sec5 = f"""
<h3>狙い</h3>
<p>構造形式の比率を<b>件数</b>と<b>延長 (km)</b>の 2 ベースで比較する。
<b>「件数では多いが延長では少ない」</b>のような構造形式 (例: 短い物揚場が多数) を識別する。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p><b>WKT (Well-Known Text)</b> とは、点・線・ポリゴンなどの幾何形状を
人間が読める文字列で表現する標準形式です。例:</p>
<pre><code>LINESTRING (132.75588 34.18490, 132.75582 34.18484)</code></pre>
<p>これを Python の <code>shapely.wkt.loads()</code> で <b>幾何オブジェクト</b>に変換し、
<code>geometry.length</code> でメートル単位の延長を計算します。
ただし<b>緯度経度のままだと度単位</b>になってしまうので、まず
<b>EPSG:6671 (平面直角座標系 第 III 系、広島県を含む)</b>に再投影します。</p>

<ul>
<li><b>入力</b>: WKT 文字列 (LINESTRING または MULTILINESTRING)</li>
<li><b>出力</b>: メートル単位の延長 (1 施設につき 1 値)</li>
<li><b>限界</b>: 線が複雑に折れていれば実際のバース長より長くなる場合あり (本データの線は単純なので問題なし)</li>
</ul>

<h3>実装</h3>
{code('''
from shapely.wkt import loads as wkt_loads
ALL["geometry"] = ALL["GIS情報"].apply(lambda s: wkt_loads(str(s)) if pd.notna(s) else None)
gdf = gpd.GeoDataFrame(ALL.dropna(subset=["geometry"]),
                       geometry="geometry", crs="EPSG:4326").to_crs("EPSG:6671")
gdf["length_m"] = gdf.geometry.length

# 件数 vs 延長 の 2 ベース
pv_kind = pd.pivot_table(ALL, index="port_category", columns="施設種類",
                          values="施設番号", aggfunc="count", fill_value=0)
pv_len = pd.pivot_table(gdf, index="port_category", columns="施設種類",
                         values="length_m", aggfunc="sum", fill_value=0) / 1000.0
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 件数構成比だけでは「小さな施設が多数」のような
構造を見落とす。延長ベースを併置することで「実際のバース総長としての貢献」が読める。</p>
{figure("assets/L33_fig2_kind_share.png", "構造形式の構成比 — 件数ベース (左) と延長ベース (右)")}
<p><b>読み取り</b>:</p>
<ul>
<li>港湾の件数では物揚場が <b>{ratios_h[1]:.0f}%</b>、岸壁 <b>{ratios_h[0]:.0f}%</b> と<b>物揚場優位</b>だが、
延長で見ると物揚場 <b>{ratios_h2[1]:.0f}%</b>、岸壁 <b>{ratios_h2[0]:.0f}%</b> と<b>岸壁の延長シェアが急増</b>。
これは岸壁 1 本が長い ({qw_h_med:.0f} m 級) ためで、
<b>「件数では物揚場が、延長では岸壁が主役」</b>という二相が読める。</li>
<li>漁港では物揚場が件数 <b>{ratios_f[1]:.0f}%</b> → 延長 <b>{ratios_f2[1]:.0f}%</b> でほぼ一致。
これは漁港の物揚場サイズが均質 (中央値 {mp_f_med:.0f}m) であることを示す。</li>
<li>港湾の<b>浮さん橋</b>は件数 <b>{ratios_h[2]:.0f}%</b> → 延長 <b>{ratios_h2[2]:.0f}%</b>。
1 本が短い (10-20m 級) ため延長シェアは低下する。</li>
<li>港湾と漁港で<b>合計延長</b>: 港湾 {total_h:.1f} km vs 漁港 {total_f:.1f} km, 比 <b>{length_ratio:.2f} 倍</b>。
H1 の「延長 2 倍以上」を <b>{'支持' if length_ratio >= 2.0 else '部分支持'}</b>。</li>
</ul>

<h3>表と読み取り</h3>
{pv_len.to_html(classes='', border=0)}
<p><b>読み取り</b>:</p>
<ul>
<li>港湾の岸壁総延長は <b>{pv_len.loc['港湾','岸壁']:.2f} km</b>。これは広島港・福山港・尾道糸崎港の
大型岸壁を中心に、商業貨物受入の主要バース総量を表す。</li>
<li>港湾の物揚場総延長は <b>{pv_len.loc['港湾','物揚場']:.2f} km</b>で岸壁を {pv_len.loc['港湾','物揚場']/max(0.01,pv_len.loc['港湾','岸壁']):.1f} 倍上回る。
件数で {int(pv_kind.loc['港湾','物揚場']/max(1,pv_kind.loc['港湾','岸壁']))} 倍の差 (425 vs 68) を考慮すると、
<b>1 本あたりの長さは岸壁が物揚場の {qw_h_med/max(0.01,mp_h_med):.1f} 倍</b>。</li>
<li>漁港の合計延長 ({total_f:.1f} km) は港湾 ({total_h:.1f} km) の {total_f/total_h*100:.0f}% で、
漁港カテゴリの「分散小型多数」性を再確認できる。</li>
</ul>
"""
sections.append(("分析 2: 構造形式の構成比と延長", sec5))


# === Section 6: 分析 3 — 港別ランキング ===
top10_html = port_agg.sort_values("total_length_km", ascending=False).head(10).to_html(index=False, classes='', border=0)
sec6 = f"""
<h3>狙い</h3>
<p>「どの港が係留容量をどれだけ持つか」を件数・延長の 2 軸で順位付け。
H4 (上位港集中) を検証し、広島県の港湾整備の地理偏在を視覚化する。</p>

<h3>手法</h3>
<p>港湾名 + カテゴリでグループ集計。施設数・総延長・構造形式種類数を取り、
延長で降順ソート。上位 15 港を 2 軸バーで比較。</p>

<h3>実装</h3>
{code('''
port_agg = gdf.groupby(["port_category", "港湾名"]).agg(
    n_facilities=("length_m", "size"),
    total_length_m=("length_m", "sum"),
    n_kinds=("施設種類", "nunique"),
).reset_index()
port_agg["total_length_km"] = port_agg["total_length_m"] / 1000
port_agg = port_agg.sort_values("total_length_km", ascending=False)
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 件数だけだと「小さい施設多数の港」が上位、延長だけだと
「大型施設少数の港」が上位になる。両方を並べることで両ベクトルの整合・不整合が読める。</p>
{figure("assets/L33_fig3_port_ranking.png", "係留施設 上位 15 港 — 件数 (左) と総延長 (右)")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>{top_all.iloc[0]['港湾名']} ({top_all.iloc[0]['port_category']})</b> が両ベクトルで首位 —
{int(top_all.iloc[0]['n_facilities'])} 件 / {top_all.iloc[0]['total_length_km']:.2f} km。
県内最大の港の規模を反映。</li>
<li><b>{top_all.iloc[1]['港湾名']} ({top_all.iloc[1]['port_category']})</b> が 2 位 —
{int(top_all.iloc[1]['n_facilities'])} 件 / {top_all.iloc[1]['total_length_km']:.2f} km。</li>
<li><b>{top_all.iloc[2]['港湾名']} ({top_all.iloc[2]['port_category']})</b> が 3 位 —
{int(top_all.iloc[2]['n_facilities'])} 件 / {top_all.iloc[2]['total_length_km']:.2f} km。</li>
<li>カテゴリ混在ランキング: 上位 15 港のうち<b>港湾 ≈ 漁港</b> が混在。
H4 の「上位集中」を検証するためには、カテゴリ別に分離して集中率を見る (下表)。</li>
<li><b>L32 ランキングとの比較</b>: 外郭ランキング上位は広島港・尾道糸崎港・倉橋・豊島だったが、
係留ランキングは {", ".join(top_all.head(4)['港湾名'].tolist())} で順序が一部変わる。
これは「外郭整備の規模 ≠ 係留整備の規模」を意味し、港の機能 (商業 vs 漁業) や開発時期で
両者の比率が変わることを示唆する (発展課題 Z2)。</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>上位 5 ({', '.join(top5_h['港湾名'].tolist())})</td>
<td>{top5_h_sum} / {int(port_agg[port_agg['port_category']=='港湾']['n_facilities'].sum())}</td>
<td>{top5_h_pct:.1f}%</td></tr>
<tr><td>漁港</td><td>上位 3 ({', '.join(top3_f['港湾名'].tolist())})</td>
<td>{top3_f_sum} / {int(port_agg[port_agg['port_category']=='漁港']['n_facilities'].sum())}</td>
<td>{top3_f_pct:.1f}%</td></tr>
</table>
<p><b>読み取り</b>: 港湾上位 5 で <b>{top5_h_pct:.0f}%</b>、漁港上位 3 で <b>{top3_f_pct:.0f}%</b> を占有。
H4 を {hypos[3]['verdict']}。整備は<b>「ごく少数の大規模港 + 多数の小規模港」</b>のべき分布で、
これは商業性 (港湾) も漁業集積度 (漁港) も同じ偏在を示す。L32 (外郭) でも同じパターンだった。</p>
"""
sections.append(("分析 3: 港別ランキングと上位集中", sec6))


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

<h3>手法</h3>
<p>WKT パース済 GeoDataFrame を EPSG:6671 (平面直角 III 系) に投影し、
<code>geopandas.plot()</code> で<b>線描画</b>。
県全域 polygon (L15 の admin_922) を背景に、カテゴリごとに色分け。
さらに 7 構造形式を別パネル (small multiples) で並べ、それぞれの偏在パターンを比較する。</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))
g_admin_pref.boundary.plot(ax=ax, color="#aaa", linewidth=0.4)
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, linewidth=2.5)
''')}

<h3>図と読み取り — 県全域</h3>
<p><b>なぜこの図か</b>: 県全域に対する係留線の分布密度を、青 (港湾) と緑 (漁港) で
直接見ることで「どこに何カテゴリが集中しているか」を一目で把握できる。
個別の港名ラベルで上位 10 港を識別。</p>
{figure("assets/L33_fig4_pref_overview.png", f"広島県 係留施設マップ — 全 {len(gdf)} 件をカテゴリで色分け、上位 10 港にラベル")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>広島湾岸 (広島市〜廿日市)</b> に港湾 (青) が密集。広島港・厳島港の大型岸壁。</li>
<li><b>福山〜尾道〜三原の東部</b>: 港湾 (青) が連続。尾道糸崎港・福山港の連鎖。</li>
<li><b>呉湾岸 + 倉橋島・江田島</b>: 漁港 (緑) が島嶼に密集。倉橋・豊島・沖浦の漁港群。</li>
<li><b>瀬戸内海島嶼の東 (大崎上島・生口島・因島)</b>: 港湾 (生口港・木江港) と漁港が混在。</li>
<li>L32 (外郭) のマップと<b>分布パターンがほぼ一致</b>。同じ港集合の同じ地理だから当然だが、
これは<b>「外郭が築かれた港には係留も整備される」</b>という強い相関の視覚化。</li>
</ul>

<h3>図と読み取り — 構造形式 small multiples</h3>
<p><b>なぜこの図か</b>: 県全域マップでは構造形式の偏在が見えない。
7 形式を別パネルにすることで「物揚場は全域、岸壁は商業港のみ」のような
形式ごとの地理パターンが分離して見える。</p>
{figure("assets/L33_fig5_kind_smallmult.png", "構造形式 7 種別 small multiples — それぞれの偏在パターン")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>岸壁 ({total_qw} 件)</b>: 広島港・福山港・尾道糸崎港・呉地区など<b>商業港に集中</b>。漁港エリアにはほぼ無い。</li>
<li><b>物揚場 ({total_mp} 件)</b>: 県全域に均等分布。<b>最も汎用的な係留構造</b>。
これが「主役」と呼ばれる所以。</li>
<li><b>浮さん橋 ({fp_h+fp_f} 件)</b>: 広島湾と尾道糸崎周辺に偏在。
旅客船発着場・海上タクシー乗り場に多い。</li>
<li><b>船揚場</b>: 漁港エリアと小規模港湾に分散。漁船陸揚げの傾斜面。</li>
<li><b>係船くい・さん橋・係船浮標</b>: 各 1〜10 件と稀少構造。
広島港と一部の港湾に偶然見られる程度。</li>
</ul>
"""
sections.append(("分析 4: 県全域マップで係留分布を観る", sec7))


# === Section 8: 分析 5 — 上位 4 港の詳細 ===
sec8 = f"""
<h3>狙い</h3>
<p>上位 4 港 ({", ".join([p for p, c in detail_ports])}) の詳細マップを並べ、
<b>各港の係留配置パターン</b>を 1 港ずつ精読する。
「少数事例を深く掘る」アプローチで、港の機能・歴史・地形が係留設計にどう反映されているかを読む。</p>

<h3>手法</h3>
<p>各港について bbox を取り、施設のみを zoom 表示。
構造形式ごとに別色で重ねることで、岸壁・物揚場・浮さん橋の位置関係が読める。</p>

<h3>実装</h3>
{code('''
top_h_ports = port_agg[port_agg["port_category"]=="港湾"].head(2)["港湾名"].tolist()
top_f_ports = port_agg[port_agg["port_category"]=="漁港"].head(2)["港湾名"].tolist()
detail_ports = [(p, "港湾") for p in top_h_ports] + [(p, "漁港") for p in top_f_ports]

fig, axes = plt.subplots(2, 2, figsize=(13, 11))
for ax, (port, cat) in zip(axes.flatten(), detail_ports):
    sub = gdf[(gdf["港湾名"] == port) & (gdf["port_category"] == cat)]
    bbox = sub.total_bounds
    ax.set_xlim(bbox[0]-pad, bbox[2]+pad)
    ax.set_ylim(bbox[1]-pad, bbox[3]+pad)
    for k in KIND_ORDER:
        sk = sub[sub["施設種類"] == k]
        sk.plot(ax=ax, color=KIND_COLOR[k], linewidth=2.4, label=k)
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 県全域マップでは個々の港の構造が潰れる。
ズーム + 構造形式色分けで「どの港がどんな係留パターンか」を 1 港単位で読める。</p>
{figure("assets/L33_fig6_top_ports_detail.png", "上位 4 港の詳細マップ — 構造形式色分け")}
<p><b>読み取り — 1 港ずつ</b>:</p>
"""
# 動的に各港の解説を生成
detail_text = "<ul>\n"
for port, cat in detail_ports:
    sub_p = gdf[(gdf["港湾名"]==port) & (gdf["port_category"]==cat)]
    n_p = len(sub_p)
    len_p = sub_p["length_m"].sum() / 1000.0
    kinds_p = sub_p["施設種類"].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_p} 件 / {len_p:.2f} km)</b>: {kinds_str}。"
    )
    # 港固有の特徴
    if port == "広島港":
        detail_text += "県内最大商業港で<b>岸壁・物揚場・浮さん橋・係船くい・さん橋・船揚場</b>の<b>多形式共存</b>。フェリー・客船・貨物が共存する複合港湾の典型。</li>\n"
    elif port == "福山港":
        detail_text += "県東部最大の商業港。<b>長大岸壁が集中</b>し鉄鋼・機械の物流バースを形成。物揚場は補助的。</li>\n"
    elif port == "尾道糸崎港":
        detail_text += "<b>細片化が著しい</b>多島港で、各島面に物揚場と短い岸壁が分散配置。</li>\n"
    elif port == "倉橋":
        detail_text += "瀬戸内屈指の漁港集積地。<b>物揚場が圧倒的多数</b>で漁船係留に特化。岸壁は中型のみ。</li>\n"
    elif port == "豊島":
        detail_text += "島嶼漁港の典型。物揚場 + 浮さん橋で島嶼集落の係留需要に対応。</li>\n"
    elif port == "沖浦":
        detail_text += "倉橋島南部の漁港。物揚場主体で外郭施設の内側に集約配置。</li>\n"
    else:
        detail_text += "中規模港の典型構成。</li>\n"
detail_text += "</ul>\n"

sec8 += detail_text
sec8 += f"""
<div class="note"><b>係留パターンの 3 類型 (発見)</b>: 上位 4 港の比較から、係留施設の組み合わせは
少なくとも以下 3 つのパターンに分類できる:
<ol>
<li><b>多形式共存型 (広島港型)</b>: 大型商業港。岸壁 + 物揚場 + 浮さん橋 + 係船くい などが共存。
複数船種 (フェリー・貨物・客船・小型船) を 1 港で受け入れる。</li>
<li><b>長大岸壁型 (福山港型)</b>: 工業港湾。長大岸壁が中心、補助的に物揚場。
鉄鋼・機械・コンテナの大型船バースに特化。</li>
<li><b>物揚場集中型 (倉橋・豊島型)</b>: 漁港。物揚場が圧倒的、岸壁は中型 1-2 本のみ。
小型漁船の水揚げ・荷役に特化。</li>
</ol>
これらは港の<b>機能</b> (商業 / 工業 / 漁業) と<b>規模</b> (大規模商業港 / 集落漁港) の積で決まる。
<b>L32 で示した「港湾は防波堤主体・漁港は防波堤+護岸」という外郭の二相設計</b>と整合する係留設計。</div>
"""
sections.append(("分析 5: 上位 4 港の詳細マップ — 係留パターンの 3 類型", sec8))


# === Section 9: 分析 6 — 延長分布 (Lorenz 等) ===
sec9 = f"""
<h3>狙い</h3>
<p>係留施設の<b>延長分布</b>を多角的に観察し、H2 (物揚場主役)・H4 (上位港集中)・H5 (大型岸壁集中) を統計的に検証する。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p>4 つの可視化手法を統合する:</p>
<ul>
<li><b>(1) ヒストグラム</b>: 岸壁延長の度数分布を港湾・漁港で重ねる。中央値線で 2 群を比較。
200 m の閾値線で大型岸壁を識別。</li>
<li><b>(2) Box plot (log 軸)</b>: 7 構造形式の延長分布を 1 枚に集約。log スケールにすることで
小型 (5 m) 〜 大型 (500 m+) の数桁レンジを同時可視化。</li>
<li><b>(3) Lorenz 曲線</b>: 累積分布。x 軸 = 港のランク (%)、y 軸 = 累積延長 (%)。
完全均等線 (対角線) からの<b>偏差</b>が集中度。経済学のジニ係数と同じ概念。</li>
<li><b>(4) 散布図 (バブル)</b>: 各港について「物揚場本数 × 平均延長」を散布。
バブルサイズで総延長を表現。1 港 = 1 バブルの構造プロファイル。</li>
</ul>

<h3>実装</h3>
{code('''
# 物揚場プロファイル
mp_per_port = gdf[gdf["施設種類"]=="物揚場"].groupby(["port_category", "港湾名"]).agg(
    n=("length_m", "size"),
    mean_m=("length_m", "mean"),
    total_m=("length_m", "sum"),
).reset_index()

# Lorenz 曲線
all_sorted = port_agg.sort_values("total_length_km", ascending=False)
all_sorted["cum_pct"] = all_sorted["total_length_km"].cumsum() / all_sorted["total_length_km"].sum() * 100
all_sorted["rank_pct"] = (np.arange(len(all_sorted)) + 1) / len(all_sorted) * 100
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 4 視点を 2x2 で並べることで、分布形状・分散・集中・港単位プロファイルを
1 枚に統合できる。これは多角的データを 1 セクションで深く掘るための定石レイアウト。</p>
{figure("assets/L33_fig7_length_distribution.png", "延長分布 4 視点 — 岸壁ヒスト、boxplot、Lorenz、物揚場散布図")}
<p><b>読み取り (左上 — 岸壁ヒストグラム)</b>:</p>
<ul>
<li>港湾岸壁の中央値 <b>{qw_h_med:.0f} m</b> vs 漁港 <b>{qw_f_med:.0f} m</b>。
港湾が漁港の {qw_h_med/max(0.1,qw_f_med):.1f} 倍。</li>
<li>港湾岸壁は<b>右裾が長い</b>: 200 m 超の長大岸壁が <b>{len(quays_large)} 本</b>。
これは外洋に面する大型商業港 (広島・福山・尾道糸崎) の岸壁。</li>
<li>漁港岸壁 ({len(qw_f)} 件) は<b>100 m 未満に集中</b>。漁港の岸壁は中型以下に限定。</li>
</ul>
<p><b>読み取り (右上 — 構造形式 Boxplot)</b>:</p>
<ul>
<li>岸壁の中央値が最大。1 本あたりの長さでは<b>岸壁が首位</b> (件数 vs 延長の主役逆転を示唆)。</li>
<li>物揚場は中央値 50-60m に集中、ばらつき小。<b>標準化された設計</b>。</li>
<li>浮さん橋は中央値 15-25m と短く、フロート構造の限界を反映 (大型化が困難)。</li>
<li>係船くい・係船浮標はサンプル数が少なく分布が不安定。</li>
</ul>
<p><b>読み取り (左下 — Lorenz 曲線)</b>:</p>
<ul>
<li>曲線は対角線から<b>大きく上に膨らむ</b> = 上位集中の証拠。</li>
<li>上位 10% の港 (= 上位 4 港) で係留線延長の <b>~{share_at_10:.0f}%</b> を保有 (具体値は図参照)。</li>
<li>上位 25% (= 上位 10 港) で <b>~{share_at_25:.0f}%</b>。</li>
<li>これは「ごく少数の主要港 + 多数の小規模港」のべき分布構造。
L32 (外郭) と<b>同じ偏在パターン</b>で、外郭整備と係留整備が連動していることを示唆。</li>
</ul>
<p><b>読み取り (右下 — 物揚場プロファイル散布図)</b>:</p>
<ul>
<li>右上 (本数多 + 平均長 大) に位置する港は<b>大型物揚場ハブ</b>: 倉橋・広島港など。</li>
<li>左下 (本数少 + 平均長 小) は<b>小規模港</b>: 多数の漁港。</li>
<li>右下 (本数多 + 平均長 小) は<b>細片型</b>: 倉橋 ({mp_per_port[(mp_per_port['port_category']=='漁港')&(mp_per_port['港湾名']=='倉橋')]['n'].iloc[0] if len(mp_per_port[(mp_per_port['port_category']=='漁港')&(mp_per_port['港湾名']=='倉橋')]) else 0:.0f} 本だが平均 {mp_per_port[(mp_per_port['port_category']=='漁港')&(mp_per_port['港湾名']=='倉橋')]['mean_m'].iloc[0] if len(mp_per_port[(mp_per_port['port_category']=='漁港')&(mp_per_port['港湾名']=='倉橋')]) else 0:.0f}m)。</li>
<li>左上 (本数少 + 平均長 大) は<b>大型 1 本型</b>: 工業港湾の少数大型物揚場。</li>
</ul>
"""
sections.append(("分析 6: 延長分布の 4 視点解析", sec9))


# === Section 10: 分析 7 — 大型岸壁マップ ===
big_quay_html = quays_large_pv.to_html(index=False, classes='', border=0)
sec10 = f"""
<h3>狙い</h3>
<p>大型岸壁 (≥ {LARGE_QUAY_THRESH:.0f} m) は商業貨物受入の<b>戦略的バース</b>。
どの港にどれだけ集中しているかを地図と表で示し、H5 を検証する。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p><b>条件抽出 + 地図プロット</b>のシンプルな手順だが、
<b>「200 m」という閾値の根拠</b>が研究的に重要:</p>
<ul>
<li><b>200 m 岸壁の意味</b>: コンテナ船・大型 RO-RO 船 (船長 150-200m) が 1 隻接岸できる長さ。
これより短いと中小型船専用、これより長いと複数船同時接岸が可能。</li>
<li><b>限界</b>: 200 m はあくまで「中型を超える」境界。<b>真の大型岸壁 (国際海上コンテナ用)</b>
は 300-400m 以上が必要だが、本データでは少数 (=広島県は中規模商業港メイン)。</li>
</ul>

<h3>実装</h3>
{code('''
LARGE_QUAY_THRESH = 200.0  # m
quays = gdf[gdf["施設種類"] == "岸壁"].copy()
quays_large = quays[quays["length_m"] >= LARGE_QUAY_THRESH].copy()
quays_large_pv = quays_large.groupby(["port_category", "港湾名"]).agg(
    n=("length_m", "size"),
    total_m=("length_m", "sum"),
    max_m=("length_m", "max"),
).reset_index().sort_values("total_m", ascending=False)
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 大型岸壁は数が少ない ({len(quays_large)} 本) ので、
<b>1 本ずつ位置と長さをラベル付き</b>で示す。中型岸壁 (100-200m) を黄色、
大型岸壁 (≥200m) を赤、その他全係留を背景灰で表示することで階層が見える。</p>
{figure("assets/L33_fig8_large_quays.png", f"大型岸壁 (≥{LARGE_QUAY_THRESH:.0f} m) の地理分布")}
<p><b>読み取り</b>:</p>
<ul>
<li>大型岸壁は <b>{len(quays_large)} 本</b>のみで、岸壁全 {len(quays)} 本の <b>{len(quays_large)/max(1,len(quays))*100:.1f}%</b>。</li>
<li>港別では<b>{', '.join(quays_large_pv['港湾名'].head(3).tolist())}</b> に集中。
H5 ({", ".join(expected_big_ports)} に集中) を <b>{hypos[4]['verdict']}</b>。</li>
<li>最長岸壁は <b>{quays_large['length_m'].max():.0f} m</b> (港: {quays_large.sort_values('length_m', ascending=False).iloc[0]['港湾名']})。
これは 1 隻の大型船が接岸し荷役する戦略バース。</li>
<li>漁港の岸壁 (12 件) はすべて 200 m 未満で、漁港カテゴリには大型岸壁ゼロ。</li>
<li>大型岸壁は<b>島嶼ではなく本土沿岸</b>に集中。これは「大型船が接岸できる水深と航路が必要」
という船舶工学的制約の表れ。</li>
</ul>

<h3>表と読み取り — 大型岸壁の港別内訳</h3>
{big_quay_html}
<p><b>読み取り</b>: 上位 3 港 ({', '.join(quays_large_pv['港湾名'].head(3).tolist())}) で
大型岸壁の {(quays_large_pv['n'].head(3).sum()/max(1,len(quays_large))*100):.0f}% を占有。
これらは広島県の<b>商業海運の戦略的拠点</b>であり、
県内貨物量・船舶往来数・経済規模との対応関係が想定される (発展課題 Z3)。</p>
"""
sections.append((f"分析 7: 大型岸壁 (≥{LARGE_QUAY_THRESH:.0f}m) の集中分析", sec10))


# === Section 11: 分析 8 — 外郭と係留の 2 層関係 ===
if L32_OK:
    near_pct_h = n_near_h / max(1, n_total_h) * 100
    near_pct_f = n_near_f / max(1, n_total_f) * 100
    near_pct_total = (n_near_h + n_near_f) / max(1, n_total_h + n_total_f) * 100
    median_dist = gdf["dist_to_outer_m"].median()
    sec11 = f"""
<h3>狙い</h3>
<p><b>L32 で示した外郭施設 (防波堤・突堤等) と本記事の係留施設の空間関係</b>を分析する。
仮説: 係留は外郭で守られた港内に配置されるはず。「裸の係留 (外郭から遠い)」はあるか?</p>

<h3>手法 (リテラシレベル解説)</h3>
<p><b>距離計算 + 閾値判定</b>:</p>
<ol>
<li><b>外郭施設を 1 つの巨大な MultiLineString に統合</b> (<code>union_all()</code>)。
これにより距離計算が高速化される (1 vs 1 の距離計算で済む)。</li>
<li><b>各係留施設の中心点から外郭への最近距離</b>を <code>geometry.distance()</code> で計算。</li>
<li><b>{PROXIMITY_M:.0f} m 以内</b>を「守られた係留」と判定。</li>
</ol>
<p><b>{PROXIMITY_M:.0f} m の根拠</b>: 港湾の典型サイズ (港口幅 100-200m + 港内奥行き 200-300m)
の概算上限。これより遠いと「同一港の外郭で守られている」とは言えない。</p>

<ul>
<li><b>入力</b>: 係留 LineString {n_total_geo} 本 + 外郭 LineString (L32 由来) </li>
<li><b>出力</b>: 各係留の <code>dist_to_outer_m</code> + <code>near_outer</code> ブール</li>
<li><b>限界</b>: 「近い」 ≠ 「守られている」。実際の防護は<b>外郭の高さ・配向・港の地形</b>に依存。
本分析は<b>位置の必要条件</b>を見ている。</li>
</ul>

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

# 外郭施設を 1 つの巨大 MultiLineString に統合 (高速化)
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"] <= 500.0  # 500m 圏
''')}

<h3>図と読み取り</h3>
<p><b>なぜこの図か</b>: 県全域 (左) と広島港ズーム (右) の 2 視点。
県全域では外郭(紫)+守られた係留(青)+裸の係留(赤)の 3 層が一目で見える。
広島港ズームで個別の港内構造を精査できる。</p>
{figure("assets/L33_fig9_outer_mooring_layers.png", "外郭 (L32) と係留 (L33) の 2 層関係マップ")}
<p><b>読み取り</b>:</p>
<ul>
<li><b>係留全体の {near_pct_total:.1f}%</b> が外郭から {PROXIMITY_M:.0f}m 圏内 (港湾 {near_pct_h:.1f}%, 漁港 {near_pct_f:.1f}%)。
H6 (80%+) を <b>{hypos[5]['verdict']}</b>。</li>
<li>係留全体の<b>外郭への最近距離 中央値: {median_dist:.0f} m</b>。
ほとんどの係留は外郭のすぐ近く (= 港内) に位置する。</li>
<li><b>裸の係留 ({(1 - near_pct_total/100)*100:.1f}%)</b>: これは外郭施設が無い小規模港、
または geom 不一致 (外郭の geom 欠損港の係留)、もしくは離れた離島の係留。
県全域マップで赤線として識別される。</li>
<li>広島港ズームでは、紫 (外郭=北防波堤・宇品防波堤等) の<b>内側</b>に、
青や色付の係留 (岸壁・物揚場・浮さん橋) が並ぶ<b>明確な階層構造</b>が見える。
これは港湾工学の教科書的設計の実データ確認。</li>
<li><b>L32 + L33 の統合的読み</b>: 41 港の港湾施設は、外側に外郭 (波浪遮断)、
内側に係留 (船舶接岸) の<b>2 層階層</b>で整備されている。
これは港湾整備の<b>「先に外郭、後に係留」</b>という時系列順序とも整合 (発展課題 Z2)。</li>
</ul>

<div class="warn"><b>解釈の注意</b>: 「{PROXIMITY_M:.0f}m 圏内」 = 「物理的に近い」だけで、
実際の防護機能は<b>外郭の高さ・港の地形・想定波高</b>に依存する。
本分析は階層配置の<b>位置検証</b>であり、機能評価ではない。
本格的な防護評価は<b>外郭高さ × 港の方位 × 想定波</b>の 3 軸で別途行う必要がある (発展課題 Z3)。</div>
"""
else:
    sec11 = f"""
<h3>狙い</h3>
<p>L32 (外郭施設) と本記事の係留施設の空間関係を分析する予定でしたが、
L32 中間 CSV が <code>lessons/assets/L32_all_facilities.csv</code> に見つからないためスキップ。</p>
<p>L32 を先に実行してから本記事を再実行することで、外郭との階層関係分析が有効化されます。</p>
"""
sections.append(("分析 8: 外郭 (L32) と係留 (L33) の 2 層関係", sec11))


# === Section 12: 仮説検証と考察 ===
sec12_table = hypos_df.to_html(index=False, classes='', border=0)

# 機能別構成比を表化
func_table_rows = []
for func, kinds in FUNC_GROUPS.items():
    total_func = int((ALL["施設種類"].isin(kinds)).sum())
    pct = total_func / max(1, n_total_all) * 100
    h_func = int(((ALL["port_category"]=="港湾") & (ALL["施設種類"].isin(kinds))).sum())
    f_func = int(((ALL["port_category"]=="漁港") & (ALL["施設種類"].isin(kinds))).sum())
    func_table_rows.append(
        f"<tr><td>{func}</td><td>{', '.join(kinds)}</td>"
        f"<td>{total_func} ({pct:.1f}%)</td>"
        f"<td>{h_func}</td><td>{f_func}</td></tr>"
    )
func_html = (
    "<table><tr><th>機能</th><th>含む構造</th><th>合計</th><th>港湾</th><th>漁港</th></tr>"
    + "".join(func_table_rows) + "</table>"
)

sec12 = f"""
<p>H1〜H6 の検証結果を 1 表で示す。</p>
{sec12_table}

<h3>機能グルーピング (4 機能で集約)</h3>
<p>7 構造形式を物理機能で 4 グループに集約すると:</p>
{func_html}
{figure("assets/L33_fig11_function_compare.png", "機能別マップ (左) + 港湾 vs 漁港 機能構成比 (右)")}
<p><b>読み取り</b>:</p>
<ul>
<li>港湾は<b>大型船バース (岸壁) + 小型船・荷役 (物揚場) + 桟橋系 (浮さん橋) の 3 機能</b>が共存。
4 機能目の「杭・浮標」は補助的に少数。</li>
<li>漁港は<b>小型船・荷役 (物揚場+船揚場) が圧倒的</b>。桟橋系は浮さん橋のみで補助。
大型船バース (岸壁) は中型 12 件のみ。</li>
<li>これは港湾と漁港の<b>機能分化</b>の明確な視覚化。</li>
</ul>

<h3>総括: 広島県の係留整備思想</h3>
<p>2 dataset から再構成した係留施設の構造分析により、以下の<b>3 つの設計思想</b>が読み取れる。</p>
<ul>
<li><b>(1) 二相設計</b>: 港湾と漁港でほぼ完全に分化した構造形式採用。
港湾は<b>多形式共存 (岸壁+物揚場+浮さん橋)</b>、
漁港は<b>物揚場主体 (件数 {int(pv_kind.loc['漁港','物揚場'])}/{n_f_total} = {int(pv_kind.loc['漁港','物揚場'])/max(1,n_f_total)*100:.0f}%)</b>。
これは<b>商業 vs 漁業</b>という機能の違いと、<b>大規模港 vs 集落漁港</b>という規模の違いの積。
L32 の「港湾は防波堤主体・漁港は防波堤+護岸の二相」と整合する。</li>
<li><b>(2) 主要港集中</b>: 港湾上位 5 + 漁港上位 3 の<b>合計 8 港</b>で
全係留施設の <b>{(top5_h_sum + top3_f_sum)/(n_h_total+n_f_total)*100:.0f}%</b>
({top5_h_sum + top3_f_sum}/{n_h_total + n_f_total}) を占有。
広島県の物流・漁業は<b>少数のハブに収束</b>する偏在構造。L32 と同パターン。</li>
<li><b>(3) 外郭の内側に係留 (階層整備)</b>: 係留の <b>{(n_near_h+n_near_f)/max(1,n_total_h+n_total_f)*100:.0f}%</b>
が外郭から {PROXIMITY_M:.0f}m 圏内に位置する。
これは「外郭→係留」の 2 層整備の実証で、
広島県の港湾は<b>「外洋からの波を外郭で遮り、その内側に船を泊める」</b>という階層的設計。</li>
</ul>

<p>本記事は<b>「係留施設は港の生産機能の核」</b>という視点を実データで裏付けた。
{n_total_all} 件の係留が、<b>2 つの法的カテゴリ</b>で 41 港の<b>船舶受入容量</b>を形成し、
<b>L32 で示した外郭施設の防護網の内側</b>に階層配置されている。
この網の幾何構造を理解することは、<b>港湾の物流・漁業・防災の総合分析</b>の出発点である。</p>

<h3>2 港の比較で見る係留思想 (=本記事タイトルの問いに答える)</h3>
<p>2 件のデータが「2 港分」ではなく「2 行政カテゴリ分」だったため、
本記事は<b>「広島港 vs 福山港」のような 2 港比較</b>ではなく、
<b>「港湾 vs 漁港」の二相比較 + 上位 4 港の精読</b>を主題とした。これにより
<b>27 + 14 = 41 港</b>の俯瞰が可能になり、研究的により価値ある分析になった。</p>
<p>具体的な係留思想の相違:</p>
<table>
<tr><th>軸</th><th>港湾 (商業港)</th><th>漁港 (漁業基地)</th></tr>
<tr><td>主役構造 (件数)</td><td>物揚場 ({ratios_h[1]:.0f}%) + 浮さん橋 ({ratios_h[2]:.0f}%)</td>
<td>物揚場 ({ratios_f[1]:.0f}%) + 浮さん橋 ({ratios_f[2]:.0f}%)</td></tr>
<tr><td>主役構造 (延長)</td><td>岸壁 ({ratios_h2[0]:.0f}%) + 物揚場 ({ratios_h2[1]:.0f}%)</td>
<td>物揚場 ({ratios_f2[1]:.0f}%) + 岸壁 ({ratios_f2[0]:.0f}%)</td></tr>
<tr><td>岸壁 1 本サイズ</td><td>中央値 {qw_h_med:.0f} m (大型化)</td><td>中央値 {qw_f_med:.0f} m (中型のみ)</td></tr>
<tr><td>大型岸壁 (≥{LARGE_QUAY_THRESH:.0f}m)</td><td>{int(quays_large_pv[quays_large_pv['port_category']=='港湾']['n'].sum() if (quays_large_pv['port_category']=='港湾').any() else 0)} 本</td><td>{int(quays_large_pv[quays_large_pv['port_category']=='漁港']['n'].sum() if (quays_large_pv['port_category']=='漁港').any() else 0)} 本</td></tr>
<tr><td>カテゴリ計延長</td><td>{total_h:.1f} km</td><td>{total_f:.1f} km</td></tr>
<tr><td>船種</td><td>大型貨物・客船・フェリー・小型船</td><td>漁船 + 小型レジャーボート</td></tr>
<tr><td>地理</td><td>湾奥・主要市街地・島嶼商業港</td><td>島嶼 + 沿岸集落</td></tr>
<tr><td>外郭から{PROXIMITY_M:.0f}m圏内率</td><td>{near_pct_h if L32_OK else 0:.1f}%</td><td>{near_pct_f if L32_OK else 0:.1f}%</td></tr>
</table>
"""
sections.append(("仮説検証と考察", sec12))


# === Section 13: 発展課題 ===
sec13 = f"""
<h3>結果 X1 → 新仮説 Y1 → 課題 Z1: 維持管理状態と岸壁劣化リスク</h3>
<ul>
<li><b>結果 X1</b>: 本記事は係留施設の<b>地理構造</b>に集中したが、
DoBoX dataset 1251/1255 にも実は<b>「維持管理情報」</b>列が存在する (本ダウンロードでは別ファイル)。</li>
<li><b>新仮説 Y1</b>: 古い港 (例: 御手洗港、明治期以来の古港) の岸壁は<b>劣化リスク</b>が高い。
特に大型岸壁 ({len(quays_large)} 本) のうち築年数 50 年超のものは早期補修が必要。
劣化マップは港の歴史と相関する。</li>
<li><b>課題 Z1</b>: 維持管理 dataset (1251-1255 系列の他リソース) を取得し、健全度・補修履歴を
本記事の地理データに空間結合。<b>「劣化スコア × 岸壁延長」</b>のマップを作成し、
補修優先岸壁を抽出する。L32 の外郭劣化と組み合わせれば<b>港全体の補修優先度</b>が出せる。</li>
</ul>

<h3>結果 X2 → 新仮説 Y2 → 課題 Z2: 外郭整備時期と係留整備時期の時系列順序</h3>
<ul>
<li><b>結果 X2</b>: 係留の {(n_near_h+n_near_f)/max(1,n_total_h+n_total_f)*100:.0f}% が外郭から {PROXIMITY_M:.0f}m 圏内に配置されており、
外郭→係留の階層整備が空間的に検証された。</li>
<li><b>新仮説 Y2</b>: 港湾整備の歴史的順序は<b>「先に外郭 (防波堤) を築き、その内側に係留 (岸壁・物揚場) を整備」</b>。
広島港の場合、明治期に北防波堤、大正期に宇品の岸壁、戦後にコンテナバース、
というように<b>外郭→係留の時系列</b>がある。</li>
<li><b>課題 Z2</b>: 港湾の整備史記録 (港湾管理者の記念誌・年報) と本データを照合し、
施設番号 (B-1-01 = 防波堤、C-1-01 = 係留) の<b>付番順</b>から整備時期を推定。
時系列マップで「先外郭・後係留」の順序を可視化する。
これは<b>港湾発達の地理史</b>研究の入口になる。</li>
</ul>

<h3>結果 X3 → 新仮説 Y3 → 課題 Z3: 大型岸壁と貨物量・船舶往来の対応</h3>
<ul>
<li><b>結果 X3</b>: 大型岸壁 ({len(quays_large)} 本) は <b>{', '.join(quays_large_pv['港湾名'].head(3).tolist())}</b> に集中。</li>
<li><b>新仮説 Y3</b>: 大型岸壁の総延長は<b>港の年間貨物取扱量 (トン)</b>と強く相関する。
広島港・福山港の貨物量データ (国交省港湾統計) と岸壁長を散布すれば、
<b>1 m あたり貨物量 (生産性指標)</b> が港間で比較できる。</li>
<li><b>課題 Z3</b>: 国交省<b>港湾統計年報</b> (公開) を取得し、
本記事の岸壁延長と港別貨物量を結合。<b>「岸壁 1 m あたりトン数」</b>を計算し、
県内で最も生産性の高い岸壁を特定。経済学の<b>資本生産性</b>の港湾版指標。</li>
</ul>

<h3>結果 X4 → 新仮説 Y4 → 課題 Z4: 浮さん橋の潮位対応の真の必要性検証</h3>
<ul>
<li><b>結果 X4</b>: 浮さん橋 {fp_h+fp_f} 件は港湾優位 ({fp_h_pct:.1f}% > 漁港 {fp_f_pct:.1f}%)、H3 を支持。</li>
<li><b>新仮説 Y4</b>: 浮さん橋の必要性は<b>潮位差 × 旅客船発着頻度</b>で決まる。
潮位差の大きい湾奥 (太田川河口、芦田川河口) で旅客船発着の多い港に集中するはず。</li>
<li><b>課題 Z4</b>: 観測情報シリーズの<b>潮位日集計 (dsid 1441)</b> から各港最寄りの潮位観測点の
日次潮位差を取得。本記事の浮さん橋位置と空間結合し、
<b>「潮位差 × 浮さん橋本数」</b>の散布図で必要性の根拠を検証する。
潮位差の小さい港にも浮さん橋がある場合は「歴史的経緯 or 旅客船需要」が要因と推察。</li>
</ul>

<h3>結果 X5 → 新仮説 Y5 → 課題 Z5: 物揚場の規格化 (中央値 {mp_h_med:.0f}m vs {mp_f_med:.0f}m) の歴史的起源</h3>
<ul>
<li><b>結果 X5</b>: 物揚場 1 本の中央値は港湾 {mp_h_med:.0f}m、漁港 {mp_f_med:.0f}m とかなり均質。</li>
<li><b>新仮説 Y5</b>: 物揚場のサイズが標準化されているのは、<b>農林水産省の漁港整備指針</b>
または<b>国土交通省の港湾施設技術基準</b>に「標準物揚場 60m」のような規格があるため。
規格を超える/下回る物揚場は<b>例外的事情</b> (例: 既存岸壁を分割改造した物揚場) を持つ。</li>
<li><b>課題 Z5</b>: 港湾施設技術基準・漁港整備指針の公開資料を読み、
<b>標準物揚場の長さ規定</b>を特定。本記事の物揚場長分布と照合し、
規格適合率を計算。例外物揚場の地理分布から<b>歴史的経緯</b>を推察する。</li>
</ul>
"""
sections.append(("発展課題", sec13))


# レンダリング
html = render_lesson(
    num=33,
    title="係留施設 2 件統合分析 — 広島県 27 港湾 + 14 漁港の船舶受入構造",
    tags=["係留", "GIS", "港湾", "漁港", "geopandas", "WKT", "岸壁", "物揚場", "外郭との階層"],
    time="40 分",
    level="リテラシ",
    data_label="DoBoX 係留施設 2 dataset (1251 港湾, 1255 漁港)",
    sections=sections,
    script_filename="L33_port_mooring.py",
)

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


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