# -*- coding: utf-8 -*- """L35 排水機場 2 件統合分析 — 広島県の 15 機場 (港湾海岸 3 + 河川 12) が支える内水浸水防御の最終線 カバー宣言: 本記事は DoBoX のシリーズ 「排水機場 港湾海岸 / 河川」 2 件 (dataset_id = 1269, 1273) を統合し、広島県内の排水機場の地理分布・ ポンプ吐出能力・建設年代を分析する研究記事である。 L32 (外郭施設) / L33 (係留施設) / L34 (臨港交通施設) は港湾の外洋・接岸・陸側を 扱った 3 層階層であった。本記事 L35 は港湾系シリーズの最終記事であるが、 扱う施設は港湾施設ではなく「内水排水施設」であり、 港湾海岸エリアの内水排水 (3 件) と 河川流域の内水排水 (12 件) の両者を扱う。 すなわち、外郭で外海を防ぎ、係留で船を繋ぎ、交通で陸へ繋いだ後の「陸側に溜まる雨水・河川氾濫水」を ポンプで強制排出する「内水浸水防御の最終線」を扱う。 研究の問い (RQ): 広島県内の排水機場 (合計 15 件) は、設置目的 (港湾海岸 / 河川)・地理分布 (沿岸低地 / 河川合流部 / 干拓地)・規模 (吐出量 m³/s) の観点でどう構成され、 内水浸水域とどう重なるのか? 各機場が守る集水区はどこか? (1) 2 dataset の構造比: 河川 12 機場 vs 港湾海岸 3 機場 (1:4 偏重)。 この偏重は内水氾濫の主因が「河川合流部の堰き止め」であることを示唆。 (2) 福山市集中: 12 河川機場のうち 7 件が福山市。 福山平野の干拓地性質と芦田川水系の合流複雑性を反映。 (3) 排水容量階層: 最大 30 m³/s (手城川・岡ノ下川) ÷ 最小 1.3 m³/s (沢) = 23倍。 2 大型機場が県全体の 46% の排水能力を担う。 (4) 建設年代と更新: 昭和40 (1965) 〜 平成28 (2016) の 51 年間で建設。 昭和期 7 件、平成期 8 件。新安川は昭和41 (第一) → 平成28 (第二) と62 年差で増設。 (5) 浸水想定区域との重なり: 各機場の周辺 1km バッファに DoBoX 河川浸水想定 (rank 別) が重なるか? 排水機場は浸水想定エリアの内側に立地するはずで、これを実データで検証する。 (6) ポンプ単機平均能力 (m³/s/台): 大型機 ≈ 7-15、小型機 ≈ 0.5-2。 新旧比較で機械の大型化 (= 1 台あたり能力増) が読めるか。 仮説 H1〜H6: H1 (主導カテゴリは河川): 河川 12/15 = 80%、港湾海岸 3/15 = 20%。 吐出量でも河川 117.5 m³/s vs 港湾海岸 12.7 m³/s = 9.2:1 比。 内水排水の主戦場は河川合流低地で、港湾海岸はサブ。 H2 (福山市偏重): 河川機場 12 件のうち福山市 = 7 件 (58%)。 福山平野は芦田川河口の干拓地で、低地が多く水が溜まりやすい構造的特性。 広島市 (3 件) より絶対数で多い。 H3 (大型機場の二極集中): 30 m³/s 機場が 2 件 (手城川・岡ノ下川) で、 これだけで県内総排水容量 130.2 m³/s の 46.1% を占有。 Lorenz 曲線で強い集中構造。 H4 (建設年代の 2 ピーク): 昭和40-50年代 (1965-80, 戦後復興後の都市化期) と 平成20年代 (2008-16, 集中豪雨対応強化期) に建設が集中。 バブル経済期 (昭和60-平成5) は新規少なく、補強・更新中心。 H5 (浸水想定との重なり): 全 15 機場のうち 13 件以上 (≥87%) が DoBoX の河川浸水想定区域 1km バッファ内に位置。 排水機場は「想定浸水エリア内の最低標高点」に置かれる工学的必然性を実証。 H6 (大型化の傾向): 平成期建設機の 1 台あたり吐出量 (m³/s/台) は 昭和期建設機より大きい。機械技術の進歩で 1 台で扱える流量が増えた。 要件 S 準拠: 1 分以内完走。データ件数が少ないので空間処理は軽量。 要件 T 準拠: 県全域マップ + 福山ズーム + 浸水想定重ね合わせ + ポンプ能力点マップ。 要件 Q 準拠: 図 9 種以上、表 8 種以上 (15 件 × 4 属性軸 (位置/年代/容量/台数))。 データ仕様 (2 件は同一 9 列スキーマ): A. 排水機場（港湾海岸） 1269: CSV, 3 行, 大竹港海岸・瀬戸田港海岸・竹原港海岸 B. 排水機場（河川） 1273: CSV, 12 行, 広島・福山・尾道の各河川合流部 共通列: 施設の名称, 事務所名, 所在地, 緯度(10進数), 経度(10進数), 海岸・河川名, 完成年月日, 総排水量(m3), ポンプ台数(台) ※ 「総排水量(m3)」列は単位表記が m3 だが、値の桁 (1.3〜30) と 排水機場の業界慣例から毎秒の吐出量 m³/s と解釈するのが妥当。 30 m³/s = 25mプール 1 杯 (約 600 m³) を 20 秒で排水できる能力。 本記事ではこの解釈で記述する。 """ 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 plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False t0 = time.time() print("=== L35 排水機場 2 件統合分析 ===", flush=True) # ============================================================================= # 0. 定数: 2 dataset_id とカテゴリ定義 # ============================================================================= DATA_DIR = ROOT / "data" / "extras" / "L35_drainage_pumps" DATA_DIR.mkdir(parents=True, exist_ok=True) ADMIN_DIR = ROOT / "data" / "extras" / "L15_admin_zones" FLOOD_SHP = ROOT / "data" / "extras" / "flood_shp" / "shinsui_souteisaidai" / "shinsui_souteisaidai.shp" TARGET_CRS = "EPSG:6671" # JGD2011 平面直角 III, m 単位 SERIES = [ (1269, "排水機場（港湾海岸）", "港湾海岸", "drainage_pump_port_coast.csv", 39770), (1273, "排水機場（河川）", "河川", "drainage_pump_river.csv", 39774), ] # カテゴリ色 CAT_COLOR = {"港湾海岸": "#0969da", "河川": "#1a7f37"} # 浸水深ランクの公式凡例 (DoBoX 河川浸水想定 rank の意味) FLOOD_RANK = { 10: "0.0–0.5m 未満", 20: "0.5–1.0m 未満", 30: "1.0–2.0m 未満", 40: "2.0–3.0m 未満", 50: "3.0–5.0m 未満", 60: "5.0–10.0m 未満", 70: "10.0–20.0m 未満", 80: "20.0m 以上", } # 青系の階調 (浸水深が深くなるほど濃く) RANK_COLORS = { 10: "#deebf7", 20: "#c6dbef", 30: "#9ecae1", 40: "#6baed6", 50: "#4292c6", 60: "#2171b5", 70: "#08519c", 80: "#08306b", } # ============================================================================= # 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"L35/{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. 列名の正規化と派生指標 # ============================================================================= print("\n[2] 派生指標の計算", flush=True) t1 = time.time() # 列名のエイリアス ALL = ALL.rename(columns={ "施設の名称": "name", "事務所名": "office", "所在地": "address", "緯度(10進数)": "lat", "経度(10進数)": "lon", "海岸・河川名": "watercourse", "完成年月日": "完成", "総排水量(m3)": "discharge_m3s", # 業界慣例で m³/s と解釈 "ポンプ台数(台)": "n_pumps", }) # 完成年を西暦に変換 def jp_to_west(s: str) -> int | None: if pd.isna(s): return None s = str(s).strip() # 「昭和41年」「平成1年」「平成28年」などを西暦に if s.startswith("昭和"): try: n = int(s.replace("昭和", "").replace("年", "")) return 1925 + n # 昭和元年=1926, n年 = 1925+n except ValueError: return None elif s.startswith("平成"): try: n = int(s.replace("平成", "").replace("年", "")) return 1988 + n # 平成元年=1989, n年 = 1988+n except ValueError: return None return None ALL["year"] = ALL["完成"].apply(jp_to_west) ALL["era"] = ALL["完成"].apply(lambda s: "昭和" if str(s).startswith("昭和") else ("平成" if str(s).startswith("平成") else "不明")) ALL["m3s_per_pump"] = ALL["discharge_m3s"] / ALL["n_pumps"] # 市町を所在地から抽出 (例: "広島市佐伯区..." → "広島市") def extract_city(addr: str) -> str: if pd.isna(addr): return "" s = str(addr) for k in ["広島市", "福山市", "尾道市", "竹原市", "大竹市", "三原市", "呉市", "東広島市", "廿日市市", "安芸高田市", "府中市", "三次市", "庄原市"]: if k in s: return k return "" ALL["city"] = ALL["address"].apply(extract_city) print(ALL[["name", "port_category", "city", "watercourse", "完成", "year", "era", "discharge_m3s", "n_pumps", "m3s_per_pump"]].to_string(index=False)) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 3. GeoDataFrame 化 (緯度経度から Point geometry) # ============================================================================= print("\n[3] GeoDataFrame 化", flush=True) t1 = time.time() ALL["geometry"] = [Point(lon, lat) for lon, lat in zip(ALL["lon"], ALL["lat"])] gdf = gpd.GeoDataFrame(ALL, geometry="geometry", crs="EPSG:4326").to_crs(TARGET_CRS) print(f" geom OK: {len(gdf)} / {len(ALL)} (位置情報完備)", flush=True) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 4. カテゴリ別集計 (件数 / 容量 / 台数) # ============================================================================= print("\n[4] カテゴリ別集計", flush=True) t1 = time.time() cat_agg = ALL.groupby("port_category").agg( n=("name", "count"), total_m3s=("discharge_m3s", "sum"), mean_m3s=("discharge_m3s", "mean"), median_m3s=("discharge_m3s", "median"), n_pumps_sum=("n_pumps", "sum"), n_pumps_mean=("n_pumps", "mean"), ).round(2) cat_agg.loc["合計"] = [ cat_agg["n"].sum(), cat_agg["total_m3s"].sum(), ALL["discharge_m3s"].mean().round(2), ALL["discharge_m3s"].median().round(2), cat_agg["n_pumps_sum"].sum(), ALL["n_pumps"].mean().round(2), ] cat_agg["n"] = cat_agg["n"].astype(int) cat_agg["n_pumps_sum"] = cat_agg["n_pumps_sum"].astype(int) cat_agg.to_csv(ASSETS / "L35_category_summary.csv", encoding="utf-8-sig") print(cat_agg.to_string()) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 5. 市町別集計 # ============================================================================= print("\n[5] 市町別集計", flush=True) t1 = time.time() city_agg = ALL.groupby(["city", "port_category"]).agg( n=("name", "count"), total_m3s=("discharge_m3s", "sum"), n_pumps=("n_pumps", "sum"), ).reset_index() # 市町だけのワイド集計 city_wide = ALL.groupby("city").agg( n=("name", "count"), total_m3s=("discharge_m3s", "sum"), n_pumps=("n_pumps", "sum"), ).round(2).sort_values("n", ascending=False) city_wide.to_csv(ASSETS / "L35_city_summary.csv", encoding="utf-8-sig") print(city_wide.to_string()) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 6. 建設年代の集計 # ============================================================================= print("\n[6] 建設年代の集計", flush=True) t1 = time.time() era_agg = ALL.groupby(["era", "port_category"]).agg( n=("name", "count"), total_m3s=("discharge_m3s", "sum"), ).reset_index() year_bins = [1960, 1970, 1980, 1990, 2000, 2010, 2020] year_labels = ["1960s", "1970s", "1980s", "1990s", "2000s", "2010s"] ALL["year_bin"] = pd.cut(ALL["year"], bins=year_bins, labels=year_labels, right=False) year_pivot = ALL.groupby(["year_bin", "port_category"], observed=True).size().unstack(fill_value=0) year_pivot.to_csv(ASSETS / "L35_year_pivot.csv", encoding="utf-8-sig") print(year_pivot.to_string()) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 7. ポンプ能力 (1 台あたり m3/s) の分析 # ============================================================================= print("\n[7] ポンプ能力分析 (m3/s/台)", flush=True) t1 = time.time() # 古い (昭和) vs 新しい (平成) の 1 台あたり吐出量比較 era_pump = ALL.groupby("era").agg( n=("name", "count"), mean_m3s_per_pump=("m3s_per_pump", "mean"), median_m3s_per_pump=("m3s_per_pump", "median"), mean_n_pumps=("n_pumps", "mean"), mean_total_m3s=("discharge_m3s", "mean"), ).round(3) era_pump.to_csv(ASSETS / "L35_era_pump_efficiency.csv", encoding="utf-8-sig") print(era_pump.to_string()) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 8. 県境 polygon の load # ============================================================================= print("\n[8] 県境の 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) # ============================================================================= # 9. 浸水想定との重なり判定 (1 km バッファ) # ============================================================================= print("\n[9] 浸水想定区域との重なり判定", flush=True) t1 = time.time() BUFFER_M = 1000.0 # 1 km バッファで「機場が守る範囲」と仮定 FLOOD_OK = FLOOD_SHP.exists() gdf["buffer_1km"] = gdf.geometry.buffer(BUFFER_M) gdf["near_flood"] = False gdf["max_rank_in_buffer"] = 0 gdf["flood_area_ha_in_buffer"] = 0.0 flood_clip = None if FLOOD_OK: print(f" flood shp ロード: {FLOOD_SHP}", flush=True) flood = gpd.read_file(FLOOD_SHP, encoding="cp932").to_crs(TARGET_CRS) print(f" flood polygons: {len(flood)} 件", flush=True) # 機場 1km バッファ内に flood polygon があるか sjoin buf_gdf = gpd.GeoDataFrame(gdf[["name", "port_category"]], geometry=gdf["buffer_1km"], crs=TARGET_CRS) # バッファ内の flood polygon を抽出 (intersects) for i, row in buf_gdf.iterrows(): sub = flood[flood.intersects(row.geometry)] if len(sub) > 0: gdf.at[i, "near_flood"] = True gdf.at[i, "max_rank_in_buffer"] = int(sub["rank"].max()) # クリップして面積算出 try: clipped = gpd.clip(sub, row.geometry) gdf.at[i, "flood_area_ha_in_buffer"] = float(clipped.geometry.area.sum() / 10000.0) except Exception: gdf.at[i, "flood_area_ha_in_buffer"] = 0.0 n_near = int(gdf["near_flood"].sum()) print(f" 浸水想定 1km 圏内: {n_near}/{len(gdf)} = {n_near/len(gdf)*100:.1f}%", flush=True) print(f" 最大ランク max: {gdf['max_rank_in_buffer'].max()} ({FLOOD_RANK.get(int(gdf['max_rank_in_buffer'].max()), '?')})", flush=True) else: print(f" WARN: {FLOOD_SHP} not found, flood overlap analysis skipped", flush=True) n_near = 0 print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 10. 中間 CSV 出力 (geometry を除いた素データ) # ============================================================================= print("\n[10] 中間 CSV 出力", flush=True) t1 = time.time() out = gdf.drop(columns=["geometry", "buffer_1km"]).copy() out.to_csv(ASSETS / "L35_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 / "L35_series.csv", index=False, encoding="utf-8-sig") print(f" saved csv ({time.time()-t1:.1f}s)", flush=True) # 値の事前整理 (HTML から参照) n_total = len(ALL) n_h = int((ALL["port_category"] == "港湾海岸").sum()) n_r = int((ALL["port_category"] == "河川").sum()) total_m3s = float(ALL["discharge_m3s"].sum()) total_m3s_h = float(ALL[ALL["port_category"] == "港湾海岸"]["discharge_m3s"].sum()) total_m3s_r = float(ALL[ALL["port_category"] == "河川"]["discharge_m3s"].sum()) total_pumps = int(ALL["n_pumps"].sum()) mean_m3s = float(ALL["discharge_m3s"].mean()) median_m3s = float(ALL["discharge_m3s"].median()) max_m3s = float(ALL["discharge_m3s"].max()) min_m3s = float(ALL["discharge_m3s"].min()) n_fukuyama = int((ALL["city"] == "福山市").sum()) n_hiroshima = int((ALL["city"] == "広島市").sum()) n_showa = int((ALL["era"] == "昭和").sum()) n_heisei = int((ALL["era"] == "平成").sum()) year_min = int(ALL["year"].min()) year_max = int(ALL["year"].max()) n_top2 = 2 # 30 m3/s 機場 (手城川・岡ノ下川) top2_share = float(ALL.nlargest(2, "discharge_m3s")["discharge_m3s"].sum() / total_m3s * 100) # ============================================================================= # 11. 図 1: 2 dataset 構造マトリクス + カード # ============================================================================= print("\n[11] 図 1: dataset 構造", flush=True) t1 = time.time() 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 1269\n排水機場（港湾海岸）", ha="center", fontsize=12, fontweight="bold", color=CAT_COLOR["港湾海岸"]) ax.text(2.5, 3.05, f"{n_h} 件 / {total_m3s_h:.1f} m³/s", ha="center", fontsize=14) ax.text(2.5, 2.3, "大竹港・瀬戸田港・竹原港の海岸\n(港の陸側に降った雨水を海へ)", ha="center", fontsize=9.5) ax.text(2.5, 1.0, "全 3 件 = 全部「港海岸」\n沿岸都市の浸水防御", ha="center", fontsize=9, 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 1273\n排水機場（河川）", ha="center", fontsize=12, fontweight="bold", color=CAT_COLOR["河川"]) ax.text(7.5, 3.05, f"{n_r} 件 / {total_m3s_r:.1f} m³/s", ha="center", fontsize=14) ax.text(7.5, 2.3, "本川 (本流) ではなく\n支流・小河川の合流部に集中", ha="center", fontsize=9.5) ax.text(7.5, 1.0, f"福山市 {n_fukuyama} 件 / 広島市 {n_hiroshima} 件\n芦田川水系の干拓地特性", ha="center", fontsize=9, color="#666") ax.text(5, 0.25, f"合計 {n_total} 件 / 排水容量 {total_m3s:.1f} m³/s / ポンプ {total_pumps} 台", ha="center", fontsize=11, fontweight="bold") # 右: 排水容量バー ax = axes[1] sub = ALL.sort_values("discharge_m3s", ascending=True) y_pos = np.arange(len(sub)) colors = [CAT_COLOR[c] for c in sub["port_category"]] ax.barh(y_pos, sub["discharge_m3s"], color=colors, edgecolor="#222", linewidth=0.5) ax.set_yticks(y_pos) ax.set_yticklabels(sub["name"], fontsize=8.5) for i, v in enumerate(sub["discharge_m3s"]): ax.text(v + 0.5, i, f"{v:.1f}", va="center", fontsize=8.5) ax.set_xlabel("排水容量 (m³/s)") ax.set_title("15 機場の排水容量 — 港湾海岸 (青) / 河川 (緑)") legend_h = [Patch(facecolor=CAT_COLOR["港湾海岸"], label="港湾海岸 (3)"), Patch(facecolor=CAT_COLOR["河川"], label="河川 (12)")] ax.legend(handles=legend_h, loc="lower right") ax.grid(axis="x", alpha=0.3) plt.tight_layout() fig.savefig(ASSETS / "L35_fig1_dataset_overview.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L35_fig1_dataset_overview.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 12. 図 2: 県全域 点マップ (緯度経度をプロット) # ============================================================================= print("\n[12] 図 2: 県全域 点マップ", flush=True) t1 = time.time() fig, ax = plt.subplots(figsize=(13, 8)) pref_diss.boundary.plot(ax=ax, color="#222", linewidth=1.0) g_admin_pref.boundary.plot(ax=ax, color="#888", linewidth=0.4, alpha=0.5) # バッファ (1km) を薄く描画 for i, row in gdf.iterrows(): buf = gpd.GeoSeries([row["buffer_1km"]], crs=TARGET_CRS) buf.plot(ax=ax, color=CAT_COLOR[row["port_category"]], alpha=0.10, edgecolor="none") # 各機場を点で描画 (サイズ = 排水容量に比例) sizes = (gdf["discharge_m3s"] / gdf["discharge_m3s"].max() * 600 + 30).values for cat, color in CAT_COLOR.items(): sub = gdf[gdf["port_category"] == cat] ax.scatter(sub.geometry.x, sub.geometry.y, s=sizes[gdf["port_category"] == cat], c=color, edgecolor="#222", linewidth=0.7, alpha=0.85, label=f"{cat} ({len(sub)})") # 主要 5 機場にラベル top5 = gdf.nlargest(5, "discharge_m3s") for _, r in top5.iterrows(): ax.annotate(f"{r['name'].replace('排水機場','').replace('排水ポンプ場','').replace('ポンプ排水場','')} ({r['discharge_m3s']:.0f})", (r.geometry.x, r.geometry.y), xytext=(8, 8), textcoords="offset points", fontsize=8.5, bbox=dict(facecolor="white", alpha=0.7, edgecolor="none", pad=1.5)) ax.set_title(f"広島県 排水機場 全 {n_total} 機場マップ — 円サイズ=排水容量, 円=1km バッファ") ax.set_axis_off() ax.legend(loc="upper left") plt.tight_layout() fig.savefig(ASSETS / "L35_fig2_pref_overview.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L35_fig2_pref_overview.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 13. 図 3: 福山市集中ズーム # ============================================================================= print("\n[13] 図 3: 福山市ズーム", flush=True) t1 = time.time() # 福山市の境界を取得 (admin_932 福山市 zip があるはず) fukuyama_zip = ADMIN_DIR / "admin_832_福山市.zip" if fukuyama_zip.exists(): g_fukuyama = load_zip_first_geo(fukuyama_zip).to_crs(TARGET_CRS) fukuyama_diss = g_fukuyama.dissolve() else: # フォールバック: 福山市の機場の bbox から fukuyama_pts = gdf[gdf["city"] == "福山市"] minx, miny, maxx, maxy = fukuyama_pts.total_bounds pad = 5000 g_fukuyama = None fukuyama_diss = None fig, axes = plt.subplots(1, 2, figsize=(15, 7)) # 左: 福山市ズーム ax = axes[0] fukuyama_pts = gdf[gdf["city"] == "福山市"].copy() fukuyama_other = gdf[gdf["city"] != "福山市"] if fukuyama_diss is not None: fukuyama_diss.boundary.plot(ax=ax, color="#222", linewidth=1.5) fukuyama_diss.plot(ax=ax, color="#fef9e7", alpha=0.5) bbox = fukuyama_diss.total_bounds else: bbox = fukuyama_pts.total_bounds pad = 5000 ax.set_xlim(bbox[0] - pad, bbox[2] + pad) ax.set_ylim(bbox[1] - pad, bbox[3] + pad) # 機場の点 (サイズ = 容量) sizes_f = (fukuyama_pts["discharge_m3s"] / 30 * 800 + 50).values for cat, color in CAT_COLOR.items(): sub = fukuyama_pts[fukuyama_pts["port_category"] == cat] if len(sub) > 0: ax.scatter(sub.geometry.x, sub.geometry.y, s=sub["discharge_m3s"]/30*800 + 50, c=color, edgecolor="#222", linewidth=0.7, alpha=0.85, label=f"{cat}") # ラベル for _, r in fukuyama_pts.iterrows(): short_name = r["name"].replace("排水機場","").replace("排水ポンプ場","").replace("ポンプ排水場","") ax.annotate(f"{short_name}\n({r['discharge_m3s']:.0f} m³/s)", (r.geometry.x, r.geometry.y), xytext=(7, 7), textcoords="offset points", fontsize=8.5, bbox=dict(facecolor="white", alpha=0.75, edgecolor="none", pad=2)) ax.set_title(f"福山市 排水機場 {n_fukuyama} 件 — 干拓地+芦田川水系の集中") ax.legend(loc="upper right") ax.set_axis_off() # 右: 広島市ズーム ax = axes[1] hiroshima_zip = ADMIN_DIR / "admin_786_広島市.zip" if hiroshima_zip.exists(): g_hiroshima = load_zip_first_geo(hiroshima_zip).to_crs(TARGET_CRS) hiroshima_diss = g_hiroshima.dissolve() hiroshima_diss.boundary.plot(ax=ax, color="#222", linewidth=1.5) hiroshima_diss.plot(ax=ax, color="#fef9e7", alpha=0.5) bbox = hiroshima_diss.total_bounds else: hiroshima_pts = gdf[gdf["city"] == "広島市"] bbox = hiroshima_pts.total_bounds ax.set_xlim(bbox[0] - pad, bbox[2] + pad) ax.set_ylim(bbox[1] - pad, bbox[3] + pad) hiroshima_pts = gdf[gdf["city"] == "広島市"].copy() for cat, color in CAT_COLOR.items(): sub = hiroshima_pts[hiroshima_pts["port_category"] == cat] if len(sub) > 0: ax.scatter(sub.geometry.x, sub.geometry.y, s=sub["discharge_m3s"]/30*800 + 50, c=color, edgecolor="#222", linewidth=0.7, alpha=0.85, label=f"{cat}") for _, r in hiroshima_pts.iterrows(): short_name = r["name"].replace("排水機場","").replace("排水ポンプ場","").replace("ポンプ排水場","") ax.annotate(f"{short_name}\n({r['discharge_m3s']:.0f} m³/s)", (r.geometry.x, r.geometry.y), xytext=(7, 7), textcoords="offset points", fontsize=8.5, bbox=dict(facecolor="white", alpha=0.75, edgecolor="none", pad=2)) ax.set_title(f"広島市 排水機場 {n_hiroshima} 件 — 太田川水系の支流分布") ax.legend(loc="upper right") ax.set_axis_off() plt.tight_layout() fig.savefig(ASSETS / "L35_fig3_city_zoom.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L35_fig3_city_zoom.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 14. 図 4: 浸水想定との重ね合わせマップ # ============================================================================= print("\n[14] 図 4: 浸水想定重ね合わせ", flush=True) t1 = time.time() if FLOOD_OK: flood = gpd.read_file(FLOOD_SHP, encoding="cp932").to_crs(TARGET_CRS) # 描画高速化: simplify でジオメトリを 100m 解像度に粗化 + ランク別 dissolve # 県全域マップの場合 100m 解像度で十分 (画面 1 px ≈ 200m) flood_simple_geom = flood.copy() flood_simple_geom["geometry"] = flood_simple_geom.geometry.simplify(100, preserve_topology=False) flood_simple = flood_simple_geom.dissolve(by="rank").reset_index() fig, axes = plt.subplots(1, 2, figsize=(15, 8)) # 左: 県全域 (粗化済みポリゴン) ax = axes[0] pref_diss.boundary.plot(ax=ax, color="#222", linewidth=0.8) # 浸水ポリゴンをランク色で for _, row_f in flood_simple.iterrows(): rank = int(row_f["rank"]) gpd.GeoSeries([row_f.geometry], crs=TARGET_CRS).plot( ax=ax, color=RANK_COLORS.get(rank, "#666"), alpha=0.65, edgecolor="none") # 機場の点 for cat, color in CAT_COLOR.items(): sub = gdf[gdf["port_category"] == cat] ax.scatter(sub.geometry.x, sub.geometry.y, s=sub["discharge_m3s"]/30*400 + 30, c=color, edgecolor="white", linewidth=1.0, alpha=0.95, label=f"{cat}", zorder=10) ax.set_title("広島県 浸水想定区域 + 排水機場 — 機場が想定エリア内にあるか?") ax.legend(loc="upper left") ax.set_axis_off() # 右: 福山ズーム (浸水詳細) ax = axes[1] fukuyama_pts = gdf[gdf["city"] == "福山市"] bbox_pts = fukuyama_pts.total_bounds pad = 8000 xmin, xmax = bbox_pts[0] - pad, bbox_pts[2] + pad ymin, ymax = bbox_pts[1] - pad, bbox_pts[3] + pad ax.set_xlim(xmin, xmax); ax.set_ylim(ymin, ymax) # ローカル flood 抽出 (cx で空間フィルタ + simplify + ランク別 dissolve) flood_local = flood.cx[xmin:xmax, ymin:ymax].copy() if len(flood_local) > 0: # 福山ズームは 30m 解像度で十分 flood_local["geometry"] = flood_local.geometry.simplify(30, preserve_topology=False) flood_local_simple = flood_local.dissolve(by="rank").reset_index() for _, row_f in flood_local_simple.iterrows(): rank = int(row_f["rank"]) gpd.GeoSeries([row_f.geometry], crs=TARGET_CRS).plot( ax=ax, color=RANK_COLORS.get(rank, "#666"), alpha=0.7, edgecolor="none") if fukuyama_diss is not None: fukuyama_diss.boundary.plot(ax=ax, color="#222", linewidth=1.2) for cat, color in CAT_COLOR.items(): sub = fukuyama_pts[fukuyama_pts["port_category"] == cat] if len(sub) > 0: ax.scatter(sub.geometry.x, sub.geometry.y, s=sub["discharge_m3s"]/30*600 + 50, c=color, edgecolor="white", linewidth=1.0, alpha=0.95, label=f"{cat}", zorder=10) for _, r in fukuyama_pts.iterrows(): short_name = r["name"].replace("排水機場","").replace("排水ポンプ場","").replace("ポンプ排水場","") ax.annotate(f"{short_name}", (r.geometry.x, r.geometry.y), xytext=(7, 7), textcoords="offset points", fontsize=8, bbox=dict(facecolor="white", alpha=0.85, edgecolor="none", pad=1.5)) # 凡例 (浸水深) if len(flood_local) > 0: legend_rank = [Patch(facecolor=RANK_COLORS.get(int(r), "#666"), label=f"rank {r}: {FLOOD_RANK.get(int(r), '?')}") for r in sorted(flood_local["rank"].unique())] else: legend_rank = [] legend_rank += [Patch(facecolor=CAT_COLOR[c], label=c) for c in CAT_COLOR] ax.legend(handles=legend_rank, loc="lower left", fontsize=7.5) ax.set_title("福山市ズーム — 浸水想定区域内の機場立地") ax.set_axis_off() plt.tight_layout() fig.savefig(ASSETS / "L35_fig4_flood_overlay.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L35_fig4_flood_overlay.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 15. 図 5: 建設年代タイムライン # ============================================================================= print("\n[15] 図 5: 建設年代タイムライン", flush=True) t1 = time.time() fig, axes = plt.subplots(2, 1, figsize=(13, 7)) # 上: 機場ごとに横棒で「年代」を表示 (ガントチャート風) ax = axes[0] sub = ALL.sort_values("year") y_pos = np.arange(len(sub)) for i, (_, r) in enumerate(sub.iterrows()): color = CAT_COLOR[r["port_category"]] ax.barh(i, 0.5, left=r["year"]-0.25, color=color, alpha=0.95, edgecolor="#222") ax.scatter(r["year"], i, s=r["discharge_m3s"]/30*200+40, c=color, edgecolor="#222", linewidth=0.8, zorder=10) ax.text(r["year"] + 1, i, f"{r['name'][:10]} ({r['discharge_m3s']:.0f}m³/s)", va="center", fontsize=8.5) ax.set_yticks(y_pos) ax.set_yticklabels([f"{r['era']}{r['完成'].replace(r['era'],'').replace('年','')}年" for _, r in sub.iterrows()], fontsize=8) ax.set_xlabel("完成年 (西暦)") ax.set_xlim(1960, 2025) ax.axvline(x=1989, color="red", linestyle="--", alpha=0.5, label="平成元年 (1989)") ax.set_title("15 機場の建設年代タイムライン — 円サイズ=排水容量") ax.grid(axis="x", alpha=0.3) ax.legend(loc="lower right") # 下: 年代区分でのヒストグラム ax = axes[1] year_pivot.plot(kind="bar", stacked=True, ax=ax, color=[CAT_COLOR.get(c, "#888") for c in year_pivot.columns]) ax.set_xlabel("建設年代") ax.set_ylabel("件数") ax.set_title("建設年代別 件数 (港湾海岸 / 河川)") ax.legend(loc="upper right", title="カテゴリ") ax.grid(axis="y", alpha=0.3) plt.xticks(rotation=0) plt.tight_layout() fig.savefig(ASSETS / "L35_fig5_era_timeline.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L35_fig5_era_timeline.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 16. 図 6: 容量の Lorenz 曲線 + 散布図 (容量 × 台数) # ============================================================================= print("\n[16] 図 6: 集中構造分析", flush=True) t1 = time.time() fig, axes = plt.subplots(2, 2, figsize=(13, 10)) # (1) 容量 boxplot (カテゴリ別) ax = axes[0, 0] data_h = ALL[ALL["port_category"] == "港湾海岸"]["discharge_m3s"].values data_r = ALL[ALL["port_category"] == "河川"]["discharge_m3s"].values bp = ax.boxplot([data_h, data_r], labels=["港湾海岸", "河川"], patch_artist=True, widths=0.5) for patch, color in zip(bp["boxes"], [CAT_COLOR["港湾海岸"], CAT_COLOR["河川"]]): patch.set_facecolor(color) patch.set_alpha(0.6) ax.scatter(np.ones(len(data_h)) + np.random.uniform(-0.05, 0.05, len(data_h)), data_h, c=CAT_COLOR["港湾海岸"], edgecolor="#222", s=70, alpha=0.85, zorder=10) ax.scatter(np.ones(len(data_r))*2 + np.random.uniform(-0.05, 0.05, len(data_r)), data_r, c=CAT_COLOR["河川"], edgecolor="#222", s=70, alpha=0.85, zorder=10) ax.set_ylabel("排水容量 (m³/s)") ax.set_title("(1) カテゴリ別 排水容量分布") ax.grid(axis="y", alpha=0.3) # (2) Lorenz 曲線 ax = axes[0, 1] sorted_disc = np.sort(ALL["discharge_m3s"].values) cum_pct = np.cumsum(sorted_disc) / sorted_disc.sum() * 100 n = len(sorted_disc) x_pct = np.arange(1, n+1) / n * 100 ax.plot([0, 100], [0, 100], "--", color="gray", label="完全均等線") ax.plot(np.concatenate([[0], x_pct]), np.concatenate([[0], cum_pct]), color=CAT_COLOR["河川"], linewidth=2.5, marker="o", label="排水容量") ax.fill_between(np.concatenate([[0], x_pct]), 0, np.concatenate([[0], cum_pct]), alpha=0.15, color=CAT_COLOR["河川"]) # Gini 係数 gini = 1 - 2 * np.trapezoid(np.concatenate([[0], cum_pct]) / 100, np.concatenate([[0], x_pct]) / 100) ax.text(5, 90, f"Gini = {gini:.3f}", fontsize=11, bbox=dict(facecolor="white", alpha=0.8)) ax.text(5, 80, f"上位 2 機場 = {top2_share:.1f}%", fontsize=10, bbox=dict(facecolor="white", alpha=0.8)) ax.set_xlabel("機場ランク (低 → 高 %)") ax.set_ylabel("累積排水容量 (%)") ax.set_title("(2) 排水容量の Lorenz 曲線") ax.legend(loc="upper left") ax.grid(alpha=0.3) # (3) 容量 vs 台数 散布図 ax = axes[1, 0] for cat, color in CAT_COLOR.items(): sub = ALL[ALL["port_category"] == cat] ax.scatter(sub["n_pumps"], sub["discharge_m3s"], c=color, edgecolor="#222", s=120, alpha=0.85, label=f"{cat}") # ラベル (主要) top4 = ALL.nlargest(4, "discharge_m3s") for _, r in top4.iterrows(): ax.annotate(r["name"][:8], (r["n_pumps"], r["discharge_m3s"]), xytext=(5, 5), textcoords="offset points", fontsize=8.5, bbox=dict(facecolor="white", alpha=0.7, edgecolor="none", pad=1.5)) ax.set_xlabel("ポンプ台数 (台)") ax.set_ylabel("排水容量 (m³/s)") ax.set_title("(3) ポンプ台数 × 排水容量 散布図") ax.legend(loc="upper left") ax.grid(alpha=0.3) # (4) 1台あたり吐出量 (era 別) ax = axes[1, 1] era_data = [] era_labels = [] for era in ["昭和", "平成"]: sub = ALL[ALL["era"] == era] if len(sub) > 0: era_data.append(sub["m3s_per_pump"].values) era_labels.append(f"{era} (n={len(sub)})") bp2 = ax.boxplot(era_data, labels=era_labels, patch_artist=True, widths=0.5) for patch, color in zip(bp2["boxes"], ["#9467bd", "#d62728"]): patch.set_facecolor(color) patch.set_alpha(0.6) for i, d in enumerate(era_data): ax.scatter(np.ones(len(d))*(i+1) + np.random.uniform(-0.05, 0.05, len(d)), d, c="#222", s=50, alpha=0.7, zorder=10) ax.set_ylabel("1 台あたり吐出量 (m³/s/台)") ax.set_title("(4) 昭和 vs 平成 — 1 台の能力比較 (大型化検証)") ax.grid(axis="y", alpha=0.3) plt.tight_layout() fig.savefig(ASSETS / "L35_fig6_distribution.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L35_fig6_distribution.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 17. 図 7: 浸水想定 1km 圏ヒートマップ + 機場別の最大ランク # ============================================================================= print("\n[17] 図 7: 浸水想定圏分析", flush=True) t1 = time.time() fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # (1) 機場ごとの最大浸水ランク + 1km 内浸水面積 ax = axes[0] sub = gdf.sort_values("flood_area_ha_in_buffer", ascending=True).copy() y_pos = np.arange(len(sub)) colors_rank = [RANK_COLORS.get(int(r), "#bbb") for r in sub["max_rank_in_buffer"]] ax.barh(y_pos, sub["flood_area_ha_in_buffer"], color=colors_rank, edgecolor="#222") for i, (v, rank) in enumerate(zip(sub["flood_area_ha_in_buffer"], sub["max_rank_in_buffer"])): rank_str = FLOOD_RANK.get(int(rank), "—") if rank > 0 else "圏外" ax.text(v + 5, i, f"{v:.0f} ha / 最大 {rank_str}", va="center", fontsize=8) ax.set_yticks(y_pos) ax.set_yticklabels([f"{r['name'][:10]} ({r['city']})" for _, r in sub.iterrows()], fontsize=8.5) ax.set_xlabel("1km バッファ内の浸水想定面積 (ha)") ax.set_title("各機場 1km 圏内の浸水想定面積 — 守備範囲の規模") ax.grid(axis="x", alpha=0.3) # (2) 容量 vs 浸水面積 散布図 (機場の能力と守備範囲の整合) ax = axes[1] mask = gdf["flood_area_ha_in_buffer"] > 0 sub_plot = gdf[mask] for cat, color in CAT_COLOR.items(): s2 = sub_plot[sub_plot["port_category"] == cat] if len(s2) > 0: ax.scatter(s2["flood_area_ha_in_buffer"], s2["discharge_m3s"], c=color, edgecolor="#222", s=120, alpha=0.85, label=f"{cat}") # ラベル for _, r in sub_plot.iterrows(): if r["discharge_m3s"] >= 5 or r["flood_area_ha_in_buffer"] >= 100: ax.annotate(r["name"][:8], (r["flood_area_ha_in_buffer"], r["discharge_m3s"]), xytext=(5, 5), textcoords="offset points", fontsize=8, bbox=dict(facecolor="white", alpha=0.7, edgecolor="none", pad=1.5)) ax.set_xlabel("1km 圏内浸水想定面積 (ha)") ax.set_ylabel("排水容量 (m³/s)") ax.set_title("守備面積 × 能力の対応関係") ax.legend(loc="lower right") ax.grid(alpha=0.3) plt.tight_layout() fig.savefig(ASSETS / "L35_fig7_flood_buffer.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L35_fig7_flood_buffer.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 18. 図 8: 事務所別管理 + 河川水系別 # ============================================================================= print("\n[18] 図 8: 事務所/水系", flush=True) t1 = time.time() # 事務所別 (上位の事務所) office_pv = ALL.groupby("office").agg( n=("name", "count"), total_m3s=("discharge_m3s", "sum"), ).sort_values("n", ascending=True) office_pv.to_csv(ASSETS / "L35_office_summary.csv", encoding="utf-8-sig") # 水系/海岸名でグルーピング (本データの「海岸・河川名」) wc_pv = ALL.groupby("watercourse").agg( n=("name", "count"), total_m3s=("discharge_m3s", "sum"), ).sort_values("total_m3s", ascending=True) wc_pv.to_csv(ASSETS / "L35_watercourse_summary.csv", encoding="utf-8-sig") fig, axes = plt.subplots(1, 2, figsize=(14, 6)) ax = axes[0] y = np.arange(len(office_pv)) ax.barh(y, office_pv["n"], color="#0969da", alpha=0.85, edgecolor="#222") for i, (n, m) in enumerate(zip(office_pv["n"], office_pv["total_m3s"])): ax.text(n + 0.1, i, f"{int(n)} 件 / {m:.1f} m³/s", va="center", fontsize=9) ax.set_yticks(y) ax.set_yticklabels(office_pv.index, fontsize=9) ax.set_xlabel("管理機場数") ax.set_title("管理事務所別 機場数") ax.grid(axis="x", alpha=0.3) ax = axes[1] y = np.arange(len(wc_pv)) ax.barh(y, wc_pv["total_m3s"], color="#1a7f37", alpha=0.85, edgecolor="#222") for i, (n, m) in enumerate(zip(wc_pv["n"], wc_pv["total_m3s"])): ax.text(m + 0.3, i, f"{int(n)} 件 / {m:.1f} m³/s", va="center", fontsize=9) ax.set_yticks(y) ax.set_yticklabels(wc_pv.index, fontsize=9) ax.set_xlabel("総排水容量 (m³/s)") ax.set_title("海岸・河川名別 排水容量") ax.grid(axis="x", alpha=0.3) plt.tight_layout() fig.savefig(ASSETS / "L35_fig8_office_watercourse.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L35_fig8_office_watercourse.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 19. 図 9: 機場プロファイル レーダー (4 軸) # ============================================================================= print("\n[19] 図 9: 機場プロファイル", flush=True) t1 = time.time() # 上位 6 機場のレーダーチャート (4 軸: 容量 / 台数 / 1台あたり / 浸水面積) top6 = gdf.nlargest(6, "discharge_m3s").copy() axes_labels = ["排水容量\n(m³/s)", "ポンプ台数\n(台)", "1台あたり\n(m³/s/台)", "1km圏\n浸水(ha)"] # 各列を 0-1 正規化 (全 15 件中の最大で割る) norm_data = pd.DataFrame({ "容量": top6["discharge_m3s"] / gdf["discharge_m3s"].max(), "台数": top6["n_pumps"] / gdf["n_pumps"].max(), "1台あたり": top6["m3s_per_pump"] / gdf["m3s_per_pump"].max(), "浸水": top6["flood_area_ha_in_buffer"] / max(1, gdf["flood_area_ha_in_buffer"].max()), }) fig, axes = plt.subplots(2, 3, figsize=(14, 9), subplot_kw=dict(polar=True)) angles = np.linspace(0, 2 * np.pi, 4, endpoint=False).tolist() angles += angles[:1] for ax, (i, row) in zip(axes.flat, top6.iterrows()): name = row["name"][:14] values = norm_data.loc[i].tolist() values += values[:1] color = CAT_COLOR[row["port_category"]] ax.plot(angles, values, color=color, linewidth=2) ax.fill(angles, values, color=color, alpha=0.25) ax.set_xticks(angles[:-1]) ax.set_xticklabels(axes_labels, fontsize=8) ax.set_ylim(0, 1) ax.set_yticks([0.25, 0.5, 0.75, 1.0]) ax.set_yticklabels(["", "0.5", "", "1.0"], fontsize=7) ax.set_title(f"{name}\n{row['city']} / {row['完成']}", fontsize=10) plt.tight_layout() fig.savefig(ASSETS / "L35_fig9_radar_top6.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L35_fig9_radar_top6.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 20. 仮説検証 # ============================================================================= print("\n[20] 仮説検証", flush=True) t1 = time.time() hypos = [] # H1: 河川主導 ratio_count = n_r / max(1, n_h) ratio_capacity = total_m3s_r / max(0.01, total_m3s_h) h1 = (n_r >= n_h * 3) and (total_m3s_r >= total_m3s_h * 5) hypos.append({ "H": "H1", "claim": "河川 12 件 (80%) で主導、容量比でも 河川:港湾海岸 ≈ 9:1", "result": f"河川 {n_r}/{n_total} = {n_r/n_total*100:.0f}%, 容量 {total_m3s_r:.1f}/{total_m3s:.1f} = {total_m3s_r/total_m3s*100:.0f}%, 件数比 {ratio_count:.1f}:1, 容量比 {ratio_capacity:.1f}:1", "verdict": "支持" if h1 else ("部分支持" if ratio_count >= 2 else "反証"), }) # H2: 福山市偏重 fukuyama_pct = n_fukuyama / n_r * 100 h2 = (fukuyama_pct >= 50) hypos.append({ "H": "H2", "claim": "河川機場 12 件のうち福山市が 50%+ を占有", "result": f"福山市 {n_fukuyama}/{n_r} = {fukuyama_pct:.0f}% (広島市 {n_hiroshima}, 尾道市 {int((ALL['city']=='尾道市').sum())}, 竹原市 {int((ALL['city']=='竹原市').sum())})", "verdict": "支持" if h2 else ("部分支持" if fukuyama_pct >= 30 else "反証"), }) # H3: 大型機場の二極集中 (上位 2 機場で 40%+) h3 = (top2_share >= 40) hypos.append({ "H": "H3", "claim": "30 m³/s 機場 2 件 (手城川・岡ノ下川) で県全体の 40%+", "result": f"上位 2 機場 {top2_share:.1f}%, Gini = {gini:.3f} (高度な集中)", "verdict": "支持" if h3 else "部分支持", }) # H4: 建設年代の 2 ピーク (1965-80, 2008-16) n_1960s_70s = int(((ALL["year"] >= 1960) & (ALL["year"] < 1985)).sum()) n_2008_16 = int(((ALL["year"] >= 2008) & (ALL["year"] <= 2016)).sum()) n_other = n_total - n_1960s_70s - n_2008_16 h4 = (n_1960s_70s >= 3) and (n_2008_16 >= 3) hypos.append({ "H": "H4", "claim": "建設に 2 ピーク: 1965-84 (戦後復興・都市化期) と 2008-16 (集中豪雨対応強化期)", "result": f"1965-84: {n_1960s_70s} 件, 2008-16: {n_2008_16} 件, その他: {n_other} 件", "verdict": "支持" if h4 else ("部分支持" if (n_1960s_70s >= 2) else "反証"), }) # H5: 浸水想定 1km 圏内の機場 ≥ 87% flood_in_pct = (gdf["near_flood"].sum() / len(gdf)) * 100 if FLOOD_OK else 0 h5 = (flood_in_pct >= 87) hypos.append({ "H": "H5", "claim": "全 15 機場のうち ≥ 87% が浸水想定 1km バッファ内に立地", "result": (f"圏内 {int(gdf['near_flood'].sum())}/{len(gdf)} = {flood_in_pct:.1f}%" if FLOOD_OK else "浸水shp 不在のため未検証"), "verdict": "支持" if h5 else ("部分支持" if flood_in_pct >= 60 else ("反証" if FLOOD_OK else "未検証")), }) # H6: 大型化 (平成 > 昭和 で 1台あたり吐出量) showa_per_pump = float(ALL[ALL["era"] == "昭和"]["m3s_per_pump"].mean()) heisei_per_pump = float(ALL[ALL["era"] == "平成"]["m3s_per_pump"].mean()) h6 = (heisei_per_pump > showa_per_pump) hypos.append({ "H": "H6", "claim": "平成期建設機の 1 台あたり吐出量は昭和期より大きい (機械大型化)", "result": f"昭和 平均 {showa_per_pump:.2f} m³/s/台, 平成 平均 {heisei_per_pump:.2f} m³/s/台", "verdict": "支持" if h6 else "反証", }) hypos_df = pd.DataFrame(hypos) hypos_df.to_csv(ASSETS / "L35_hypothesis_results.csv", index=False, encoding="utf-8-sig") with (ASSETS / "L35_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) # ============================================================================= # 21. HTML 生成 # ============================================================================= print("\n[21] HTML 生成", flush=True) t1 = time.time() def dl(name): return f'{name}' sections = [] # === Section 1: 学習目標と問い === sec1 = f"""
カバー宣言: 本記事は DoBoX のシリーズ 「排水機場 港湾海岸 / 河川」 2 件 (dataset_id = 1269, 1273) を統合し、広島県の排水機場の地理分布・ポンプ吐出能力・ 建設年代・浸水想定との重なりを分析する研究記事です。 本シリーズは港湾系シリーズ (L32-L34) の最終記事ですが、扱う施設は 港湾施設ではなく「内水排水施設」であり、 港湾海岸 ({n_h} 件) と河川 ({n_r} 件) の 2 つの設置目的を持ちます。

シリーズ構造判定: (a) 設置目的による 2 分割型 + (b) 完全相補

このシリーズの 2 件は、「排水機場の管轄」で 2 分されている相補型構造です。 L32-L34 の 「港湾 / 漁港」分割と異なり、「港湾海岸 / 河川」という 内水氾濫の主因に基づく分類です。

L32-L34 (港湾施設の 3 層) と同じ「港」を見る記事ではない点に注意。 L32-L34 は港湾の外洋・接岸・陸側の 3 層を扱ったが、 L35 は陸側にさらに溜まった雨水・氾濫水をポンプで汲み出す施設を扱う。 これは「内水浸水防御の最終線」であり、外郭・係留・交通とは目的が根本的に異なる。

独自用語の定義 (本記事内で固定)

用語定義
排水機場集水区の雨水・地下水・河川氾濫水をポンプで強制排水する施設。 内水氾濫を防ぐ最後の砦。本データに {n_total} 件。 建屋 (ポンプ室) + ポンプ {total_pumps} 台 + 吐出管 + 吸込水路 + 吐出水路 + 樋門 で構成される複合施設。
内水氾濫河川本流が氾濫しなくても、支流・側溝・下水管路の容量を超えた雨水が 低地に溜まる現象。本流の水位が高くて支流の水を受け入れられない (=堰き止め) ときに発生。 排水機場はこの内水を強制排水する。外水氾濫 (本流の決壊・越水) とは別物で、 排水機場は内水専用の防御施設。
吐出量 (m³/s)ポンプが 1 秒間に排水できる水量。本データの「総排水量(m3)」列は 表記上 m³ だが、値の桁 (1.3〜30) と排水機場の業界慣例から毎秒の吐出量 m³/sと解釈。 30 m³/s = 25 m プール 1 杯 (600 m³) を 20 秒で空にする能力。 1 m³/s = 約 1 トン/秒 = 大型ダンプカー 1 台/秒分の水量。
干拓地もとは海・湖だった場所を堤防で囲って水を抜き、陸地化した土地。 広島県では福山市芦田川河口の福山平野が代表例で、江戸期からの干拓により 広大な低地 (海抜 0 m 前後) が形成。排水機場が必須な土地条件で、 本記事で福山市に機場が集中する理由。
強制排水自然流下で水が引かない場合に、ポンプで動力的に水を持ち上げて排水すること。 河川本流の水位が支流より高い時 (= 自然流下不能)、 あるいは満潮時に港湾海岸が排水できない時に作動する。 電力 + 設備 + 維持管理が必要で、「最後の手段」として整備される。
守備範囲 (1 km バッファ)本記事独自指標。各機場が排水を担う集水区を 正確に特定するには排水計画書が必要だが、 本記事では機場から半径 1 km 圏を「物理的に近い守備候補範囲」と仮定し、 この圏内の浸水想定面積で機場の仕事量を測る。これは概算指標であり、 正確な集水区は別途取得が必要 (発展課題)。
浸水想定区域DoBoX が公開する河川浸水想定 Shapefile (39 河川統合版)。 1000 年確率規模降雨時の浸水域とその深さ (rank 10〜80 = 0.5m 未満〜20m 以上) を示す。 排水機場は本来この想定エリア内の低標高点に置かれるべきで、 機場立地と浸水想定の重なりは工学的整合性の指標。

研究の問い (RQ)

広島県の排水機場 ({n_total} 件) は、設置目的 (港湾海岸 {n_h} / 河川 {n_r})・ 地理分布 (沿岸低地 / 河川合流部 / 干拓地)・規模 (吐出量 m³/s)・建設年代の観点で どう構成され、内水浸水想定とどう重なるのか?

  1. 2 dataset の構造比は? 河川主導が想定されるが、容量比でも本当に河川偏重か?
  2. 市町別の偏在 (福山市・広島市) は何を反映しているか? 干拓地・水系の特性?
  3. 排水容量はどう分布するか? 上位機場への集中度は?
  4. 建設年代の傾向は? 戦後の都市化期・近年の集中豪雨対応で集中ピークがあるか?
  5. 機場は浸水想定区域内に立地しているか? 位置の工学的整合性?
  6. 1 台あたり吐出量 (機械の大型化) は時代と共に増えているか?

仮説 H1〜H6

到達点

2 dataset を「港湾海岸 + 河川」の設置目的別の 2 カテゴリとして読み解き、 {n_total} 件の排水機場の地理分布・容量分布・建設年代・浸水想定との重なりを統合分析する。 これにより、広島県の内水浸水防御戦略が 「河川合流部での強制排水を主、港湾海岸でのサブ排水を従」とする構成であり、 かつ福山平野の干拓地に重点投資されていることを実データで裏付ける。 さらに、上位 2 機場が県全体の半分近くを担う強い集中構造と、 機械の大型化トレンドを統計的に示す。

本記事のスコープ外: 本記事は機場の立地構造と容量に集中する。 ポンプの消費電力・運転履歴・故障率、 集水区の正確な範囲と排水負荷曲線、 洪水時の運転実績などは本記事では扱わず発展課題とする。 また「総排水量(m3)」の単位については業界慣例の m³/s と解釈するが、 これは管理者の運用記録 (運転日誌等) で確認すべき (発展課題 Z6)。
""" sections.append(("学習目標と問い", sec1)) # === Section 2: 使用データ === sec2_table = "" for r in series_df.itertuples(): n_rows = int(ALL[ALL["dsid"] == r.dsid].shape[0]) sec2_table += ( f"" f"" f"" f"" ) sec2_table += "
dsidカテゴリ件数形式列数DoBoX
{r.dsid}{r.category}{n_rows}CSV (UTF-8 BOM)9#{r.dsid}
" # 列スキーマ表 schema = [ ("施設の名称", "機場の固有名称 (例: 御幸排水機場、岡ノ下川排水機場)"), ("事務所名", "管理事務所 (例: 西部建設事務所、東部建設事務所、各支所)"), ("所在地", "市町名 + 番地 (例: 福山市東手城町二丁目)"), ("緯度(10進数)", "緯度 (10進数表記)"), ("経度(10進数)", "経度 (10進数表記)"), ("海岸・河川名", "排水先または接続水系の名称 (例: 岡ノ下川、大竹港海岸)"), ("完成年月日", "和暦の年 (例: 昭和41年、平成28年)"), ("総排水量(m3)", "業界慣例で m³/s と解釈。1.3〜30 の範囲"), ("ポンプ台数(台)", "ポンプの基数 (1〜4 台)"), ] schema_html = "" + "".join( f"" for c, m in schema ) + "
列名意味
{c}{m}
" # 全 15 件の生データ表 all_table = ALL[["name", "port_category", "city", "watercourse", "完成", "year", "discharge_m3s", "n_pumps", "m3s_per_pump"]].copy() all_table = all_table.rename(columns={ "name": "機場名", "port_category": "カテゴリ", "city": "市町", "watercourse": "海岸/河川", "完成": "完成", "year": "西暦", "discharge_m3s": "容量m³/s", "n_pumps": "台数", "m3s_per_pump": "台/台m³/s", }) all_table["容量m³/s"] = all_table["容量m³/s"].round(1) all_table["台/台m³/s"] = all_table["台/台m³/s"].round(2) all_table_html = all_table.to_html(index=False, classes="", border=0) sec2 = f"""

本記事が使用する 2 dataset_id の一覧。両者は同一スキーマ (9 列) で、 データセット名の括弧内 (港湾海岸 / 河川) で設置目的が分かれます。

{sec2_table}

列スキーマ (両 dataset 共通)

{schema_html}

L32-L34 シリーズとの重要な相違

軸L32 外郭 / L33 係留 / L34 交通L35 排水機場 (本記事)
件数合計 1,520 件{n_total} 件 (極小)
地理分布瀬戸内海岸の港湾区域港湾海岸 + 河川合流部 (内陸も含む)
主指標延長 m / 面積 m²吐出量 m³/s + ポンプ台数
カテゴリ軸港湾 / 漁港 (運営主体)港湾海岸 / 河川 (排水先)
GISWKT POLYGON / LINESTRING緯度経度のみ (Point)
分析の主役3 層空間関係容量分布・建設年代・浸水想定との重なり

15 機場 全件リスト (本記事のデータ全部)

{all_table_html}

データ品質メモ

""" 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'
  • {r.label} (dsid={r.dsid}, {r.category}): ' f'直 DL ' f'/ DoBoX
    • ' ) sec3 = f"""

    生データ (DoBoX 直リンク)

    中間データ (本記事生成 CSV)

    図 (本記事生成 PNG)

    再現用 Python スクリプト

    L35_drainage_pumps.py を取得して プロジェクトルートで py -X utf8 lessons/L35_drainage_pumps.py を実行。 データが無ければ自動取得します (DoBoX 浸水想定 shp が data/extras/flood_shp/ にあれば浸水重ね分析も自動有効化)。

    """ sections.append(("ダウンロード", sec3)) # === Section 4: 分析 1 — 2 dataset の構造 === sec4 = f"""

    狙い

    「港湾海岸」と「河川」の 2 dataset が、件数・容量・地理分布でどう分化しているかを 1 枚の絵で示す。河川主導 (H1) の量的検証。

    手法 (簡潔に)

    2 dataset を縦結合し、port_category 列 (港湾海岸 / 河川) で 件数・総容量・台数を集計。15 件全部を 1 枚の横棒グラフで容量順に並べ、 カテゴリ色で塗り分けて視覚的に二分構造を示す。

    実装

    {code(''' SERIES = [ (1269, "排水機場（港湾海岸）", "港湾海岸", "drainage_pump_port_coast.csv", 39770), (1273, "排水機場（河川）", "河川", "drainage_pump_river.csv", 39774), ] 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 dfs.append(df) ALL = pd.concat(dfs, ignore_index=True) # 全 15 行 cat_agg = ALL.groupby("port_category").agg( n=("施設の名称", "count"), total_m3s=("総排水量(m3)", "sum"), n_pumps_sum=("ポンプ台数(台)", "sum"), ) ''')}

    図と読み取り

    なぜこの図か: 2 dataset の規模を文字 + バーで同時に伝える。 左 (カード) は「カテゴリ別の件数+容量+概要」をテキストで、 右 (バー) は「15 機場の容量を一望」できる横棒。両方を 1 枚に置くことで 「カテゴリの偏重 (H1)」と「容量分布の偏在 (H3)」が同時に読める。

    {figure("assets/L35_fig1_dataset_overview.png", "2 dataset の構造概観 — カードビュー (左) と 15 機場の排水容量バー (右)")}

    読み取り:

    表と読み取り (カテゴリ別集計)

    {cat_agg.to_html(classes='', border=0)}

    読み取り: 件数の偏重 (河川 80%) より容量の偏重 (河川 90%) の方が強い。 これは河川機場が大型機を持ちやすいことを反映 (合流部での大流量を扱うため)。 平均値で見ると 河川 {ALL[ALL['port_category']=='河川']['discharge_m3s'].mean():.1f} m³/s vs 港湾海岸 {ALL[ALL['port_category']=='港湾海岸']['discharge_m3s'].mean():.1f} m³/s で約 {ALL[ALL['port_category']=='河川']['discharge_m3s'].mean()/ALL[ALL['port_category']=='港湾海岸']['discharge_m3s'].mean():.1f} 倍。

    """ sections.append(("分析 1: 2 dataset の構造を可視化", sec4)) # === Section 5: 分析 2 — 県全域マップで地理分布 === sec5 = f"""

    狙い

    件数や容量だけでは見えない地理偏在を、県全域の点マップで一望する。 広島県のどの市町に排水機場が分布し、どこが手薄かを可視化。

    手法 (リテラシレベル解説)

    緯度経度 → Point geometry → matplotlib scatter の 3 ステップ:

    1. 各機場の緯度・経度から shapely.Point を作る。
    2. EPSG:4326 (WGS84) → EPSG:6671 (平面直角 III 系、m 単位) に投影変換。 これで距離計算とバッファ生成が m 単位で正確になる。
    3. matplotlib scatter で点を描画、サイズ = 排水容量に比例 (バブル散布)。 1 km バッファも薄塗りで重ね、「機場が守る範囲の目安」を視覚化。

    実装

    {code(''' from shapely.geometry import Point ALL["geometry"] = [Point(lon, lat) for lon, lat in zip(ALL["lon"], ALL["lat"])] gdf = gpd.GeoDataFrame(ALL, geometry="geometry", crs="EPSG:4326").to_crs("EPSG:6671") gdf["buffer_1km"] = gdf.geometry.buffer(1000.0) # 1 km バッファ fig, ax = plt.subplots(figsize=(13, 8)) pref_diss.boundary.plot(ax=ax, color="#222", linewidth=1.0) # サイズ = 容量に比例 sizes = gdf["discharge_m3s"] / gdf["discharge_m3s"].max() * 600 + 30 for cat, color in CAT_COLOR.items(): sub = gdf[gdf["port_category"] == cat] ax.scatter(sub.geometry.x, sub.geometry.y, s=sizes[gdf.port_category==cat], c=color, edgecolor="#222", alpha=0.85) ''')}

    図と読み取り

    なぜこの図か: 県全域に対する 15 機場の分布密度を、青 (港湾海岸) と緑 (河川) で 直接見ることで「どこにどの目的の機場が集中しているか」を一目で把握できる。 バッファ円で守備範囲の目安も同時に示される。

    {figure("assets/L35_fig2_pref_overview.png", f"広島県 排水機場 全 {n_total} 件マップ — 円サイズ=排水容量、円=1km バッファ")}

    読み取り:

    """ sections.append(("分析 2: 県全域マップで地理偏在を観る", sec5)) # === Section 6: 分析 3 — 福山市・広島市の集中ズーム === sec6 = f"""

    狙い

    2 大集中地 (福山市 {n_fukuyama} 件 + 広島市 {n_hiroshima} 件) を市単位でズームして、 機場の分布パターンと集中の理由を地形・水系から読み解く。 H2 (福山市偏重) の地理的検証。

    手法

    市町境界 Shapefile (admin_832 福山市 / admin_786 広島市) に 機場点を重ね、サイズ = 容量のバブル + ラベルを配置。 2 市を並列に並べることで分布形状の違いを比較できる。

    実装

    {code(''' fukuyama_zip = ADMIN_DIR / "admin_832_福山市.zip" g_fukuyama = load_zip_first_geo(fukuyama_zip).to_crs("EPSG:6671") fukuyama_diss = g_fukuyama.dissolve() # 市境 1 ポリゴン fig, axes = plt.subplots(1, 2, figsize=(15, 7)) # 福山市 ax = axes[0] fukuyama_diss.plot(ax=ax, color="#fef9e7", alpha=0.5) fukuyama_pts = gdf[gdf["city"] == "福山市"] ax.scatter(fukuyama_pts.geometry.x, fukuyama_pts.geometry.y, s=fukuyama_pts["discharge_m3s"]/30*800 + 50, c=...) ''')}

    図と読み取り

    なぜこの図か: 県全域マップでは見えない市町内の分布を市境ポリゴン上に拡大し、 同じ市内でも機場が集中するエリアと空白エリアが分かる。 2 市を並べる small multiples 配置で対比が直観的。

    {figure("assets/L35_fig3_city_zoom.png", "福山市 (左) と広島市 (右) の機場ズーム — サイズ=容量")}

    読み取り:

    """ sections.append(("分析 3: 福山市・広島市ズームで集中構造を読む", sec6)) # === Section 7: 分析 4 — 浸水想定との重ね合わせ === near_flood_n = int(gdf["near_flood"].sum()) sec7 = f"""

    狙い

    機場の立地は浸水想定区域内にあるか? 工学的には機場は浸水域の最低標高点に 置かれるべきで、機場立地と浸水想定の重なりは工学的整合性の指標。 H5 の検証。

    手法 (リテラシレベル解説)

    空間結合 + バッファ判定:

    1. DoBoX 河川浸水想定区域 Shapefile (39 河川統合版、613 ポリゴン) をロード。 rank 列が浸水深 (10=0.5m未満 〜 80=20m以上)。
    2. 各機場に 1 km バッファを生成 (geometry.buffer(1000))。
    3. バッファと浸水想定ポリゴンが交差するか gpd.GeoDataFrame.intersects() で判定。 1 件以上交差すれば「浸水想定圏内立地」と判定。
    4. 各機場の 1km 圏内の浸水想定面積 (ha) を gpd.clip() で算出 (= 守備規模の指標)。

    実装

    {code(''' flood = gpd.read_file("data/extras/flood_shp/shinsui_souteisaidai/shinsui_souteisaidai.shp", encoding="cp932").to_crs("EPSG:6671") gdf["buffer_1km"] = gdf.geometry.buffer(1000.0) for i, row in gdf.iterrows(): sub = flood[flood.intersects(row["buffer_1km"])] if len(sub) > 0: gdf.at[i, "near_flood"] = True gdf.at[i, "max_rank_in_buffer"] = sub["rank"].max() clipped = gpd.clip(sub, row["buffer_1km"]) gdf.at[i, "flood_area_ha_in_buffer"] = clipped.geometry.area.sum() / 10000 ''')}

    図と読み取り — 県全域 + 福山ズーム

    なぜこの図か: 浸水想定 (青階調) と機場 (色マーカー) を同一マップに重ねることで、 「機場が浸水想定エリア内に立地するか」が直視できる。 県全域 (左) で全体像、福山ズーム (右) で詳細を見る。

    {figure("assets/L35_fig4_flood_overlay.png", "浸水想定 + 排水機場 — 機場が想定エリア内に立地するかの視覚検証")}

    読み取り:

    図と読み取り — 機場別 1km 圏内浸水面積ランキング

    なぜこの図か: 機場ごとの「守るべき面積」を 1km 圏内の浸水想定面積で代替指標化し、 バーで順位付け。機場容量と守備面積の整合を後段で検証する基礎データ。

    {figure("assets/L35_fig7_flood_buffer.png", "各機場 1km 圏内の浸水想定面積 (左) と容量×面積散布 (右)")}

    読み取り (左バー):

    読み取り (右散布図):

    解釈の注意: 1 km バッファは近接の目安であり、 実際の集水区は機場の排水計画書で個別定義される。 正確な守備関係を知るには、機場ごとの排水区域図を取得する必要がある (発展課題 Z3)。
    """ sections.append(("分析 4: 浸水想定との重なり — 機場立地の工学的整合性", sec7)) # === Section 8: 分析 5 — 容量分布と集中構造 === sec8 = f"""

    狙い

    排水機場の容量分布を多角的に観察し、 カテゴリ別の標準サイズ・集中度・台数との関係・新旧の機械大型化を統計的に検証する。 H3 (上位 2 機場集中) と H6 (機械大型化) の検証。

    手法 (リテラシレベル解説)

    4 つの可視化手法を 2x2 で統合する:

    実装

    {code(''' # (2) Lorenz 曲線と Gini 係数 sorted_disc = np.sort(ALL["discharge_m3s"].values) cum_pct = np.cumsum(sorted_disc) / sorted_disc.sum() * 100 n = len(sorted_disc) x_pct = np.arange(1, n+1) / n * 100 gini = 1 - 2 * np.trapezoid(np.concatenate([[0], cum_pct]) / 100, np.concatenate([[0], x_pct]) / 100) # Gini = 0 (均等) 〜 1 (1 機独占)。容量で 0.4-0.5 なら強い集中。 # (4) 1 台あたり吐出量 ALL["m3s_per_pump"] = ALL["discharge_m3s"] / ALL["n_pumps"] era_pump = ALL.groupby("era")["m3s_per_pump"].mean() ''')}

    図と読み取り

    なぜこの図か: 4 視点を 2x2 で並べることで、分布形状・集中度・関係性・時系列変化を 1 枚に統合できる。

    {figure("assets/L35_fig6_distribution.png", "容量分布 4 視点 — boxplot、Lorenz、容量×台数、新旧 1台あたり")}

    読み取り (左上 — カテゴリ別 boxplot):

    読み取り (右上 — Lorenz 曲線):

    読み取り (左下 — 容量 × 台数 散布図):

    読み取り (右下 — 昭和 vs 平成 1台あたり吐出量):

    """ sections.append(("分析 5: 容量分布と集中構造 — Gini 係数 + 機械大型化検証", sec8)) # === Section 9: 分析 6 — 建設年代タイムライン === sec9 = f"""

    狙い

    機場の建設年代を時系列で並べ、整備の歴史を読み解く。 昭和40 ({year_min}) 〜 平成28 ({year_max}) の {year_max-year_min} 年間で 15 件が建設された。 H4 (2 ピーク仮説: 1965-84 戦後復興期 + 2008-16 集中豪雨対応期) の検証。

    手法

    和暦 (例: 昭和41年) を西暦 (1966) に変換し、横軸=西暦のガントチャート風タイムラインを描く。 バブルサイズで容量、色でカテゴリを表現。下段に10 年区分のヒストグラムを並べて 集中時期を統計的に確認。

    実装

    {code(''' def jp_to_west(s: str) -> int | None: s = str(s).strip() if s.startswith("昭和"): n = int(s.replace("昭和", "").replace("年", "")) return 1925 + n # 昭和元年=1926, n年 = 1925+n elif s.startswith("平成"): n = int(s.replace("平成", "").replace("年", "")) return 1988 + n # 平成元年=1989, n年 = 1988+n return None ALL["year"] = ALL["完成"].apply(jp_to_west) ALL["era"] = ALL["完成"].apply(lambda s: "昭和" if str(s).startswith("昭和") else "平成") # 10 年区分のヒストグラム year_bins = [1960, 1970, 1980, 1990, 2000, 2010, 2020] year_labels = ["1960s", "1970s", "1980s", "1990s", "2000s", "2010s"] ALL["year_bin"] = pd.cut(ALL["year"], bins=year_bins, labels=year_labels, right=False) year_pivot = ALL.groupby(["year_bin", "port_category"], observed=True).size().unstack(fill_value=0) ''')}

    図と読み取り

    なぜこの図か: タイムライン (上) で各機場の建設順序と容量、ヒストグラム (下) で集中時期を 2 段に並べることで「いつ・誰が・どれくらい大きい機場を建てたか」が読める。

    {figure("assets/L35_fig5_era_timeline.png", "建設年代タイムライン (上) と 10 年区分ヒストグラム (下)")}

    読み取り:

    """ sections.append(("分析 6: 建設年代タイムライン — 整備の歴史と 2 ピーク仮説", sec9)) # === Section 10: 分析 7 — 事務所別管理 + 海岸/河川名別 === sec10 = f"""

    狙い

    機場の管轄行政組織 (建設事務所) と排水先水系 (海岸・河川名) で集計し、 誰がどこを守っているかの組織構造を読む。

    手法

    2 軸クロス集計とバー比較。事務所軸は管理体制 (人的責任)、 水系軸は地理的責任を反映する。

    実装

    {code(''' office_pv = ALL.groupby("office").agg( n=("name", "count"), total_m3s=("discharge_m3s", "sum"), ).sort_values("n", ascending=True) wc_pv = ALL.groupby("watercourse").agg( n=("name", "count"), total_m3s=("discharge_m3s", "sum"), ).sort_values("total_m3s", ascending=True) ''')}

    図と読み取り

    なぜこの図か: 同じ 15 機場を 2 つの集計軸で並列に見ることで 「組織責任 (誰)」と「地理責任 (どこ)」の両面が読める。

    {figure("assets/L35_fig8_office_watercourse.png", "管理事務所別 (左) と 海岸/河川名別 (右) — 組織と地理の責任分担")}

    読み取り (左バー — 事務所別):

    読み取り (右バー — 水系別):

    """ sections.append(("分析 7: 事務所別管理と水系別配置", sec10)) # === Section 11: 分析 8 — 上位機場プロファイル レーダー === sec11 = f"""

    狙い

    容量上位 6 機場の4 軸プロファイルをレーダーチャートで並べ、 各機場の特性 (バランス型 vs 偏った型) を比較する。

    手法 (リテラシレベル解説)

    レーダーチャート (Radar Chart):

    実装

    {code(''' top6 = gdf.nlargest(6, "discharge_m3s") norm_data = pd.DataFrame({{ "容量": top6["discharge_m3s"] / gdf["discharge_m3s"].max(), "台数": top6["n_pumps"] / gdf["n_pumps"].max(), "1台あたり": top6["m3s_per_pump"] / gdf["m3s_per_pump"].max(), "浸水": top6["flood_area_ha_in_buffer"] / gdf["flood_area_ha_in_buffer"].max(), }}) fig, axes = plt.subplots(2, 3, figsize=(14, 9), subplot_kw=dict(polar=True)) angles = np.linspace(0, 2 * np.pi, 4, endpoint=False).tolist() angles += angles[:1] for ax, (i, row) in zip(axes.flat, top6.iterrows()): values = norm_data.loc[i].tolist() + [norm_data.loc[i].iloc[0]] ax.plot(angles, values, color=color, linewidth=2) ax.fill(angles, values, color=color, alpha=0.25) ''')}

    図と読み取り

    なぜこの図か: 上位 6 機場のプロファイルを同じスケールの 6 つのレーダーで 並べることで、機場ごとの個性が一目で比較できる。

    {figure("assets/L35_fig9_radar_top6.png", "上位 6 機場の 4 軸プロファイル — 容量/台数/効率/守備面積")}

    読み取り:

    """ sections.append(("分析 8: 上位 6 機場のプロファイル比較 (レーダー)", sec11)) # === Section 12: 仮説検証と考察 === sec12_table = hypos_df.to_html(index=False, classes='', border=0) sec12 = f"""

    H1〜H6 の検証結果を 1 表で示す。全 6 仮説が支持された。

    {sec12_table}

    総括: 広島県の内水浸水防御戦略

    2 dataset から再構成した排水機場 {n_total} 件の構造分析により、以下の5 つの設計思想が読み取れる。

    本記事は「排水機場は内水浸水防御の最後の砦」という視点を実データで裏付けた。 {n_total} 件の機場が、2 つの設置目的 (港湾海岸 / 河川) で広島県の {total_m3s:.0f} m³/s の強制排水容量を形成し、 福山平野・広島市の干拓地と河川合流低地に集中投資されている。 この網はL8 (多重浸水想定) で扱った浸水域の内側に意図的に置かれた防御線であり、 浸水想定 → 機場立地 → 排水という多重防御の最終層を成す。

    L32 / L33 / L34 / L35 シリーズ比較表

    軸L32 外郭L33 係留L34 交通L35 排水機場 (本記事)
    件数8421,224296{n_total}
    カテゴリ軸港湾/漁港港湾/漁港港湾/漁港港湾海岸/河川
    主指標延長 m延長 m面積 m²吐出量 m³/s + ポンプ台数
    地理範囲瀬戸内海岸瀬戸内海岸瀬戸内海岸沿岸 + 内陸河川合流部
    主役立地防波堤 (海面)岸壁 (汀線)道路 (陸地)低地・干拓地
    機能波防御船接岸陸海接続内水排水
    関係する浸水津波——河川浸水想定
    集中度上位 3 港 44%上位 3 港 46%上位 3 港 50%上位 2 機 {top2_share:.0f}%
    整備時期明治-現代大正-現代戦後-現代1965-2016 (51年間)

    L32-L34 が「港湾」の階層、L35 は「内水排水」の独立施設。 本シリーズで広島県の「外水→内水」防御の全層を扱い終えた。

    """ sections.append(("仮説検証と考察", sec12)) # === Section 13: 発展課題 === sec13 = f"""

    結果 X1 → 新仮説 Y1 → 課題 Z1: 集水区の正確な範囲特定

    結果 X2 → 新仮説 Y2 → 課題 Z2: 機場運転実績との照合

    結果 X3 → 新仮説 Y3 → 課題 Z3: 内水ハザードマップとの統合

    結果 X4 → 新仮説 Y4 → 課題 Z4: 機場立地の標高・地形分析

    結果 X5 → 新仮説 Y5 → 課題 Z5: 福山平野の干拓地史と機場整備

    結果 X6 → 新仮説 Y6 → 課題 Z6: 「総排水量(m3)」単位の確認

    """ sections.append(("発展課題", sec13)) # レンダリング html = render_lesson( num=35, title=f"排水機場 2 件統合分析 — 広島県 {n_total} 機場 ({n_h} 港湾海岸 + {n_r} 河川) が支える内水浸水防御の最終線", tags=["排水機場", "内水氾濫", "ポンプ", "GIS", "geopandas", "浸水想定", "干拓地", "強制排水"], time="40 分", level="リテラシ", data_label="DoBoX 排水機場 2 dataset (1269 港湾海岸, 1273 河川)", sections=sections, script_filename="L35_drainage_pumps.py", ) (LESSONS / "L35_drainage_pumps.html").write_text(html, encoding="utf-8") print(f" saved L35_drainage_pumps.html ({time.time()-t1:.1f}s)", flush=True) print(f"\n=== 完了: 全所要 {time.time()-t0:.1f}s ===", flush=True)