# -*- coding: utf-8 -*- """L45 ため池浸水想定区域 2 件統合分析 — 広島県内の防災重点ため池決壊リスクを「規模属性 × 浸水範囲」で読み解く カバー宣言: 本記事は DoBoX のシリーズ「ため池」系 2 dataset (id = 62, 63) を 厳密に統合し、広島県内の防災重点ため池決壊浸水想定の構造を分析する。 「ため池」とは農業用水確保のために谷や低地を堰き止めた小規模貯水池で、 日本には約 15 万箇所、その多くが江戸期以前から続く農業遺構である。 広島県は全国有数のため池県 (兵庫・広島・香川が御三家) で、 瀬戸内型の少雨気候下で水田農業を支えてきた基幹インフラ。 しかし 2018 年西日本豪雨でため池決壊が相次ぎ (県内 32 箇所損傷)、 従来「人里離れた静かな農業遺構」だったものが急に「下流に脅威を撒く 決壊リスク源」として再認識された。これを受け 2019 年に 「農業用ため池の管理及び保全に関する法律」が施行され、 決壊時に下流被害が想定される池は「防災重点農業用ため池」 (=「特定農業用ため池」) として県知事が指定するようになった。 本記事はその指定の結果として作られた 2 dataset を扱う: - dataset 62: ため池基本情報 (CSV, 6,754 件) 広島県内の特定農業用ため池の属性表 (堤高/貯水量/緯度経度/所管市町/指定年月日) - dataset 63: ため池浸水想定区域情報 (Shapefile, 7,223 polygons / 6,730 unique 池) ため池ごとに、決壊時の下流浸水想定エリアをポリゴンで記述 この 2 件は「同じため池群の異なる断面」であり、 ため池番号 (CSV) と FIELD001 (SHP) を結合キーに統合分析できる。 本記事は L4-L11 (河川浸水) ・L44 (高潮浸水) と厳密に独立。 ため池決壊浸水と河川浸水・高潮浸水は水文機構が異なる別災害として扱う。 L07 (浸水×インフラ統合) で dataset 62 を「インフラ施設」の 1 として既参照しているが、 本記事はため池そのものを主役に据え、Shapefile 浸水想定区域を主幾何データとする。 研究の問い (主 RQ): 広島県の特定農業用ため池 6,754 件と、その決壊時浸水想定区域 6,730 件は、 ため池の規模属性 (堤高・貯水量) と浸水範囲・市町分布でどう関係し、 「県のため池決壊リスクの地理学」はどう描けるか？ 仮説 H1〜H6: H1 (カバレッジ完全性): 基本情報 (CSV) に登録された 6,754 池のうち、 浸水想定 (SHP) を持たないため池が少数 (< 1%) 存在する。 存在するなら、その非カバレッジ位置には地理的偏りがある (H1 強支持の場合) か、ランダムに散らばる (反証の場合) か。 H2 (規模 ↔ 浸水面積の正相関): 堤高 (m) や貯水量 (千m³) が大きい池ほど、 決壊時の下流浸水想定面積 (km²) も大きい。Pearson 相関 r > 0.6 と予測。 ただし堤高より貯水量の方が強い相関を示す (=決壊時に流出する水量で決まる)。 H3 (中山間集中): ため池数の上位 5 市町は東広島・福山・庄原・三次等の 内陸〜中山間地域で、沿岸都市部 (広島市・呉市・尾道市) よりはるかに多い。 これは「ため池=農業遺構」の歴史的分布と一致する。 H4 (連鎖重複の地理学): 個別ため池の浸水想定面積を単純合計すると、 全体の unary_union 面積より大幅に大きい (重複比率 > 30%)。 重複は「上流ため池の浸水想定下流に下流ため池の浸水想定が連鎖する」 谷地形で発生し、同じ住宅街が複数のため池決壊で繰り返し被災想定される 連鎖リスク構造を示す。 H5 (指定の波): 特定農業用ため池の指定日は 2020-2023 (R2-R5) に集中し、 R3.5.31 (令和 3 年 5 月 31 日) が最大ピークになっている。 これは農業用ため池法の経過措置期限と整合し、 一斉指定が行政手続き的に行われたことを示す。 H6 (規模分布の対数性): ため池の貯水量分布は対数的 (右に長い裾)。 中央値は数千 m³ 級だが上位は 100 万 m³ 級まで存在し、 貯水量で 3 桁以上のレンジを持つ。これは「小池の分散」と 「大池の少数」が共存する典型的な農業遺構分布である。 要件 S 準拠: - 重い空間処理 (143MB SHP の dissolve, unary_union ~87 秒) は `lessons/_l45_build_cache.py` で事前計算し、 `data/extras/L45_pond_inundation/_cache/` に GPKG/CSV キャッシュ。 - 本スクリプトはキャッシュ読込のみで 1 分以内完走。 要件 T 準拠: ため池位置点マップ + 浸水想定ポリゴンマップ + 市町コロプレス + 福山ズームマップ + small multiples (上位 4 市町並置) + 規模散布。 要件 Q 準拠: 図 9 / 表 9 (ため池数 / 浸水面積 / 規模分布 / 重複 / 指定日 / 市町別 など多角度)。 データ仕様: - dataset 62 (CSV): 6,754 行 × 12 列, 1.1MB 列: ため池番号(キー), ため池名称, 特定農業用ため池の指定 (令和年月日), 緯度, 経度, 堤高（m）, 貯水量（千m3）, 所管市町, PDF1-4 (浸水想定 PDF へのリンク) - dataset 63 (Shapefile): 7,223 polygon, 143MB の SHP 列: IDD, BCODE=130, SCODE=130, HCODE=130, ため池ID(=ため池名称), 市町村名(全 None), FIELD001(=ため池番号 9 桁), geometry (Polygon, EPSG:2445) - CRS: EPSG:2445 (旧 JGD2000 平面 III) → 内部処理は EPSG:6671 (現 JGD2011 平面 III) m 単位 - 結合キー: ため池番号 (CSV) ↔ FIELD001 (SHP), 完全結合 6,730 / 部分結合 24 """ from __future__ import annotations import sys, time, subprocess from pathlib import Path sys.path.insert(0, str(Path(__file__).parent)) from _common import ROOT, ASSETS, LESSONS, render_lesson, code, figure import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.patches import Patch from matplotlib.lines import Line2D import geopandas as gpd plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False t_all = time.time() print("=== L45 ため池浸水想定区域 2 件統合分析 ===", flush=True) # ============================================================================= # 0. 定数・データパス # ============================================================================= TARGET_CRS = "EPSG:6671" DATA_DIR = ROOT / "data" / "extras" / "L45_pond_inundation" CACHE = DATA_DIR / "_cache" # 2 dataset メタ DATASETS = { 62: {"jp": "ため池基本情報", "rid": 90, "format": "CSV", "size": "1.1 MB", "rows": 6754, "告示": "2024-update"}, 63: {"jp": "ため池浸水想定区域情報_Shapefile", "rid": 98562, "format": "Shapefile", "size": "81.2 MB (ZIP) / 143 MB (展開後)", "rows": 7223, "告示": "2025-01-22"}, } # CITY_CD → 市町名 (L44 と同じ広島県 都市計画区域行政区域コード) CITY_NAME = { 101: "広島市中区", 102: "広島市東区", 103: "広島市南区", 104: "広島市西区", 105: "広島市安佐南区", 106: "広島市安佐北区", 107: "広島市安芸区", 108: "広島市佐伯区", 202: "呉市", 203: "竹原市", 204: "三原市", 205: "尾道市", 207: "福山市", 208: "府中市", 209: "三次市", 210: "庄原市", 211: "大竹市", 212: "東広島市", 213: "廿日市市", 214: "安芸高田市", 215: "江田島市", 302: "府中町", 304: "海田町", 307: "熊野町", 309: "坂町", 369: "安芸太田町", 462: "世羅町", } # 所管市町 (CSV 表記) → 代表 CITY_CD のマップ (本記事内マッピング) SHOKAN_TO_CITYCD = { "広島市": 105, # 安佐南区を代表 (中山間部のため池が多い) "呉市": 202, "竹原市": 203, "三原市": 204, "尾道市": 205, "福山市": 207, "府中市": 208, "三次市": 209, "庄原市": 210, "大竹市": 211, "東広島市": 212, "廿日市市": 213, "安芸高田市": 214, "江田島市": 215, "安芸郡府中町": 302, "安芸郡海田町": 304, "安芸郡熊野町": 307, "熊野町": 307, "山県郡安芸太田町": 369, "山県郡北広島町": 462, "世羅郡世羅町": 462, "豊田郡大崎上島町": 203, # 竹原市と同区域 "神石郡神石高原町": 208, # 府中市と同区域 } # 中山間 vs 沿岸 の分類 (本記事独自定義) COASTAL_CITIES = {"広島市", "呉市", "竹原市", "三原市", "尾道市", "福山市", "大竹市", "廿日市市", "江田島市", "豊田郡大崎上島町", "安芸郡府中町", "安芸郡海田町", "安芸郡熊野町", "熊野町"} INLAND_CITIES = {"府中市", "三次市", "庄原市", "東広島市", "安芸高田市", "山県郡安芸太田町", "山県郡北広島町", "世羅郡世羅町", "神石郡神石高原町"} # ============================================================================= # 1. キャッシュ確保 # ============================================================================= print("\n[1] キャッシュ確保", flush=True) t1 = time.time() CACHE_FILES = [ CACHE / "diss_per_pond.gpkg", CACHE / "pond_inund_merge.csv", CACHE / "inund_union.gpkg", CACHE / "admin_diss.gpkg", CACHE / "per_city_stats.csv", ] if not all(p.exists() for p in CACHE_FILES): print(" キャッシュ未作成 → _l45_build_cache.py を実行 (~2 分, 初回のみ)", flush=True) subprocess.run([sys.executable, "-X", "utf8", str(LESSONS / "_l45_build_cache.py")], cwd=ROOT, check=True) else: print(" キャッシュ全 5 ファイル OK", flush=True) # 結合済み属性表 merge = pd.read_csv(CACHE / "pond_inund_merge.csv", dtype={"ため池番号": str, "FIELD001": str}) # Boolean 復元 (CSV は文字列で読まれることがある) merge["has_inund"] = merge["area_km2"].notna() # ため池ごと dissolve 済 ポリゴン (浸水想定区域) diss = gpd.read_file(CACHE / "diss_per_pond.gpkg").to_crs(TARGET_CRS) # 県全体 unary_union (面積二重カウント排除用) union_gdf = gpd.read_file(CACHE / "inund_union.gpkg").to_crs(TARGET_CRS) total_union_area = float(union_gdf["area_km2"].iloc[0]) # 都市計画区域行政区域 (L15 由来) admin_diss = gpd.read_file(CACHE / "admin_diss.gpkg").to_crs(TARGET_CRS) # 市町別統計 city_stats = pd.read_csv(CACHE / "per_city_stats.csv", encoding="utf-8-sig") print(f" merge {len(merge)} rows, diss {len(diss)} polys, " f"admin {len(admin_diss)} cities, city_stats {len(city_stats)} 市町, " f"{time.time()-t1:.1f}s", flush=True) # ============================================================================= # 2. 全体サマリ (= 表 1) # ============================================================================= print("\n[2] 全体サマリ", flush=True) t1 = time.time() n_total = len(merge) n_with = int(merge["has_inund"].sum()) n_without = n_total - n_with sum_indiv_km2 = float(merge["area_km2"].sum()) # 個別合計 (重複あり) overlap_km2 = sum_indiv_km2 - total_union_area overlap_ratio = overlap_km2 / sum_indiv_km2 * 100 if sum_indiv_km2 > 0 else 0.0 n_cities = merge["所管市町"].nunique() mean_dam = float(merge["堤高_m"].mean()) median_dam = float(merge["堤高_m"].median()) max_dam = float(merge["堤高_m"].max()) mean_cap = float(merge["貯水量_千m3"].mean()) median_cap = float(merge["貯水量_千m3"].median()) max_cap = float(merge["貯水量_千m3"].max()) T_overall = pd.DataFrame([ ("ため池総数 (CSV)", f"{n_total:,} 池"), ("浸水想定 SHP あり", f"{n_with:,} 池 ({n_with/n_total*100:.1f}%)"), ("浸水想定 SHP なし", f"{n_without} 池 ({n_without/n_total*100:.1f}%)"), ("所管市町数", f"{n_cities} 市町"), ("個別浸水面積の合計 (重複あり)", f"{sum_indiv_km2:.2f} km²"), ("浸水想定区域の真の面積 (union)", f"{total_union_area:.2f} km²"), ("重複面積", f"{overlap_km2:.2f} km² ({overlap_ratio:.1f}%)"), ("平均堤高", f"{mean_dam:.2f} m (中央値 {median_dam:.1f} m, 最大 {max_dam:.1f} m)"), ("平均貯水量", f"{mean_cap:.2f} 千 m³ (中央値 {median_cap:.2f}, 最大 {max_cap:,.0f})"), ], columns=["指標", "値"]) T_overall.to_csv(ASSETS / "L45_overall.csv", index=False, encoding="utf-8-sig") print(T_overall.to_string(index=False)) print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 3. 浸水想定 SHP 無しため池 24 件の特定 # ============================================================================= print("\n[3] 浸水想定なし 24 件の特定", flush=True) t1 = time.time() no_shp = merge[~merge["has_inund"]].copy().reset_index(drop=True) no_shp_show = no_shp[["ため池番号", "ため池名称", "所管市町", "堤高_m", "貯水量_千m3", "特定農業用ため池の指定"]].copy() no_shp_show.columns = ["ため池番号", "ため池名称", "所管市町", "堤高 m", "貯水量 千m³", "特定指定日"] no_shp_show.to_csv(ASSETS / "L45_no_inund.csv", index=False, encoding="utf-8-sig") print(f" 浸水想定なし {len(no_shp)} 件 / 所管市町: {no_shp['所管市町'].value_counts().to_dict()}") print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 4. 規模属性 (堤高・貯水量) ↔ 浸水面積 の関係 # ============================================================================= print("\n[4] 規模属性 ↔ 浸水面積", flush=True) t1 = time.time() m2 = merge[merge["has_inund"]].copy() # 0 値のログ計算回避 m2 = m2[(m2["堤高_m"] > 0) & (m2["貯水量_千m3"] > 0) & (m2["area_km2"] > 0)].copy() m2["log_dam"] = np.log10(m2["堤高_m"]) m2["log_cap"] = np.log10(m2["貯水量_千m3"]) m2["log_area"] = np.log10(m2["area_km2"]) corr_dam_area = m2[["堤高_m", "area_km2"]].corr().iloc[0, 1] corr_cap_area = m2[["貯水量_千m3", "area_km2"]].corr().iloc[0, 1] corr_logdam_logarea = m2[["log_dam", "log_area"]].corr().iloc[0, 1] corr_logcap_logarea = m2[["log_cap", "log_area"]].corr().iloc[0, 1] T_corr = pd.DataFrame([ ("堤高 (生値) ↔ 浸水面積", f"{corr_dam_area:.3f}"), ("貯水量 (生値) ↔ 浸水面積", f"{corr_cap_area:.3f}"), ("log10 堤高 ↔ log10 面積", f"{corr_logdam_logarea:.3f}"), ("log10 貯水量 ↔ log10 面積", f"{corr_logcap_logarea:.3f}"), ], columns=["相関の組", "Pearson r"]) T_corr.to_csv(ASSETS / "L45_correlation.csv", index=False, encoding="utf-8-sig") print(T_corr.to_string(index=False)) print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 5. 特定農業用ため池 指定日の分布 # ============================================================================= print("\n[5] 指定日分布", flush=True) t1 = time.time() # 令和年月日を西暦に変換 (R2 = 2020, R3 = 2021, ...) def reiwa_to_date(s): if pd.isna(s) or not isinstance(s, str): return None s = s.strip() # 例: R3.5.31, R2.12.28 if not s.startswith("R"): return None try: body = s[1:] parts = body.split(".") if len(parts) != 3: return None r_year = int(parts[0]) month = int(parts[1]) day = int(parts[2]) # 令和元年 = 2019 ad = 2018 + r_year return f"{ad:04d}-{month:02d}-{day:02d}" except Exception: return None merge["指定日_西暦"] = merge["特定農業用ため池の指定"].apply(reiwa_to_date) desig = merge["指定日_西暦"].dropna() desig_count = desig.value_counts().sort_index() T_desig = pd.DataFrame({"指定日": desig_count.index, "ため池数": desig_count.values}) T_desig.to_csv(ASSETS / "L45_designation_dates.csv", index=False, encoding="utf-8-sig") print(f" ユニーク指定日数: {len(desig_count)}") print(T_desig.to_string(index=False)) # 月単位サマリ desig_month = pd.to_datetime(desig).dt.to_period("M") month_counts = desig_month.value_counts().sort_index() print(f"\n 月別累計 (上位 5): {month_counts.head(5).to_dict()}") print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 6. 市町別 集計表 (= 表 6) # ============================================================================= print("\n[6] 市町別集計", flush=True) t1 = time.time() # 中山間/沿岸 列を追加 def categorize(name): if name in INLAND_CITIES: return "中山間" if name in COASTAL_CITIES: return "沿岸" return "その他" city_stats["地理分類"] = city_stats["所管市町"].map(categorize) # 浸水想定密度 = 合計浸水想定 / ため池数 city_stats["池あたり浸水_km2"] = (city_stats["合計浸水想定_km2"] / city_stats["ため池数"]).round(4) # rounded city_stats_view = city_stats.copy() for c in ["平均堤高_m", "最大堤高_m", "平均貯水量_千m3"]: city_stats_view[c] = city_stats_view[c].round(2) for c in ["合計貯水量_千m3"]: city_stats_view[c] = city_stats_view[c].round(0).astype(int) for c in ["合計浸水想定_km2", "平均浸水想定_km2"]: city_stats_view[c] = city_stats_view[c].round(3) city_stats_view.to_csv(ASSETS / "L45_per_city.csv", index=False, encoding="utf-8-sig") print(city_stats_view.head(10).to_string(index=False)) print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 7. Top 浸水想定面積ため池 (= 表 7) # ============================================================================= print("\n[7] Top 浸水想定 ため池", flush=True) t1 = time.time() top_inund = merge.sort_values("area_km2", ascending=False).head(15).copy() top_inund_view = top_inund[["ため池名称", "所管市町", "堤高_m", "貯水量_千m3", "area_km2", "特定農業用ため池の指定"]].copy() top_inund_view.columns = ["ため池名称", "所管市町", "堤高 m", "貯水量 千m³", "浸水想定面積 km²", "特定指定日"] top_inund_view["堤高 m"] = top_inund_view["堤高 m"].round(1) top_inund_view["貯水量 千m³"] = top_inund_view["貯水量 千m³"].round(1) top_inund_view["浸水想定面積 km²"] = top_inund_view["浸水想定面積 km²"].round(3) top_inund_view = top_inund_view.reset_index(drop=True) top_inund_view.index = top_inund_view.index + 1 top_inund_view = top_inund_view.reset_index().rename(columns={"index": "順位"}) top_inund_view.to_csv(ASSETS / "L45_top_inund.csv", index=False, encoding="utf-8-sig") print(top_inund_view.to_string(index=False)) print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 8. 図の生成 (9 枚) # ============================================================================= print("\n[8] 図の生成", flush=True) t1 = time.time() def save_fig(name, dpi=120): p = ASSETS / name plt.savefig(p, dpi=dpi, bbox_inches="tight", facecolor="white") plt.close('all') return p # ----- 図 1: ため池所管市町別 ため池数 (横棒) ----- fig, ax = plt.subplots(figsize=(9, 7.5)) cs = city_stats.sort_values("ため池数", ascending=True) colors = cs["地理分類"].map({"中山間": "#388e3c", "沿岸": "#1976d2", "その他": "#9e9e9e"}) ax.barh(cs["所管市町"], cs["ため池数"], color=colors, edgecolor="#333", linewidth=0.4) for i, (n, v) in enumerate(zip(cs["所管市町"], cs["ため池数"])): ax.text(v + 15, i, f"{v:,}", va="center", fontsize=9) ax.set_xlabel("ため池数 (件)", fontsize=11) ax.set_title("図1: 所管市町別 特定農業用ため池数 (n=6,754)\n緑=中山間 / 青=沿岸", fontsize=12) handles = [Patch(facecolor="#388e3c", label="中山間 (内陸)"), Patch(facecolor="#1976d2", label="沿岸"), Patch(facecolor="#9e9e9e", label="その他")] ax.legend(handles=handles, loc="lower right", fontsize=9) ax.grid(axis="x", alpha=0.3) plt.tight_layout() save_fig("L45_fig1_pond_count_by_city.png") # ----- 図 2: 規模分布 (堤高・貯水量) ヒストグラム log scale ----- fig, axes = plt.subplots(1, 2, figsize=(13, 5)) ax = axes[0] ax.hist(merge["堤高_m"].dropna(), bins=40, color="#5e35b1", edgecolor="#333", linewidth=0.3) ax.axvline(median_dam, color="#cf222e", linestyle="--", label=f"中央値 {median_dam:.1f} m") ax.axvline(mean_dam, color="#0969da", linestyle=":", label=f"平均 {mean_dam:.2f} m") ax.set_xlabel("堤高 (m)", fontsize=11) ax.set_ylabel("ため池数", fontsize=11) ax.set_title(f"堤高分布 (n={n_total:,})", fontsize=11) ax.legend(fontsize=9) ax.grid(alpha=0.3) ax = axes[1] cap_pos = merge["貯水量_千m3"][merge["貯水量_千m3"] > 0].dropna() ax.hist(cap_pos, bins=np.logspace(np.log10(cap_pos.min()), np.log10(cap_pos.max()), 40), color="#00838f", edgecolor="#333", linewidth=0.3) ax.set_xscale("log") ax.axvline(median_cap, color="#cf222e", linestyle="--", label=f"中央値 {median_cap:.2f} 千m³") ax.axvline(mean_cap, color="#0969da", linestyle=":", label=f"平均 {mean_cap:.2f}") ax.set_xlabel("貯水量 (千 m³, log10 軸)", fontsize=11) ax.set_ylabel("ため池数", fontsize=11) ax.set_title("貯水量分布 (log scale)", fontsize=11) ax.legend(fontsize=9) ax.grid(alpha=0.3, which="both") fig.suptitle("図2: ため池規模分布 — 堤高は正規様, 貯水量は対数様 (右に長い裾)", fontsize=12) plt.tight_layout() save_fig("L45_fig2_size_dist.png") # ----- 図 3: 規模 vs 浸水面積 散布図 (log-log) ----- fig, axes = plt.subplots(1, 2, figsize=(13, 5.5)) ax = axes[0] ax.scatter(m2["堤高_m"], m2["area_km2"], s=8, alpha=0.4, c="#5e35b1", edgecolor="none") ax.set_xscale("log"); ax.set_yscale("log") ax.set_xlabel("堤高 (m, log)", fontsize=11) ax.set_ylabel("浸水想定面積 (km², log)", fontsize=11) ax.set_title(f"堤高 vs 浸水面積\nPearson r (生値) = {corr_dam_area:.3f}, " f"log-log r = {corr_logdam_logarea:.3f}", fontsize=11) ax.grid(alpha=0.3, which="both") ax = axes[1] ax.scatter(m2["貯水量_千m3"], m2["area_km2"], s=8, alpha=0.4, c="#00838f", edgecolor="none") # 回帰線 (log-log で OLS) from numpy.polynomial.polynomial import polyfit x_l = m2["log_cap"].values y_l = m2["log_area"].values b0, b1 = polyfit(x_l, y_l, 1) xr = np.linspace(x_l.min(), x_l.max(), 50) ax.plot(10**xr, 10**(b0 + b1*xr), "-", color="#cf222e", linewidth=2, label=f"OLS log: y = {b0:.2f} + {b1:.2f}x") ax.set_xscale("log"); ax.set_yscale("log") ax.set_xlabel("貯水量 (千m³, log)", fontsize=11) ax.set_ylabel("浸水想定面積 (km², log)", fontsize=11) ax.set_title(f"貯水量 vs 浸水面積\nPearson r (生値) = {corr_cap_area:.3f}, " f"log-log r = {corr_logcap_logarea:.3f}", fontsize=11) ax.legend(fontsize=9, loc="lower right") ax.grid(alpha=0.3, which="both") fig.suptitle("図3: ため池規模 ↔ 決壊時下流浸水想定面積 の log-log スケーリング", fontsize=12) plt.tight_layout() save_fig("L45_fig3_size_vs_area.png") # ----- 図 4: 県全域 ため池位置点マップ (規模で色分け) ----- fig, ax = plt.subplots(figsize=(13, 7)) admin_diss.plot(ax=ax, color="#f8f8f8", edgecolor="#aaa", linewidth=0.4) # 緯度経度 → EPSG:6671 に投影変換 (簡易: 緯度経度 4326 を直接変換) pond_pts = gpd.GeoDataFrame( merge, geometry=gpd.points_from_xy(merge["経度"], merge["緯度"]), crs="EPSG:4326" ).to_crs(TARGET_CRS) # size = log10 貯水量 にマッピング sizes = np.clip(np.log10(merge["貯水量_千m3"].clip(lower=0.1)) + 2, 1, 50) * 5 # 色 = 浸水想定の有無 colors = pond_pts["has_inund"].map({True: "#cf222e", False: "#0969da"}).fillna("#aaa") ax.scatter(pond_pts.geometry.x, pond_pts.geometry.y, s=sizes, c=colors, alpha=0.55, edgecolor="none") ax.set_title(f"図4: 広島県 特定農業用ため池 位置マップ (n={n_total:,})\n" f"円サイズ ∝ log10 貯水量, 色: 赤=浸水想定あり ({n_with:,}), " f"青=なし ({n_without})", fontsize=11) ax.set_xlabel("EPSG:6671 X (m)") ax.set_ylabel("EPSG:6671 Y (m)") ax.set_aspect("equal") handles = [Line2D([0], [0], marker='o', color='w', markerfacecolor='#cf222e', markersize=10, label=f'浸水想定あり ({n_with:,} 池)'), Line2D([0], [0], marker='o', color='w', markerfacecolor='#0969da', markersize=10, label=f'浸水想定なし ({n_without} 池)')] ax.legend(handles=handles, loc="lower right", fontsize=9) plt.tight_layout() save_fig("L45_fig4_pond_points_map.png") # ----- 図 5: 浸水想定区域 ポリゴンマップ (県全域) ----- fig, ax = plt.subplots(figsize=(13, 7)) admin_diss.plot(ax=ax, color="#f8f8f8", edgecolor="#aaa", linewidth=0.4) # 浸水想定 全 6730 ポリゴンを薄赤で diss.plot(ax=ax, color="#e64a19", edgecolor="none", alpha=0.55) ax.set_title(f"図5: ため池決壊時 浸水想定区域マップ (6,730 池)\n" f"個別合計 = {sum_indiv_km2:.1f} km², 真の面積 (union) = " f"{total_union_area:.1f} km², 重複 = {overlap_km2:.1f} km² " f"({overlap_ratio:.0f}%)", fontsize=11) ax.set_xlabel("EPSG:6671 X (m)") ax.set_ylabel("EPSG:6671 Y (m)") ax.set_aspect("equal") plt.tight_layout() save_fig("L45_fig5_inund_polygon_map.png") # ----- 図 6: 上位 4 市町 small multiples (ため池 + 浸水) ----- fig, axes = plt.subplots(2, 2, figsize=(14, 10)) top4_cities = city_stats.head(4)["所管市町"].tolist() # 各市町の bbox を計算 (CSV 緯度経度ベース) for ax, city in zip(axes.flat, top4_cities): # 該当する CITY_CD citycd = SHOKAN_TO_CITYCD.get(city) sub_pond_csv = merge[merge["所管市町"] == city] sub_pond_pts = pond_pts[pond_pts["所管市町"] == city] sub_diss = diss[diss["FIELD001"].isin(sub_pond_csv["ため池番号"].dropna())] if len(sub_diss) == 0: continue # bbox (浸水ポリゴン基準, 1km マージン) xmin, ymin, xmax, ymax = sub_diss.total_bounds margin = 1500 xmin -= margin; ymin -= margin; xmax += margin; ymax += margin # 背景行政 admin_diss.plot(ax=ax, color="#f8f8f8", edgecolor="#aaa", linewidth=0.4) # 浸水想定 sub_diss.plot(ax=ax, color="#e64a19", edgecolor="none", alpha=0.55) # ため池点 ax.scatter(sub_pond_pts.geometry.x, sub_pond_pts.geometry.y, s=8, c="#1976d2", alpha=0.7, edgecolor="none") ax.set_xlim(xmin, xmax); ax.set_ylim(ymin, ymax) ax.set_aspect("equal") ax.set_xticks([]); ax.set_yticks([]) ax.set_title(f"{city}: ため池 {len(sub_pond_csv)} 件, " f"浸水想定 {sub_diss['area_km2'].sum() if len(sub_diss) else 0:.1f} km²", fontsize=10) fig.suptitle("図6: ため池数上位 4 市町 ズームマップ (赤=浸水想定 / 青点=ため池位置)", fontsize=13) plt.tight_layout() save_fig("L45_fig6_top4_cities_zoom.png") # ----- 図 7: 福山市の浸水想定 詳細マップ (服部大池ほか) ----- fig, ax = plt.subplots(figsize=(13, 8)) fukuyama_csv = merge[merge["所管市町"] == "福山市"] fukuyama_pond_pts = pond_pts[pond_pts["所管市町"] == "福山市"] fukuyama_diss = diss[diss["FIELD001"].isin(fukuyama_csv["ため池番号"].dropna())] xmin, ymin, xmax, ymax = fukuyama_diss.total_bounds margin = 2000 xmin -= margin; ymin -= margin; xmax += margin; ymax += margin admin_diss.plot(ax=ax, color="#f8f8f8", edgecolor="#aaa", linewidth=0.4) # 浸水想定 (色は面積 log で濃淡) fk_with_area = fukuyama_diss.copy() fk_with_area["area_km2"] = fk_with_area.geometry.area / 1e6 fk_with_area.plot(ax=ax, column="area_km2", cmap="OrRd", edgecolor="none", alpha=0.65, legend=True, legend_kwds={"label": "浸水想定面積 (km²)", "shrink": 0.6}) # ため池点 ax.scatter(fukuyama_pond_pts.geometry.x, fukuyama_pond_pts.geometry.y, s=10, c="#0d47a1", alpha=0.8, edgecolor="none") # Top3 池の名前 top3_fk = fukuyama_csv.sort_values("area_km2", ascending=False).head(3) top3_fk_pts = pond_pts[pond_pts["ため池番号"].isin(top3_fk["ため池番号"])] for _, p_csv in top3_fk.iterrows(): p_pt = top3_fk_pts[top3_fk_pts["ため池番号"] == p_csv["ため池番号"]] if len(p_pt) == 0: continue cen = p_pt.geometry.iloc[0] ax.annotate(f"{p_csv['ため池名称']}\n{p_csv['area_km2']:.2f} km²", (cen.x, cen.y), fontsize=9, ha="center", weight="bold", xytext=(10, 10), textcoords="offset points", bbox=dict(boxstyle="round,pad=0.3", fc="white", alpha=0.9, ec="#333")) ax.set_xlim(xmin, xmax); ax.set_ylim(ymin, ymax) ax.set_aspect("equal") ax.set_title(f"図7: 福山市 ため池決壊浸水想定 詳細マップ " f"({len(fukuyama_csv)} 池, 合計 {fukuyama_csv['area_km2'].sum():.1f} km²)\n" f"最大池: 服部大池 ({fukuyama_csv['area_km2'].max():.2f} km², " f"福山市内で最大)", fontsize=11) ax.set_xlabel("EPSG:6671 X (m)") ax.set_ylabel("EPSG:6671 Y (m)") plt.tight_layout() save_fig("L45_fig7_fukuyama_detail.png") # ----- 図 8: 指定日時系列ヒストグラム ----- fig, ax = plt.subplots(figsize=(11, 5.2)) desig_dt = pd.to_datetime(desig) # 月単位に集計 months = desig_dt.dt.to_period("M").value_counts().sort_index() xs = [p.to_timestamp() for p in months.index] ax.bar(xs, months.values, width=20, color="#388e3c", edgecolor="#333", linewidth=0.3) # Top 3 月 にラベル top_months = months.sort_values(ascending=False).head(3) for p, n in top_months.items(): ax.annotate(f"{p.strftime('%Y-%m')}\n{n} 池", (p.to_timestamp(), n), fontsize=9, ha="center", va="bottom") ax.set_xlabel("特定農業用ため池 指定日 (月単位)", fontsize=11) ax.set_ylabel("指定ため池数", fontsize=11) ax.set_title(f"図8: 特定農業用ため池 指定の波 (n={int(desig_count.sum())} / " f"未指定 {int((merge['指定日_西暦'].isna()).sum())} を除く)\n" f"R3.5.31 (令和 3 年 5 月末日 = 経過措置期限) が突出", fontsize=12) ax.grid(axis="y", alpha=0.3) plt.tight_layout() save_fig("L45_fig8_designation_timeline.png") # ----- 図 9: 市町別 ため池密度 vs 浸水想定密度 (散布) ----- fig, ax = plt.subplots(figsize=(10, 7)) plot_df = city_stats.copy() plot_df["地理色"] = plot_df["地理分類"].map( {"中山間": "#388e3c", "沿岸": "#1976d2", "その他": "#9e9e9e"}) ax.scatter(plot_df["ため池数"], plot_df["合計浸水想定_km2"], s=120, c=plot_df["地理色"], edgecolor="#333", alpha=0.85) # 全 市町ラベル for _, r in plot_df.iterrows(): ax.annotate(r["所管市町"], (r["ため池数"], r["合計浸水想定_km2"]), fontsize=9, ha="left", xytext=(5, 3), textcoords="offset points") # 1:1 参照線 (= 1 池あたり 1 km² の場合) maxn = plot_df["ため池数"].max() ax.plot([0, maxn], [0, maxn * 0.05], ":", color="#aaa", label="1 池 = 0.05 km² 参考線") ax.set_xlabel("ため池数 (件)", fontsize=11) ax.set_ylabel("合計浸水想定面積 (km²)", fontsize=11) ax.set_title("図9: 市町別 ため池数 vs 浸水想定面積 (中山間/沿岸 色分け)\n" "東広島・福山が突出し、両者で県全体の約 35% を占める", fontsize=11) handles = [Patch(facecolor="#388e3c", label="中山間"), Patch(facecolor="#1976d2", label="沿岸"), Patch(facecolor="#9e9e9e", label="その他"), Line2D([0], [0], color="#aaa", linestyle=":", label="1 池 = 0.05 km² 参考線")] ax.legend(handles=handles, loc="lower right", fontsize=9) ax.grid(alpha=0.3) plt.tight_layout() save_fig("L45_fig9_density_scatter.png") print(f" 9 図 saved ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 9. 仮説検証 # ============================================================================= print("\n[9] 仮説検証", flush=True) t1 = time.time() H_results = [] # H1 カバレッジ完全性 nc_ratio = n_without / n_total * 100 nc_cities = no_shp["所管市町"].value_counts() # 全件が「中山間市町セット (INLAND_CITIES)」に集中するか nc_inland = sum(1 for c in nc_cities.index if c in INLAND_CITIES) all_inland = (nc_inland == len(nc_cities)) H1_supp = ("強支持 (中山間市町に集中, 地理偏りあり)" if (n_without > 0 and all_inland and n_without/n_total < 0.01) else ("支持 (少数だが多市町に散在)" if (n_without > 0 and n_without/n_total < 0.01) else "反証")) nc_breakdown = ", ".join(f"{k}={v}" for k, v in nc_cities.items()) H_results.append({ "仮説": "H1 カバレッジ完全性", "想定": "浸水想定なしため池が < 1% で、地理的偏り (中山間市町集中) がある", "実測": f"{n_without}/{n_total} ({nc_ratio:.2f}%) が浸水想定なし, " f"内訳: {nc_breakdown} (全件が中山間 {nc_inland}/{len(nc_cities)} 市町)", "判定": H1_supp, }) # H2 規模 ↔ 浸水面積 正相関 H2_strong = (corr_cap_area > corr_dam_area) and (corr_cap_area > 0.6) H2_supp = ("強支持" if (H2_strong and corr_cap_area > 0.7) else ("支持" if H2_strong else "部分支持")) H_results.append({ "仮説": "H2 規模 ↔ 浸水面積 正相関", "想定": "Pearson r > 0.6, 貯水量 > 堤高", "実測": f"r(堤高)={corr_dam_area:.3f}, r(貯水量)={corr_cap_area:.3f}, " f"log-log: r(貯水量)={corr_logcap_logarea:.3f}", "判定": H2_supp, }) # H3 中山間集中 top5 = city_stats.head(5)["所管市町"].tolist() inland_in_top5 = sum(1 for c in top5 if c in INLAND_CITIES) H3_supp = ("強支持" if inland_in_top5 >= 3 else ("支持" if inland_in_top5 >= 2 else "反証")) H_results.append({ "仮説": "H3 中山間集中", "想定": "ため池数 上位 5 市町のうち 3 件以上が中山間 (内陸)", "実測": f"上位 5 = {top5}, うち中山間 = {inland_in_top5} 件", "判定": H3_supp, }) # H4 連鎖重複 H4_supp = ("強支持" if overlap_ratio > 30 else ("支持" if overlap_ratio > 20 else "反証")) H_results.append({ "仮説": "H4 連鎖重複の地理学", "想定": "個別合計 - union = 30% 以上の重複面積", "実測": f"個別合計={sum_indiv_km2:.1f} km², union={total_union_area:.1f} km², " f"overlap={overlap_km2:.1f} km² ({overlap_ratio:.1f}%)", "判定": H4_supp, }) # H5 指定日のピーク R3.5.31 top_date = desig_count.idxmax() top_date_count = int(desig_count.max()) total_designated = int(desig_count.sum()) top_share = top_date_count / total_designated * 100 H5_supp = ("強支持" if (top_date == "2021-05-31" and top_share > 30) else ("支持" if top_share > 25 else "部分支持")) H_results.append({ "仮説": "H5 指定の波 R3.5.31 ピーク", "想定": "R3.5.31 (2021-05-31) が単一最大ピークで > 30%", "実測": f"top_date={top_date} ({top_date_count} 池, {top_share:.1f}% / 総指定数 {total_designated})", "判定": H5_supp, }) # H6 規模分布の対数性 log_range = np.log10(max_cap) - np.log10(merge["貯水量_千m3"][merge["貯水量_千m3"]>0].min()) H6_supp = ("強支持" if log_range >= 5 else ("支持" if log_range >= 3 else "反証")) H_results.append({ "仮説": "H6 規模分布の対数性", "想定": "貯水量レンジが log10 で 3 桁以上", "実測": f"min={merge['貯水量_千m3'][merge['貯水量_千m3']>0].min():.4f}, " f"max={max_cap:,.1f} 千 m³, log10 range={log_range:.2f}", "判定": H6_supp, }) df_hypo = pd.DataFrame(H_results) df_hypo.to_csv(ASSETS / "L45_hypothesis.csv", index=False, encoding="utf-8-sig") print(df_hypo.to_string(index=False)) print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 10. HTML レンダリング # ============================================================================= print("\n[10] HTML 生成", flush=True) t1 = time.time() def df_to_html(df, max_rows=None): if max_rows and len(df) > max_rows: df = df.head(max_rows) return df.to_html(index=False, escape=False, classes="data") # データ仕様表 T_dataspec = pd.DataFrame([ {"dataset_id": 62, "シリーズ": "ため池基本情報", "形式": "CSV (Shift_JIS, BOM なし)", "ZIP/サイズ": "1.1 MB", "件数": "6,754 行 × 12 列", "更新": "2024 年版 (随時更新)", "本記事での役割": "ため池属性 (堤高/貯水量/緯度経度/所管市町/特定指定日)"}, {"dataset_id": 63, "シリーズ": "ため池浸水想定区域情報_Shapefile", "形式": "Shapefile (.shp/.shx/.dbf/.prj/.sbn/.sbx)", "ZIP/サイズ": "81.2 MB (ZIP) / 143 MB (展開後 .shp 単体)", "件数": "7,223 polygon / 6,730 unique ため池", "更新": "2025-01-22", "本記事での役割": "ため池ごとの決壊時下流浸水想定エリア (主幾何データ)"}, ]) T_dataspec.to_csv(ASSETS / "L45_dataspec.csv", index=False, encoding="utf-8-sig") # 結合キー解説表 T_joinkey = pd.DataFrame([ {"段階": "STEP 0", "内容": "tameike_basic.csv を読込", "サイズ/件数": "6,754 行 × 12 列"}, {"段階": "STEP 1", "内容": "ため池浸水想定区域.shp を読込 → EPSG:6671 投影変換", "サイズ/件数": "7,223 polygon (元の CRS = EPSG:2445)"}, {"段階": "STEP 2", "内容": "30 件の不正幾何 (自己交差等) を buffer(0) で修復", "サイズ/件数": "30 polygon を修復 → 7,223 polygon (有効)"}, {"段階": "STEP 3", "内容": "FIELD001 (= ため池番号) で dissolve", "サイズ/件数": "7,223 → 6,730 ため池ポリゴン (1 池 = 1 MultiPolygon)"}, {"段階": "STEP 4", "内容": "merge: CSV.ため池番号 LEFT JOIN SHP.FIELD001", "サイズ/件数": "6,754 行 (うち 6,730 結合成功 + 24 結合無し)"}, {"段階": "STEP 5", "内容": "全 7,223 ポリゴン unary_union で「県全体 真の浸水想定面積」", "サイズ/件数": f"単一 MultiPolygon, 面積 = {total_union_area:.2f} km²"}, ]) T_joinkey.to_csv(ASSETS / "L45_joinkey_steps.csv", index=False, encoding="utf-8-sig") sections = [] # (1) 学習目標と問い sections.append(("学習目標と問い", f"""

本レッスンでは、広島県オープンデータ DoBoX が公開する ため池系 2 dataset (id = 62, 63) を統合し、 広島県内の「特定農業用ため池の決壊リスクと浸水想定の構造」を分析する。

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

主 RQ (研究の問い)

広島県の特定農業用ため池 6,754 件と、その決壊時浸水想定区域 6,730 件は、 ため池の規模属性 (堤高・貯水量) と浸水範囲・市町分布でどう関係し、 「県のため池決壊リスクの地理学」はどう描けるか？

仮説 H1〜H6 (本記事で検証する)

到達点

2 dataset を ため池番号 をキーに完全結合し、規模属性と浸水想定面積を クロスして 6 仮説を量的に検証する。学習者は (a) Shapefile + CSV の異種統合、 (b) dissolve(by=...) による属性集約と unary_union による 重複排除の使い分け、(c) ため池決壊という農業×防災の交差問題のデータ表現 を体得する。

厳密性の宣言: 本記事は L4-L11 (河川浸水) ・L44 (高潮浸水) と 厳密に独立。河川浸水 (流域降雨) ・高潮浸水 (海面異常上昇) ・ ため池決壊浸水 (人工堤体破壊) は水文機構が別の災害として扱う。 L07 (浸水×インフラ統合) で dataset 62 を「インフラ施設」の 1 として既参照しているが、 本記事はため池そのものを主役に据え、Shapefile 浸水想定区域を主幾何データに用いる 新しい視点である。
""")) # (2) 使用データ sections.append(("使用データ", f"""

本記事は DoBoX の「ため池」関連 2 dataset を統合する:

表 (1) — 2 dataset の仕様サマリ

{df_to_html(T_dataspec)}

この表から読み取れること: (1) 形式が完全に異なる (CSV vs Shapefile) ため、統合には結合キー設計が必須。 (2) サイズは Shapefile が 100 倍以上 (1 MB → 143 MB) で、 ポリゴン幾何が支配的。CSV は属性のみだが、地理空間処理を支える「場所情報」と 「規模指標」の両方を持つ。 (3) dataset 62 のキー = ため池番号 (9 桁数字) は dataset 63 の FIELD001 列と完全一致するため、属性結合 (merge) でデータ統合できる。 (4) dataset 63 の Shapefile は 1 池に複数 polygon (浸水範囲が分断された場合) を 持つことがあり、dissolve(by="FIELD001") で 7,223 → 6,730 に集約する。

表 (2) — 結合 STEP の段階表 (要件 K Before/After)

{df_to_html(T_joinkey)}

この表から読み取れること: 入力 6,754 行 (CSV) が STEP 5 まで 通って union 単一 polygon (面積 {total_union_area:.2f} km²) に至るまでの 段階的サイズ変化を示す。STEP 3 で 7,223 → 6,730 に減るのは 「1 池の浸水範囲が分断された複数 polygon」が dissolve で結合されるため、 STEP 4 で 6,754 - 6,730 = 24 件が結合無しで残るのは 「CSV にあって SHP に対応 polygon 無いため池」=分析 1 で詳細特定する。

形式の詳細

""")) # (3) ダウンロード sections.append(("ダウンロード（再現用データ・中間データ・図）", f"""

HTML から直接以下を取得できる:

(A) DoBoX 直リンク (2 dataset)

datasetカタログresource_download (直 DL)形式サイズ
62 ため池基本情報 dataset 62 直 DL (rid 90) CSV1.1 MB
63 ため池浸水想定区域情報_Shapefile dataset 63 直 DL (rid 98562) Shapefile (ZIP)81.2 MB

(B) PowerShell 一括取得 (再現用)

cd "2026 DoBoX 教材"
mkdir -Force data\\extras\\L45_pond_inundation
iwr "https://hiroshima-dobox.jp/resource_download/90"    -OutFile "data\\extras\\tameike_basic.csv"
iwr "https://hiroshima-dobox.jp/resource_download/98562" -OutFile "data\\extras\\L45_pond_inundation\\tameike_inundation_shp.zip"
# ZIP 展開
Expand-Archive -Force "data\\extras\\L45_pond_inundation\\tameike_inundation_shp.zip" "data\\extras\\L45_pond_inundation\\shp"
# キャッシュビルド (~2 分, 初回のみ)
py -X utf8 lessons\\_l45_build_cache.py
# 本記事スクリプト (~1 分)
py -X utf8 lessons\\L45_pond_inundation.py

(C) 中間データ・図 (本記事生成物の直リンク)

""")) # (4) 分析 1: カバレッジ完全性 — 浸水想定なし 24 件の特定 sections.append(("分析 1: カバレッジ完全性 — CSV と SHP の重ね合わせ", f"""

狙い (RQ ↔ 仮説 H1)

2 dataset のカバレッジ差を量的に評価する。 CSV (62) は 6,754 池、SHP (63) は 7,223 polygon (= 6,730 unique 池)。 カバレッジ差 24 件は「CSV にあって SHP 無し」のため池であり、 地理的偏りがあれば「データ整備の優先度」に関する制度的事実が データから直接読める。

手法 — pandas merge LEFT JOIN による属性結合

pandas の DataFrame.merge(how='left') で CSV を主表、SHP の polygon area を従属表として結合する。 結合キー は CSV.ため池番号 = SHP.FIELD001。 LEFT JOIN を選んだ理由は「CSV 側を主張として残し、 SHP に対応 polygon が無い池は area_km2 が NaN になる」ことで 非カバレッジ件数がそのまま得られるため。

結合キーの設計 (要件 O STEP 分け):

実装コード

{code(''' import pandas as pd, geopandas as gpd from pathlib import Path CACHE = Path("data/extras/L45_pond_inundation/_cache") # (1) ため池基本情報 (CSV) を読込, 数値列は to_numeric df = pd.read_csv("data/extras/tameike_basic.csv", dtype=str) df.columns = [c.strip() for c in df.columns] df["ため池番号"] = df["ため池番号"].str.strip() df["堤高_m"] = pd.to_numeric(df["堤高（m）"], errors="coerce") df["貯水量_千m3"] = pd.to_numeric(df["貯水量（千m3）"], errors="coerce") # (2) 浸水想定 SHP は事前 dissolve 済 GPKG (= ため池ごと 1 polygon) を読込 diss = gpd.read_file(CACHE / "diss_per_pond.gpkg") # diss["FIELD001"] が結合キー, diss["area_km2"] が浸水想定面積 # (3) LEFT JOIN: 結合無しは area_km2 = NaN として残る merge = df.merge( diss[["FIELD001", "area_km2"]], left_on="ため池番号", right_on="FIELD001", how="left", ) merge["has_inund"] = merge["area_km2"].notna() # (4) 浸水想定なし のため池を抽出 no_shp = merge[~merge["has_inund"]] print(f"非カバレッジ {{len(no_shp)}} 件 / 所管市町: {{no_shp['所管市町'].value_counts().to_dict()}}") ''')}

結果の図

図4 を選んだ理由: 6,754 池の地理分布を一目で把握するには 点マップが最適。色 = 浸水想定の有無、サイズ = log10 貯水量 の 2 軸を 1 図で示すことで、「どこに大きい池があるか」「どこの池が 非カバレッジか」を同時に確認できる。棒グラフでは地理が見えない。

{figure('assets/L45_fig4_pond_points_map.png', '図4: 広島県 特定農業用ため池 位置マップ')}

図4 から読み取れること:

結果の表 (1) — 全体サマリ

{df_to_html(T_overall)}

表 (1) から読み取れること:

結果の表 (2) — 浸水想定なし 24 件の内訳

{df_to_html(no_shp_show.head(15))}

(表は上位 15 件のみ表示, 全 24 件は L45_no_inund.csv 参照)

表 (2) から読み取れること:

""")) # (5) 分析 2: 規模属性 ↔ 浸水面積 のスケーリング sections.append(("分析 2: 規模属性 ↔ 浸水面積 のスケーリング", f"""

狙い (RQ ↔ 仮説 H2)

「ため池の規模 (堤高・貯水量) が大きいほど決壊時の下流浸水面積も大きい」 という直感的予測を、Pearson 相関と log-log 散布で検証する。 堤高と貯水量のどちらが浸水面積をより強く支配するかも比較する。

手法 — log-log 散布 + OLS 回帰

ため池の規模指標 (堤高 m, 貯水量 千 m³) は桁数のレンジが広い (堤高 0.2-33 m = 約 2 桁、貯水量 0.002-1,053 千 m³ = 約 6 桁)。 よって生値の散布図では大池の点だけが目立ち、小池の構造が潰れる。 log-log 軸 ((x, y) の両方を log10 化) で表示すると べき乗則 (power law) の傾きが直線として読める。

Pearson 相関 r とは (要件 B/J): 2 変数 (x, y) の線形関係の強さを -1〜+1 で示す。 +1 = 完全な正の直線、0 = 無相関、-1 = 完全な負の直線。 入力: 同じ長さの 2 列の数値、出力: 1 つの相関値。 限界: 非線形関係 (放物線等) は捉えられない。代替: Spearman 順位相関、 log-log 化での Pearson (本記事はこれを併用)。

実装コード

{code(''' import numpy as np import pandas as pd # 結合済 merge 表から「浸水想定あり」だけ抽出 (NaN 除外) m2 = merge[merge["has_inund"]].copy() m2 = m2[(m2["堤高_m"] > 0) & (m2["貯水量_千m3"] > 0) & (m2["area_km2"] > 0)].copy() # 生値 Pearson 相関 corr_dam_area = m2[["堤高_m", "area_km2"]].corr().iloc[0, 1] corr_cap_area = m2[["貯水量_千m3", "area_km2"]].corr().iloc[0, 1] # log10 で再相関 (べき乗則の検証) m2["log_dam"] = np.log10(m2["堤高_m"]) m2["log_cap"] = np.log10(m2["貯水量_千m3"]) m2["log_area"] = np.log10(m2["area_km2"]) log_corr_cap_area = m2[["log_cap", "log_area"]].corr().iloc[0, 1] # log-log OLS 回帰: log(area) = b0 + b1 * log(cap) b0, b1 = np.polynomial.polynomial.polyfit(m2["log_cap"], m2["log_area"], 1) print(f"log-log 回帰: log(area) = {{b0:.2f}} + {{b1:.2f}} * log(cap)") # b1 ≈ 0.6 → power law: area ∝ cap^0.6 (亜線形スケーリング) ''')}

結果の図 (1) — 規模分布 (要件 L 次元・サイズの混同回避)

図2 を選んだ理由: 規模属性の分布形状を見るのが先。 分布が偏っていれば次の散布図の解釈も変わる (= 散布点が密集する場所が偏る)。 左パネル = 堤高 (生値ヒストグラム) は右に裾を引く正規様。 右パネル = 貯水量 (log10 軸ヒストグラム) は対数正規様で、log scale で見て 初めて分布の中身が読める。

{figure('assets/L45_fig2_size_dist.png', '図2: ため池規模分布 — 堤高は正規様, 貯水量は対数様')}

図2 から読み取れること:

結果の図 (2) — 規模 vs 浸水面積 散布

図3 を選んだ理由: 2 変数の関係を見るには散布図が王道。 log-log 軸を採用するのは「規模も面積も対数で広く分布する」ため。 log-log 軸でべき乗則が直線になることを利用すれば、傾きから スケーリング指数 (= area ∝ cap^β) が読める。

{figure('assets/L45_fig3_size_vs_area.png', '図3: ため池規模 ↔ 決壊時下流浸水想定面積')}

図3 から読み取れること:

結果の表 — 4 種の Pearson r

{df_to_html(T_corr)}

表から読み取れること: log-log 化すると貯水量との r が {corr_cap_area:.3f} → {corr_logcap_logarea:.3f} に上昇するが、 堤高はあまり変わらない ({corr_dam_area:.3f} → {corr_logdam_logarea:.3f})。 これは「貯水量 vs 浸水面積はべき乗則的、堤高 vs 浸水面積は弱い線形」 という構造の違いを示す。水量が直接浸水範囲を決め、堤体の高さは間接的。

""")) # (6) 分析 3: 中山間集中 vs 沿岸 — 市町別ランキング sections.append(("分析 3: 中山間集中 vs 沿岸 — 市町別ランキング", f"""

狙い (RQ ↔ 仮説 H3)

ため池の地理分布を市町別に集計し、中山間 (内陸) 市町と 沿岸市町でどう違うかを定量化する。歴史的にため池は農業遺構なので 内陸の水田地帯に集中するはずである。

手法 — 市町別 集計 + 中山間/沿岸 分類

CSV.所管市町 列で groupby して ため池数・平均堤高・合計貯水量・合計浸水想定 km² を集計する。 さらに各市町を本記事独自定義の中山間 (内陸 9 市町) / 沿岸 (海岸線あり 14 市町) / その他 に分類して カテゴリ別の差を可視化する。

実装コード

{code(''' # 中山間/沿岸の独自分類 (本記事内マッピング) INLAND_CITIES = {"府中市", "三次市", "庄原市", "東広島市", "安芸高田市", "山県郡安芸太田町", "山県郡北広島町", "世羅郡世羅町", "神石郡神石高原町"} COASTAL_CITIES = {"広島市", "呉市", "竹原市", "三原市", "尾道市", "福山市", "大竹市", "廿日市市", "江田島市", "豊田郡大崎上島町", "安芸郡府中町", "安芸郡海田町", "安芸郡熊野町"} # 市町別 集計 (groupby + agg) g = merge.groupby("所管市町").agg( ため池数 =("ため池番号", "count"), 平均堤高_m =("堤高_m", "mean"), 合計貯水量_千m3 =("貯水量_千m3", "sum"), 合計浸水想定_km2 =("area_km2", "sum"), ).reset_index().sort_values("ため池数", ascending=False) g["地理分類"] = g["所管市町"].map( lambda x: "中山間" if x in INLAND_CITIES else ("沿岸" if x in COASTAL_CITIES else "その他")) ''')}

結果の図 (1) — 市町別 ため池数 横棒

図1 を選んだ理由: 23 市町のラベルは長い文字列を含むので horizontal bar が読みやすい。色で中山間/沿岸を区別することで 「上位は中山間中心」が一目で分かる。pie chart や treemap だと比較的順位が分かりにくい。

{figure('assets/L45_fig1_pond_count_by_city.png', '図1: 所管市町別 特定農業用ため池数')}

図1 から読み取れること:

結果の図 (2) — ため池数 vs 浸水想定 散布

図9 を選んだ理由: 「池数」と「浸水想定面積」の関係を見る 2 次元散布。1 池あたりの平均浸水想定面積 (傾き) が中山間と沿岸で違うかを見る。

{figure('assets/L45_fig9_density_scatter.png', '図9: 市町別 ため池数 vs 浸水想定面積')}

図9 から読み取れること:

結果の表 — 市町別 集計 (上位 10 + その他)

{df_to_html(city_stats_view)}

表から読み取れること:

""")) # (7) 分析 4: 連鎖重複の地理学 sections.append(("分析 4: 連鎖重複の地理学 — unary_union vs 個別合計", f"""

狙い (RQ ↔ 仮説 H4)

本記事の核心的発見。各ため池の浸水想定面積を単純に足すと {sum_indiv_km2:.1f} km² だが、unary_union で重複を排除すると {total_union_area:.1f} km²。差分 {overlap_km2:.1f} km² ({overlap_ratio:.1f}%) は 「同じ地点が複数のため池決壊で繰り返し被災想定される」連鎖重複である。 これがどんな地形パターンで発生するかを地図で可視化する。

手法 — shapely.unary_union と差分計算

shapely.unary_union は複数の polygon を1 つの MultiPolygon にマージする (重なりは 1 度だけ計算)。これと個別 polygon の単純合計面積の差が 重複面積になる。

unary_union vs dissolve の違い (要件 J ツール化視点): 本分析は「県全体 1 つの MultiPolygon」が欲しいので unary_union。

実装コード

{code(''' import shapely import geopandas as gpd # 全 7,223 ポリゴン (ため池ごと dissolve 前) を読込 raw = gpd.read_file("data/extras/L45_pond_inundation/shp/ため池浸水想定区域.shp", encoding="cp932").to_crs("EPSG:6671") # 不正幾何を buffer(0) で修復 (30 件) geoms = raw.geometry.values inv = ~shapely.is_valid(geoms) new = [g.buffer(0) if isi else g for g, isi in zip(geoms, inv)] raw["geometry"] = new # unary_union で全 polygon を結合 (~87 秒, 重い処理 → キャッシュ化) union_geom = shapely.unary_union(raw.geometry.values) union_km2 = union_geom.area / 1e6 print(f"県全体 真の浸水想定面積 = {{union_km2:.2f}} km²") # 個別合計 vs union: 差が連鎖重複 individual_km2 = (raw.geometry.area / 1e6).sum() overlap_km2 = individual_km2 - union_km2 print(f"個別合計 = {{individual_km2:.2f}} km², 重複 = {{overlap_km2:.2f}} km² " f"({{overlap_km2/individual_km2*100:.1f}}%)") ''')}

結果の図 (1) — 県全域 浸水想定マップ

図5 を選んだ理由: 全 6,730 池の浸水想定 polygon を 1 枚の地図に重ね、 「線状に集中している場所」を見つけるため。線状集中=複数ため池の 浸水想定が同じ谷に並んで連鎖している証拠。

{figure('assets/L45_fig5_inund_polygon_map.png', '図5: ため池決壊時 浸水想定区域マップ (全 6,730 池)')}

図5 から読み取れること:

結果の図 (2) — 上位 4 市町ズーム small multiples

図6 を選んだ理由: ため池数上位 4 市町 (東広島・福山・庄原・三次) を 同じレイアウトで並べることで市町ごとの分布パターンの違いが直感的に分かる。 small multiples は「同じレンズで複数地域を比較する」ための定番構図。

{figure('assets/L45_fig6_top4_cities_zoom.png', '図6: ため池数上位 4 市町 ズームマップ')}

図6 から読み取れること:

結果の図 (3) — 福山市 詳細マップ

図7 を選んだ理由: 福山平野は少数大池パターンの典型なので、 個別池の名前と浸水想定範囲を 1:1 で対応させて見るのに最適。 カラーマップは「浸水想定面積で濃淡」に塗ることで大池の規模感が掴める。

{figure('assets/L45_fig7_fukuyama_detail.png', '図7: 福山市 ため池決壊浸水想定 詳細マップ')}

図7 から読み取れること:

""")) # (8) 分析 5: 指定の波 — 行政の制度的タイミング sections.append(("分析 5: 指定の波 — 行政の制度的タイミング", f"""

狙い (RQ ↔ 仮説 H5)

「特定農業用ため池」の指定日 (CSV 列「特定農業用ため池の指定」) は令和年月日で記述されている。これを西暦に変換して時系列ヒストグラムを 作り、制度的なピークがいつ発生したかを可視化する。

手法 — 令和年号パーサ + 月単位集計

CSV の指定日は "R3.5.31" のような形式。令和元年 = 2019 として 正規表現で年月日を抽出し、pandas Timestamp に変換する。 月単位で集計してヒストグラムを描き、Top 月をラベル付けする。

実装コード

{code(''' import pandas as pd def reiwa_to_date(s): """R3.5.31 → 2021-05-31""" if pd.isna(s) or not s.startswith("R"): return None parts = s[1:].split(".") if len(parts) != 3: return None r_year, mo, day = map(int, parts) return f"{{2018 + r_year:04d}}-{{mo:02d}}-{{day:02d}}" merge["指定日_西暦"] = merge["特定農業用ため池の指定"].apply(reiwa_to_date) desig = pd.to_datetime(merge["指定日_西暦"].dropna()) # 月単位集計 months = desig.dt.to_period("M").value_counts().sort_index() top_month = months.idxmax() print(f"最大月 = {{top_month}} ({{int(months.max())}} 池)") ''')}

結果の図

図8 を選んだ理由: 月単位棒グラフで時系列のピークと 急増・収束のパターンを見るのに適する。Top 3 月にラベル付けすることで 重要な制度的時点を強調する。

{figure('assets/L45_fig8_designation_timeline.png', '図8: 特定農業用ため池 指定の波')}

図8 から読み取れること:

結果の表 — 指定日別 ため池数 (Top 10)

{df_to_html(T_desig.sort_values('ため池数', ascending=False).head(10))}

表から読み取れること: Top 6 はすべて令和 2-3 年で、それぞれが 法的に意味のある「区切り日」(月末・年度末・経過措置期限) に対応している。 これは指定が個別の池の決壊リスク評価ベースではなく、 行政手続きベースで進められたことを示す。

""")) # (9) 分析 6: Top 浸水想定 ため池 — 規模の代表例 sections.append(("分析 6: Top 浸水想定ため池 — 規模の代表例", f"""

狙い (RQ ↔ 仮説 H2 補強)

浸水想定面積が大きい上位 15 池を抽出し、規模属性 (堤高・貯水量) との対応 を表で確認する。仮説 H2 (規模 ↔ 面積 正相関) の個別事例として、 最上位の池がどんな歴史的背景を持つかを定性的に補強する。

手法 — sort_values で Top N

結合済 merge 表で area_km2 列の降順ソートし上位 15 件を抽出。 特定指定日も合わせて表示することで、「最大級の池がいつ指定されたか」 の情報も得る。

実装コード

{code(''' top_inund = merge.sort_values("area_km2", ascending=False).head(15) print(top_inund[["ため池名称", "所管市町", "堤高_m", "貯水量_千m3", "area_km2"]].to_string()) ''')}

結果の表 — Top 15 浸水想定ため池

{df_to_html(top_inund_view)}

表から読み取れること:

""")) # (10) 仮説検証と考察 sections.append(("仮説検証と考察", f"""

6 仮説の検証結果

{df_to_html(df_hypo)}

主要発見の整理

2 dataset 相互関係の構造発見

本記事の最重要発見は、CSV (62) と Shapefile (63) が「同じため池群の異なる断面」 として完全結合できる一方、カバレッジと表現の質が大きく異なる点である。

2 件は「ため池の identity (CSV) と impact (SHP)」として相補的であり、 結合することで初めて「規模 ↔ 浸水面積」のスケーリングと 「連鎖重複の地理学」という研究的問いに答えられる。 単独の dataset では得られない知見が、結合してはじめて見える典型例。

さらに 2018 年西日本豪雨後の制度進化が R3.5.31 のピークとして データに刻まれており、「災害 → 法律 → 指定 → データ整備」という 社会的時間軸が 2 dataset の更新タイミング差 (CSV は随時更新、SHP は 2025-01-22 と最近) としても観察できる。

本研究の限界

""")) # (11) 発展課題 sections.append(("発展課題", f"""

本記事の結果から論理的に導かれる新仮説と、それを検証する具体的課題:

課題 1: 連鎖決壊シミュレーションの数値モデル化

課題 2: 浸水想定区域内の住民・建物の空間結合

課題 3: ため池決壊 × 河川浸水 × 高潮 「3 機構ハザード重ね」

課題 4: 三次市・庄原市 24 件 非カバレッジ池の現地調査

課題 5: ため池規模分布のべき乗則 vs 対数正規 統計検定

""")) html = render_lesson( num=45, title="L45 ため池浸水想定区域 2 件統合分析 — 広島県の防災重点ため池決壊リスクを「規模属性 × 浸水範囲」で読み解く", tags=["ため池", "決壊シナリオ", "防災重点ため池", "西日本豪雨", "農業用ため池法", "geopandas", "Shapefile", "中山間"], time="60〜90 分", level="リテラシ〜中級 (geopandas/CSV-Shapefile 統合/log-log スケーリング)", data_label="DoBoX dataset 62, 63 (ため池基本情報 + ため池浸水想定区域情報) — CSV 1.1 MB + Shapefile 81 MB", sections=sections, script_filename="L45_pond_inundation.py", ) out_html = LESSONS / "L45_pond_inundation.html" out_html.write_text(html, encoding="utf-8") print(f" HTML written: {out_html} ({out_html.stat().st_size:,} B), {time.time()-t1:.1f}s", flush=True) print(f"\n=== L45 done in {time.time()-t_all:.1f}s ===", flush=True)