# -*- coding: utf-8 -*-
"""L88 M4 産業物語 — 福山JFE・大竹コンビナート・呉海軍工廠・工業専用地域

カバー宣言:
  本記事は DoBoX の 用途地域 (#794-946, #932) / 港湾施設 (#1250-1252) /
  漁港施設 (#1254-1255) / 高潮浸水想定 max (#45) / 橋梁基本情報 (#11) /
  行政界 (L15 cache) / 海岸保全 (#1253) を統合し、
  「<b>広島県の産業物語 — 戦前 - 戦後の産業立地が織りなす都市
  の物語</b>」 を 6 章構成で物語る<b>テーマ統合 (M系) 記事</b>である。
  単一データの分析ではなく、 「<b>戦前 - 戦後 - 沿岸 - 防災</b>」 の時間
  + 空間軸で、 広島県の産業地理を物語る。

  M系 (テーマ統合) の意味:
    M1 (L85) = 太田川水系を 1 水系深掘りで物語化 (河川治水の物語)
    M2 (L86) = 埋蔵文化財を時間軸で物語化 (古代地理の物語)
    M3 (L87) = 観光関連 4 dataset を空間軸 (沿岸-島嶼-しまなみ) で物語化
    M4 (本記事) = 工業専用地域 + 港湾 + 防災 を 4 産業拠点 (福山 / 大竹 /
                  呉 / 尾道) の物語軸で再構成 (産業地理物語)

研究の問い (主 RQ):
  「広島県の産業物語 — 戦前 - 戦後の産業立地が織りなす都市の物語とは
  何か?」 戦前 (旧海軍工廠呉) ・戦前化学 (大竹コンビナート) ・戦後高度
  成長 (福山JFE) ・尾道造船 という 4 産業拠点は、 工業専用地域・港湾
  施設・橋梁・浸水リスクの量と質でどう描けるか?  「臨海重工業」 という
  広島県の<b>産業 DNA</b>はデータでどう浮かぶか?

仮説 H1〜H5:
  H1 (臨海重工業集中仮説): 工業専用地域 (YOTO_CD=12) 全 18 市町中、
       <b>沿岸 4 産業拠点 (福山・呉・大竹・尾道)</b>の合計面積比率は
       <b>50% 以上</b>。 広島県の重工業帯は<b>瀬戸内沿岸に偏在</b>する。
       本仮説は L17 用途地域データの市町別 dissolve で量化する。

  H2 (福山JFE圧倒仮説): 工業専用地域 (YOTO_CD=12) の<b>市町別最大</b>
       は<b>福山市</b>であり、 県全体の<b>20% 以上</b>を占める。 戦後
       高度成長期 (1965 年福山製鉄所開設) が広島県工業地理を再編した
       量的証拠。 また<b>福山港の係留施設</b>は産業港として上位に位置
       する。

  H3 (港湾規模 vs 工業面積 相関仮説): 各市町の<b>工業専用地域面積</b>と
       <b>港湾係留施設件数</b>は<b>正の相関</b> (Pearson r ≥ 0.5) を
       示す。 大規模工業地帯 = 港湾物流の集積地、 という<b>臨海工業
       立地原理</b>を示す。

  H4 (沿岸工業 vs 高潮交差仮説): 工業専用地域ポリゴンのうち、 <b>高潮
       浸水想定エリアと交差</b>する面積比は<b>40% 以上</b>。 沿岸埋立地
       中心の工業地は<b>気候災害の最前線</b>に位置する。

  H5 (産業港 vs 漁港 二極化仮説): 港湾施設 (港湾) と 漁港施設 (漁港)
       は地理的に<b>明確に分離</b>される。 港湾の<b>係留施設件数 上位
       5 港</b>と漁港の<b>外郭施設件数 上位 5 港</b>の<b>市町重複</b>は
       <b>1 件以下</b>。 「重工業港 vs 漁業港」 という機能分化が空間
       的に成立。

要件 S 準拠 (1 分以内完走):
  - L17 cache (21 zip × 平均 1MB) の dissolve 集計のみ
  - 港湾 CSV は軽量 (合計 < 3000 行)
  - 高潮 shp は 7 行のみ (SINSUIRANK 別 dissolved)
  - 津波 shp は 1.25M メッシュだが bbox フィルタ (各産業拠点 5km 圏内
    のみ) で実用速度 + sjoin
  - 想定完走時間 30 - 90 秒

要件 T 準拠 (位置情報あり = 地図必須):
  - 工業専用地域 県全域 主題図 (市町別色分け)
  - 4 産業拠点ズーム × 4 (福山 / 呉 / 大竹 / 尾道)
  - 工業専用地域 vs 高潮浸水想定 重ね合わせマップ
  - 工業専用地域 vs 港湾係留施設 重ね合わせマップ
  - 産業 4 拠点 + 橋梁 small multiples

要件 Q 準拠: 図 8+ / 表 11+ (市町 × 用途 × 港湾 × 防災 多角度)

X系 / L17 / L21 / L32-L34 / L52-L55 / L66 との重複回避 (徹底):
  L17 (用途地域 21 市町統合) = 13 用途地域 + 3 類型 + k-means プロファイル
                              + 工業専用地域は 1 章扱いのみ
  L32 (港湾外郭) = 防波堤・防砂堤 480 件の単独深掘り (港湾名別ランキング)
  L33 (港湾係留) = 物揚場・浮さん橋 802 件単独深掘り
  L34 (臨港交通) = 道路・駐車場 278 件単独深掘り
  L21 (宅地開発 1,038 件) = 期間・面積・年別の単独深掘り
  L52-L55 (盛土系) = 盛土申請・届出・工区の単独深掘り
  L66 (橋梁 4,203 件) = 架設年度・延長・幅員の単独深掘り
  L44 (高潮浸水想定) = SINSUIRANK 別深さ分布の単独深掘り

  L88 (本記事) は<b>これら全てと厳密に重複しない切り口</b>:
    (i)   <b>工業専用地域 (YOTO_CD=12) のみ</b>を取り出して、 4 産業拠点
          の物語軸で<b>再分類</b>。 L17 は 13 用途地域並列、 L88 は
          工業専用 1 用途のみで深掘り。
    (ii)  <b>4 産業拠点 (福山JFE / 呉海軍工廠 / 大竹コンビナート /
          尾道造船)</b>の<b>歴史軸</b>で語る。 L17 は時系列を扱わない。
    (iii) <b>港湾係留施設 × 工業専用地域</b>のクロス相関 (H3) は本記事が
          初。 L33 は港湾名別ランキング、 L17 は用途地域単独。
    (iv)  <b>工業専用地域 × 高潮浸水想定 (#1149)</b> 交差量化は本記事が
          初。 L44 は浸水想定単独、 L17 は災害扱わず。
    (v)   <b>産業港 vs 漁港 機能分化</b> (H5) の空間検証は本記事が初。
          L32-L34 は港湾全体、 漁港との比較なし。
    (vi)  既存記事は<b>シリーズ網羅 / 多角集計 / 単独深掘り</b>に集中、
          L88 は<b>歴史物語軸</b>でデータを再編する切り口。

実行:
    cd "2026 DoBoX 教材"
    py -X utf8 lessons/L88_M4_industry_story.py
"""
from __future__ import annotations

import sys
import time
import zipfile
import io
import warnings
from pathlib import Path
from html import escape

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

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

warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)

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

t_all = time.time()
print("=== L88 M4 産業物語 — 福山JFE・大竹コンビナート・"
      "呉海軍工廠・工業専用地域 ===", flush=True)


# =============================================================================
# 0. 定数・パス
# =============================================================================
TARGET_CRS = "EPSG:6671"  # JGD2011 平面直角第 III 系 (m)

DATA_EX = ROOT / "data" / "extras"
DATA_DIR = DATA_EX / "L88_M4_industry_story"
DATA_DIR.mkdir(parents=True, exist_ok=True)

USEZONE_DIR = DATA_EX / "L17_use_zones"
ADMIN_DIR = DATA_EX / "L15_admin_zones"
PORT_OUTER = DATA_EX / "L32_port_breakwaters" / "harbor_outer_facility.csv"
FISH_OUTER = DATA_EX / "L32_port_breakwaters" / "fishing_port_outer_facility.csv"
PORT_MOORING = DATA_EX / "L33_port_mooring" / "harbor_mooring_facility.csv"
FISH_MOORING = DATA_EX / "L33_port_mooring" / "fishing_port_mooring_facility.csv"
PORT_TRAFFIC = DATA_EX / "L34_port_traffic" / "harbor_traffic_facility.csv"
STORM_SHP = DATA_EX / "L44_storm_surge" / "max_20210803" / "takashio_max.shp"
TSUNAMI_SHP = (DATA_EX / "tsunami_extracted"
               / "340006_tsunami_inundation_assumption_map_20251203"
               / "浸水メッシュ.shp")
BRIDGE_CSV = DATA_EX / "bridge_basic.csv"
COAST_CSV = DATA_EX / "L64_coast_protection" / "340006_coastal_equipment_20230509.csv"

# 用途地域 21 市町 (L17 と同じ)
CITY_DEFS = [
    (794, "広島市",     True,  "政令市"),
    (804, "呉市",       True,  "市"),
    (811, "竹原市",     True,  "市"),
    (821, "三原市",     True,  "市"),
    (831, "尾道市",     True,  "市"),
    (837, "福山市",     True,  "市"),
    (847, "府中市",     False, "市"),
    (855, "三次市",     False, "市"),
    (861, "庄原市",     False, "市"),
    (867, "大竹市",     True,  "市"),
    (875, "東広島市",   True,  "市"),
    (885, "廿日市市",   True,  "市"),
    (893, "安芸高田市", False, "市"),
    (899, "江田島市",   True,  "離島自治体"),
    (904, "府中町",     False, "町"),
    (910, "海田町",     True,  "町"),
    (915, "熊野町",     False, "町"),
    (921, "坂町",       True,  "町"),
    (940, "北広島町",   False, "町"),
    (946, "世羅町",     False, "町"),
]

# 4 産業拠点 (本記事独自) — 名称・略歴・lat/lon中心
INDUSTRIAL_HUBS = {
    "福山JFE":   {
        "city": "福山市", "lat": 34.490, "lon": 133.395,
        "story": "戦後高度成長期 1965 年に開設した福山製鉄所 "
                 "(現 JFEスチール西日本製鉄所福山地区)。 "
                 "鉄を作るために臨海埋立地を造成し、 "
                 "高炉 4 基 + 圧延・厚板を擁する世界級鉄鋼コンビナート",
        "era": "戦後高度成長 (1965)",
        "color": "#cf222e",
    },
    "大竹コンビナート": {
        "city": "大竹市", "lat": 34.230, "lon": 132.230,
        "story": "戦前から続く化学コンビナート。 "
                 "三井化学・三菱レイヨンなどの石油化学・繊維産業が "
                 "瀬戸内沿岸の狭い平地に密集する。 ナフサ輸入港として "
                 "大竹港が機能",
        "era": "戦前化学 (1920s-)",
        "color": "#9a6700",
    },
    "呉海軍工廠":   {
        "city": "呉市", "lat": 34.232, "lon": 132.566,
        "story": "戦前は東洋一の海軍工廠 (1903 - 1945)。 "
                 "大和・武蔵級戦艦を建造。 戦後は IHI 呉造船所・神戸製鋼 "
                 "として民需転換。 軍需 → 造船 → 鉄鋼 の DNA 継承",
        "era": "戦前軍需 (1903) → 戦後造船・鉄鋼",
        "color": "#0969da",
    },
    "尾道造船":     {
        "city": "尾道市", "lat": 34.405, "lon": 133.185,
        "story": "戦前から続く尾道造船・常石造船・神原汽船などの "
                 "造船・海運集積地。 瀬戸内最大の造船クラスタの一翼を "
                 "担い、 戦後の輸出造船景気で拡大",
        "era": "戦前造船 → 戦後輸出造船",
        "color": "#8250df",
    },
}

# 沿岸 4 産業拠点市町 (H1 検証用)
HUB_CITIES = ["福山市", "呉市", "大竹市", "尾道市"]

# 産業地区 zoom 用 bbox (lat_min, lat_max, lon_min, lon_max) - 拠点 ±5km 程度
HUB_BBOX = {
    "福山JFE":   (34.43, 34.55, 133.32, 133.45),
    "大竹コンビナート": (34.18, 34.27, 132.18, 132.27),
    "呉海軍工廠":   (34.18, 34.28, 132.50, 132.62),
    "尾道造船":   (34.36, 34.43, 133.10, 133.25),
}


YOTO_TARGET = 12  # 工業専用地域
YOTO_NAMES = {
    1: "1低専", 2: "2低専", 3: "1中高専", 4: "2中高専",
    5: "1住居", 6: "2住居", 7: "準住居",
    8: "近商",  9: "商業",
    10: "準工", 11: "工業", 12: "工業専用",
    13: "田園住居",
}


# =============================================================================
# 1. データ確認
# =============================================================================
print("\n[1] データ確認", flush=True)
t1 = time.time()

required = [
    (USEZONE_DIR, "L17 用途地域 21 zip"),
    (PORT_OUTER, "L32 港湾外郭"),
    (PORT_MOORING, "L33 港湾係留"),
    (PORT_TRAFFIC, "L34 臨港交通"),
    (FISH_OUTER, "漁港外郭"),
    (FISH_MOORING, "漁港係留"),
    (STORM_SHP, "L44 高潮浸水想定 max"),
    (BRIDGE_CSV, "橋梁 #1062"),
    (COAST_CSV, "海岸保全 #1085"),
    (ADMIN_DIR, "L15 行政界"),
]
for p, label in required:
    if not p.exists():
        raise FileNotFoundError(f"{label}: {p}")
print(f"  全 10 系統データ確認 OK ({time.time()-t1:.2f}s)", flush=True)


# =============================================================================
# 2. 用途地域 GeoJSON 統合 + 工業専用地域抽出
# =============================================================================
print("\n[2] 用途地域 21 市町統合 + 工業専用地域抽出", flush=True)
t2 = time.time()


def load_geojson_zip(zip_path: Path) -> gpd.GeoDataFrame:
    with zipfile.ZipFile(zip_path) as zf:
        gjs = [n for n in zf.namelist() if n.endswith(".geojson")][0]
        with zf.open(gjs) as f:
            return gpd.read_file(io.BytesIO(f.read()))


frames = []
for dsid, name, coastal, ctype in CITY_DEFS:
    z = USEZONE_DIR / f"use_zone_{dsid}_{name}.zip"
    g = load_geojson_zip(z)
    g["source_city"] = name
    g["source_dsid"] = dsid
    g["coastal"] = coastal
    frames.append(g)
zone = gpd.GeoDataFrame(pd.concat(frames, ignore_index=True),
                        geometry="geometry", crs=frames[0].crs).to_crs(TARGET_CRS)
zone["area_km2"] = zone.geometry.area / 1e6

# 工業専用地域 (YOTO_CD=12) のみ抽出
indzone = zone[zone["YOTO_CD"] == YOTO_TARGET].copy()
N_INDZONE = len(indzone)
A_INDZONE = float(indzone["area_km2"].sum())

# 市町別 工業専用地域 集計
city_indarea = (indzone.groupby("source_city")["area_km2"].sum()
                .sort_values(ascending=False))

# 4 産業拠点市町の合計
hub_area = float(city_indarea.reindex(HUB_CITIES).fillna(0).sum())
hub_share_pct = hub_area / max(A_INDZONE, 1e-9) * 100

# 市町別 全用途地域 (背景)
zone_diss = zone.dissolve(by="source_city").reset_index()

print(f"  用途地域 全 {len(zone):,} ポリゴン", flush=True)
print(f"  工業専用地域 (YOTO_CD={YOTO_TARGET}): "
      f"{N_INDZONE:,} ポリゴン / 総面積 {A_INDZONE:.2f} km²", flush=True)
print(f"  4 拠点市町 (福山+呉+大竹+尾道) 合計 "
      f"{hub_area:.2f} km² ({hub_share_pct:.1f}%)", flush=True)
print(f"  ({time.time()-t2:.2f}s)", flush=True)


# =============================================================================
# 3. 港湾施設 + 漁港施設 統合
# =============================================================================
print("\n[3] 港湾 + 漁港施設 統合", flush=True)
t3 = time.time()


def load_port(csv_path, source_label, kind):
    df = pd.read_csv(csv_path, encoding="utf-8-sig")
    df["lat"] = pd.to_numeric(df["開始位置緯度"], errors="coerce")
    df["lon"] = pd.to_numeric(df["開始位置経度"], errors="coerce")
    df = df.dropna(subset=["lat", "lon"]).copy()
    df["source"] = source_label
    df["kind"] = kind  # "港湾" or "漁港"
    return df[["source", "kind", "港湾名", "施設分類", "施設種類",
               "施設名称", "市区町村１", "lat", "lon"]]


port_outer = load_port(PORT_OUTER, "港湾_外郭", "港湾")
port_moor = load_port(PORT_MOORING, "港湾_係留", "港湾")
port_traf = load_port(PORT_TRAFFIC, "港湾_臨港交通", "港湾")
fish_outer = load_port(FISH_OUTER, "漁港_外郭", "漁港")
fish_moor = load_port(FISH_MOORING, "漁港_係留", "漁港")

ports_all = pd.concat([port_outer, port_moor, port_traf,
                       fish_outer, fish_moor], ignore_index=True)
N_PORT_ALL = len(ports_all)
N_PORT_HARBOR = int((ports_all["kind"] == "港湾").sum())
N_PORT_FISH = int((ports_all["kind"] == "漁港").sum())

ports_gdf = gpd.GeoDataFrame(
    ports_all,
    geometry=[Point(lon, lat) for lon, lat in zip(ports_all["lon"],
                                                    ports_all["lat"])],
    crs="EPSG:4326").to_crs(TARGET_CRS)
ports_gdf["x_m"] = ports_gdf.geometry.x
ports_gdf["y_m"] = ports_gdf.geometry.y

print(f"  港湾 (3 系統 = 外郭+係留+臨港交通) "
      f"{N_PORT_HARBOR:,} 件", flush=True)
print(f"  漁港 (2 系統 = 外郭+係留) {N_PORT_FISH:,} 件", flush=True)
print(f"  統合: {N_PORT_ALL:,} 件", flush=True)
print(f"  ({time.time()-t3:.2f}s)", flush=True)


# =============================================================================
# 4. 港湾名別 ランキング (H5 検証用)
# =============================================================================
print("\n[4] 港湾名別 施設件数ランキング (H5)", flush=True)
t4 = time.time()

# 港湾 (重工業港) — 係留施設 件数
harbor_moor_rank = (port_moor.groupby("港湾名").size()
                    .reset_index(name="係留施設数")
                    .sort_values("係留施設数", ascending=False))

# 漁港 — 外郭施設 件数 (漁港のメイン施設は外郭= 防波堤)
fishport_outer_rank = (fish_outer.groupby("港湾名").size()
                       .reset_index(name="外郭施設数")
                       .sort_values("外郭施設数", ascending=False))

top5_harbor = harbor_moor_rank.head(5)["港湾名"].tolist()
top5_fishport = fishport_outer_rank.head(5)["港湾名"].tolist()
overlap_top5 = set(top5_harbor) & set(top5_fishport)
N_OVERLAP = len(overlap_top5)

print(f"  港湾係留 Top5: {top5_harbor}", flush=True)
print(f"  漁港外郭 Top5: {top5_fishport}", flush=True)
print(f"  重複: {N_OVERLAP} 件 ({list(overlap_top5)})", flush=True)


# =============================================================================
# 5. 工業専用地域 vs 高潮浸水想定 交差 (H4)
# =============================================================================
print("\n[5] 工業専用地域 vs 高潮浸水想定 交差", flush=True)
t5 = time.time()

storm_gdf = gpd.read_file(STORM_SHP).to_crs(TARGET_CRS)
storm_union = storm_gdf.geometry.union_all()
indzone["storm_overlap_m2"] = indzone.geometry.intersection(storm_union).area
indzone["storm_overlap_pct"] = (
    indzone["storm_overlap_m2"] / indzone.geometry.area.clip(lower=1) * 100)
indzone["in_storm"] = (indzone["storm_overlap_m2"] > 1).astype(int)
storm_overlap_total_m2 = float(indzone["storm_overlap_m2"].sum())
storm_overlap_total_km2 = storm_overlap_total_m2 / 1e6
storm_pct = storm_overlap_total_km2 / max(A_INDZONE, 1e-9) * 100

print(f"  高潮浸水想定 {len(storm_gdf)} ランク → 統合済", flush=True)
print(f"  工業専用 ∩ 高潮 = {storm_overlap_total_km2:.2f} km² "
      f"({storm_pct:.1f}%)", flush=True)
print(f"  ({time.time()-t5:.2f}s)", flush=True)


# =============================================================================
# 6. 4 産業拠点 zoom 用フィルタ
# =============================================================================
print("\n[6] 4 産業拠点 zoom 集計", flush=True)
t6 = time.time()


def in_bbox(lat, lon, bbox):
    """bbox = (lat_min, lat_max, lon_min, lon_max)"""
    return bbox[0] <= lat <= bbox[1] and bbox[2] <= lon <= bbox[3]


hub_summary = []
for hub_name, info in INDUSTRIAL_HUBS.items():
    bbox = HUB_BBOX[hub_name]
    # 該当市町の工業専用地域面積
    city = info["city"]
    ind_area_city = float(city_indarea.get(city, 0.0))
    # 拠点 bbox 内の港湾施設件数
    sub_harbor = ports_gdf[
        (ports_gdf["kind"] == "港湾") &
        (ports_gdf.geometry.bounds["minx"].notna())
    ]
    # bbox は lat/lon なので元データで判定
    sub = ports_all[
        (ports_all["kind"] == "港湾") &
        (ports_all["lat"].between(bbox[0], bbox[1])) &
        (ports_all["lon"].between(bbox[2], bbox[3]))
    ]
    n_port = len(sub)
    sub_fish = ports_all[
        (ports_all["kind"] == "漁港") &
        (ports_all["lat"].between(bbox[0], bbox[1])) &
        (ports_all["lon"].between(bbox[2], bbox[3]))
    ]
    n_fish = len(sub_fish)
    hub_summary.append({
        "拠点": hub_name,
        "市町": city,
        "歴史": info["era"],
        "工業専用km2": round(ind_area_city, 2),
        "拠点bbox港湾件数": n_port,
        "拠点bbox漁港件数": n_fish,
    })
T_hub_summary = pd.DataFrame(hub_summary)

print(f"  拠点別サマリー完成", flush=True)
print(T_hub_summary.to_string(index=False), flush=True)
print(f"  ({time.time()-t6:.2f}s)", flush=True)


# =============================================================================
# 7. 行政界 (背景)
# =============================================================================
print("\n[7] 行政界 (L15 cache) 読込", flush=True)
t7 = time.time()

ADMIN_MAP = [
    ("広島市", 786), ("呉市", 797), ("竹原市", 807),
    ("三原市", 814), ("尾道市", 824), ("福山市", 832),
    ("府中市", 840), ("三次市", 850), ("庄原市", 856),
    ("大竹市", 862), ("東広島市", 868), ("廿日市市", 878),
    ("安芸高田市", 888), ("江田島市", 894), ("府中町", 900),
    ("海田町", 905), ("熊野町", 911), ("坂町", 916),
    ("北広島町", 935), ("世羅町", 941),
]
admin_frames = []
for cname, dsid in ADMIN_MAP:
    z = ADMIN_DIR / f"admin_{dsid}_{cname}.zip"
    if not z.exists():
        continue
    g = load_geojson_zip(z)
    g["city"] = cname
    admin_frames.append(g)
admin_all = gpd.GeoDataFrame(
    pd.concat(admin_frames, ignore_index=True),
    geometry="geometry", crs=admin_frames[0].crs).to_crs(TARGET_CRS)
admin_diss = admin_all.dissolve(by="city", as_index=False)
print(f"  行政界 {len(admin_diss)} 市町 dissolved "
      f"({time.time()-t7:.2f}s)", flush=True)


# =============================================================================
# 8. 工業専用地域 × 港湾係留 相関 (H3)
# =============================================================================
print("\n[8] 工業専用地域面積 × 港湾係留件数 相関 (H3)", flush=True)
t8 = time.time()

# 港湾CSV の「市区町村１」 列はすべて NaN なので、 緯度経度から
# 行政界に sjoin して市町名を付与する (要件 R: geopandas で実空間結合)
port_moor_gdf = gpd.GeoDataFrame(
    port_moor,
    geometry=[Point(lon, lat) for lon, lat in zip(port_moor["lon"],
                                                    port_moor["lat"])],
    crs="EPSG:4326").to_crs(TARGET_CRS)
# admin_diss は city 列を持つ
port_moor_with_city = gpd.sjoin(
    port_moor_gdf[["geometry"]], admin_diss[["city", "geometry"]],
    how="left", predicate="within"
)
# 市町別件数
city_harbor_moor = (port_moor_with_city["city"].value_counts()
                    .reset_index())
city_harbor_moor.columns = ["市町", "係留件数"]

# 市町別 工業専用面積
city_ind_df = city_indarea.reset_index()
city_ind_df.columns = ["市町", "工業専用km2"]

# join
corr_df = city_ind_df.merge(
    city_harbor_moor, on="市町", how="outer").fillna(0)
corr_df["工業専用km2"] = corr_df["工業専用km2"].astype(float)
corr_df["係留件数"] = corr_df["係留件数"].astype(int)

# 相関係数 (沿岸市町のみで計算 — 内陸は両方ゼロで意味薄)
mask_coastal = corr_df["市町"].isin(
    [c[1] for c in CITY_DEFS if c[2]])  # coastal=True
corr_coastal = corr_df[mask_coastal].copy()
if len(corr_coastal) >= 3 and corr_coastal["係留件数"].std() > 0:
    r_pearson = float(corr_coastal["工業専用km2"].corr(corr_coastal["係留件数"]))
else:
    r_pearson = 0.0

# 福山市・呉市・大竹市・尾道市の係留件数も別途確認用
hub_harbor_count = {
    h: int(city_harbor_moor[city_harbor_moor["市町"]
                             == INDUSTRIAL_HUBS[h]["city"]]["係留件数"].sum())
    for h in INDUSTRIAL_HUBS
}

print(f"  port_moor sjoin 結果: 市町割当 "
      f"{port_moor_with_city['city'].notna().sum()} / "
      f"{len(port_moor_with_city)}", flush=True)
print(f"  港湾係留 拠点別: {hub_harbor_count}", flush=True)
print(f"  沿岸 {len(corr_coastal)} 市町 工業専用 vs 係留件数 "
      f"Pearson r = {r_pearson:.3f}", flush=True)


# =============================================================================
# 9. 仮説検証
# =============================================================================
print("\n[9] 仮説検証 (5 仮説)", flush=True)
t9 = time.time()

# H1: 沿岸 4 拠点市町の工業専用地域シェア ≥ 50%
h1_value = hub_share_pct
h1_ok = h1_value >= 50.0

# H2: 福山市の工業専用地域シェア ≥ 20% (県全体)
fukuyama_area = float(city_indarea.get("福山市", 0.0))
fukuyama_share_pct = fukuyama_area / max(A_INDZONE, 1e-9) * 100
fukuyama_is_max = (city_indarea.idxmax() == "福山市")
h2_value = fukuyama_share_pct
h2_ok = (fukuyama_is_max and h2_value >= 20.0)

# H3: Pearson r ≥ 0.5
h3_value = r_pearson
h3_ok = r_pearson >= 0.5

# H4: 工業専用 ∩ 高潮 比率 ≥ 40%
h4_value = storm_pct
h4_ok = h4_value >= 40.0

# H5: 港湾係留 Top5 vs 漁港外郭 Top5 重複 ≤ 1
h5_value = N_OVERLAP
h5_ok = N_OVERLAP <= 1

T_hypothesis = pd.DataFrame([
    {"仮説": "H1 (臨海重工業集中)",
     "予測": "沿岸 4 拠点市町 工業専用地域シェア ≥ 50%",
     "実測": f"{hub_area:.2f} km² / {A_INDZONE:.2f} km² = {h1_value:.1f}%",
     "判定": "支持" if h1_ok else "反証"},
    {"仮説": "H2 (福山JFE圧倒)",
     "予測": "福山市が工業専用地域市町別最大 + 県シェア ≥ 20%",
     "実測": (f"福山市 {fukuyama_area:.2f} km² "
              f"({h2_value:.1f}%) "
              f"/ {('最大' if fukuyama_is_max else '最大ではない')}"),
     "判定": "支持" if h2_ok else "反証"},
    {"仮説": "H3 (港湾規模 vs 工業面積 相関)",
     "予測": "沿岸市町 工業専用 vs 係留件数 Pearson r ≥ 0.5",
     "実測": f"r = {h3_value:.3f}",
     "判定": "支持" if h3_ok else "反証"},
    {"仮説": "H4 (沿岸工業 vs 高潮交差)",
     "予測": "工業専用 ∩ 高潮浸水想定 ≥ 40%",
     "実測": f"{storm_overlap_total_km2:.2f} km² / "
            f"{A_INDZONE:.2f} km² = {h4_value:.1f}%",
     "判定": "支持" if h4_ok else "反証"},
    {"仮説": "H5 (産業港 vs 漁港 二極化)",
     "予測": "港湾係留 Top5 vs 漁港外郭 Top5 重複 ≤ 1",
     "実測": f"重複 = {h5_value} 港 ({list(overlap_top5)})",
     "判定": "支持" if h5_ok else "反証"},
])

n_supported = int((T_hypothesis["判定"] == "支持").sum())
print(f"  支持: {n_supported}/5", flush=True)
print(T_hypothesis.to_string(index=False), flush=True)
print(f"  ({time.time()-t9:.2f}s)", flush=True)


# =============================================================================
# 10. 集計表 (11 種以上)
# =============================================================================
print("\n[10] 集計表 11+ 種", flush=True)
t10 = time.time()

# (1) 市町別 工業専用地域 ランキング (top 10)
T_city_ind_rank = (city_indarea.head(10).reset_index()
                   .rename(columns={"source_city": "市町",
                                      "area_km2": "工業専用km2"}))
T_city_ind_rank["工業専用km2"] = T_city_ind_rank["工業専用km2"].round(2)
T_city_ind_rank["シェア%"] = (T_city_ind_rank["工業専用km2"]
                              / max(A_INDZONE, 1e-9) * 100).round(1)

# (2) 用途地域 13 類型 全県集計 (背景)
T_zone_type = (zone.groupby("YOTO_CD")["area_km2"].sum()
               .reset_index()
               .rename(columns={"YOTO_CD": "YOTO_CD",
                                "area_km2": "面積km2"}))
T_zone_type["用途地域名"] = T_zone_type["YOTO_CD"].map(YOTO_NAMES)
T_zone_type["シェア%"] = (T_zone_type["面積km2"]
                          / T_zone_type["面積km2"].sum() * 100).round(2)
T_zone_type["面積km2"] = T_zone_type["面積km2"].round(2)
T_zone_type = T_zone_type[["YOTO_CD", "用途地域名", "面積km2", "シェア%"]]

# (3) 港湾名別 係留施設件数 Top10 (重工業港候補)
T_harbor_top = harbor_moor_rank.head(10)

# (4) 漁港名別 外郭施設件数 Top10 (漁業港候補)
T_fish_top = fishport_outer_rank.head(10)

# (5) 4 拠点 サマリー (上で作成済 T_hub_summary)
# (6) 工業専用 vs 高潮 拠点別交差
hub_storm = []
for hub_name, info in INDUSTRIAL_HUBS.items():
    city = info["city"]
    sub = indzone[indzone["source_city"] == city]
    if len(sub) == 0:
        hub_storm.append({"拠点": hub_name, "市町": city,
                          "工業専用km2": 0.0,
                          "高潮交差km2": 0.0, "高潮交差%": 0.0})
        continue
    a_total = float(sub["area_km2"].sum())
    a_storm = float(sub["storm_overlap_m2"].sum() / 1e6)
    pct = a_storm / max(a_total, 1e-9) * 100
    hub_storm.append({"拠点": hub_name, "市町": city,
                      "工業専用km2": round(a_total, 2),
                      "高潮交差km2": round(a_storm, 2),
                      "高潮交差%": round(pct, 1)})
T_hub_storm = pd.DataFrame(hub_storm)

# (7) 4 拠点 港湾係留 詳細 (拠点市町別)
# 漁港 (fish_outer) も「市区町村１」 が NaN なので sjoin で市町付与
fish_outer_gdf = gpd.GeoDataFrame(
    fish_outer,
    geometry=[Point(lon, lat) for lon, lat in zip(fish_outer["lon"],
                                                    fish_outer["lat"])],
    crs="EPSG:4326").to_crs(TARGET_CRS)
fish_outer_with_city = gpd.sjoin(
    fish_outer_gdf[["港湾名", "geometry"]], admin_diss[["city", "geometry"]],
    how="left", predicate="within"
)
fish_city_count = (fish_outer_with_city["city"].value_counts()
                   .to_dict())

T_hub_port = []
for hub_name, info in INDUSTRIAL_HUBS.items():
    city = info["city"]
    n_port = hub_harbor_count.get(hub_name, 0)
    n_fish = int(fish_city_count.get(city, 0))
    # 港湾名 (係留 max) — sjoin 結果に港湾名があるか確認
    pmw = port_moor_with_city.copy()
    pmw["港湾名"] = port_moor.reset_index()["港湾名"]
    pmw_city = pmw[pmw["city"] == city]
    if len(pmw_city) > 0:
        max_port = pmw_city["港湾名"].mode().iloc[0]
    else:
        max_port = "—"
    T_hub_port.append({
        "拠点": hub_name,
        "市町": city,
        "港湾係留件数": n_port,
        "漁港外郭件数": n_fish,
        "港湾名 (係留 max)": max_port,
    })
T_hub_port = pd.DataFrame(T_hub_port)

# (8) 1 件追跡 (要件 K)
sample_poly = indzone.iloc[indzone["area_km2"].argmax()]
T_track = pd.DataFrame([
    {"段階": "0. 元 GeoJSON",
     "値": f"YOTO_CD={int(sample_poly['YOTO_CD'])} "
          f"({YOTO_NAMES[int(sample_poly['YOTO_CD'])]}) "
          f"in {sample_poly['source_city']}"},
    {"段階": "1. 投影 EPSG:6671",
     "値": f"area = {sample_poly['area_km2']:.4f} km²"},
    {"段階": "2. 抽出 (YOTO_CD=12)",
     "値": "工業専用地域フィルタ後の最大ポリゴン"},
    {"段階": "3. 重心 (lat, lon)",
     "値": (f"{sample_poly.geometry.centroid.coords[0]}"
            if sample_poly.geometry.centroid else "—")},
    {"段階": "4. 高潮交差 m²",
     "値": f"{sample_poly['storm_overlap_m2']:,.0f} m²"},
    {"段階": "5. 高潮交差比率",
     "値": f"{sample_poly['storm_overlap_pct']:.1f}%"},
    {"段階": "6. in_storm flag",
     "値": f"{int(sample_poly['in_storm'])}"},
])

# (9) 既存記事比較 (重複ゼロ宣言)
T_compare = pd.DataFrame([
    {"記事": "L17 (用途地域 21 統合)",
     "主軸": "13 用途地域 + 3 類型 + k-means プロファイル",
     "主結果": "市町別ベクトル → 4 クラスタ",
     "L88 との重複": "L88 は工業専用 (YOTO_CD=12) のみ抽出"},
    {"記事": "L21 (宅地開発 1,038 件)",
     "主軸": "事業期間・面積・年別 単独深掘り",
     "主結果": "中央値 12 ヶ月、対数正規",
     "L88 との重複": "L88 は宅地開発扱わず、産業地のみ"},
    {"記事": "L32 (港湾外郭 480 件)",
     "主軸": "防波堤・防砂堤の港湾名別",
     "主結果": "広島港 61 件 max",
     "L88 との重複": "L88 は産業港の機能分化のみ"},
    {"記事": "L33 (港湾係留 802 件)",
     "主軸": "物揚場・浮さん橋ランキング",
     "主結果": "広島港 125 件 max",
     "L88 との重複": "L88 は工業専用地域との相関のみ"},
    {"記事": "L34 (臨港交通 278 件)",
     "主軸": "道路・駐車場ランキング",
     "主結果": "—",
     "L88 との重複": "L88 は臨港交通扱わず"},
    {"記事": "L44 (高潮浸水想定 7 ランク)",
     "主軸": "SINSUIRANK 別深さ分布",
     "主結果": "—",
     "L88 との重複": "L88 は工業専用との交差のみ"},
    {"記事": "L66 (橋梁 4,203 件)",
     "主軸": "架設年度・延長・幅員",
     "主結果": "中央値 1978 年",
     "L88 との重複": "L88 は産業拠点周辺の橋梁のみ参照"},
    {"記事": "L52-L55 (盛土系)",
     "主軸": "盛土申請・届出・工区",
     "主結果": "—",
     "L88 との重複": "L88 は盛土扱わず、産業地のみ"},
    {"記事": "L88 (本記事, 物語)",
     "主軸": "工業専用 + 港湾 + 防災 を 4 産業拠点で物語化",
     "主結果": "5 仮説検証 (臨海工業 + 福山JFE + 港湾相関 + 高潮 + 機能分化)",
     "L88 との重複": "—"},
])

# (10) 13 用途地域 × 沿岸/内陸 ピボット (背景比較)
zone_grouped = zone.copy()
zone_grouped["coast_label"] = zone_grouped["coastal"].map(
    {True: "沿岸", False: "内陸"})
T_zone_coast = (zone_grouped.groupby(["YOTO_CD", "coast_label"])["area_km2"]
                .sum().unstack(fill_value=0).round(2))
T_zone_coast["用途地域名"] = T_zone_coast.index.map(YOTO_NAMES)
T_zone_coast = T_zone_coast.reset_index()
T_zone_coast = T_zone_coast[["YOTO_CD", "用途地域名", "沿岸", "内陸"]]

# (11) 仮説検証総合 (T_hypothesis 既に作成済)
# (12) 工業専用地域 × 高潮 拠点別 (T_hub_storm)

print(f"  集計表 12+ 完成 ({time.time()-t10:.2f}s)", flush=True)


# =============================================================================
# 11. 中間 CSV 保存
# =============================================================================
print("\n[11] 中間 CSV 保存", flush=True)
t11 = time.time()

T_city_ind_rank.to_csv(ASSETS / "L88_city_industrial_rank.csv",
                        index=False, encoding="utf-8-sig")
T_zone_type.to_csv(ASSETS / "L88_zone_type_summary.csv",
                    index=False, encoding="utf-8-sig")
T_harbor_top.to_csv(ASSETS / "L88_harbor_mooring_top.csv",
                     index=False, encoding="utf-8-sig")
T_fish_top.to_csv(ASSETS / "L88_fishport_outer_top.csv",
                   index=False, encoding="utf-8-sig")
T_hub_summary.to_csv(ASSETS / "L88_hub_summary.csv",
                      index=False, encoding="utf-8-sig")
T_hub_storm.to_csv(ASSETS / "L88_hub_storm_overlap.csv",
                    index=False, encoding="utf-8-sig")
T_hub_port.to_csv(ASSETS / "L88_hub_port.csv",
                   index=False, encoding="utf-8-sig")
T_track.to_csv(ASSETS / "L88_track_one_polygon.csv",
                index=False, encoding="utf-8-sig")
T_compare.to_csv(ASSETS / "L88_compare_existing.csv",
                  index=False, encoding="utf-8-sig")
T_zone_coast.to_csv(ASSETS / "L88_zone_coast_split.csv",
                     index=False, encoding="utf-8-sig")
T_hypothesis.to_csv(ASSETS / "L88_hypothesis.csv",
                     index=False, encoding="utf-8-sig")
corr_df.round(2).to_csv(ASSETS / "L88_correlation_industrial_port.csv",
                         index=False, encoding="utf-8-sig")

# 工業専用地域ポリゴン (簡易表) を CSV へ (geometry除外)
ind_export = indzone.drop(columns="geometry").copy()
ind_export["area_km2"] = ind_export["area_km2"].round(4)
ind_export["storm_overlap_m2"] = ind_export["storm_overlap_m2"].round(0)
ind_export["storm_overlap_pct"] = ind_export["storm_overlap_pct"].round(2)
ind_export[["source_city", "YOTO_CD", "area_km2",
             "storm_overlap_m2", "storm_overlap_pct",
             "in_storm"]].to_csv(
    ASSETS / "L88_industrial_zones.csv",
    index=False, encoding="utf-8-sig")

print(f"  CSV 13 件 保存完了 ({time.time()-t11:.2f}s)", flush=True)


# =============================================================================
# 12. 図 8 枚
# =============================================================================
print("\n[12] 図作成", flush=True)
t12 = time.time()


def _format_admin_bg(ax):
    admin_diss.plot(ax=ax, facecolor="#fafafa",
                     edgecolor="#888", linewidth=0.4)


# --- Fig 1: 工業専用地域 県全域 主題図 (市町別色分け) ---
fig, ax = plt.subplots(1, 1, figsize=(15, 10))
_format_admin_bg(ax)

# 工業専用地域を市町で塗り分け
hub_color_map = {
    "福山市": "#cf222e", "呉市": "#0969da",
    "大竹市": "#9a6700", "尾道市": "#8250df",
    "広島市": "#1f883d",
}
for city in indzone["source_city"].unique():
    sub = indzone[indzone["source_city"] == city]
    color = hub_color_map.get(city, "#888")
    sub.plot(ax=ax, facecolor=color, alpha=0.85,
             edgecolor="white", linewidth=0.3)

# 4 拠点中心点を ★ で
for hub_name, info in INDUSTRIAL_HUBS.items():
    pt_4326 = gpd.GeoSeries([Point(info["lon"], info["lat"])],
                              crs="EPSG:4326").to_crs(TARGET_CRS).iloc[0]
    ax.scatter([pt_4326.x], [pt_4326.y], s=500, marker="*",
                facecolor=info["color"], edgecolor="black",
                linewidth=1.5, zorder=10)
    ax.annotate(f"  {hub_name}",
                 (pt_4326.x, pt_4326.y),
                 xytext=(15, -3), textcoords="offset points",
                 fontsize=11, fontweight="bold",
                 bbox=dict(boxstyle="round,pad=0.3",
                           facecolor=info["color"],
                           edgecolor="black", alpha=0.85),
                 color="white")

ax.set_aspect("equal")
ax.set_xticks([]); ax.set_yticks([])
ax.set_title(
    f"Fig 1: 広島県の工業専用地域 (YOTO_CD=12) — "
    f"市町別 {N_INDZONE:,} ポリゴン / {A_INDZONE:.1f} km²\n"
    f"4 産業拠点 (福山JFE / 呉海軍工廠 / 大竹コンビナート / 尾道造船) を ★ で表示",
    fontsize=12, pad=10
)
legend_elems = []
for city, color in hub_color_map.items():
    a = float(city_indarea.get(city, 0))
    legend_elems.append(Patch(facecolor=color, alpha=0.85,
                               label=f"{city} ({a:.2f} km²)"))
legend_elems.append(Patch(facecolor="#888", alpha=0.85,
                           label="その他市町"))
for hub_name, info in INDUSTRIAL_HUBS.items():
    legend_elems.append(Line2D(
        [0], [0], marker="*", color="white",
        markerfacecolor=info["color"], markeredgecolor="black",
        markersize=14, label=f"★ {hub_name}"
    ))
ax.legend(handles=legend_elems, loc="upper right", fontsize=9,
          framealpha=0.95)
plt.tight_layout()
plt.savefig(ASSETS / "L88_fig01_industrial_zones_overview.png", dpi=130)
plt.close("all")
print(f"  Fig 1 完成 ({time.time()-t12:.1f}s)", flush=True)


# --- Fig 2: 4 産業拠点 zoom (small multiples 2×2) ---
fig, axes = plt.subplots(2, 2, figsize=(16, 13))
axes = axes.flatten()
for ax, (hub_name, info) in zip(axes, INDUSTRIAL_HUBS.items()):
    bbox = HUB_BBOX[hub_name]
    # bbox to projected
    bbox_4326 = shapely.geometry.box(bbox[2], bbox[0], bbox[3], bbox[1])
    bbox_proj = (gpd.GeoSeries([bbox_4326], crs="EPSG:4326")
                 .to_crs(TARGET_CRS).iloc[0])
    xmin, ymin, xmax, ymax = bbox_proj.bounds

    # 行政界
    admin_diss.plot(ax=ax, facecolor="#fafafa",
                     edgecolor="#888", linewidth=0.4)

    # 工業専用地域
    sub_ind = indzone[indzone["source_city"] == info["city"]]
    sub_ind.plot(ax=ax, facecolor=info["color"], alpha=0.55,
                  edgecolor="white", linewidth=0.3)

    # 港湾施設 (港湾系のみ)
    sub_port = ports_gdf[
        (ports_gdf["kind"] == "港湾") &
        (ports_gdf.geometry.x >= xmin - 2000) &
        (ports_gdf.geometry.x <= xmax + 2000) &
        (ports_gdf.geometry.y >= ymin - 2000) &
        (ports_gdf.geometry.y <= ymax + 2000)
    ]
    if len(sub_port) > 0:
        sub_port.plot(ax=ax, color="#0a4d68", markersize=22,
                       marker="o", alpha=0.7, edgecolor="white",
                       linewidth=0.4)

    # 拠点中心 ★
    pt_4326 = gpd.GeoSeries([Point(info["lon"], info["lat"])],
                              crs="EPSG:4326").to_crs(TARGET_CRS).iloc[0]
    ax.scatter([pt_4326.x], [pt_4326.y], s=400, marker="*",
                facecolor="#ffd92f", edgecolor="black",
                linewidth=1.5, zorder=10)

    ax.set_xlim(xmin - 1500, xmax + 1500)
    ax.set_ylim(ymin - 1500, ymax + 1500)
    ax.set_aspect("equal")
    ax.set_xticks([]); ax.set_yticks([])
    a = float(city_indarea.get(info["city"], 0))
    ax.set_title(f"{hub_name} ({info['era']})\n"
                  f"{info['city']} 工業専用 {a:.2f} km²"
                  f" / 港湾施設 {len(sub_port)} 件",
                  fontsize=11)

plt.suptitle("Fig 2: 4 産業拠点 zoom — 工業専用地域 + 港湾施設 重ね合わせ",
             fontsize=14, y=1.00)
plt.tight_layout()
plt.savefig(ASSETS / "L88_fig02_hub_zoom.png", dpi=130)
plt.close("all")
print(f"  Fig 2 完成 ({time.time()-t12:.1f}s)", flush=True)


# --- Fig 3: 市町別 工業専用地域 面積ランキング ---
fig, axes = plt.subplots(1, 2, figsize=(17, 7))

# 左: 横棒
ax0 = axes[0]
top_show = T_city_ind_rank.head(10).iloc[::-1]
colors = ["#cf222e" if c in HUB_CITIES else "#888"
           for c in top_show["市町"]]
ax0.barh(top_show["市町"], top_show["工業専用km2"],
          color=colors, edgecolor="white")
for i, (c, a, s) in enumerate(zip(top_show["市町"],
                                    top_show["工業専用km2"],
                                    top_show["シェア%"])):
    ax0.text(a + 0.1, i, f"{a:.2f} km² ({s:.1f}%)",
              va="center", fontsize=10)
ax0.set_xlabel("工業専用地域 面積 (km²)")
ax0.set_title(f"市町別 工業専用地域 ランキング Top10 "
              f"(赤=4 拠点市町)\n4 拠点合計 = {hub_area:.2f} km² "
              f"({hub_share_pct:.1f}%)", fontsize=11)
ax0.grid(True, alpha=0.3, axis="x")

# 右: 沿岸 vs 内陸 13 用途地域 stacked
ax1 = axes[1]
zone_show = T_zone_coast.copy()
zone_show["合計"] = zone_show["沿岸"] + zone_show["内陸"]
zone_show = zone_show.sort_values("YOTO_CD")
y = np.arange(len(zone_show))
ax1.barh(y, zone_show["沿岸"], color="#0969da", label="沿岸市町",
          edgecolor="white")
ax1.barh(y, zone_show["内陸"], left=zone_show["沿岸"],
          color="#1f883d", label="内陸市町", edgecolor="white")
ax1.set_yticks(y)
ax1.set_yticklabels([f"{r['YOTO_CD']}={r['用途地域名']}"
                      for _, r in zone_show.iterrows()],
                     fontsize=9)
ax1.set_xlabel("面積 km²")
ax1.set_title("13 用途地域 沿岸 vs 内陸 比較\n"
              "(YOTO_CD=12 工業専用は 沿岸 100% に近い)",
              fontsize=11)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3, axis="x")

plt.suptitle("Fig 3: 第1章 — 工業専用地域の市町別ランキング + "
             "用途地域 沿岸/内陸 二極化", fontsize=13, y=1.00)
plt.tight_layout()
plt.savefig(ASSETS / "L88_fig03_city_ranking.png", dpi=130)
plt.close("all")
print(f"  Fig 3 完成 ({time.time()-t12:.1f}s)", flush=True)


# --- Fig 4: 工業専用 vs 高潮浸水想定 重ね合わせ ---
fig, axes = plt.subplots(1, 2, figsize=(17, 9))

ax0 = axes[0]
_format_admin_bg(ax0)
storm_gdf.plot(ax=ax0, facecolor="#0969da", alpha=0.30,
                edgecolor="none")
indzone.plot(ax=ax0, facecolor="#cf222e", alpha=0.85,
              edgecolor="white", linewidth=0.3)
# 拠点 ★
for hub_name, info in INDUSTRIAL_HUBS.items():
    pt_4326 = gpd.GeoSeries([Point(info["lon"], info["lat"])],
                              crs="EPSG:4326").to_crs(TARGET_CRS).iloc[0]
    ax0.scatter([pt_4326.x], [pt_4326.y], s=400, marker="*",
                 facecolor="#ffd92f", edgecolor="black",
                 linewidth=1.5, zorder=10)
ax0.set_aspect("equal")
ax0.set_xticks([]); ax0.set_yticks([])
ax0.set_title(f"工業専用地域 (赤) vs 高潮浸水想定 (青)\n"
              f"交差面積 {storm_overlap_total_km2:.2f} km² "
              f"({h4_value:.1f}%)", fontsize=11)
ax0.legend(handles=[
    Patch(facecolor="#cf222e", alpha=0.85,
           label=f"工業専用地域 ({A_INDZONE:.1f} km²)"),
    Patch(facecolor="#0969da", alpha=0.30,
           label=f"高潮浸水想定 max"),
    Line2D([0], [0], marker="*", color="white",
            markerfacecolor="#ffd92f",
            markeredgecolor="black", markersize=14,
            label="4 産業拠点中心"),
], loc="upper right", fontsize=10, framealpha=0.95)

# 右: 拠点別 高潮交差比率 棒
ax1 = axes[1]
hub_storm_show = T_hub_storm.copy()
y = np.arange(len(hub_storm_show))
colors_b = [INDUSTRIAL_HUBS[h]["color"] for h in hub_storm_show["拠点"]]
ax1.barh(y, hub_storm_show["高潮交差%"],
          color=colors_b, alpha=0.85, edgecolor="white")
for i, (h, p, a) in enumerate(zip(hub_storm_show["拠点"],
                                    hub_storm_show["高潮交差%"],
                                    hub_storm_show["高潮交差km2"])):
    ax1.text(p + 1, i, f"{p:.1f}% ({a:.2f} km²)",
              va="center", fontsize=10)
ax1.axvline(40, color="black", linestyle="--", linewidth=1.5,
             label="40% 閾値 (H4)")
ax1.set_yticks(y); ax1.set_yticklabels(hub_storm_show["拠点"])
ax1.set_xlabel("高潮浸水想定との交差比率 (%)")
ax1.set_title("拠点別 工業専用 vs 高潮 交差%", fontsize=11)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3, axis="x")

plt.suptitle("Fig 4: 第5章 — 産業と防災の交差 (工業専用地域 vs 高潮浸水想定)",
             fontsize=13, y=1.00)
plt.tight_layout()
plt.savefig(ASSETS / "L88_fig04_storm_overlay.png", dpi=130)
plt.close("all")
print(f"  Fig 4 完成 ({time.time()-t12:.1f}s)", flush=True)


# --- Fig 5: 港湾係留 vs 漁港外郭 ランキング比較 (H5) ---
fig, axes = plt.subplots(1, 2, figsize=(17, 7))

ax0 = axes[0]
top_h = T_harbor_top.head(10).iloc[::-1]
ax0.barh(top_h["港湾名"], top_h["係留施設数"],
          color="#0a4d68", edgecolor="white")
for i, (n, v) in enumerate(zip(top_h["港湾名"],
                                top_h["係留施設数"])):
    ax0.text(v + 1, i, f"{int(v)}", va="center", fontsize=10)
ax0.set_xlabel("係留施設件数")
ax0.set_title(f"港湾_係留施設 港湾名別 Top10\n"
              f"(産業港候補)", fontsize=11)
ax0.grid(True, alpha=0.3, axis="x")

ax1 = axes[1]
top_f = T_fish_top.head(10).iloc[::-1]
ax1.barh(top_f["港湾名"], top_f["外郭施設数"],
          color="#1f883d", edgecolor="white")
for i, (n, v) in enumerate(zip(top_f["港湾名"],
                                top_f["外郭施設数"])):
    ax1.text(v + 0.5, i, f"{int(v)}", va="center", fontsize=10)
ax1.set_xlabel("外郭施設件数")
ax1.set_title(f"漁港_外郭施設 港湾名別 Top10\n"
              f"(漁業港候補)\nTop5 重複: {N_OVERLAP} 件",
              fontsize=11)
ax1.grid(True, alpha=0.3, axis="x")

plt.suptitle("Fig 5: 第4章 — 産業港 vs 漁港 機能分化 (H5 検証)",
             fontsize=13, y=1.00)
plt.tight_layout()
plt.savefig(ASSETS / "L88_fig05_port_vs_fishport.png", dpi=130)
plt.close("all")
print(f"  Fig 5 完成 ({time.time()-t12:.1f}s)", flush=True)


# --- Fig 6: 工業専用面積 × 港湾係留件数 散布 (H3) ---
fig, ax = plt.subplots(1, 1, figsize=(13, 8))
sub = corr_df[mask_coastal].copy().sort_values("工業専用km2",
                                                ascending=False)
# プロット (沿岸市町)
for _, r in sub.iterrows():
    ax.scatter(r["工業専用km2"], r["係留件数"],
                s=180, alpha=0.85,
                c=hub_color_map.get(r["市町"], "#888"),
                edgecolor="white", linewidth=0.8)
    ax.annotate(r["市町"], (r["工業専用km2"], r["係留件数"]),
                 xytext=(8, 6), textcoords="offset points",
                 fontsize=10,
                 bbox=dict(boxstyle="round,pad=0.3",
                           facecolor="#fff3cd",
                           edgecolor="#cf222e", alpha=0.85))

# 内陸市町
sub_inland = corr_df[~mask_coastal].copy()
for _, r in sub_inland.iterrows():
    ax.scatter(r["工業専用km2"], r["係留件数"],
                s=70, alpha=0.5, c="#888",
                edgecolor="white", linewidth=0.5)

# 線形回帰線 (沿岸のみ)
if len(sub) >= 3:
    slope, intercept = np.polyfit(sub["工業専用km2"],
                                    sub["係留件数"], 1)
    xs = np.linspace(0, sub["工業専用km2"].max() * 1.05, 50)
    ys = slope * xs + intercept
    ax.plot(xs, ys, color="black", linestyle="--", linewidth=1.5,
             label=f"線形回帰 (沿岸) y = {slope:.1f}x + {intercept:.1f}")

ax.set_xlabel("工業専用地域 面積 (km²)")
ax.set_ylabel("港湾_係留施設 件数")
ax.set_title(f"Fig 6: 第3章 — 工業専用地域 vs 港湾係留 件数 相関 (H3)\n"
             f"沿岸 {len(sub)} 市町 Pearson r = {h3_value:.3f} "
             f"({'支持' if h3_ok else '反証'})",
             fontsize=12, pad=8)
ax.legend(fontsize=10, framealpha=0.95)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(ASSETS / "L88_fig06_correlation.png", dpi=130)
plt.close("all")
print(f"  Fig 6 完成 ({time.time()-t12:.1f}s)", flush=True)


# --- Fig 7: 13 用途地域 構成円グラフ + 工業専用シェア時系列風 ---
fig, axes = plt.subplots(1, 2, figsize=(17, 7))

# 左: 13 用途地域 円グラフ (全県)
ax0 = axes[0]
zone_show = T_zone_type.sort_values("面積km2", ascending=False)
labels = [f"{r['YOTO_CD']} {r['用途地域名']} "
          f"({r['シェア%']:.1f}%)"
          for _, r in zone_show.iterrows()]
colors_pie = []
for _, r in zone_show.iterrows():
    cd = int(r["YOTO_CD"])
    if cd == 12:
        colors_pie.append("#cf222e")
    elif cd in (10, 11):
        colors_pie.append("#9a6700")
    elif cd in (8, 9):
        colors_pie.append("#fdae61")
    elif cd == 13:
        colors_pie.append("#ffd92f")
    else:
        colors_pie.append("#1f883d")
ax0.pie(zone_show["面積km2"], labels=labels, colors=colors_pie,
         autopct=lambda p: f"{p:.1f}%" if p > 5 else "",
         pctdistance=0.78, textprops={"fontsize": 9})
ax0.set_title(f"13 用途地域 全県 面積構成\n"
              f"(工業専用赤={zone_show[zone_show['YOTO_CD']==12]['シェア%'].iloc[0]:.1f}%)",
              fontsize=11)

# 右: 4 拠点別 工業専用面積 縦棒 (歴史軸ストーリー)
ax1 = axes[1]
hubs = list(INDUSTRIAL_HUBS.keys())
hub_areas = [float(city_indarea.get(INDUSTRIAL_HUBS[h]["city"], 0))
              for h in hubs]
colors_h = [INDUSTRIAL_HUBS[h]["color"] for h in hubs]
bars = ax1.bar(hubs, hub_areas, color=colors_h, alpha=0.85,
                edgecolor="white")
for b, a, h in zip(bars, hub_areas, hubs):
    ax1.text(b.get_x() + b.get_width()/2, b.get_height() + 0.1,
              f"{a:.2f} km²\n{INDUSTRIAL_HUBS[h]['era']}",
              ha="center", fontsize=9)
ax1.set_ylabel("工業専用地域 面積 (km²)")
ax1.set_title("4 産業拠点 — 歴史軸 vs 現在面積", fontsize=11)
ax1.grid(True, alpha=0.3, axis="y")
plt.setp(ax1.get_xticklabels(), rotation=15)

plt.suptitle("Fig 7: 第6章 — 用途地域 全県構成 と 4 拠点 歴史軸",
             fontsize=13, y=1.00)
plt.tight_layout()
plt.savefig(ASSETS / "L88_fig07_pie_history.png", dpi=130)
plt.close("all")
print(f"  Fig 7 完成 ({time.time()-t12:.1f}s)", flush=True)


# --- Fig 8: 仮説検証総合ダッシュボード ---
fig, axes = plt.subplots(2, 2, figsize=(16, 11))

# (1) 5 仮説検証 横棒
ax0 = axes[0, 0]
hyp_short = ["H1 臨海重工業集中", "H2 福山JFE圧倒",
              "H3 港湾相関", "H4 高潮交差", "H5 機能分化"]
hyp_actual = [h1_value, h2_value, h3_value, h4_value, h5_value]
hyp_threshold_disp = ["≥50%", "≥20%", "≥0.5", "≥40%", "≤1"]
hyp_pass = [h1_ok, h2_ok, h3_ok, h4_ok, h5_ok]
y = np.arange(len(hyp_short))
colors = ["#1f883d" if p else "#cf222e" for p in hyp_pass]
# 値のスケールが異なるので相対化 (各仮説 0-1 範囲)
display_vals = [
    h1_value / 50,    # 50% を 1
    h2_value / 20,    # 20% を 1
    h3_value / 0.5,   # 0.5 を 1
    h4_value / 40,    # 40% を 1
    1 - h5_value / 5, # 5 を 0 (反転)
]
ax0.barh(y, display_vals, color=colors, alpha=0.7, edgecolor="white")
ax0.axvline(1.0, color="black", linestyle=":", linewidth=1.5,
             label="閾値ライン")
for i, (a, t, p, name) in enumerate(zip(hyp_actual, hyp_threshold_disp,
                                          hyp_pass, hyp_short)):
    mark = "支持" if p else "反証"
    raw = (f"{a:.1f}%" if i in (0, 1, 3) else
           f"r={a:.2f}" if i == 2 else
           f"{int(a)}件")
    ax0.text(display_vals[i] + 0.05, i,
              f"{raw} [{t}] → {mark}", va="center", fontsize=9)
ax0.set_yticks(y); ax0.set_yticklabels(hyp_short)
ax0.set_xlabel("実測 / 閾値 (相対値, 1.0 = 閾値)")
ax0.set_title(f"5 仮説 検証結果 (支持 {n_supported}/5)", fontsize=11)
ax0.legend(fontsize=9)
ax0.grid(True, alpha=0.3, axis="x")

# (2) 4 拠点別 工業専用 面積
ax1 = axes[0, 1]
hubs = list(INDUSTRIAL_HUBS.keys())
hub_a = [float(city_indarea.get(INDUSTRIAL_HUBS[h]["city"], 0))
          for h in hubs]
hub_storm_a = [r["高潮交差km2"] for _, r in T_hub_storm.iterrows()]
x = np.arange(len(hubs))
ax1.bar(x - 0.2, hub_a, width=0.4, color="#cf222e",
         label="工業専用 km²", edgecolor="white")
ax1.bar(x + 0.2, hub_storm_a, width=0.4, color="#0969da",
         label="高潮交差 km²", edgecolor="white")
ax1.set_xticks(x); ax1.set_xticklabels(hubs, rotation=15)
ax1.set_ylabel("面積 km²")
ax1.set_title("4 拠点別 工業専用 vs 高潮交差", fontsize=11)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3, axis="y")

# (3) 港湾係留 Top5 vs 漁港外郭 Top5
ax2 = axes[1, 0]
labels_h = T_harbor_top.head(5)["港湾名"].tolist()
vals_h = T_harbor_top.head(5)["係留施設数"].tolist()
labels_f = T_fish_top.head(5)["港湾名"].tolist()
vals_f = T_fish_top.head(5)["外郭施設数"].tolist()
y_h = np.arange(5)
ax2.barh(y_h - 0.2, vals_h, height=0.4, color="#0a4d68",
          label="港湾_係留 Top5", edgecolor="white")
ax2.barh(y_h + 0.2, vals_f, height=0.4, color="#1f883d",
          label="漁港_外郭 Top5", edgecolor="white")
ax2.set_yticks(y_h)
ax2.set_yticklabels([f"{i+1}位" for i in range(5)])
for i, (lh, vh) in enumerate(zip(labels_h, vals_h)):
    ax2.text(vh + 1, i - 0.2, f" {lh} ({vh})",
              va="center", fontsize=8)
for i, (lf, vf) in enumerate(zip(labels_f, vals_f)):
    ax2.text(vf + 1, i + 0.2, f" {lf} ({vf})",
              va="center", fontsize=8)
ax2.set_xlabel("施設件数")
ax2.set_title(f"H5 — 産業港 vs 漁港 機能分化 (Top5 重複 {N_OVERLAP})",
              fontsize=11)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3, axis="x")

# (4) 沿岸 vs 内陸 13 用途地域 比較
ax3 = axes[1, 1]
zs = T_zone_coast.sort_values("YOTO_CD")
y = np.arange(len(zs))
ax3.barh(y, zs["沿岸"], color="#0969da",
          label="沿岸市町", edgecolor="white")
ax3.barh(y, zs["内陸"], left=zs["沿岸"],
          color="#1f883d", label="内陸市町", edgecolor="white")
ax3.set_yticks(y)
ax3.set_yticklabels([f"{r['YOTO_CD']}={r['用途地域名']}"
                      for _, r in zs.iterrows()], fontsize=8)
ax3.set_xlabel("面積 km²")
ax3.set_title("13 用途地域 沿岸/内陸 比較", fontsize=11)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3, axis="x")

plt.suptitle("Fig 8: 仮説検証総合ダッシュボード — "
             "産業地理 4 視点", fontsize=13, y=1.00)
plt.tight_layout()
plt.savefig(ASSETS / "L88_fig08_dashboard.png", dpi=130)
plt.close("all")
print(f"  Fig 8 完成 ({time.time()-t12:.1f}s)", flush=True)

print(f"  全 8 図完成 ({time.time()-t12:.1f}s)", flush=True)


# =============================================================================
# 13. HTML レンダー
# =============================================================================
print("\n[13] HTML レンダー", flush=True)
t13 = time.time()


def df_to_html_table(df, max_rows=None):
    if max_rows:
        df = df.head(max_rows)
    cols = "".join(f"<th>{escape(str(c))}</th>" for c in df.columns)
    rows_html = ""
    for _, row in df.iterrows():
        cells = ""
        for v in row:
            if isinstance(v, float):
                if pd.isna(v):
                    s = ""
                elif abs(v) >= 1000:
                    s = f"{v:,.0f}"
                else:
                    s = f"{v:.2f}".rstrip("0").rstrip(".")
                    if s == "" or s == "-":
                        s = "0"
            else:
                s = "" if pd.isna(v) else str(v)
            cells += f"<td>{escape(s)}</td>"
        rows_html += f"<tr>{cells}</tr>"
    return f"<table><tr>{cols}</tr>{rows_html}</table>"


sections = []


# --- 1. 学習目標と問い ---
intro_html = f"""
<p>1903 年、 呉に東洋一の海軍工廠が開設された。 戦艦大和を生んだその巨大な
鉄の街は、 1945 年の終戦で解体され、 戦後は IHI 呉造船・神戸製鋼として
民需に転換した。 その北東 70 km の福山では、 1965 年に川崎製鉄
(現 JFE スチール) が高度成長の象徴として鉄を吹き始めた。 さらに南西、
大竹の狭い平地には戦前から三井化学のコンビナートが煙突を並べ、 尾道では
向島の入江で常石・尾道造船が世界中の船を造っている。</p>

<p><b>広島県の産業地理は、 戦前 - 戦後の連続性を持って瀬戸内海岸線に
書き込まれている</b>。 海軍工廠跡は造船所に、 戦前化学コンビナートは
石油精製に、 戦後高度成長の埋立地は世界級鉄鋼コンビナートに。 「臨海重
工業」 という広島県の<b>産業 DNA</b> は、 都市計画法の用途地域指定 (工業専用
地域 YOTO_CD=12) として、 そして港湾施設として、 オープンデータの中で
今もそのまま読み取れる。</p>

<p>本記事は、 DoBoX の <b>用途地域 21 dataset (#794-946, #932)</b> +
<b>港湾施設 (#1250-1252) {N_PORT_HARBOR:,} 件</b> +
<b>漁港施設 (#1254-1255) {N_PORT_FISH:,} 件</b> +
<b>高潮浸水想定 max (#45) {len(storm_gdf)} ランク</b> +
<b>橋梁 (#11) {len(pd.read_csv(BRIDGE_CSV, encoding='utf-8-sig')):,} 件</b>
+ <b>海岸保全 (#1253)</b> + <b>L15 行政界</b> を統合し、 「<b>広島県の産業物語 —
戦前 - 戦後の産業立地が織りなす都市の物語</b>」 を 6 章の物語として読み
解く<b>テーマ統合 (M系) 記事</b>である。</p>

<h3>独自用語の定義 (この記事だけで使う)</h3>
<ul class="kv">
  <li><b>4 産業拠点</b>: 本記事独自分類。 広島県の臨海重工業を 4 つの
    歴史的拠点で代表させる: 福山JFE (戦後高度成長 1965)、 大竹コンビナート
    (戦前化学 1920s-)、 呉海軍工廠 (戦前軍需 1903)、 尾道造船 (戦前造船)。
    厳密な行政区域分けではなく、 産業史の物語的拠点。</li>
  <li><b>工業専用地域</b>: 都市計画法・建築基準法上の用途地域 13 種類のうち
    最も工業特化が強い類型 (YOTO_CD=12)。 住宅・学校・病院などの建築が
    禁止され、 大規模工場・コンビナートのみが認められる。 広島県では
    {N_INDZONE:,} ポリゴン / {A_INDZONE:.2f} km² が指定されている。</li>
  <li><b>コンビナート</b>: 石油・化学・鉄鋼などの装置産業が、 原料・中間品・
    エネルギーを物理パイプラインで結合した一体型工場群。 大竹市 (三井化学
    系)、 福山市 (JFE 系) が代表例。</li>
  <li><b>海軍工廠</b>: 戦前の旧帝国海軍が艦艇・兵器を直接生産・整備した
    国営工場。 呉海軍工廠 (1903-1945) は東洋最大級の海軍生産拠点で、
    戦艦大和を生んだ。 戦後は IHI・神戸製鋼などに継承され、 民需造船・
    鉄鋼の集積地となった。</li>
  <li><b>産業港</b>: 鉄鋼・化学・造船など重工業の物流を担う港。 港湾施設の
    係留施設件数 (物揚場・浮さん橋・岸壁) で代理する。 広島県では
    広島港・尾道糸崎港・蒲刈港が上位、 福山港 ({port_moor[port_moor['港湾名']=='福山港'].shape[0]} 件)
    が産業港として続く。</li>
  <li><b>漁業港</b>: 漁業活動を支える港。 漁港施設 (倉橋・豊島・沖浦など)
    で代理する。 「産業港」 とは機能が異なり、 水産庁所管の漁港法に基づき
    指定される。</li>
  <li><b>臨海重工業帯</b>: 本記事独自概念。 沿岸 4 産業拠点市町 (福山・呉・
    大竹・尾道) を統合した、 広島県の主要工業地帯。 H1 で量化検証する。</li>
  <li><b>沿岸工業 vs 高潮交差</b>: 工業専用地域ポリゴンと高潮浸水想定エリア
    の地理的重なり。 「工業生産 (経済的価値) と 気候災害リスク (浸水) の
    同居」 という瀬戸内臨海工業地の構造的特徴を量化する独自概念。</li>
  <li><b>機能分化</b>: 港湾係留施設 Top5 と漁港外郭施設 Top5 の市町重複が
    少ないこと。 広島県の沿岸は<b>重工業港 (係留集積)</b> と<b>漁業港
    (外郭集積)</b> に空間的に二極化していることを示す。 H5 で検証。</li>
  <li><b>産業 DNA</b>: 本記事独自概念。 戦前 - 戦後を貫く広島県の臨海
    重工業の連続性。 海軍工廠 → 造船・鉄鋼、 戦前化学 → 石油精製、
    戦後高度成長 → 鉄鋼コンビナート、 という時間連続性を、 用途地域
    指定が現代まで保存している。</li>
</ul>

<h3>研究の問い (主 RQ)</h3>
<p><b>広島県の産業物語 — 戦前 - 戦後の産業立地が織りなす都市の物語とは
何か?  戦前 (旧海軍工廠呉) ・戦前化学 (大竹コンビナート) ・戦後高度成長
(福山JFE) ・尾道造船 という 4 産業拠点は、 工業専用地域・港湾施設・
浸水リスクの量と質でどう描けるか?  「臨海重工業」 という広島県の産業 DNA
はデータでどう浮かぶか?</b></p>

<h3>仮説 (H1〜H5)</h3>
<ul>
  <li><b>H1 (臨海重工業集中仮説)</b>: 工業専用地域 (YOTO_CD=12) のうち、
    沿岸 4 産業拠点市町 (福山・呉・大竹・尾道) の合計面積比率は
    <b>50% 以上</b>。 広島県の重工業帯は瀬戸内沿岸に偏在する。</li>
  <li><b>H2 (福山JFE圧倒仮説)</b>: 工業専用地域の市町別最大は<b>福山市</b>
    であり、 県全体の<b>20% 以上</b>を占める。 戦後高度成長期 (1965 年
    福山製鉄所開設) が広島県工業地理を再編した量的証拠。</li>
  <li><b>H3 (港湾規模 vs 工業面積 相関仮説)</b>: 各市町の工業専用地域面積
    と港湾係留施設件数は<b>正の相関</b> (Pearson r ≥ 0.5)。
    大規模工業地帯 = 港湾物流集積、という臨海工業立地原理。</li>
  <li><b>H4 (沿岸工業 vs 高潮交差仮説)</b>: 工業専用地域ポリゴンの
    <b>高潮浸水想定エリアとの交差面積比は 40% 以上</b>。 沿岸埋立地の
    工業地は気候災害の最前線。</li>
  <li><b>H5 (産業港 vs 漁港 二極化仮説)</b>: 港湾係留施設 Top5 と漁港外郭
    施設 Top5 の<b>市町重複は 1 件以下</b>。 「重工業港 vs 漁業港」 の
    機能分化が空間的に成立。</li>
</ul>

<h3>到達点</h3>
<ul>
  <li>4 産業拠点 (福山JFE・呉海軍工廠・大竹コンビナート・尾道造船) という
    本記事独自の<b>歴史軸</b>で、 工業専用地域 + 港湾施設を統合 GDF として
    1 枚の地図で読める。</li>
  <li>沿岸 4 拠点市町の工業専用シェア ({hub_share_pct:.1f}%) と福山市の
    シェア ({fukuyama_share_pct:.1f}%) を量化し、 広島県の臨海重工業
    集中構造を示す。</li>
  <li>4 拠点 zoom (small multiples) で、 産業地と港湾施設の地理的
    結合を可視化する。</li>
  <li>工業専用地域 vs 高潮浸水想定の重ね合わせで、 「経済的価値とリスク」
    の交差 ({h4_value:.1f}% が高潮想定圏内) を量化する。</li>
  <li>港湾係留 vs 漁港外郭 ランキングで、 産業港と漁業港の<b>機能分化</b>
    (Top5 重複 {N_OVERLAP} 件) を検証する。</li>
</ul>

<div class="note">
<b>L17 / L21 / L32-L34 / L52-L55 / L66 との明確な切り分け</b>: 既存の
L 系記事はそれぞれ「13 用途地域並列」 「宅地開発単独」 「港湾外郭・係留・
臨港交通の単独深掘り」 「橋梁単独」 を扱った。 L88 はこれらと完全に独立
した切り口で、 <b>「工業専用 1 用途 × 港湾相関 × 高潮交差 × 4 拠点歴史
軸」</b> という新しい統合視点で物語を組む。 同じデータが切り口を変えれば
全く異なる物語を語る — それが M系 (テーマ統合) の本質である。
</div>
"""
sections.append(("学習目標と問い", intro_html))


# --- 2. 使用データ ---
data_html = f"""
<p>本記事は産業物語を語るため、 <b>主軸 1 系統 (用途地域 21 dataset) +
従属 7 系統 = 計 8 系統</b>を組み合わせて使う。 主軸の用途地域から
工業専用 (YOTO_CD=12) のみを抽出し、 港湾・漁港・橋梁・浸水想定・海岸
保全を従属参照する。</p>

<table>
  <tr><th>論題</th><th>dataset</th><th>件数</th><th>用途</th></tr>
  <tr><td><b>都市計画区域情報_用途地域 21 件</b></td>
      <td><a href="https://hiroshima-dobox.jp/datasets/932" target="_blank">DoBoX #932 (県全域)</a> +
          21 市町別 (#794-946)</td>
      <td>{len(zone):,} ポリゴン</td>
      <td>本記事の<b>主軸</b>。 工業専用 {N_INDZONE:,} 件を抽出</td></tr>
  <tr><td>港湾施設_外郭施設</td>
      <td><a href="https://hiroshima-dobox.jp/datasets/1250" target="_blank">DoBoX #1250</a></td>
      <td>{len(port_outer):,}</td>
      <td>防波堤・防砂堤 (港湾)</td></tr>
  <tr><td>港湾施設_係留施設</td>
      <td><a href="https://hiroshima-dobox.jp/datasets/1251" target="_blank">DoBoX #1251</a></td>
      <td>{len(port_moor):,}</td>
      <td>物揚場・浮さん橋・岸壁 (産業港代理指標, H3)</td></tr>
  <tr><td>港湾施設_臨港交通施設</td>
      <td><a href="https://hiroshima-dobox.jp/datasets/1252" target="_blank">DoBoX #1252</a></td>
      <td>{len(port_traf):,}</td>
      <td>道路・駐車場 (港湾の陸側)</td></tr>
  <tr><td>漁港施設_外郭施設</td>
      <td><a href="https://hiroshima-dobox.jp/datasets/1254" target="_blank">DoBoX #1254</a></td>
      <td>{len(fish_outer):,}</td>
      <td>漁港の防波堤等 (漁業港代理指標, H5)</td></tr>
  <tr><td>漁港施設_係留施設</td>
      <td><a href="https://hiroshima-dobox.jp/datasets/1255" target="_blank">DoBoX #1255</a></td>
      <td>{len(fish_moor):,}</td>
      <td>漁港の物揚場等</td></tr>
  <tr><td>高潮浸水想定区域 max</td>
      <td><a href="https://hiroshima-dobox.jp/datasets/45" target="_blank">DoBoX #45</a></td>
      <td>{len(storm_gdf)} ランク</td>
      <td>第5章 産業 vs 高潮 交差判定 (H4)</td></tr>
  <tr><td>橋梁基本情報</td>
      <td><a href="https://hiroshima-dobox.jp/datasets/11" target="_blank">DoBoX #11</a></td>
      <td>4,203 件</td>
      <td>背景参照 (拠点周辺の橋梁集積)</td></tr>
  <tr><td>(従属) 海岸保全施設基本情報</td>
      <td><a href="https://hiroshima-dobox.jp/datasets/1253" target="_blank">DoBoX #1253</a></td>
      <td>1,646 件</td>
      <td>沿岸線推定 (背景)</td></tr>
  <tr><td>(派生) L15 行政界 21 市町</td>
      <td>L15 cache (admin_*.zip)</td>
      <td>20</td>
      <td>背景マップ + dissolve</td></tr>
</table>

<h3>派生指標 (本記事独自)</h3>
<ul>
  <li><b>area_km2</b>: 工業専用地域ポリゴン面積 (EPSG:6671 で正確計算)</li>
  <li><b>storm_overlap_m2 / pct</b>: 高潮浸水想定との交差面積 / 比率</li>
  <li><b>in_storm</b>: 高潮交差フラグ (0/1)</li>
  <li><b>HUB_BBOX</b>: 4 産業拠点 ±5km bbox (zoom 用)</li>
  <li><b>r_pearson</b>: 沿岸市町 工業専用面積 vs 港湾係留件数 相関係数</li>
</ul>

<h3>1 件追跡 (要件 K — Before/After)</h3>
{df_to_html_table(T_track)}

<h3>この表から読み取れること</h3>
<ul>
  <li>工業専用地域の最大ポリゴン (1 件) が、 元 GeoJSON → 投影 → 抽出 →
    重心 → 高潮交差判定 という<b>7 段階</b>を経て統合される過程を追跡。</li>
  <li>1 ポリゴンから 7 つの派生情報 (面積・座標・高潮交差量・比率・
    フラグ) を導出することで、 「同じデータでも切り口を変えれば多角的
    物語が語れる」 ことを示す。</li>
</ul>
"""
sections.append(("使用データ", data_html))


# --- 3. ダウンロード ---
dl_html = f"""
<h3>DoBoX 元データ (直リンク)</h3>
<div>
  <a href="https://hiroshima-dobox.jp/datasets/932" target="_blank"
     style="display:inline-block;padding:8px 14px;margin:4px;
     background:#cf222e;color:white;border-radius:6px;text-decoration:none;
     font-weight:bold;">▶ #932 用途地域 県全域</a>
  <a href="https://hiroshima-dobox.jp/datasets/837" target="_blank"
     style="display:inline-block;padding:6px 10px;margin:3px;
     background:#9a6700;color:white;border-radius:4px;text-decoration:none;
     font-size:13px;">#837 福山市</a>
  <a href="https://hiroshima-dobox.jp/datasets/804" target="_blank"
     style="display:inline-block;padding:6px 10px;margin:3px;
     background:#0969da;color:white;border-radius:4px;text-decoration:none;
     font-size:13px;">#804 呉市</a>
  <a href="https://hiroshima-dobox.jp/datasets/867" target="_blank"
     style="display:inline-block;padding:6px 10px;margin:3px;
     background:#1f883d;color:white;border-radius:4px;text-decoration:none;
     font-size:13px;">#867 大竹市</a>
  <a href="https://hiroshima-dobox.jp/datasets/831" target="_blank"
     style="display:inline-block;padding:6px 10px;margin:3px;
     background:#8250df;color:white;border-radius:4px;text-decoration:none;
     font-size:13px;">#831 尾道市</a>
  <a href="https://hiroshima-dobox.jp/datasets/45" target="_blank"
     style="display:inline-block;padding:6px 10px;margin:3px;
     background:#bf3989;color:white;border-radius:4px;text-decoration:none;
     font-size:13px;">#45 高潮浸水想定 max</a>
  <a href="https://hiroshima-dobox.jp/datasets/11" target="_blank"
     style="display:inline-block;padding:6px 10px;margin:3px;
     background:#0a4d68;color:white;border-radius:4px;text-decoration:none;
     font-size:13px;">#11 橋梁</a>
  <a href="https://hiroshima-dobox.jp/datasets/1251" target="_blank"
     style="display:inline-block;padding:6px 10px;margin:3px;
     background:#0a4d68;color:white;border-radius:4px;text-decoration:none;
     font-size:13px;">#1251 港湾係留</a>
  <a href="https://hiroshima-dobox.jp/datasets/1253" target="_blank"
     style="display:inline-block;padding:6px 10px;margin:3px;
     background:#666;color:white;border-radius:4px;text-decoration:none;
     font-size:13px;">#1253 海岸保全</a>
</div>

<h3>本記事が生成した中間データ (再現用 — 直リンク)</h3>
<ul>
  <li><a href="assets/L88_industrial_zones.csv" download>L88_industrial_zones.csv</a>
    — 工業専用地域 {N_INDZONE:,} ポリゴン + 高潮交差量</li>
  <li><a href="assets/L88_city_industrial_rank.csv" download>L88_city_industrial_rank.csv</a>
    — 市町別 工業専用地域 ランキング</li>
  <li><a href="assets/L88_zone_type_summary.csv" download>L88_zone_type_summary.csv</a>
    — 13 用途地域 全県集計</li>
  <li><a href="assets/L88_harbor_mooring_top.csv" download>L88_harbor_mooring_top.csv</a>
    — 港湾_係留 港湾名別 ランキング</li>
  <li><a href="assets/L88_fishport_outer_top.csv" download>L88_fishport_outer_top.csv</a>
    — 漁港_外郭 港湾名別 ランキング</li>
  <li><a href="assets/L88_hub_summary.csv" download>L88_hub_summary.csv</a>
    — 4 拠点サマリー</li>
  <li><a href="assets/L88_hub_storm_overlap.csv" download>L88_hub_storm_overlap.csv</a>
    — 4 拠点別 高潮交差量</li>
  <li><a href="assets/L88_hub_port.csv" download>L88_hub_port.csv</a>
    — 4 拠点別 港湾施設件数</li>
  <li><a href="assets/L88_correlation_industrial_port.csv" download>L88_correlation_industrial_port.csv</a>
    — 工業専用 × 港湾係留 散布データ (H3)</li>
  <li><a href="assets/L88_zone_coast_split.csv" download>L88_zone_coast_split.csv</a>
    — 13 用途地域 沿岸/内陸 比較</li>
  <li><a href="assets/L88_track_one_polygon.csv" download>L88_track_one_polygon.csv</a>
    — 1 件追跡 (Before/After)</li>
  <li><a href="assets/L88_compare_existing.csv" download>L88_compare_existing.csv</a>
    — 既存 L 系記事との重複ゼロ宣言</li>
  <li><a href="assets/L88_hypothesis.csv" download>L88_hypothesis.csv</a>
    — 仮説検証 5 行</li>
</ul>

<h3>図 8 枚 (PNG, 直 DL 可)</h3>
<ul>
  <li><a href="assets/L88_fig01_industrial_zones_overview.png" download>fig01 工業専用 県全域</a></li>
  <li><a href="assets/L88_fig02_hub_zoom.png" download>fig02 4拠点 zoom (2×2)</a></li>
  <li><a href="assets/L88_fig03_city_ranking.png" download>fig03 市町ランキング+用途沿岸内陸</a></li>
  <li><a href="assets/L88_fig04_storm_overlay.png" download>fig04 工業 vs 高潮</a></li>
  <li><a href="assets/L88_fig05_port_vs_fishport.png" download>fig05 港湾 vs 漁港</a></li>
  <li><a href="assets/L88_fig06_correlation.png" download>fig06 工業 vs 港湾 相関</a></li>
  <li><a href="assets/L88_fig07_pie_history.png" download>fig07 13用途円+4拠点歴史</a></li>
  <li><a href="assets/L88_fig08_dashboard.png" download>fig08 仮説検証ダッシュボード</a></li>
</ul>

<h3>再現スクリプト</h3>
<ul>
  <li><a href="L88_M4_industry_story.py" download>L88_M4_industry_story.py</a>
    — 本記事の全処理</li>
</ul>

<h3>実行</h3>
<pre><code>cd "2026 DoBoX 教材"
py -X utf8 lessons/L88_M4_industry_story.py</code></pre>
<p class="tnote">用途地域 21 zip は L17 cache (<code>data/extras/L17_use_zones/</code>)
を再利用、 港湾は L32-L34 cache、 漁港は同じくキャッシュ、 高潮は L44 cache、
橋梁・海岸保全・行政界も既存 cache を参照。 追加 DL なしで即実行可能。</p>
"""
sections.append(("ダウンロード", dl_html))


# --- 4. 第1章: 工業専用地域の地理 ---
ch1_purpose = f"""
<p><b>狙い</b>: 広島県の重工業帯はどこか?  まず<b>工業専用地域 (YOTO_CD=12)</b>
{N_INDZONE:,} ポリゴン / {A_INDZONE:.2f} km² の<b>市町別分布</b>を見て、
沿岸 4 拠点市町 (福山・呉・大竹・尾道) の集中度を量化する。 これが
広島県産業地理の<b>地理的基盤</b>である。</p>

<p><b>手法の要点</b>: L17 で読み込んだ用途地域 21 市町 GeoJSON のうち
<b>YOTO_CD=12 のポリゴンのみフィルタ</b>し、 EPSG:6671 平面直角座標で
正確な面積計算 (m²) を行う。 さらに市町別に <code>groupby</code> で集計し、
ランキング Top10 を作成。 沿岸 4 拠点市町の合計面積比 (H1) を検証する。</p>

<p><b>用途地域とは (リテラシ説明)</b>: 都市計画法・建築基準法に基づき、
都市計画区域内の土地を <b>13 種類</b>に分類する制度。 住居系 (1-7) /
商業系 (8-9) / 工業系 (10-12) / 田園系 (13)。 工業系のうち最も特化が
強いのが<b>工業専用地域 (YOTO_CD=12)</b> で、 住宅・学校・病院などの
建築が禁止され、 大規模工場のみが認められる。</p>
"""

ch1_code = code('''
import zipfile, io
import geopandas as gpd
import pandas as pd

CITY_DEFS = [
    (794, "広島市", True), (804, "呉市", True), (837, "福山市", True),
    (867, "大竹市", True), (831, "尾道市", True), # ... 21 市町
]

# 21 市町統合
frames = []
for dsid, name, coastal in CITY_DEFS:
    z = f"data/extras/L17_use_zones/use_zone_{dsid}_{name}.zip"
    with zipfile.ZipFile(z) as zf:
        gjs = [n for n in zf.namelist() if n.endswith(".geojson")][0]
        with zf.open(gjs) as f:
            g = gpd.read_file(io.BytesIO(f.read()))
    g["source_city"] = name
    g["coastal"] = coastal
    frames.append(g)

zone = gpd.GeoDataFrame(pd.concat(frames, ignore_index=True),
                        geometry="geometry", crs=frames[0].crs)
zone = zone.to_crs("EPSG:6671")  # 平面直角第 III 系 (m)

# 工業専用地域のみ抽出 (YOTO_CD=12)
indzone = zone[zone["YOTO_CD"] == 12].copy()
indzone["area_km2"] = indzone.geometry.area / 1e6

# 市町別ランキング
city_indarea = (indzone.groupby("source_city")["area_km2"]
                .sum().sort_values(ascending=False))
print(city_indarea.head(10))
# 例:
# 福山市     6.57 km²
# 広島市     7.42 km²
# 呉市       3.66 km²
# ...

# H1 検証 — 4 拠点合計シェア
HUB_CITIES = ["福山市", "呉市", "大竹市", "尾道市"]
hub_total = float(city_indarea.reindex(HUB_CITIES).fillna(0).sum())
hub_share = hub_total / city_indarea.sum() * 100
print(f"4 拠点合計 = {hub_total:.2f} km² ({hub_share:.1f}%)")
''')

ch1_html = (
    "<h3>狙い・手法</h3>" + ch1_purpose +
    "<h3>実装</h3>" + ch1_code +
    "<h3>結果図 — 工業専用地域 県全域 主題図</h3>" +
    "<p><b>なぜこの図か</b>: 「広島県の重工業帯はどこか」 を一目で伝えるには、 "
    "<b>工業専用地域ポリゴン</b>を市町別に色分けした<b>主題図 (choropleth)</b>"
    "が最強。 4 拠点中心 ★ も併記し、 産業史と地理を結ぶ。</p>" +
    figure("assets/L88_fig01_industrial_zones_overview.png",
           f"Fig 1: 工業専用地域 (YOTO_CD=12) 県全域 — "
           f"{N_INDZONE:,} ポリゴン / {A_INDZONE:.2f} km²") +
    f"""<h3>この図から読み取れること</h3>
<ul>
  <li>工業専用地域は<b>瀬戸内沿岸</b>に集中する。 内陸の三次・庄原・
    安芸高田などにはほぼ無い (内陸は工業地域 YOTO_CD=11 までで止まる)。</li>
  <li>赤 (福山市): 福山市東部の臨海埋立地に大規模ポリゴンが連続。
    JFE スチール福山地区を中心とする戦後高度成長コンビナート。</li>
  <li>青 (呉市): 呉湾岸に複数の中規模ポリゴン。 旧海軍工廠跡地が
    現代の工業専用地域として継承されている。</li>
  <li>茶 (大竹市): 山口県境付近の細長いコンビナート帯。 三井化学・
    三菱レイヨンなどの戦前化学産業の継承。</li>
  <li>紫 (尾道市): 向島・三原境付近の造船クラスタ。
    常石造船・尾道造船の集積。</li>
  <li>緑 (広島市): 南区・佐伯区などに分散。 マツダ・三菱重工広島
    製作所などの戦後企業が立地。</li>
  <li>4 拠点合計 = <b>{hub_area:.2f} km² ({hub_share_pct:.1f}%)</b>。
    H1 (≥50% 仮説) は<b>{"支持" if h1_ok else "反証"}</b>。</li>
</ul>""" +
    f"""<h3>市町別 工業専用地域 ランキング Top10</h3>
{df_to_html_table(T_city_ind_rank)}

<h3>この表から読み取れること</h3>
<ul>
  <li>1 位 = <b>{T_city_ind_rank.iloc[0]['市町']}</b>
    ({T_city_ind_rank.iloc[0]['工業専用km2']} km²,
    シェア {T_city_ind_rank.iloc[0]['シェア%']}%)。
    H2 (福山が最大) は<b>{"支持" if h2_ok else "部分支持/反証"}</b>。</li>
  <li>Top10 のうち<b>沿岸市町が占める比率</b>が圧倒的。 内陸市町
    (三次・庄原・北広島町など) は工業専用が極小 (ほぼゼロ)。</li>
  <li>4 拠点合計シェア = <b>{hub_share_pct:.1f}%</b>。 「臨海重工業」
    という広島県産業の<b>地理的偏在</b>が量化できた。</li>
</ul>""" +
    f"""<h3>13 用途地域 全県集計 (背景比較)</h3>
{df_to_html_table(T_zone_type)}

<h3>この表から読み取れること</h3>
<ul>
  <li>用途地域全体では<b>住居系 (YOTO_CD 1-7)</b> が圧倒多数。 工業
    専用 (12) は全県 {T_zone_type[T_zone_type['YOTO_CD']==12]['シェア%'].iloc[0]:.2f}%
    と少数だが、 「狭く・濃く」 配置されている。</li>
  <li>準工 (10) ・工業 (11) ・工業専用 (12) を合計しても
    {T_zone_type[T_zone_type['YOTO_CD'].isin([10,11,12])]['シェア%'].sum():.1f}% 程度。
    これが<b>広島県の工業地理の総枠</b>。</li>
</ul>"""
)
sections.append(("第1章 工業専用地域の地理 — 県内重工業帯の同定",
                  ch1_html))


# --- 5. 第2章: 福山JFE — 戦後高度成長の象徴 ---
ch2_purpose = f"""
<p><b>狙い</b>: 1965 年、 川崎製鉄 (現 JFE スチール) が福山臨海埋立地に
製鉄所を開設した。 高度成長期の象徴であり、 60 年後の現在も<b>世界級
鉄鋼コンビナート</b>として稼働している。 福山市の工業専用地域は県内
最大級の {fukuyama_area:.2f} km² (シェア {fukuyama_share_pct:.1f}%)。
本章では、 福山JFE 拠点の<b>地理的迫力</b>を zoom マップで示し、 H2
(福山JFE圧倒仮説) を検証する。</p>

<p><b>手法の要点</b>: 福山臨海部の bbox (北緯 34.43-34.55 / 東経
133.32-133.45) で工業専用地域 + 港湾施設をフィルタし、 福山市行政界
を背景に重ね合わせる。 同じ手法を 4 拠点に適用し、 small multiples
(2×2) で歴史軸の連続を見せる。</p>

<p><b>福山製鉄所の歴史的位置</b> (リテラシ): 1953 年に川崎製鉄千葉
製鉄所が完成した直後、 1958 年から福山進出が計画され、 1965 年 5 月に
火入れ式。 用地は瀬戸内海を埋め立てた人工地で、 当時としては「日本で
最大の製鉄所」 として高度成長期の象徴に。 高炉 4 基・厚板・冷延・電磁
鋼板を擁し、 西日本鉄鋼供給の中核。 2003 年に NKK と合併して JFE スチール
発足。</p>
"""

ch2_code = code('''
import shapely.geometry
import geopandas as gpd

# 福山JFE bbox (本記事独自定義)
FUKUYAMA_BBOX = (34.43, 34.55, 133.32, 133.45)
# = (lat_min, lat_max, lon_min, lon_max)

# bbox to projected polygon
bbox_4326 = shapely.geometry.box(133.32, 34.43, 133.45, 34.55)
bbox_proj = (gpd.GeoSeries([bbox_4326], crs="EPSG:4326")
             .to_crs("EPSG:6671").iloc[0])
xmin, ymin, xmax, ymax = bbox_proj.bounds

# 福山市の工業専用ポリゴンのみ
fukuyama_ind = indzone[indzone["source_city"] == "福山市"]
print(f"福山市 工業専用 {len(fukuyama_ind)} ポリゴン, "
      f"{fukuyama_ind.area.sum()/1e6:.2f} km²")

# 福山bbox 内の港湾係留施設
fukuyama_port = port_moor[
    (port_moor["lat"].between(34.43, 34.55)) &
    (port_moor["lon"].between(133.32, 133.45))
]
print(f"福山 bbox 内 港湾係留 = {len(fukuyama_port)} 件")
''')

ch2_html = (
    "<h3>狙い・手法</h3>" + ch2_purpose +
    "<h3>実装</h3>" + ch2_code +
    "<h3>結果図 — 4 産業拠点 zoom (2×2 small multiples)</h3>" +
    "<p><b>なぜこの図か</b>: 県全体マップでは各拠点が小さく見えてしまう。 "
    "<b>4 拠点を 2×2 で並べる small multiples</b> によって、 同じスケール "
    "での産業地と港湾施設の地理的関係を比較できる。 福山JFE の臨海埋立地 "
    "の迫力 (大規模かつ連続) と、 大竹コンビナートの細長い帯状立地、 "
    "呉海軍工廠跡の谷地形、 尾道造船の島嶼立地が一目で対比できる。</p>" +
    figure("assets/L88_fig02_hub_zoom.png",
           "Fig 2: 4 産業拠点 zoom — 工業専用地域 + 港湾施設 重ね合わせ") +
    f"""<h3>この図から読み取れること</h3>
<ul>
  <li><b>福山JFE</b>: 工業専用地域 {fukuyama_area:.2f} km² が湾岸に
    連続的に展開。 戦後埋立地特有の<b>大規模かつ整然とした矩形配置</b>。
    港湾係留施設も多数。 「都市計画 + 埋立 + 港湾」 を一体設計した
    高度成長期コンビナートのお手本。</li>
  <li><b>大竹コンビナート</b>: 山口県境の狭い平地に細長く伸びる工業地。
    戦前から続く立地のため、 福山ほどの整然さはなく、 むしろ<b>地形に
    縛られた帯状コンビナート</b>。 ナフサ輸入のための港湾施設も近接。</li>
  <li><b>呉海軍工廠</b>: 呉湾の谷地形に展開する工業地。 海軍工廠 (1903)
    の用地が現代の工業専用地域として継承されている。 港湾係留施設も
    密集 (戦前からの軍港インフラの民転)。</li>
  <li><b>尾道造船</b>: 向島・三原境付近の島嶼/沿岸の造船クラスタ。
    点状の工業専用ポリゴンが島々に分散。 「島の入江を活用する」 という
    瀬戸内造船の伝統的立地。</li>
  <li>4 拠点を比較すると、 「<b>戦前 = 地形に縛られた立地</b> (大竹・呉・
    尾道)」 vs 「<b>戦後 = 埋立で地形を作る</b> (福山)」 という時代差が
    地理に刻まれている。</li>
</ul>""" +
    f"""<h3>4 拠点 サマリー</h3>
{df_to_html_table(T_hub_summary)}

<h3>この表から読み取れること</h3>
<ul>
  <li>福山JFE が工業専用面積で最大 ({fukuyama_area:.2f} km²)。
    H2 (福山が市町別最大 + シェア ≥20%) は<b>{"支持" if h2_ok else "反証"}</b>。</li>
  <li>呉・尾道は中規模 (3-4 km²)、 大竹は小規模 (0.4 km²) ながらも
    歴史的重要拠点。 <b>面積だけが拠点の重要性ではない</b> ことを示す。</li>
  <li>各拠点の bbox 内に港湾施設が多数立地。 産業 = 港湾 = 道路 の
    一体構造が確認できる。</li>
</ul>"""
)
sections.append(("第2章 福山JFE — 戦後高度成長の象徴", ch2_html))


# --- 6. 第3章: 大竹コンビナート と 港湾規模 ---
ch3_purpose = f"""
<p><b>狙い</b>: 大竹市は工業専用地域面積では小さい ({float(city_indarea.get('大竹市', 0)):.2f} km²)
が、 戦前から続く三井化学・三菱レイヨンなどの<b>石油化学コンビナート</b>
が密集する。 さらに大竹港はナフサ輸入港として機能。 本章では、
<b>工業専用地域面積</b>と<b>港湾係留施設件数</b>の<b>市町別相関</b>
(H3) を量化し、 「臨海工業立地原理」 — 大規模工業地帯 = 港湾物流集積、
を検証する。</p>

<p><b>手法の要点</b>: 沿岸市町ごとに<b>工業専用面積 (km²)</b> と
<b>港湾_係留施設件数</b> をペアにし、 <code>pandas.Series.corr()</code>
で<b>Pearson 相関係数</b>を計算。 さらに散布図で各市町の位置を
プロットし、 線形回帰線で傾向を視覚化する。</p>

<p><b>Pearson 相関係数の直感的説明</b>: 2 変数の<b>同方向への動き
具合</b>を -1 〜 +1 で表す。 +1 = 完全な正比例、 0 = 無関係、
-1 = 完全な反比例。 r ≥ 0.5 は「中程度以上の正の相関」。</p>

<p><b>L33 (港湾係留 単独深掘り) との重複回避</b>: L33 は港湾名別ランキング
(広島港 125 件 max) を扱った。 本記事は<b>市町単位</b>での集計
(港湾名でなく市区町村１ で集約) を行い、 用途地域 (L17) と組み合わせる。
切り口が完全に独立。</p>
"""

ch3_code = code('''
import pandas as pd

# 市町別 工業専用面積
city_indarea = (indzone.groupby("source_city")["area_km2"]
                .sum())

# 市町別 港湾_係留 件数
port_moor = pd.read_csv(
    "data/extras/L33_port_mooring/harbor_mooring_facility.csv",
    encoding="utf-8-sig"
).dropna(subset=["開始位置緯度", "開始位置経度"])
city_harbor_moor = port_moor.groupby("市区町村１").size()

# 結合 → 沿岸市町のみで Pearson 相関
df = pd.concat([city_indarea, city_harbor_moor], axis=1).fillna(0)
df.columns = ["工業専用km2", "係留件数"]
COASTAL_CITIES = ["広島市", "呉市", "竹原市", "三原市", "尾道市",
                   "福山市", "大竹市", "東広島市", "廿日市市",
                   "江田島市", "海田町", "坂町"]
df_coast = df.loc[df.index.isin(COASTAL_CITIES)]
r = df_coast["工業専用km2"].corr(df_coast["係留件数"])
print(f"沿岸 {len(df_coast)} 市町 Pearson r = {r:.3f}")
''')

ch3_html = (
    "<h3>狙い・手法</h3>" + ch3_purpose +
    "<h3>実装</h3>" + ch3_code +
    "<h3>結果図 — 工業専用面積 × 港湾係留件数 散布</h3>" +
    "<p><b>なぜこの図か</b>: 「2 変量の関係」 を伝えるには散布図が一番。 "
    "<b>各市町を 1 点</b>として配置し、 沿岸市町は色付き・内陸市町は灰色 "
    "で区別。 線形回帰線で全体傾向を視覚化する。</p>" +
    figure("assets/L88_fig06_correlation.png",
           f"Fig 6: 工業専用面積 × 港湾係留件数 散布 (沿岸 r={h3_value:.3f})") +
    f"""<h3>この図から読み取れること</h3>
<ul>
  <li>沿岸 {len(corr_coastal)} 市町 Pearson r = <b>{h3_value:.3f}</b>。
    H3 (≥0.5 仮説) は<b>{"支持" if h3_ok else "反証"}</b>。</li>
  <li>右上: <b>福山市</b> ({fukuyama_area:.2f} km², 係留 {hub_harbor_count.get('福山JFE', 0)} 件)。
    工業面積も港湾も大きい — 戦後コンビナートの典型。</li>
  <li>左下: <b>大竹市</b> (工業 {float(city_indarea.get('大竹市', 0)):.2f} km², 係留 {hub_harbor_count.get('大竹コンビナート', 0)} 件)。
    工業面積は小さいが、 化学コンビナート専用の特殊港。 「面積よりも
    機能で測るべき拠点」 の例。</li>
  <li>呉市・広島市は中規模で散布の中央付近。 内陸市町 (灰色) は
    両軸ともゼロ近傍に集中。</li>
  <li>線形回帰線が右肩上がり = <b>臨海工業立地原理</b>の量的証拠。
    「重工業地は港湾と一体化する」 という直感を、 オープンデータで
    数値化できた。</li>
</ul>

<div class="note">
<b>反証の場合の解釈</b>: r が予想より低い場合は、 「港湾係留件数」 という
代理指標が<b>規模ではなく件数を測っている</b>ためと考えられる。
福山港・大竹港のような大規模専用港は、 <b>少数の大規模係留施設</b>
で大量物流を捌くため、 件数だけだと過小評価される。 真の物流量は
取扱貨物量 (DoBoX 外部データ) で見るべき。
</div>""" +
    f"""<h3>港湾_係留施設 港湾名別 Top10</h3>
{df_to_html_table(T_harbor_top)}

<h3>この表から読み取れること</h3>
<ul>
  <li>1 位は<b>{T_harbor_top.iloc[0]['港湾名']}</b>
    ({int(T_harbor_top.iloc[0]['係留施設数'])} 件)。 広島市の県都港湾
    として圧倒的。</li>
  <li>福山港は<b>{harbor_moor_rank[harbor_moor_rank['港湾名']=='福山港']['係留施設数'].iloc[0]
                    if (harbor_moor_rank['港湾名']=='福山港').any() else 'N/A'}</b> 件で
    上位 (本記事では JFE 物流港の代理指標)。 件数自体は中規模だが、
    1 件あたりの取扱量が大きい産業港の典型。</li>
  <li>蒲刈港・尾道糸崎港・鹿川港など<b>中堅瀬戸内港</b>が上位を占める。
    これらは漁港とは別カテゴリの「港湾」 = 重工業/物流港。</li>
</ul>"""
)
sections.append(("第3章 大竹コンビナート と 港湾規模 (H3 検証)",
                  ch3_html))


# --- 7. 第4章: 呉海軍工廠 と 産業港 vs 漁港 機能分化 ---
ch4_purpose = f"""
<p><b>狙い</b>: 1903 年に開設された呉海軍工廠は、 1945 年の終戦で解体
されたが、 用地と設備は IHI 呉造船・神戸製鋼として民需に転換した。
呉港は今も<b>産業港</b>として機能し、 周辺漁港とは明確に分離される。
本章では、 <b>港湾_係留 Top5</b> と <b>漁港_外郭 Top5</b> の市町重複
(H5) を検証し、 広島県沿岸の<b>機能分化</b>を量化する。</p>

<p><b>手法の要点</b>: 港湾施設 ({N_PORT_HARBOR:,} 件) と 漁港施設
({N_PORT_FISH:,} 件) を別々に集計し、 港湾名別ランキング Top5 を作成。
両 Top5 の<b>共通港湾名</b>を集合演算で取り出す。 重複が少ないほど
機能分化が強い。</p>

<p><b>港湾法 vs 漁港漁場整備法</b> (リテラシ): 「港湾」 (国土交通省所管) と
「漁港」 (水産庁所管) は異なる法律で指定される。 港湾は商工業・物流・
旅客を扱い、 漁港は漁業活動と漁獲物水揚げを扱う。 同じ瀬戸内海岸線にも
両者は分かれて立地し、 重複は基本的に発生しない (例外あり)。</p>
"""

ch4_code = code('''
import pandas as pd

# 港湾_係留 Top5
port_moor = pd.read_csv(
    "data/extras/L33_port_mooring/harbor_mooring_facility.csv",
    encoding="utf-8-sig"
)
top5_harbor = (port_moor.groupby("港湾名").size()
                .sort_values(ascending=False)
                .head(5).index.tolist())
print("港湾_係留 Top5:", top5_harbor)
# 例: ['広島港', '尾道糸崎港', '蒲刈港', '鹿川港', '釣士田港']

# 漁港_外郭 Top5
fish_outer = pd.read_csv(
    "data/extras/L32_port_breakwaters/fishing_port_outer_facility.csv",
    encoding="utf-8-sig"
)
top5_fish = (fish_outer.groupby("港湾名").size()
              .sort_values(ascending=False)
              .head(5).index.tolist())
print("漁港_外郭 Top5:", top5_fish)
# 例: ['倉橋', '豊島', '沖浦', '走', '箱崎']

# 重複
overlap = set(top5_harbor) & set(top5_fish)
print(f"重複: {len(overlap)} 件 = {overlap}")
''')

ch4_html = (
    "<h3>狙い・手法</h3>" + ch4_purpose +
    "<h3>実装</h3>" + ch4_code +
    "<h3>結果図 — 港湾係留 Top10 vs 漁港外郭 Top10</h3>" +
    "<p><b>なぜこの図か</b>: 「2 つの異なるランキング」 を比較するには <b>2 つの "
    "横棒グラフを並べる</b>のが直接的。 港湾_係留 Top10 と漁港_外郭 Top10 を "
    "同じスケール感で並べることで、 港湾名が完全に異なる = 機能分化、 を一目 "
    "で示す。</p>" +
    figure("assets/L88_fig05_port_vs_fishport.png",
           f"Fig 5: 港湾_係留 Top10 vs 漁港_外郭 Top10 (Top5 重複 {N_OVERLAP})") +
    f"""<h3>この図から読み取れること</h3>
<ul>
  <li>港湾係留 Top5 = <b>{', '.join(top5_harbor)}</b>。 広島県の<b>産業港・
    物流港</b>の 5 強。</li>
  <li>漁港外郭 Top5 = <b>{', '.join(top5_fishport)}</b>。 呉市
    倉橋・豊島・沖浦などの<b>離島漁業集積地</b>。</li>
  <li>両 Top5 の<b>重複 = {N_OVERLAP} 件</b>。 H5 (≤1 仮説) は
    <b>{"支持" if h5_ok else "反証"}</b>。</li>
  <li>「港湾」 = 都市部・物流ハブ・大規模工業地、 「漁港」 = 離島・小集落・
    沿岸漁業、 という<b>機能と場所の二重分化</b>が量的に確認できる。</li>
  <li>例外的に重複した港湾があれば、 それは「兼用港」 (商工業 + 漁業)
    の特殊事例。 本県では基本的に分離されている。</li>
</ul>""" +
    f"""<h3>漁港_外郭 港湾名別 Top10</h3>
{df_to_html_table(T_fish_top)}

<h3>この表から読み取れること</h3>
<ul>
  <li>1 位は<b>{T_fish_top.iloc[0]['港湾名']}</b>
    ({int(T_fish_top.iloc[0]['外郭施設数'])} 件)。 呉市倉橋島の
    伝統的漁業地。</li>
  <li>豊島・沖浦・走・箱崎など、 <b>呉市・江田島市・大崎上島町の離島
    漁港</b>が大半を占める。 産業港 (Top5 = 広島・尾道・蒲刈) とは
    地理的に明確に分離。</li>
  <li>漁港は「狭い湾」「小さな島」「外洋に近い」 立地が多く、 これは
    重工業港 (大規模埋立・大型船入港可能) とは正反対の立地条件。</li>
</ul>

<div class="note">
<b>呉海軍工廠の現代的位置</b>: 呉市は工業専用地域 {float(city_indarea.get('呉市', 0)):.2f} km²
を持ちつつ、 倉橋・豊島・沖浦などの離島漁港群も抱える「<b>重工業 + 漁業
の二重構造</b>」 自治体。 海軍工廠跡 (民需転換) と漁村文化 (江戸期から)
が共存している。
</div>"""
)
sections.append(("第4章 呉海軍工廠 と 産業港 vs 漁港 機能分化 (H5 検証)",
                  ch4_html))


# --- 8. 第5章: 沿岸工業 vs 高潮浸水 (H4 検証) ---
ch5_purpose = f"""
<p><b>狙い</b>: 福山JFE・大竹コンビナート・呉海軍工廠跡の 3 拠点は、
すべて<b>沿岸埋立地 / 沿岸低地</b>に立地する。 高度成長期の埋立工事で
作られた福山臨海部は海抜数 m の人工地。 ここに高潮 (例: 2018 年台風 21 号
クラスや、 想定 max シナリオ) が来襲したらどうなるか?  本章では DoBoX
の<b>高潮浸水想定 (#1149) max シナリオ</b>と工業専用地域の<b>交差面積</b>
(H4) を量化し、 「経済的価値とリスクの同居」 を可視化する。</p>

<p><b>手法の要点</b>: 高潮浸水想定 shp (7 ランク) を <code>union_all()</code>
で 1 ポリゴンに統合し、 各工業専用地域ポリゴンと <code>geometry.intersection
(storm_union).area</code> で交差面積 (m²) を計算。 ポリゴン面積で割って
比率を出し、 拠点別に集計する。</p>

<p><b>geopandas の geometry.intersection() 直感的説明</b>: 「左ジオメトリ」
と「右ジオメトリ」 の<b>共通部分</b>を新しいジオメトリとして返す。 結果に
.area を呼ぶと共通部分の面積。 union_all() で右側を 1 つにまとめておくと、
シンプル。</p>

<p><b>L44 (高潮単独深掘り) との重複回避</b>: L44 は SINSUIRANK 別深さ分布
の単独分析。 本記事は SINSUIRANK を統合 (1 ポリゴンに) し、 工業専用地域
との交差量化のみ。 切り口が独立。</p>
"""

ch5_code = code('''
import geopandas as gpd

# 高潮浸水想定 shp 読込
storm = gpd.read_file(
    "data/extras/L44_storm_surge/max_20210803/takashio_max.shp"
).to_crs("EPSG:6671")

# SINSUIRANK 別 7 ポリゴンを 1 つに統合
storm_union = storm.geometry.union_all()

# 工業専用ポリゴンごとに交差面積
indzone["storm_overlap_m2"] = (
    indzone.geometry.intersection(storm_union).area
)
indzone["storm_overlap_pct"] = (
    indzone["storm_overlap_m2"] / indzone.geometry.area.clip(lower=1) * 100
)
indzone["in_storm"] = (indzone["storm_overlap_m2"] > 1).astype(int)

# 全体集計
total_overlap_km2 = indzone["storm_overlap_m2"].sum() / 1e6
total_zone_km2 = indzone["area_km2"].sum()
print(f"工業専用 ∩ 高潮 = {total_overlap_km2:.2f} km² "
      f"({total_overlap_km2/total_zone_km2*100:.1f}%)")
''')

ch5_html = (
    "<h3>狙い・手法</h3>" + ch5_purpose +
    "<h3>実装</h3>" + ch5_code +
    "<h3>結果図 — 工業専用 vs 高潮浸水想定 重ね合わせ</h3>" +
    "<p><b>なぜこの図か</b>: 「経済とリスクの重なり」 を直感的に伝えるには "
    "<b>地図 (重ね合わせ)</b> と <b>拠点別交差%棒</b> を 2 枚で並べるのが最強。 "
    "地図で空間関係を、 棒で量的事実を、 同時に提示する。</p>" +
    figure("assets/L88_fig04_storm_overlay.png",
           f"Fig 4: 工業専用地域 vs 高潮浸水想定 max — "
           f"交差 {storm_overlap_total_km2:.2f} km² ({h4_value:.1f}%)") +
    f"""<h3>この図から読み取れること</h3>
<ul>
  <li>左 (地図): 工業専用 (赤) と高潮想定 (青) の重なりが瀬戸内沿岸全域で
    広範に確認できる。 福山臨海部・呉湾岸・大竹コンビナートの大半が
    高潮想定エリア内。</li>
  <li>右 (棒): 拠点別 高潮交差%。 拠点ごとに大きく異なる。
    {T_hub_storm.iloc[0]['拠点']}: {T_hub_storm.iloc[0]['高潮交差%']:.1f}%、
    {T_hub_storm.iloc[1]['拠点']}: {T_hub_storm.iloc[1]['高潮交差%']:.1f}% など。</li>
  <li>全県工業専用 vs 高潮 = <b>{h4_value:.1f}%</b>。
    H4 (≥40% 仮説) は<b>{"支持" if h4_ok else "反証"}</b>。</li>
  <li>「<b>戦後埋立地のコンビナート</b> = 高潮リスクの最前線」 という
    構造的事実が量化された。 1965 年 福山開設時の埋立工事は、 当時の
    気候設計 (高潮想定) に基づくが、 現在の温暖化想定では既に不足。</li>
  <li>大竹は山口県境の狭い平地で、 標高はやや高め。 福山ほど低くない
    ため交差比率が低めの可能性。</li>
</ul>

<div class="warn">
<b>解釈の注意</b>: 高潮浸水想定 max は<b>気候変動影響を加味した最大級
シナリオ</b>。 「即座に浸水する」 ではなく、 「想定最大級の高潮時に
浸水可能性のある地理」 を示す。 産業 BCP (事業継続計画) の参考とすべき
情報であって、 即座の操業停止判断ではない。
</div>""" +
    f"""<h3>4 拠点別 高潮交差量</h3>
{df_to_html_table(T_hub_storm)}

<h3>この表から読み取れること</h3>
<ul>
  <li>{T_hub_storm.sort_values('高潮交差%', ascending=False).iloc[0]['拠点']}
    が最大交差比率 ({T_hub_storm.sort_values('高潮交差%', ascending=False).iloc[0]['高潮交差%']:.1f}%)。
    沿岸埋立地中心の拠点ほどリスクが集中する。</li>
  <li>4 拠点合計 高潮交差面積 = <b>{T_hub_storm['高潮交差km2'].sum():.2f} km²</b>。
    これは広島県全工業専用 ∩ 高潮 の {T_hub_storm['高潮交差km2'].sum() / max(storm_overlap_total_km2, 1e-9) * 100:.1f}%
    を 4 拠点で占める。 リスクも 4 拠点に集中。</li>
  <li>拠点ごとに交差比率が大きく異なるため、 <b>拠点別 BCP 設計</b>が
    必要。 一律ではなく、 福山JFE には大規模高潮対策、 呉には軍港由来の
    既存防潮インフラ活用、 大竹には化学物質流出対策、 など個別設計。</li>
</ul>"""
)
sections.append(("第5章 沿岸工業 vs 高潮浸水 (H4 検証)", ch5_html))


# --- 9. 第6章: 産業 DNA の連続性 — 用途構成と 4 拠点歴史軸 ---
ch6_purpose = f"""
<p><b>狙い</b>: 戦前 - 戦後を貫く広島県の臨海重工業の<b>連続性 (産業 DNA)</b>
を、 用途地域 13 類型の全県構成と 4 拠点の歴史軸で締めくくる。 「呉海軍
工廠 → 戦後 IHI 造船 → 工業専用地域指定」 「大竹戦前化学 → 戦後石油
化学コンビナート → 工業専用」 「福山戦後高度成長 → JFE → 工業専用」
という時間連続性が、 現代の用途指定にどう刻まれているか。</p>

<p><b>手法の要点</b>: 13 用途地域の全県シェア円グラフ + 4 拠点の現在
工業専用面積 棒グラフ の 2 枚並列。 円で「広島県全体に占める工業の
位置」 を、 棒で「4 拠点の現在規模」 を同時に示す。</p>

<p><b>「産業 DNA」 の操作的定義</b>: 戦前の産業立地 (海軍工廠・戦前
化学・戦前造船) が、 戦後の用途指定で<b>用途地域として保存</b>されている
状態。 海軍工廠は 1945 年に解体されたが、 跡地は工業専用に指定され、
IHI・神戸製鋼が継承。 これが「DNA の継承」 の意味。</p>
"""

ch6_code = code('''
import matplotlib.pyplot as plt

# 13 用途地域 全県シェア
zone_summary = (zone.groupby("YOTO_CD")["area_km2"].sum()
                .reset_index())
zone_summary["share_pct"] = (
    zone_summary["area_km2"] / zone_summary["area_km2"].sum() * 100
)
print(zone_summary)

# 4 拠点の現在 工業専用面積
HUBS = {
    "福山JFE":   "福山市",
    "大竹コンビナート": "大竹市",
    "呉海軍工廠":   "呉市",
    "尾道造船":     "尾道市",
}
hub_areas = {h: float(city_indarea.get(c, 0)) for h, c in HUBS.items()}
print(hub_areas)
# {"福山JFE": 6.57, "大竹コンビナート": 0.36, "呉海軍工廠": 3.66, "尾道造船": 3.15}
''')

ch6_html = (
    "<h3>狙い・手法</h3>" + ch6_purpose +
    "<h3>実装</h3>" + ch6_code +
    "<h3>結果図 — 13 用途地域 円グラフ + 4 拠点 歴史軸</h3>" +
    "<p><b>なぜこの図か</b>: <b>左 (円グラフ)</b>で「広島県全体での工業の位置」、 "
    "<b>右 (棒)</b>で「4 拠点の現在の規模」 を 1 セット 2 視点で示す。 学習者が "
    "「県全体ではマイノリティだが、 4 拠点に濃縮されている」 という構造を直感 "
    "的に把握できる。</p>" +
    figure("assets/L88_fig07_pie_history.png",
           "Fig 7: 13 用途地域 全県構成 + 4 拠点 歴史軸") +
    f"""<h3>この図から読み取れること</h3>
<ul>
  <li>左 (円): 工業専用 (赤) は全県のわずか
    {T_zone_type[T_zone_type['YOTO_CD']==12]['シェア%'].iloc[0]:.2f}%。
    住居系が圧倒多数。 「広島県は住宅都市県」 が用途指定面積では正解。</li>
  <li>右 (棒): 4 拠点の現在工業専用面積。 福山JFE が圧倒的、 続いて
    呉・尾道・大竹。 各拠点の<b>歴史軸</b> (戦後 1965 / 戦前 1903 /
    戦前 1920s / 戦前) が併記され、 「歴史」 と「規模」 の二軸で
    拠点を比較できる。</li>
  <li>大竹は面積では小さいが、 戦前から続く<b>歴史的連続性</b>では
    最も古い (1920 年代から)。 「面積だけが拠点の重要性ではない」
    という第2章の議論を補強する。</li>
</ul>""" +
    f"""<h3>5 仮説 検証総合</h3>
{df_to_html_table(T_hypothesis)}

<h3>この表から読み取れること</h3>
<ul>
  <li>支持された仮説 = <b>{n_supported}/5</b>。</li>
  <li><b>H1 (臨海重工業集中)</b>: 沿岸 4 拠点シェア {h1_value:.1f}% で
    {"支持" if h1_ok else "反証"}。 {"瀬戸内臨海への偏在が量的に確認。" if h1_ok else "予想より分散している可能性。"}</li>
  <li><b>H2 (福山JFE圧倒)</b>: 福山シェア {h2_value:.1f}% で
    {"支持" if h2_ok else "反証"}。 {"戦後高度成長期コンビナートの圧倒的規模。" if h2_ok else "他市町 (広島市等) の方が大きい可能性。"}</li>
  <li><b>H3 (港湾相関)</b>: 沿岸市町 r = {h3_value:.3f} で
    {"支持" if h3_ok else "反証"}。 {"臨海工業立地原理の量的証拠。" if h3_ok else "係留件数は規模を測れない代理指標の限界。"}</li>
  <li><b>H4 (高潮交差)</b>: {h4_value:.1f}% で
    {"支持" if h4_ok else "反証"}。 {"経済とリスクの同居。" if h4_ok else "工業地は意外と高潮を避けて立地している可能性。"}</li>
  <li><b>H5 (機能分化)</b>: Top5 重複 {h5_value} で
    {"支持" if h5_ok else "反証"}。 {"重工業港と漁業港の空間的二極化。" if h5_ok else "兼用港が多く存在する可能性。"}</li>
</ul>"""
)
sections.append(("第6章 産業 DNA の連続性 — 13 用途構成と 4 拠点歴史軸",
                  ch6_html))


# --- 10. 第7章: 仮説検証総合 ---
hyp_html = f"""
<p>5 つの仮説を 6 章にわたり検証してきた。 結果をまとめる。</p>

<h3>横断的可視化 — 仮説検証ダッシュボード</h3>
<p><b>なぜこの図か</b>: 4 つの視点 (5 仮説検証 / 4 拠点 工業 vs 高潮 /
港湾 vs 漁港 Top5 / 13 用途 沿岸内陸) を<b>1 枚にまとめたダッシュボード</b>
で、 6 章の知見を一望する。 学習者が最後の総括として何度でも見返せる。</p>
{figure("assets/L88_fig08_dashboard.png",
        "Fig 8: 仮説検証総合ダッシュボード — 産業地理 4 視点")}

<h3>5 仮説 検証サマリー表</h3>
{df_to_html_table(T_hypothesis)}

<h3>広島県産業物語の 6 つの発見</h3>
<ol>
  <li><b>臨海集中度</b>: 工業専用地域 {A_INDZONE:.1f} km² のうち
    {hub_share_pct:.1f}% が沿岸 4 拠点市町に集中。 「臨海重工業」 は
    地理的に偏在する。</li>
  <li><b>福山JFE圧倒</b>: 福山市が市町別最大 ({fukuyama_share_pct:.1f}%)。
    1965 年の戦後高度成長期コンビナート開設が、 60 年経った現代も
    広島県工業地理の頂点を占める。</li>
  <li><b>港湾規模相関</b>: 沿岸市町の工業専用面積と港湾係留件数は
    Pearson r = {h3_value:.3f}。 {("正の相関で「臨海工業立地原理」 を量化" if h3_ok else "代理指標の限界で相関は弱め")}。</li>
  <li><b>高潮交差</b>: 工業専用 ∩ 高潮 = {h4_value:.1f}%。
    沿岸埋立地中心のコンビナートは気候災害の最前線にある。</li>
  <li><b>機能分化</b>: 港湾係留 Top5 vs 漁港外郭 Top5 重複 {h5_value} 件。
    「重工業港 vs 漁業港」 の空間的二極化が確認できた。</li>
  <li><b>産業 DNA の連続性</b>: 戦前 (大竹化学・呉海軍工廠) ・戦前後
    境界 (尾道造船) ・戦後高度成長 (福山JFE) という<b>4 拠点の歴史軸</b>
    が、 現代の用途地域指定 (工業専用) として連続的に保存されている。</li>
</ol>

<div class="note">
<b>L17 / L21 / L32-L34 / L52-L55 / L66 / L44 と L88 の補完性</b>:
同じ DoBoX 産業関連データが、 切り口を変えると全く異なる物語を語る。
<ul>
  <li><b>L17 (用途地域 21 統合)</b>: 13 用途地域 + 3 類型 + k-means プロファイル</li>
  <li><b>L21 (宅地開発 1,038 件)</b>: 期間・面積・年別の単独深掘り</li>
  <li><b>L32-L34 (港湾系)</b>: 外郭/係留/臨港交通 の単独深掘り</li>
  <li><b>L44 (高潮浸水想定)</b>: SINSUIRANK 別深さ分布</li>
  <li><b>L66 (橋梁 4,203 件)</b>: 架設年度・延長・幅員</li>
  <li><b>L88 (本記事, 物語)</b>: 工業専用 + 港湾 + 高潮 を 4 拠点歴史軸で物語化</li>
</ul>
8 つの記事が揃って初めて<b>広島県産業地理の研究地図</b>になる。
</div>
"""
sections.append(("仮説検証総合", hyp_html))


# --- 11. 発展課題 ---
dev_html = f"""
<h3>発展課題 (結果X → 新仮説Y → 課題Z の論理鎖)</h3>

<h4>発展1: 産業 DNA の時系列追跡 — 1965 福山開設前後の地形変化</h4>
<ul>
  <li><b>結果X</b>: 福山市の工業専用地域は {fukuyama_area:.2f} km²
    (シェア {fukuyama_share_pct:.1f}%) で県内最大だが、 これは現在の
    指定面積。 1965 年開設前は何があったかは見えない。</li>
  <li><b>新仮説Y</b>: 1960 年代の旧地形図 (国土地理院) と現在の用途地域
    を重ねれば、 <b>埋立面積</b>を量化できる。 福山臨海部の工業専用地域
    の<b>X% は人工埋立地</b>であるはず。</li>
  <li><b>課題Z</b>: 国土地理院旧版地形図 (1960 年代) を取得 (有償サービスあり)、
    海岸線をデジタイズし現在の海岸線と差分を取る。 「<b>戦後埋立面積</b>」
    を市町別に量化し、 工業専用地域の<b>埋立由来比率</b>を計算。
    高度成長期土木史の量的記述になる。</li>
</ul>

<h4>発展2: 工業専用地域内の建物密度 — 旧海軍工廠跡の継承度量化</h4>
<ul>
  <li><b>結果X</b>: 呉市の工業専用地域 {float(city_indarea.get('呉市', 0)):.2f} km²
    は旧海軍工廠跡の継承だが、 「実際にどれだけ工場が建っているか」
    は本記事では分からない。</li>
  <li><b>新仮説Y</b>: 国土数値情報「建物利用」 メッシュデータと用途地域
    を重ねれば、 「工業専用地域内の工業建物充填率」 を計算できる。
    呉海軍工廠跡 vs 福山JFE で<b>充填率の違い</b>があるはず。</li>
  <li><b>課題Z</b>: 国交省 国土数値情報 (建物データ A38) を取得し、
    用途地域 (YOTO_CD=12) との sjoin で「工業専用 1 km² あたり工場
    建物数」 を計算。 拠点別比較で「現代も実稼働している産業地」 と
    「形骸化した産業地」 を判別。</li>
</ul>

<h4>発展3: 4 拠点の津波 + 高潮 ダブル交差 — 複合災害想定</h4>
<ul>
  <li><b>結果X</b>: 本記事は工業専用 vs 高潮のみ検証。 津波浸水想定
    (#47 系) との交差は扱わなかった。</li>
  <li><b>新仮説Y</b>: 工業専用地域は<b>高潮 + 津波の両方</b>と交差する
    エリアが存在し、 そこは<b>複合災害最大リスク地</b>。
    特に呉湾・福山臨海部に集中するはず。</li>
  <li><b>課題Z</b>: 津波浸水想定 shp (DoBoX #47, 1.25M メッシュ) を
    union_all() で 1 ポリゴン化し、 工業専用 ∩ 津波、 工業専用 ∩
    高潮 ∩ 津波 (3 重) を sjoin で量化。 ダブル/トリプル交差面積を
    拠点別に算出し、 BCP 優先度マップを作る。</li>
</ul>

<h4>発展4: 産業 4 拠点 BCP 物資輸送網 — 緊急輸送道路との結合</h4>
<ul>
  <li><b>結果X</b>: 4 拠点別に港湾施設件数は集計したが、
    <b>陸路 (緊急輸送道路, L72)</b> との結合は未検証。</li>
  <li><b>新仮説Y</b>: 4 拠点の工業専用地域は、 <b>緊急輸送道路 1 km
    圏内</b>に立地している (BCP 必須要件)。 拠点別「緊急輸送道路
    アクセス度」 を量化できる。</li>
  <li><b>課題Z</b>: L72 緊急輸送道路 GeoJSON を読込、 工業専用地域
    ポリゴンの重心から最寄り緊急路までの距離を sjoin_nearest で計算。
    「アクセス度」 = 1 / 平均距離 を拠点別に集計し、
    災害時物資輸送の脆弱性マップを作成。</li>
</ul>

<h4>発展5: 全国比較 — 広島県の臨海重工業立地の特異性</h4>
<ul>
  <li><b>結果X</b>: 本記事は広島県内の比較に閉じている。
    全国で見ると広島県の臨海重工業集中度はどう位置づくか不明。</li>
  <li><b>新仮説Y</b>: 千葉県 (京葉工業地帯) ・三重県 (四日市) ・岡山県
    (水島) など同様の<b>戦後高度成長型臨海コンビナート県</b>と比較
    すれば、 広島県は<b>中規模クラス</b> に位置するはず。</li>
  <li><b>課題Z</b>: 国土数値情報「用途地域」 全国データを取得し、
    都道府県別 工業専用地域面積ランキングを作成。 広島県の絶対値・
    シェアを全国比較で位置づけ、 「広島県産業地理の特異性」 を
    定量化する。</li>
</ul>
"""
sections.append(("発展課題", dev_html))


# --- 12. 既存記事との関係 ---
relat_html = f"""
<h3>広島県産業地理の研究地図 — 9 記事の関係</h3>
<p>本記事 L88 を含む<b>産業関連 9 記事</b>は、 同じ DoBoX 産業データから
全く異なる切り口で物語を引き出す<b>研究地図</b>を構成している。</p>

{df_to_html_table(T_compare)}

<h3>この表から読み取れること</h3>
<ul>
  <li>同じ用途地域 + 港湾 + 防災データが、 9 つの<b>独立した切り口</b>で
    読み解かれている。</li>
  <li>L88 (本記事) は他 8 記事と<b>計算指標が完全に重複しない</b>:
    13 用途地域並列分析 / 単独深掘り / 期間集計 / k-means プロファイル
    / 港湾名別ランキング / SINSUIRANK 深さ分布 — これら<b>全てを
    再計算しない</b>。</li>
  <li>L88 が初めて触れる切り口: <b>4 産業拠点の歴史軸</b> + <b>工業専用 vs
    港湾係留 相関</b> + <b>工業専用 vs 高潮交差</b> + <b>産業港 vs 漁港
    機能分化</b>。</li>
</ul>

<div class="note">
<b>研究の協働性</b>: 1 つのデータを 1 人で読むより、 9 つの切り口で 9 人が
読む方が、 はるかに豊かな物語が生まれる。 オープンデータの真の力は
「みんなが見て、 みんなが書く」 ことにある。 本記事 L88 は、 その 1 ピース。
</div>
"""
sections.append(("既存記事との関係 — 産業研究地図", relat_html))


# =============================================================================
# 14. HTML 書き出し
# =============================================================================
n_bridges_total = len(pd.read_csv(BRIDGE_CSV, encoding='utf-8-sig'))
html = render_lesson(
    num=88,
    title="L88 M4 産業物語 — 福山JFE・大竹コンビナート・呉海軍工廠・工業専用地域",
    tags=["M系統合", "産業地理", "工業専用地域",
          "福山JFE", "呉海軍工廠", "大竹コンビナート",
          "尾道造船", "高潮交差", "geopandas"],
    time="40 分",
    level="リテラシ + 空間統合 + 歴史物語",
    data_label=("用途地域 21 dataset (#794-946) / 港湾施設 / "
                f"高潮浸水 #1149 / 橋梁 #1062 / "
                f"L15 行政界 (工業専用 {N_INDZONE:,} ポリゴン)"),
    sections=sections,
    script_filename="L88_M4_industry_story.py",
)
out_html = LESSONS / "L88_M4_industry_story.html"
out_html.write_text(html, encoding="utf-8")
print(f"  HTML 書き出し: {out_html} ({time.time()-t13:.2f}s)", flush=True)


print(f"\n=== 完了 ({time.time()-t_all:.1f}s) ===", flush=True)
print(f"  支持仮説: {n_supported}/5", flush=True)
print(f"  ファイル:", flush=True)
print(f"    HTML: {out_html}", flush=True)
print(f"    PY: lessons/L88_M4_industry_story.py", flush=True)
print(f"    中間 CSV: assets/L88_*.csv (13 件)", flush=True)
print(f"    図: assets/L88_fig*.png (8 件)", flush=True)