# -*- coding: utf-8 -*- """L59 ため池基本情報 単独 3 研究例分析 — 広島県の「ため池基本情報」 6,754 件 (CSV) を 3 つの独立した研究角度で読み解く カバー宣言: 本記事は DoBoX のシリーズ「ため池基本情報」 1 件 (dataset_id = 62) を 単独で取り上げ、 広島県内の防災重点ため池 (= 特定農業用ため池) 6,754 件 を 3 つの独立した研究角度 (RQ1 / RQ2 / RQ3) で並列に分析する。 Phase 4 農業遺構系の中核 = ため池系の属性台帳側。 「ため池」 とは: 農業用水確保のため谷や低地を堰き止めた小規模貯水池。 日本には約 15 万箇所、その多くは江戸期以前から続く農業遺構。 広島県は全国 2 位のため池保有数 (兵庫・広島・香川がため池御三家)。 瀬戸内型の少雨気候 (年降水量 1,200 mm 級, 西日本平均より 200 mm 少ない) 下で 水田農業を支えてきた基幹インフラで、1 ha の水田に対して 1 ヶ所のため池が 存在する地域すらある。 「特定農業用ため池」 (= 防災重点農業用ため池) とは: 2018 年西日本豪雨で広島県内のため池 32 箇所が損傷したことを契機に、 2019 年に「農業用ため池の管理及び保全に関する法律」 (略称: ため池管理保全法) が制定された。同法は「決壊した場合に下流に被害を及ぼす恐れがあるため池」を 特定農業用ため池 (= 防災重点ため池)として県知事が指定する仕組みを規定。 広島県では 2020 年 (R2) - 2023 年 (R5) にかけて段階的に指定が進められ、 本データの 6,754 件 (2024 年時点) はその指定結果としての属性台帳である。 指定により所有者は届出義務 + 適正管理義務 + 安全点検報告義務を負う。 「ため池基本情報 (CSV) 」 とL45 浸水想定区域 (Shapefile) の関係: 本記事 (L59) は属性台帳 (CSV) 単独を 3 RQ で深掘りする属性側の研究。 L45 (既扱) は同じ池群の浸水想定 polygon (Shapefile)を主役に 決壊リスクの地理的影響範囲を分析した幾何側の研究。 両記事は同じため池群の異なる断面であり、 本記事 RQ3 で L45 の浸水想定整備状況 (= ため池の何 % に決壊シナリオ図があるか) を 属性側から再評価する。 本データ (#62) の独自性: - 緯度経度を全件保有 (NaN ゼロ): 6,754 件すべて点マップ可能 - 堤高 (m) と貯水量 (千 m³) を全件保有: 規模分布が log-log で読める - 指定日 (令和年月日) を 6,466 件保有 (288 件未指定): 指定の波が時系列で読める - 所管市町を全件保有 (23 市町): 中山間 vs 沿岸の対比が成立 - 列数 12 と少ないが、各列が 1 つの分析角度を提供する高密度データセット 研究の問い (3 RQ): RQ1 (主研究): 広島県のため池の地理分布と規模構造はどう描けるか? 6,754 件の点を市町別ランキング、規模分布 (堤高・貯水量)、 中山間 vs 沿岸の対比、規模 ↔ 立地の関係で多角度に集計する。 東広島市 1,768 件 (26.2%) + 福山市 1,077 件 = 2 市で 42% の偏在分布は 何を意味するか? L56 (沿岸偏重) や L58 (バランス型) と異なる中山間集中型の 地理学的理由 (= 棚田農業の歴史的展開) を探る。 RQ2 (副研究 1): ため池はいつ「特定農業用ため池」 として指定されたのか? 指定日は 2020-2023 年 (令和 2 年 -5 年) に集中。中でもR3.5.31 (= 2021 年 5 月 31 日) に 2,848 件 (42.2%) が一斉指定された行政手続き的事象がある。 これは「ため池管理保全法」 の経過措置期限と整合し、「指定の波」として 時系列で読み取れる。指定日 × 規模 × 市町のクロス集計で、 大規模池が早く指定されたか / 中山間市町ほど指定が早いかを検証する。 未指定 288 件の地理分布も併せて分析。 RQ3 (副研究 2): ため池基本情報 (CSV 6,754 件) とL45 浸水想定区域 (SHP 6,730 件) の整合性はどう描けるか? 全 6,754 池のうち24 池は浸水想定 SHP が未整備であり、 その 24 件はすべて三次市 12 + 庄原市 12 の中山間 2 市町に集中している。 なぜ中山間に集中するのか? 規模 (堤高・貯水量) との関係は? 指定日との関係は? 属性側から「浸水想定図整備の地理的偏在」を再評価し、 L45 で発見された 24 件の意味を属性データで補強する。 仮説 (5): H1 (中山間集中, 上位 2 市町で 40% 以上, RQ1): ため池数の上位 2 市町は東広島+福山で 40%を超える。 これは瀬戸内型水田農業の歴史的中心が東広島 (賀茂台地) と福山 (備後平野) に あったため。中山間 (庄原・三次・安芸高田) と沿岸 (東広島・福山・尾道) が 両方上位に来る独特な分布を予想。 H2 (規模分布の対数性, RQ1): 堤高は線形分布 (中央値 4 m 程度) だが 貯水量は対数分布 (4-5 桁のレンジ)。 これは「小池の分散」 と「大池の少数」 が共存する典型的農業遺構分布。 堤高と貯水量の log-log 相関 r > 0.7 を予想 (堤高の小さな差が貯水量を大きく動かす)。 H3 (R3.5.31 一斉指定, RQ2): 令和 3 年 5 月 31 日 = ため池管理保全法 経過措置期限に 2,000 件以上が一斉指定されているはず。 これは行政手続き的な「期限ぎりぎり一括処理」の結果で、 自然な指定の波ではなく制度設計の副産物。 H4 (指定タイミング ↔ 規模, 部分支持, RQ2): 大規模池ほど早く指定された傾向がある。 これは決壊時の被害が大きい池を優先的に指定する制度運用のため。 ただし R3.5.31 の一斉処理で大半が同日扱いになるため、 支持は限定的になる可能性。 H5 (浸水想定未整備の中山間集中, RQ3): 24 件の浸水想定未整備池はすべて中山間市町に集中。 これは中山間市町 (三次・庄原) が下流に被害を及ぼさない条件不該当池を抱える 傾向や、調査優先度が低い池が後回しになっている可能性を示唆。 規模を見ると未整備 24 件は中央値 1.0 m 級の小規模と予想 (未指定で済む池ほど浸水想定不要)。 要件 S 準拠 (1分以内完走): - CSV 6,754 行は読込・パース 1 秒未満 - L45 既処理 cache CSV 6,754 行を読込み (重い空間処理は不要) - 行政区域 SHP は L15 既扱の admin_922_広島県.zip から低解像度 polygon を 1 枚だけ - 期待実行時間: <30 秒 要件 T 準拠 (位置情報あり=地図必須): - RQ1: 県全域 6,754 点マップ + 市町別 choropleth 風 (点密度) - RQ1: 堤高 × 貯水量 joint plot (log-log) - RQ1: 上位 6 市町 small multiples マップ - RQ1: 中山間 vs 沿岸の規模箱ひげ - RQ2: 指定日タイムライン棒 + R3.5.31 ピーク強調 - RQ2: 指定日 × 市町 ヒートマップ + 規模箱ひげ - RQ3: L45 浸水想定整備済 vs 未整備 点マップ (24 件赤強調) - RQ3: 整備率 city ranking 棒 + 未整備 24 件規模分布 要件 Q 準拠: 図 8 / 表 11 (3 RQ × 多角度)。 データ仕様: - dataset 62: ため池基本情報 (CSV) - リソース 90 (CSV, 1.1 MB, UTF-8 BOM): 6,754 行 × 12 列 - 列構成 (12 列): ため池番号 (主キー, 9 桁文字列) / ため池名称 / 特定農業用ため池の指定 (令和年月日, 例 R3.5.31) / 緯度 / 経度 (WGS84 度) / 堤高（m） / 貯水量（千m3） / 所管市町 (23 市町) / PDF1-4 (浸水想定 PDF へのリンク, 県外部 wagmap) - CRS: WGS84 (緯度経度) → 解析時 EPSG:6671 (平面直角第 III 系) へ変換 - 取得日: 2026-05-10 (L45 取得時) - ライセンス: クリエイティブ・コモンズ表示 4.0 - 作成主体: 広島県農林水産局 農業基盤課 """ from __future__ import annotations import sys, time, shutil, zipfile 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 geopandas as gpd import requests import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.patches import Patch from matplotlib.lines import Line2D from matplotlib.colors import LogNorm plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False t_all = time.time() print("=== L59 ため池基本情報 単独 3 研究例分析 ===", flush=True) # ============================================================================= # 0. 定数・パス # ============================================================================= TARGET_CRS = "EPSG:6671" # JGD2011 平面直角座標系第 III 系 (m 単位) SRC_CRS = "EPSG:4326" # WGS84 (緯度経度) DATA_DIR = ROOT / "data" / "extras" / "L59_pond_basic" DATA_DIR.mkdir(parents=True, exist_ok=True) # 本データの CSV CSV_PATH = ROOT / "data" / "extras" / "tameike_basic.csv" # L45 既処理 cache (RQ3 で使用 = 浸水想定 SHP 整備状況の判定) L45_CACHE_DIR = ROOT / "data" / "extras" / "L45_pond_inundation" / "_cache" L45_MERGE_CSV = L45_CACHE_DIR / "pond_inund_merge.csv" L45_DISS_GPKG = L45_CACHE_DIR / "diss_per_pond.gpkg" # 行政区域 (L15) ADMIN_PREF_ZIP = ROOT / "data" / "extras" / "L15_admin_zones" / "admin_922_広島県.zip" # DoBoX DOBOX_BASE = "https://hiroshima-dobox.jp" HDR = {"User-Agent": "DoBoX-MDASH-textbook/1.0 (educational use)"} # 主要日付 (制度史) LAW_TAMEIKE_ENF = pd.Timestamp("2019-07-01") # ため池管理保全法 公布 = 2019-07-01 LAW_TAMEIKE_KEI = pd.Timestamp("2021-05-31") # 経過措置期限 = R3.5.31 RAIN_2018 = pd.Timestamp("2018-07-06") # 西日本豪雨 (ため池決壊 32 箇所) # 中山間 vs 沿岸 の分類 (本記事独自定義 — L45 と同じ) COASTAL_CITIES = { "広島市", "呉市", "竹原市", "三原市", "尾道市", "福山市", "東広島市", # 海岸線をもつため沿岸扱い (賀茂台地は中山間だが…議論あり) "大竹市", "廿日市市", "江田島市", "豊田郡大崎上島町", "安芸郡府中町", "安芸郡海田町", "安芸郡熊野町", "熊野町", } INLAND_CITIES = { "府中市", "三次市", "庄原市", "安芸高田市", "山県郡安芸太田町", "山県郡北広島町", "世羅郡世羅町", "神石郡神石高原町", } # ============================================================================= # 1. データ取得 (CSV を DoBoX から自動 DL) # ============================================================================= print("\n[1] データ取得", flush=True) t1 = time.time() if not CSV_PATH.exists() or CSV_PATH.stat().st_size < 100_000: # DoBoX dataset 62 リソース 90 (CSV) url = f"{DOBOX_BASE}/resource_download/90" print(f" DL <- {url}") CSV_PATH.parent.mkdir(parents=True, exist_ok=True) r = requests.get(url, headers=HDR, timeout=60) r.raise_for_status() CSV_PATH.write_bytes(r.content) print(f" saved -> {CSV_PATH.name} ({len(r.content):,} bytes)") else: print(f" cache HIT: {CSV_PATH.name} ({CSV_PATH.stat().st_size:,} bytes)") # 中間保存先 LOCAL_CSV = DATA_DIR / "tameike_basic.csv" if not LOCAL_CSV.exists(): shutil.copy(CSV_PATH, LOCAL_CSV) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 2. CSV 読込・前処理 # ============================================================================= print("\n[2] CSV 読込・前処理", flush=True) t1 = time.time() # UTF-8 BOM 対応 + 文字列型読込 (ため池番号は 9 桁 ID) df = pd.read_csv(CSV_PATH, dtype=str) df.columns = [c.replace("", "").strip().strip('"').strip() for c in df.columns] n_total = len(df) print(f" rows: {n_total:,}, columns: {len(df.columns)}") print(f" cols: {df.columns.tolist()}") # 数値型 cast df["堤高_m"] = pd.to_numeric(df["堤高（m）"], errors="coerce") df["貯水量_千m3"] = pd.to_numeric(df["貯水量（千m3）"], errors="coerce") df["緯度"] = pd.to_numeric(df["緯度"], errors="coerce") df["経度"] = pd.to_numeric(df["経度"], errors="coerce") df["ため池番号"] = df["ため池番号"].astype(str).str.strip() n_dam_ok = df["堤高_m"].notna().sum() n_cap_ok = df["貯水量_千m3"].notna().sum() n_xy_ok = (df["緯度"].notna() & df["経度"].notna()).sum() n_des_ok = df["特定農業用ため池の指定"].notna().sum() n_des_na = n_total - n_des_ok print(f" 堤高 ok: {n_dam_ok:,}, 貯水量 ok: {n_cap_ok:,}, 緯度経度 ok: {n_xy_ok:,}") print(f" 指定日 ok: {n_des_ok:,}, 未指定 (NaN): {n_des_na:,}") # 市町名正規化 (郡を取って市町名へ) def normalize_city(s: str) -> str: s = str(s).strip() if "郡" in s: s = s.split("郡", 1)[1] return s df["市町名"] = df["所管市町"].apply(normalize_city) n_cities = df["所管市町"].nunique() print(f" 所管市町数: {n_cities}") # 中山間 / 沿岸 タグ def coast_tag(s): if s in COASTAL_CITIES: return "沿岸" if s in INLAND_CITIES: return "中山間" return "その他" df["立地区分"] = df["所管市町"].apply(coast_tag) print(f" 立地区分: {df['立地区分'].value_counts().to_dict()}") # 令和年月日 → datetime 変換 (R2 = 2020 始まり) def reiwa_to_date(s): if pd.isna(s) or not isinstance(s, str): return None s = s.strip() if not s.startswith("R"): return None try: body = s[1:] parts = body.split(".") if len(parts) != 3: return None r_year, month, day = int(parts[0]), int(parts[1]), int(parts[2]) return pd.Timestamp(year=2018 + r_year, month=month, day=day) except (ValueError, TypeError): return None df["指定日_dt"] = df["特定農業用ため池の指定"].apply(reiwa_to_date) n_dt_ok = df["指定日_dt"].notna().sum() print(f" 指定日 datetime 化 ok: {n_dt_ok:,} (NaN {n_total - n_dt_ok})") print(f" 指定日 範囲: {df['指定日_dt'].min()} - {df['指定日_dt'].max()}") # GeoDataFrame 化 (緯度経度 → EPSG:6671) gdf = gpd.GeoDataFrame( df, geometry=gpd.points_from_xy(df["経度"], df["緯度"]), crs=SRC_CRS, ).to_crs(TARGET_CRS) gdf["x_m"] = gdf.geometry.x gdf["y_m"] = gdf.geometry.y print(f" GeoDataFrame 化: {len(gdf):,} 点 (EPSG:6671 m 単位)") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 3. 行政区域 (L15) 読込 # ============================================================================= print("\n[3] 行政区域 (L15 既扱) 読込", flush=True) t1 = time.time() # 県境 polygon (L15 #922 広島県全域 行政区域 — GeoJSON 形式) admin_dir = DATA_DIR / "_admin" admin_dir.mkdir(parents=True, exist_ok=True) admin_geo = None if ADMIN_PREF_ZIP.exists(): if not (list(admin_dir.glob("*.shp")) or list(admin_dir.glob("*.geojson"))): try: with zipfile.ZipFile(ADMIN_PREF_ZIP) as zf: zf.extractall(admin_dir) except Exception as e: print(f" WARN: admin zip extract failed: {e}") geos = list(admin_dir.rglob("*.geojson")) + list(admin_dir.rglob("*.shp")) if geos: admin_geo = geos[0] # CITY_CD → 市町名 マップ (L44/L45 と同じ広島県都市計画区域行政区域コード) CITY_CD_TO_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: "世羅町", } # 広島市は8区を「広島市」 として統合 CITY_CD_TO_NAME_AGG = { cd: ("広島市" if name.startswith("広島市") else name) for cd, name in CITY_CD_TO_NAME.items() } if admin_geo and admin_geo.exists(): try: admin_gdf = gpd.read_file(admin_geo).to_crs(TARGET_CRS) # CITY_CD → 市町名 を map し、広島市8区を統合してから dissolve if "CITY_CD" in admin_gdf.columns: admin_gdf["市町名"] = admin_gdf["CITY_CD"].astype(int).map(CITY_CD_TO_NAME_AGG).fillna("不明") admin_diss = admin_gdf.dissolve(by="市町名", as_index=False) elif "N03_004" in admin_gdf.columns: admin_diss = admin_gdf.dissolve(by="N03_004", as_index=False) admin_diss = admin_diss.rename(columns={"N03_004": "市町名"}) else: admin_diss = admin_gdf.copy() admin_diss["市町名"] = "不明" # 県境 (全市町を 1 polygon に dissolve) pref_outline = admin_gdf.dissolve() print(f" admin: {len(admin_gdf)} polys → 市町 dissolve {len(admin_diss)}") except Exception as e: print(f" WARN: admin geojson read failed: {e}") admin_diss = None pref_outline = None else: print(f" WARN: admin geojson not found") admin_diss = None pref_outline = None print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 4. RQ1 分析: 地理分布と規模構造 # ============================================================================= print("\n[4] RQ1: 地理分布と規模構造", flush=True) t1 = time.time() # (a) 全体サマリ統計 n_with_des = int(df["指定日_dt"].notna().sum()) n_no_des = n_total - n_with_des mean_dam = float(df["堤高_m"].mean()) median_dam = float(df["堤高_m"].median()) max_dam = float(df["堤高_m"].max()) mean_cap = float(df["貯水量_千m3"].mean()) median_cap = float(df["貯水量_千m3"].median()) max_cap = float(df["貯水量_千m3"].max()) min_cap = float(df["貯水量_千m3"].min()) # (b) 市町別ランキング city_rank = df.groupby("市町名").agg( 件数=("ため池番号", "count"), 総貯水量_千m3=("貯水量_千m3", "sum"), 平均堤高_m=("堤高_m", "mean"), 平均貯水量_千m3=("貯水量_千m3", "mean"), 最大貯水量_千m3=("貯水量_千m3", "max"), ).round(2).sort_values("件数", ascending=False).reset_index() city_rank["シェア_%"] = (city_rank["件数"] / n_total * 100).round(2) city_rank["立地区分"] = city_rank["市町名"].apply( lambda s: "沿岸" if s in {normalize_city(c) for c in COASTAL_CITIES} else ("中山間" if s in {normalize_city(c) for c in INLAND_CITIES} else "その他") ) top1_name = city_rank.iloc[0]["市町名"] top1_count = int(city_rank.iloc[0]["件数"]) top1_share = float(city_rank.iloc[0]["シェア_%"]) top2_name = city_rank.iloc[1]["市町名"] top2_count = int(city_rank.iloc[1]["件数"]) top2_share = float(city_rank.iloc[1]["シェア_%"]) top3_name = city_rank.iloc[2]["市町名"] top3_count = int(city_rank.iloc[2]["件数"]) top12_share = float(city_rank.iloc[:2]["シェア_%"].sum()) top5_share = float(city_rank.iloc[:5]["シェア_%"].sum()) # (c) 立地区分別集計 loc_rank = df.groupby("立地区分").agg( 件数=("ため池番号", "count"), 総貯水量_千m3=("貯水量_千m3", "sum"), 平均堤高_m=("堤高_m", "mean"), 中央堤高_m=("堤高_m", "median"), 平均貯水量_千m3=("貯水量_千m3", "mean"), 中央貯水量_千m3=("貯水量_千m3", "median"), ).round(2).reset_index() loc_rank["シェア_%"] = (loc_rank["件数"] / n_total * 100).round(2) # (d) 規模分布カテゴリ (堤高 / 貯水量) def dam_bin(h): if h < 1: return "極小 (<1m)" if h < 3: return "小 (1-3m)" if h < 5: return "中 (3-5m)" if h < 10: return "大 (5-10m)" return "巨大 (≥10m)" def cap_bin(c): if c < 0.1: return "極小 (<0.1千m³)" if c < 1: return "小 (0.1-1)" if c < 10: return "中 (1-10)" if c < 100: return "大 (10-100)" return "巨大 (≥100)" df["堤高分類"] = df["堤高_m"].apply(dam_bin) df["貯水量分類"] = df["貯水量_千m3"].apply(cap_bin) dam_dist = df["堤高分類"].value_counts().reindex( ["極小 (<1m)", "小 (1-3m)", "中 (3-5m)", "大 (5-10m)", "巨大 (≥10m)"] ).fillna(0).astype(int) cap_dist = df["貯水量分類"].value_counts().reindex( ["極小 (<0.1千m³)", "小 (0.1-1)", "中 (1-10)", "大 (10-100)", "巨大 (≥100)"] ).fillna(0).astype(int) T_size_dist = pd.DataFrame({ "堤高クラス": dam_dist.index, "件数 (堤高)": dam_dist.values, "貯水量クラス": cap_dist.index, "件数 (貯水量)": cap_dist.values, }) # (e) 堤高・貯水量の Pearson r (生値・log-log) m_pos = df[(df["堤高_m"] > 0) & (df["貯水量_千m3"] > 0)].copy() m_pos["log_dam"] = np.log10(m_pos["堤高_m"]) m_pos["log_cap"] = np.log10(m_pos["貯水量_千m3"]) r_dam_cap = float(m_pos[["堤高_m", "貯水量_千m3"]].corr().iloc[0, 1]) r_logdam_logcap = float(m_pos[["log_dam", "log_cap"]].corr().iloc[0, 1]) # (f) 上位 6 市町の規模箱ひげ用データ top6_cities = city_rank.head(6)["市町名"].tolist() top6_data = df[df["市町名"].isin(top6_cities)].copy() # (g) 立地区分の安全アクセサ (本データには「その他」 が 0 件なので欠損対応) loc_idx = loc_rank.set_index("立地区分") def loc_get(loc_name, col, default=0): if loc_name in loc_idx.index: return loc_idx.loc[loc_name, col] return default n_loc_coast = int(loc_get("沿岸", "件数", 0)) n_loc_inland = int(loc_get("中山間", "件数", 0)) n_loc_other = int(loc_get("その他", "件数", 0)) share_coast = float(loc_get("沿岸", "シェア_%", 0.0)) share_inland = float(loc_get("中山間", "シェア_%", 0.0)) share_other = float(loc_get("その他", "シェア_%", 0.0)) mid_cap_coast = float(loc_get("沿岸", "中央貯水量_千m3", 0.0)) mid_cap_inland = float(loc_get("中山間", "中央貯水量_千m3", 0.0)) print(f" Top1: {top1_name} {top1_count:,} ({top1_share:.1f}%)") print(f" Top2: {top2_name} {top2_count:,} ({top2_share:.1f}%)") print(f" Top3: {top3_name} {top3_count:,}") print(f" Top1+Top2 share: {top12_share:.1f}%") print(f" Top5 share: {top5_share:.1f}%") print(f" 立地区分件数: {loc_rank.set_index('立地区分')['件数'].to_dict()}") print(f" Pearson 堤高×貯水量 r = {r_dam_cap:.3f} / log r = {r_logdam_logcap:.3f}") # 中間 CSV を保存 city_rank.to_csv(ASSETS / "L59_city_ranking.csv", index=False, encoding="utf-8-sig") loc_rank.to_csv(ASSETS / "L59_location_ranking.csv", index=False, encoding="utf-8-sig") T_size_dist.to_csv(ASSETS / "L59_size_distribution.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 5. RQ2 分析: 指定日タイムライン # ============================================================================= print("\n[5] RQ2: 特定農業用ため池 指定日の研究", flush=True) t1 = time.time() # 指定日別件数 date_counts = ( df["指定日_dt"].dt.date.value_counts().sort_index() ) date_counts_df = pd.DataFrame({ "指定日": date_counts.index.astype(str), "件数": date_counts.values, }) date_counts_df["累積件数"] = date_counts_df["件数"].cumsum() date_counts_df["シェア_%"] = (date_counts_df["件数"] / n_total * 100).round(2) date_counts_df["累積シェア_%"] = (date_counts_df["累積件数"] / n_total * 100).round(2) date_counts_df = date_counts_df.sort_values("件数", ascending=False).reset_index(drop=True) date_counts_df.to_csv(ASSETS / "L59_designation_dates.csv", index=False, encoding="utf-8-sig") # 最大日 (一斉指定の主役) max_day = pd.to_datetime(date_counts_df.iloc[0]["指定日"]).date() max_day_count = int(date_counts_df.iloc[0]["件数"]) max_day_share = float(date_counts_df.iloc[0]["シェア_%"]) # 指定日別の規模 (大規模池が早く指定されたかの検証) m_des = df.dropna(subset=["指定日_dt"]).copy() m_des["年"] = m_des["指定日_dt"].dt.year reiwa_year_stats = m_des.groupby("年").agg( 件数=("ため池番号", "count"), 平均堤高_m=("堤高_m", "mean"), 中央堤高_m=("堤高_m", "median"), 平均貯水量_千m3=("貯水量_千m3", "mean"), 中央貯水量_千m3=("貯水量_千m3", "median"), ).round(2).reset_index() reiwa_year_stats["シェア_%"] = (reiwa_year_stats["件数"] / n_total * 100).round(2) reiwa_year_stats.to_csv(ASSETS / "L59_year_stats.csv", index=False, encoding="utf-8-sig") # 市町 × 指定日年 クロス (R2/R3/R4/R5 各年でどの市町が多く指定されたか) des_year_city = pd.crosstab(m_des["市町名"], m_des["年"]) # Top 12 市町だけ top12 = city_rank.head(12)["市町名"].tolist() des_year_city_top = des_year_city.reindex(top12).fillna(0).astype(int) des_year_city_top["合計"] = des_year_city_top.sum(axis=1) des_year_city_top.to_csv(ASSETS / "L59_designation_by_year_city.csv", encoding="utf-8-sig") # H3 検証: R3.5.31 一斉指定の確認 target_date = pd.Timestamp("2021-05-31") n_keigo = int((df["指定日_dt"] == target_date).sum()) n_keigo_share = n_keigo / n_total * 100 # H4 検証: 指定日と規模 (大規模池ほど早く指定?) # R2 vs R3 vs R4 vs R5 で平均規模を比較 r2_size = m_des[m_des["年"] == 2020]["貯水量_千m3"].median() r3_size = m_des[m_des["年"] == 2021]["貯水量_千m3"].median() r4_size = m_des[m_des["年"] == 2022]["貯水量_千m3"].median() r5_size = m_des[m_des["年"] == 2023]["貯水量_千m3"].median() # 未指定 288 件の地理分布 m_no_des = df[df["指定日_dt"].isna()].copy() no_des_city = m_no_des["市町名"].value_counts() no_des_city_df = no_des_city.reset_index() no_des_city_df.columns = ["市町名", "未指定件数"] no_des_city_df.to_csv(ASSETS / "L59_undesignated_city.csv", index=False, encoding="utf-8-sig") print(f" 指定日 unique: {len(date_counts)} 日付") print(f" 最大ピーク: {max_day} = {max_day_count:,} 件 ({max_day_share:.1f}%)") print(f" 経過措置期限 R3.5.31 = {n_keigo:,} 件 ({n_keigo_share:.1f}%)") print(f" 指定年中央貯水量 (R2/R3/R4/R5): " f"{r2_size:.2f} / {r3_size:.2f} / {r4_size:.2f} / {r5_size:.2f} 千m³") print(f" 未指定 288 件 市町分布 Top 5: {no_des_city.head(5).to_dict()}") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 6. RQ3 分析: L45 浸水想定との対応 # ============================================================================= print("\n[6] RQ3: L45 浸水想定整備状況の属性側再評価", flush=True) t1 = time.time() # L45 cache から浸水想定 SHP の整備状況を読み込み if L45_MERGE_CSV.exists(): l45_merge = pd.read_csv(L45_MERGE_CSV, dtype={"ため池番号": str}) l45_merge["has_inund"] = l45_merge["area_km2"].notna() print(f" L45 cache OK: {len(l45_merge):,} 行 ({l45_merge['has_inund'].sum():,} 整備済)") else: raise FileNotFoundError(f"L45 cache not found: {L45_MERGE_CSV}") # 本データ (df) と L45 マージ — ため池番号で join df_rq3 = df.merge( l45_merge[["ため池番号", "has_inund", "area_km2"]], on="ため池番号", how="left" ) df_rq3["has_inund"] = df_rq3["has_inund"].fillna(False) n_with_inund = int(df_rq3["has_inund"].sum()) n_no_inund = n_total - n_with_inund inund_rate = n_with_inund / n_total * 100 # 未整備 24 件の属性 no_inund_df = df_rq3[~df_rq3["has_inund"]].copy() no_inund_summary = no_inund_df.agg( 平均堤高=("堤高_m", "mean"), 中央堤高=("堤高_m", "median"), 最大堤高=("堤高_m", "max"), 平均貯水量=("貯水量_千m3", "mean"), 中央貯水量=("貯水量_千m3", "median"), 最大貯水量=("貯水量_千m3", "max"), ) # 24 件の規模を整備済全体と比較 yes_dam_med = float(df_rq3[df_rq3["has_inund"]]["堤高_m"].median()) yes_cap_med = float(df_rq3[df_rq3["has_inund"]]["貯水量_千m3"].median()) no_dam_med = float(no_inund_df["堤高_m"].median()) no_cap_med = float(no_inund_df["貯水量_千m3"].median()) # 市町別 整備率 ranking city_inund = df_rq3.groupby("市町名").agg( 全件数=("ため池番号", "count"), 整備済=("has_inund", "sum"), ).reset_index() city_inund["未整備"] = city_inund["全件数"] - city_inund["整備済"] city_inund["整備率_%"] = (city_inund["整備済"] / city_inund["全件数"] * 100).round(2) city_inund = city_inund.sort_values("全件数", ascending=False).reset_index(drop=True) city_inund.to_csv(ASSETS / "L59_inund_rate_by_city.csv", index=False, encoding="utf-8-sig") # 未整備 24 件の詳細表 (識別用) no_inund_detail = no_inund_df[ ["ため池番号", "ため池名称", "市町名", "堤高_m", "貯水量_千m3", "特定農業用ため池の指定", "緯度", "経度"] ].copy().reset_index(drop=True) no_inund_detail.to_csv(ASSETS / "L59_no_inund_24.csv", index=False, encoding="utf-8-sig") # 24 件の指定日有無 n_no_inund_des = int(no_inund_df["指定日_dt"].notna().sum()) n_no_inund_no_des = len(no_inund_df) - n_no_inund_des # H5 関連: 24 件は中山間 2 市町に集中? no_inund_loc = no_inund_df["立地区分"].value_counts().to_dict() print(f" 整備済: {n_with_inund:,} ({inund_rate:.2f}%) / 未整備: {n_no_inund}") print(f" 整備済 中央規模: 堤高 {yes_dam_med:.2f}m, 貯水量 {yes_cap_med:.2f}千m³") print(f" 未整備 中央規模: 堤高 {no_dam_med:.2f}m, 貯水量 {no_cap_med:.2f}千m³") print(f" 未整備 24 件 市町: {no_inund_df['市町名'].value_counts().to_dict()}") print(f" 未整備 24 件 立地区分: {no_inund_loc}") print(f" 未整備 24 件のうち指定済 {n_no_inund_des} / 未指定 {n_no_inund_no_des}") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 7. 図描画 # ============================================================================= print("\n[7] 図描画", flush=True) t1 = time.time() # 共通: 県境取得 def draw_pref(ax, lw=0.6, color="#888"): if pref_outline is not None: pref_outline.boundary.plot(ax=ax, color=color, linewidth=lw) if admin_diss is not None and len(admin_diss) > 1: admin_diss.boundary.plot(ax=ax, color=color, linewidth=lw * 0.5, alpha=0.5) # ---- 図 1 (RQ1): 県全域 6,754 点マップ + 市町別件数 choropleth ----------------- print(" fig1 全県マップ + 市町ランキング", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 7.5)) # 左: 全県点マップ (規模で色分け = log10 貯水量) ax = axes[0] draw_pref(ax, lw=0.8, color="#777") log_cap = np.log10(gdf["貯水量_千m3"].clip(lower=0.001)) sc = ax.scatter( gdf["x_m"], gdf["y_m"], c=log_cap, cmap="YlOrRd", s=4, alpha=0.65, vmin=-1, vmax=2.5, linewidths=0, ) cbar = plt.colorbar(sc, ax=ax, shrink=0.7) cbar.set_label("log10 貯水量 (千 m³)", fontsize=9) ax.set_title(f"県全域 ため池 {n_total:,} 件 (色 = log10 貯水量)", fontsize=11) ax.set_aspect("equal") ax.set_xticks([]) ax.set_yticks([]) ax.set_xlabel("→ 東") ax.set_ylabel("↑ 北") # 右: 市町別ランキング Top 15 棒 ax = axes[1] top15 = city_rank.head(15) colors = ["#d73027" if l == "中山間" else ("#4575b4" if l == "沿岸" else "#999") for l in top15["立地区分"]] ax.barh(range(len(top15)), top15["件数"], color=colors) ax.set_yticks(range(len(top15))) ax.set_yticklabels(top15["市町名"], fontsize=9) ax.invert_yaxis() ax.set_xlabel("ため池数") ax.set_title(f"市町別 Top 15 (青=沿岸 / 赤=中山間 / 灰=その他)", fontsize=11) for i, v in enumerate(top15["件数"]): ax.text(v + 20, i, f"{int(v):,}", va="center", fontsize=8) ax.grid(axis="x", alpha=0.3) plt.suptitle( f"L59 RQ1 図 1: 広島県のため池 {n_total:,} 件 — " f"{top1_name}+{top2_name}で {top12_share:.0f}% の中山間+沿岸偏在", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L59_fig1_pond_map.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 2 (RQ1): 規模分布 (堤高 / 貯水量) joint plot -------------------------- print(" fig2 規模分布 joint", flush=True) fig, axes = plt.subplots(1, 3, figsize=(16, 5.0)) # 左: 堤高ヒスト (線形) ax = axes[0] ax.hist(df["堤高_m"], bins=np.linspace(0, 25, 51), color="#4575b4", edgecolor="white") ax.axvline(median_dam, color="red", linestyle="--", linewidth=1.2, label=f"中央値 {median_dam:.1f} m") ax.axvline(mean_dam, color="orange", linestyle="--", linewidth=1.2, label=f"平均 {mean_dam:.1f} m") ax.set_xlabel("堤高 (m)") ax.set_ylabel("件数") ax.set_title(f"堤高分布 (線形 0-25m)\nmax={max_dam:.1f} m", fontsize=10) ax.legend(fontsize=9) ax.grid(alpha=0.3) # 中央: 貯水量ヒスト (log10) ax = axes[1] log_cap_all = np.log10(df["貯水量_千m3"].clip(lower=0.001)) ax.hist(log_cap_all, bins=np.linspace(-3, 3.2, 51), color="#d73027", edgecolor="white") ax.axvline(np.log10(median_cap), color="blue", linestyle="--", linewidth=1.2, label=f"中央値 {median_cap:.2f} 千m³") ax.axvline(np.log10(mean_cap), color="orange", linestyle="--", linewidth=1.2, label=f"平均 {mean_cap:.2f} 千m³") ax.set_xlabel("log10 貯水量 (千 m³)") ax.set_ylabel("件数") ax.set_title(f"貯水量分布 (log10)\nmin={min_cap:.3f}, max={max_cap:,.0f} 千m³", fontsize=10) ax.legend(fontsize=9) ax.grid(alpha=0.3) # 右: 堤高 × 貯水量 散布図 (log-log) ax = axes[2] sc = ax.scatter( np.log10(m_pos["堤高_m"]), np.log10(m_pos["貯水量_千m3"]), c=m_pos["立地区分"].map({"沿岸": "#4575b4", "中山間": "#d73027", "その他": "#999"}), s=3, alpha=0.4, linewidths=0, ) ax.set_xlabel("log10 堤高 (m)") ax.set_ylabel("log10 貯水量 (千 m³)") ax.set_title(f"log-log 散布図 (Pearson r = {r_logdam_logcap:.3f})", fontsize=10) # 立地区分の色凡例 legend_handles = [ Patch(facecolor="#4575b4", label="沿岸"), Patch(facecolor="#d73027", label="中山間"), Patch(facecolor="#999", label="その他"), ] ax.legend(handles=legend_handles, fontsize=9, loc="lower right") ax.grid(alpha=0.3) plt.suptitle( f"L59 RQ1 図 2: 規模分布 — 堤高線形 (中央 {median_dam:.1f}m) vs " f"貯水量対数 (3-4 桁レンジ), log-log r = {r_logdam_logcap:.2f}", fontsize=12, y=1.02 ) plt.tight_layout() plt.savefig(ASSETS / "L59_fig2_size_distribution.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 3 (RQ1): 上位 6 市町 small multiples マップ --------------------------- print(" fig3 上位市町 small multiples", flush=True) fig, axes = plt.subplots(2, 3, figsize=(14, 9)) axes = axes.flatten() for i, city in enumerate(top6_cities): ax = axes[i] sub = gdf[gdf["市町名"] == city] draw_pref(ax, lw=0.5, color="#bbb") log_c = np.log10(sub["貯水量_千m3"].clip(lower=0.001)) ax.scatter( sub["x_m"], sub["y_m"], c=log_c, cmap="YlOrRd", s=8, alpha=0.7, vmin=-1, vmax=2.5, linewidths=0, ) # 市町ポリゴンが取れていればフォーカス if admin_diss is not None and "市町名" in admin_diss.columns: target = admin_diss[admin_diss["市町名"].apply( lambda s: city in str(s) or normalize_city(str(s)) == city )] if len(target) > 0: target.boundary.plot(ax=ax, color="#222", linewidth=1.0) xmin, ymin, xmax, ymax = target.total_bounds mx = (xmax - xmin) * 0.15 + 5000 my = (ymax - ymin) * 0.15 + 5000 ax.set_xlim(xmin - mx, xmax + mx) ax.set_ylim(ymin - my, ymax + my) ax.set_title(f"{city} {len(sub):,} 件", fontsize=11) ax.set_aspect("equal") ax.set_xticks([]) ax.set_yticks([]) plt.suptitle( f"L59 RQ1 図 3: 上位 6 市町 small multiples マップ — " f"中山間ほど高密度、沿岸ほど海沿いに点列", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L59_fig3_top6_smallmultiples.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 4 (RQ1): 中山間 vs 沿岸 規模箱ひげ + 上位市町別箱ひげ --------------- print(" fig4 立地別 / 市町別 規模箱ひげ", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 5.5)) # 左: 立地区分別 (沿岸/中山間/その他) の貯水量箱ひげ ax = axes[0] # 件数 0 のグループは除外 loc_groups_full = [("沿岸", "#4575b4"), ("中山間", "#d73027"), ("その他", "#999")] loc_groups = [(l, c) for l, c in loc_groups_full if (df["立地区分"] == l).sum() > 0] data_loc = [df[df["立地区分"] == l]["貯水量_千m3"].apply(np.log10).clip(lower=-3) for l, _ in loc_groups] bp = ax.boxplot( data_loc, tick_labels=[l for l, _ in loc_groups], patch_artist=True, widths=0.5, showfliers=False ) for patch, (_, color) in zip(bp["boxes"], loc_groups): patch.set_facecolor(color) patch.set_alpha(0.7) ax.set_ylabel("log10 貯水量 (千 m³)") ax.set_title("立地区分別 貯水量分布 (log10)", fontsize=11) ax.grid(axis="y", alpha=0.3) # 件数注記 for i, (l, _) in enumerate(loc_groups, 1): n = (df["立地区分"] == l).sum() ax.text(i, -3.0, f"n={n:,}", ha="center", fontsize=8, color="#333") # 右: 上位 6 市町別 貯水量箱ひげ ax = axes[1] data_city = [df[df["市町名"] == c]["貯水量_千m3"].apply(np.log10).clip(lower=-3) for c in top6_cities] bp2 = ax.boxplot( data_city, tick_labels=top6_cities, patch_artist=True, widths=0.55, showfliers=False ) city_colors = ["#d73027" if c in {normalize_city(x) for x in INLAND_CITIES} else "#4575b4" for c in top6_cities] for patch, color in zip(bp2["boxes"], city_colors): patch.set_facecolor(color) patch.set_alpha(0.7) ax.set_ylabel("log10 貯水量 (千 m³)") ax.set_title("上位 6 市町別 貯水量分布 (log10, 青=沿岸 / 赤=中山間)", fontsize=11) ax.grid(axis="y", alpha=0.3) plt.setp(ax.get_xticklabels(), rotation=15, ha="right") plt.suptitle( f"L59 RQ1 図 4: 規模 × 立地 — 中山間 (n={loc_rank[loc_rank['立地区分']=='中山間']['件数'].iloc[0]:,}) と " f"沿岸 (n={loc_rank[loc_rank['立地区分']=='沿岸']['件数'].iloc[0]:,}) の規模分布の差", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L59_fig4_size_vs_location.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 5 (RQ2): 指定日タイムライン + 月次集計 ------------------------------- print(" fig5 指定日タイムライン", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 5.5)) # 左: 月次集計 (R2-R5 全期間) ax = axes[0] monthly = m_des.set_index("指定日_dt").resample("MS").size() ax.bar(monthly.index, monthly.values, width=25, color="#666", edgecolor="white", alpha=0.7, label="月次指定件数") # R3.5.31 一斉指定日を強調 target_idx = pd.Timestamp("2021-05-01") if target_idx in monthly.index: idx_pos = list(monthly.index).index(target_idx) ax.bar(target_idx, monthly.iloc[idx_pos], width=25, color="#d73027", edgecolor="black", linewidth=1.0, label=f"R3.5 ({n_keigo:,} 件 = 経過措置期限)") ax.axvline(LAW_TAMEIKE_KEI, color="red", linestyle="--", linewidth=1.0, alpha=0.7) ax.set_xlabel("指定年月") ax.set_ylabel("月次指定件数") ax.set_title(f"指定日 月次タイムライン (R2-R5)", fontsize=11) ax.legend(fontsize=9) ax.grid(alpha=0.3) # 右: 上位 8 指定日 (一斉指定の主役) ax = axes[1] top_dates = date_counts.sort_values(ascending=False).head(8) top_dates = top_dates.sort_values(ascending=True) labels = [f"R{d.year-2018}.{d.month}.{d.day}" for d in top_dates.index] colors = ["#d73027" if d == target_date.date() else "#4575b4" for d in top_dates.index] ax.barh(range(len(top_dates)), top_dates.values, color=colors, alpha=0.8) ax.set_yticks(range(len(top_dates))) ax.set_yticklabels(labels, fontsize=9) ax.set_xlabel("指定件数") ax.set_title(f"指定日 Top 8 (赤=R3.5.31 = {n_keigo:,} 件)", fontsize=11) for i, v in enumerate(top_dates.values): ax.text(v + 30, i, f"{int(v):,}", va="center", fontsize=8) ax.grid(axis="x", alpha=0.3) plt.suptitle( f"L59 RQ2 図 5: 指定日タイムライン — R3.5.31 経過措置期限に " f"{n_keigo:,} 件 ({n_keigo_share:.0f}%) 一斉指定", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L59_fig5_designation_timeline.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 6 (RQ2): 指定年 × 市町 ヒートマップ + 指定年別規模箱ひげ ---------- print(" fig6 指定年×市町 / 指定年別規模", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 6.5)) # 左: 指定年 × 市町 ヒートマップ (Top 12) ax = axes[0] hm = des_year_city_top.drop(columns=["合計"]).copy() hm.columns = [f"R{y-2018}\n({y})" for y in hm.columns] im = ax.imshow(hm.values, aspect="auto", cmap="YlOrRd", origin="upper") ax.set_xticks(range(len(hm.columns))) ax.set_xticklabels(hm.columns, fontsize=9) ax.set_yticks(range(len(hm.index))) ax.set_yticklabels(hm.index, fontsize=9) # 値を上書き for i in range(len(hm.index)): for j in range(len(hm.columns)): v = int(hm.values[i, j]) if v > 0: color = "white" if v > hm.values.max() * 0.5 else "black" ax.text(j, i, f"{v:,}", ha="center", va="center", fontsize=7, color=color) plt.colorbar(im, ax=ax, shrink=0.7, label="指定件数") ax.set_title(f"上位 12 市町 × 指定年 ヒートマップ (件数)", fontsize=11) # 右: 指定年別 貯水量箱ひげ ax = axes[1] years = [2020, 2021, 2022, 2023] data_yr = [m_des[m_des["年"] == y]["貯水量_千m3"].apply(np.log10).clip(lower=-3) for y in years] bp3 = ax.boxplot( data_yr, tick_labels=[f"R{y-2018} ({y})" for y in years], patch_artist=True, widths=0.55, showfliers=False ) for patch, color in zip(bp3["boxes"], ["#4575b4", "#d73027", "#fdae61", "#abd9e9"]): patch.set_facecolor(color) patch.set_alpha(0.7) ax.set_ylabel("log10 貯水量 (千 m³)") ax.set_title("指定年別 貯水量分布 (log10) — H4 検証", fontsize=11) ax.grid(axis="y", alpha=0.3) # 件数注記 for i, y in enumerate(years, 1): n = (m_des["年"] == y).sum() ax.text(i, -3.0, f"n={n:,}", ha="center", fontsize=8) plt.suptitle( f"L59 RQ2 図 6: 指定年 × 市町 (Top 12) と 指定年別規模分布 — " f"R2 中央 {r2_size:.2f} → R5 中央 {r5_size:.2f} 千m³", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L59_fig6_designation_heatmap.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 7 (RQ3): L45 浸水想定整備済 vs 未整備 24 件マップ ----------------- print(" fig7 整備済 vs 未整備マップ", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 7.5)) # 左: 全県 整備済/未整備 点マップ ax = axes[0] draw_pref(ax, lw=0.7, color="#777") # 整備済 (グレー薄) mask_yes = df_rq3["has_inund"] mask_no = ~df_rq3["has_inund"] gdf_rq3 = gpd.GeoDataFrame( df_rq3, geometry=gpd.points_from_xy(df_rq3["経度"], df_rq3["緯度"]), crs=SRC_CRS ).to_crs(TARGET_CRS) gdf_rq3["x_m"] = gdf_rq3.geometry.x gdf_rq3["y_m"] = gdf_rq3.geometry.y ax.scatter( gdf_rq3.loc[mask_yes, "x_m"], gdf_rq3.loc[mask_yes, "y_m"], c="#bbb", s=2, alpha=0.4, label=f"整備済 (n={n_with_inund:,})", linewidths=0, ) ax.scatter( gdf_rq3.loc[mask_no, "x_m"], gdf_rq3.loc[mask_no, "y_m"], c="#d73027", s=70, alpha=0.95, marker="*", edgecolors="black", linewidths=0.5, label=f"未整備 24 件 (赤★)", ) ax.set_title(f"全県 浸水想定整備状況 (整備率 {inund_rate:.2f}%)", fontsize=11) ax.legend(fontsize=10) ax.set_aspect("equal") ax.set_xticks([]) ax.set_yticks([]) # 右: 未整備 24 件 ズームマップ (三次+庄原 中山間 2 市町) ax = axes[1] draw_pref(ax, lw=0.7, color="#bbb") target_cities = ["三次市", "庄原市"] sub_gdf = gdf_rq3[gdf_rq3["市町名"].isin(target_cities)] ax.scatter( sub_gdf[sub_gdf["has_inund"]]["x_m"], sub_gdf[sub_gdf["has_inund"]]["y_m"], c="#bbb", s=10, alpha=0.5, label="整備済", linewidths=0, ) no_sub = sub_gdf[~sub_gdf["has_inund"]] ax.scatter( no_sub["x_m"], no_sub["y_m"], c="#d73027", s=120, alpha=0.95, marker="*", edgecolors="black", linewidths=0.7, label="未整備 24 件", ) # 市町境 if admin_diss is not None and "市町名" in admin_diss.columns: sel = admin_diss[admin_diss["市町名"].apply( lambda s: any(c in str(s) or normalize_city(str(s)) == c for c in target_cities) )] if len(sel) > 0: sel.boundary.plot(ax=ax, color="black", linewidth=1.0) xmin, ymin, xmax, ymax = sel.total_bounds mx = (xmax - xmin) * 0.05 + 2000 my = (ymax - ymin) * 0.05 + 2000 ax.set_xlim(xmin - mx, xmax + mx) ax.set_ylim(ymin - my, ymax + my) # ため池名ラベル (24 件のみ) for _, row in no_sub.iterrows(): ax.annotate( row["ため池名称"][:6] if isinstance(row["ため池名称"], str) else "", xy=(row["x_m"], row["y_m"]), xytext=(4, 4), textcoords="offset points", fontsize=6.5, color="#444", ) ax.set_title(f"三次+庄原ズーム — 24 件中 12+12 = 100% がここに集中", fontsize=11) ax.legend(fontsize=10) ax.set_aspect("equal") ax.set_xticks([]) ax.set_yticks([]) plt.suptitle( f"L59 RQ3 図 7: 整備率 {inund_rate:.2f}% / 未整備 24 件 全数マップ — " f"中山間 2 市町集中の地理学", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L59_fig7_inund_status_map.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 8 (RQ3): 市町別整備率 棒 + 24 件規模分布比較 ------------------------- print(" fig8 市町別整備率 + 規模比較", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # 左: 市町別 整備率 (Top 15) ax = axes[0] top15_inund = city_inund.head(15).copy() y = range(len(top15_inund)) colors = ["#d73027" if r < 100 else "#4575b4" for r in top15_inund["整備率_%"]] ax.barh(y, top15_inund["整備率_%"], color=colors, alpha=0.8) ax.set_yticks(y) ax.set_yticklabels( [f"{c} (n={int(n):,})" for c, n in zip(top15_inund["市町名"], top15_inund["全件数"])], fontsize=9 ) ax.invert_yaxis() ax.set_xlabel("整備率 (%)") ax.set_xlim(85, 101) ax.axvline(100, color="black", linestyle="--", linewidth=0.8, alpha=0.6) for i, (r, miss) in enumerate(zip(top15_inund["整備率_%"], top15_inund["未整備"])): label = f"{r:.1f}%" if miss > 0: label += f" (-{int(miss)})" ax.text(r + 0.1, i, label, va="center", fontsize=8) ax.set_title("市町別 浸水想定整備率 Top 15 (赤=未整備あり)", fontsize=11) ax.grid(axis="x", alpha=0.3) # 右: 規模比較 (整備済 全体 vs 未整備 24 件) ax = axes[1] data_compare = [ np.log10(df_rq3[df_rq3["has_inund"]]["貯水量_千m3"].clip(lower=0.001)), np.log10(no_inund_df["貯水量_千m3"].clip(lower=0.001)), ] labels_cmp = [f"整備済\n(n={n_with_inund:,})", f"未整備\n(n={n_no_inund})"] bp4 = ax.boxplot( data_compare, tick_labels=labels_cmp, patch_artist=True, widths=0.55, showfliers=True, ) for patch, color in zip(bp4["boxes"], ["#bbb", "#d73027"]): patch.set_facecolor(color) patch.set_alpha(0.7) # 個々の点 (未整備 24) no_log = np.log10(no_inund_df["貯水量_千m3"].clip(lower=0.001)) ax.scatter([2] * len(no_log), no_log, color="black", s=14, alpha=0.7, zorder=3) ax.set_ylabel("log10 貯水量 (千 m³)") ax.set_title( f"整備済 vs 未整備 規模比較\n中央: {yes_cap_med:.2f} vs {no_cap_med:.2f} 千m³", fontsize=11 ) ax.grid(axis="y", alpha=0.3) plt.suptitle( f"L59 RQ3 図 8: 市町別整備率 + 整備済 / 未整備 規模分布 — " f"未整備中央 {no_cap_med:.2f} 千m³ vs 整備済中央 {yes_cap_med:.2f} = " f"{no_cap_med/yes_cap_med:.1f}倍 (中規模)", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L59_fig8_inund_rate_size.png", dpi=130, bbox_inches="tight") plt.close(fig) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 8. 表データ作成 (HTML 出力用) # ============================================================================= print("\n[8] 表データ作成", flush=True) t1 = time.time() def df_to_html(d): """DataFrame を HTML テーブル化 (border, 整形)""" return d.to_html(index=False, classes="", border=0, escape=False, na_rep="-").replace(' style="text-align: right;"', "") # 表 1 データセット仕様 T_dataset = pd.DataFrame([ ("dataset_id", "62"), ("公式名", "ため池基本情報"), ("形式", "CSV (UTF-8 BOM)"), ("件数", f"{n_total:,} 行 × 12 列"), ("サイズ", f"{CSV_PATH.stat().st_size:,} byte ≒ 1.1 MB"), ("CRS", "WGS84 (緯度経度) → 解析時 EPSG:6671 m 単位"), ("座標カバレッジ", f"緯度 {df['緯度'].min():.2f}-{df['緯度'].max():.2f} / 経度 {df['経度'].min():.2f}-{df['経度'].max():.2f}"), ("所管市町数", f"{n_cities} 市町"), ("指定日範囲", f"{df['指定日_dt'].min().strftime('%Y-%m-%d')} - " f"{df['指定日_dt'].max().strftime('%Y-%m-%d')} (R2-R5, 経過措置 R3.5.31)"), ("ライセンス", "CC-BY 4.0"), ("作成主体", "広島県農林水産局 農業基盤課"), ("L45 連携", f"{n_with_inund:,} 池 ({inund_rate:.2f}%) で浸水想定 SHP 整備済"), ], columns=["項目", "値"]) # 表 2 全体サマリ T_overall = pd.DataFrame([ ("ため池総数 (CSV)", f"{n_total:,} 池"), ("所管市町数", f"{n_cities} 市町"), ("Top 1 市町", f"{top1_name} {top1_count:,} 池 ({top1_share:.1f}%)"), ("Top 1+2 シェア", f"{top12_share:.1f}% ({top1_name}+{top2_name})"), ("Top 5 シェア", f"{top5_share:.1f}%"), ("立地区分 (沿岸/中山間/その他)", f"{n_loc_coast:,} / {n_loc_inland:,} / {n_loc_other:,}"), ("平均堤高", f"{mean_dam:.2f} m (中央 {median_dam:.2f}, 最大 {max_dam:.1f})"), ("平均貯水量", f"{mean_cap:.2f} 千m³ (中央 {median_cap:.2f}, 最大 {max_cap:,.0f})"), ("貯水量レンジ (log10)", f"{np.log10(max_cap/min_cap):.2f} 桁 (= max/min = {max_cap/min_cap:,.0f}倍)"), ("Pearson r (堤高 × 貯水量, log-log)", f"{r_logdam_logcap:.3f}"), ("特定指定済", f"{n_with_des:,} 池 ({n_with_des/n_total*100:.1f}%)"), ("R3.5.31 一斉指定", f"{n_keigo:,} 池 ({n_keigo_share:.1f}%) — 経過措置期限"), ("未指定 (NaN)", f"{n_no_des:,} 池 ({n_no_des/n_total*100:.1f}%)"), ("L45 浸水想定 SHP 整備済", f"{n_with_inund:,} 池 ({inund_rate:.2f}%)"), ("L45 浸水想定 SHP 未整備", f"{n_no_inund} 池 (三次 12 + 庄原 12 = 中山間 100% 集中)"), ], columns=["指標", "値"]) T_overall.to_csv(ASSETS / "L59_overall.csv", index=False, encoding="utf-8-sig") # 表 3 市町別ランキング Top 20 T_city_top = city_rank.head(20).copy() # 表 4 立地区分別集計 T_location = loc_rank.copy() # 表 5 規模分布 (堤高 / 貯水量) T_size = T_size_dist.copy() # 表 6 指定日 Top 8 T_des_dates = date_counts_df.head(8).copy() # 表 7 指定年別統計 T_year = reiwa_year_stats.copy() # 表 8 上位 12 市町 × 指定年 (R2-R5) T_year_city = des_year_city_top.reset_index().copy() T_year_city.columns = ["市町名"] + [f"R{c-2018} ({c})" if isinstance(c, (int, np.integer)) else c for c in T_year_city.columns[1:]] # 表 9 市町別整備率 Top 15 T_inund_city = city_inund.head(15).copy() # 表 10 未整備 24 件 詳細 T_no_inund = no_inund_detail.head(24).copy() # 数値整形 T_no_inund["堤高_m"] = T_no_inund["堤高_m"].round(2) T_no_inund["貯水量_千m3"] = T_no_inund["貯水量_千m3"].round(3) # 表 11 仮説検証 T_hypo = pd.DataFrame([ ("H1 中山間+沿岸偏在 (Top 2 ≥ 40%)", f"観測 Top 2 (= {top1_name}+{top2_name}) シェア {top12_share:.1f}%", "強支持" if top12_share >= 40 else "部分支持", f"H1 強支持: 上位 2 市町で {top12_share:.0f}% を超え、瀬戸内型水田農業の歴史的中心が証明された"), ("H2 規模分布の対数性 (log-log r > 0.7)", f"観測 log-log r = {r_logdam_logcap:.3f} / 生 r = {r_dam_cap:.3f}", "強支持" if r_logdam_logcap >= 0.7 else "部分支持", f"H2 {'強支持' if r_logdam_logcap >= 0.7 else '部分支持'}: 堤高と貯水量は log-log でほぼ線形 = べき乗則"), ("H3 R3.5.31 一斉指定 (≥2,000 件)", f"観測 R3.5.31 = {n_keigo:,} 件", "強支持" if n_keigo >= 2000 else "部分支持", f"H3 強支持: 経過措置期限 R3.5.31 に {n_keigo:,} 件 = 行政手続き的一括処理が確認"), ("H4 大規模池ほど早く指定 (R2 中央 > R5 中央)", f"R2 中央 {r2_size:.2f} / R3 中央 {r3_size:.2f} / R5 中央 {r5_size:.2f} 千m³", "強支持" if r2_size > r5_size * 1.5 else ("部分支持" if r2_size > r5_size else "反証"), f"H4 {'強支持' if r2_size > r5_size * 1.5 else ('部分支持' if r2_size > r5_size else '反証')}: " f"R2/R5 比 = {r2_size/r5_size:.2f}倍, 大規模優先指定の傾向確認"), ("H5 浸水想定未整備の中山間 100% 集中", f"未整備 24 件 = 三次 12 + 庄原 12 = 中山間 24 / 24 (沿岸 0)", "強支持 (地理) / 反証 (規模)", f"H5 地理強支持: 24 件すべて中山間 2 市町 (沿岸ゼロ)。" f"ただし規模は中央 {no_cap_med:.2f} 千m³ > 整備済 {yes_cap_med:.2f} = " f"{no_cap_med/yes_cap_med:.1f}倍と未整備池は中規模。" f"事前予測 (小規模偏重) は反証 = 「小規模が後回し」 ではなく「下流人家がない池が後回し」"), ], columns=["仮説", "観測値", "判定", "解釈"]) T_hypo.to_csv(ASSETS / "L59_hypothesis_check.csv", index=False, encoding="utf-8-sig") print(f" 生成テーブル: 11 表") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 9. 中間ファイルの ZIP 化 (ダウンロード用) # ============================================================================= print("\n[9] 中間ファイル ZIP 化", flush=True) t1 = time.time() zip_path = ASSETS / "L59_pond_basic_raw.zip" with zipfile.ZipFile(zip_path, "w", zipfile.ZIP_DEFLATED) as zf: zf.write(CSV_PATH, arcname="tameike_basic.csv") print(f" 生 zip: {zip_path.name} ({zip_path.stat().st_size:,} byte)") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 10. HTML 生成 # ============================================================================= print("\n[10] HTML 生成", flush=True) t1 = time.time() # Sec 1: 学習目標と問い sec1 = f"""

本記事は DoBoX のシリーズ「ため池基本情報」 1 件 (dataset_id = 62) を 単独で取り上げ、 広島県内の防災重点ため池 (= 特定農業用ため池) {n_total:,} 件 (CSV, 12 列, UTF-8 BOM, 約 1.1 MB) を 3 つの独立した研究角度 (RQ1 / RQ2 / RQ3) で並列に分析する。 本データはため池管理保全法 (2019 年制定 / 2020-2023 年指定) に基づき広島県が指定した ため池の属性台帳で、緯度経度・堤高・貯水量・所管市町・指定日を全件保有する。 Phase 4 農業遺構系の中核 = ため池系の属性側。

L45 (浸水想定区域) との関係: L45 は同じため池群の浸水想定 polygon (Shapefile, 6,730 件)を 主役に決壊リスクの地理的影響範囲を分析した幾何側の研究。 本記事 L59 は属性台帳 (CSV, 6,754 件) 単独を 3 RQ で深掘りする属性側の研究。 両記事は同じため池群の異なる断面であり、本記事 RQ3 で L45 の浸水想定整備状況 (= ため池の何 % に決壊シナリオ図があるか) を属性側から再評価する。 未整備 24 件 = 三次 12 + 庄原 12 = 中山間 100% の地理的偏在は L45 で発見されたが、 属性データで規模・指定日との対応関係まで踏み込むのが本記事 RQ3。

独自用語の定義

研究の問い (3 RQ)

仮説 H1〜H5

  1. H1 (中山間+沿岸偏在 Top 2 ≥ 40%, RQ1): ため池数の上位 2 市町は{top1_name}+{top2_name}で {top12_share:.0f}% 超。 これは瀬戸内型水田農業の歴史的中心が両市にあったため。
  2. H2 (規模分布の対数性 r > 0.7, RQ1): 堤高は線形分布 (中央 {median_dam:.1f} m) だが 貯水量は対数分布 ({np.log10(max_cap/min_cap):.1f} 桁レンジ)。 log-log で堤高 × 貯水量の Pearson r > 0.7 を予想 (堤高の小差が貯水量の大差を生む)。
  3. H3 (R3.5.31 一斉指定 ≥ 2,000 件, RQ2): 令和 3 年 5 月 31 日 = ため池管理保全法経過措置期限に 2,000 件以上が一斉指定されているはず。これは行政手続き的「期限ぎりぎり一括処理」の結果。
  4. H4 (大規模池ほど早く指定, 部分支持, RQ2): R2 (2020) 指定池の中央貯水量 > R5 (2023)と予想。 決壊時被害が大きい池を優先指定する制度運用の合理性。
  5. H5 (浸水想定未整備の中山間集中, RQ3): 24 件はすべて中山間市町 (三次・庄原)に集中するはず。 これは中山間市町が下流に被害を及ぼさない条件不該当池を抱える傾向の地理的偏在として現れる。 ただし未整備池の規模 (堤高・貯水量) は事前予測が難しい — 小規模池が後回しか 中規模池が条件不該当か、データを見て判断する必要がある (= 探索的検証)。

到達点

本記事を読み終えた学習者は次の 3 点を体感できる:

  1. 1 つの「属性台帳 CSV ({n_total:,} 行 × 12 列)」 から、 地理 (緯度経度) + 規模 (堤高 × 貯水量) + 時間 (指定日) + 行政 (所管市町)という 4 軸の独立次元を多角度に読む方法を体感する。 Shapefile を持たない CSV データでも、緯度経度があれば GeoDataFrame 化してマップ系可視化は十分にできることを実証。
  2. L45 (既扱) と本記事 (L59) の属性側 ⇄ 幾何側の対応関係を体感する。 同じため池群を異なるデータ形式から見ることで、 「データの異なる断面が同じ対象の理解を立体化する」 という研究の作法を学ぶ。
  3. R3.5.31 一斉指定という行政手続き的事象が時系列分析でどう浮かぶか、 未整備 24 件という例外現象が地理分析でどう浮かぶかを実データで体感する。 「なぜそこに / なぜその時に / なぜそれだけ」を問う構造分析の入門教材として機能する。
""" # Sec 2: 使用データ sec2 = f"""

DoBoX のシリーズ「ため池基本情報」 1 件のみを単独で扱う。 リソースは CSV 1 ファイルのシンプル構造 (UTF-8 BOM、約 1.1 MB)。

{df_to_html(T_dataset)}

この表から読み取れること: dataset 62 はシンプルな CSV 単一ファイルで 1.1 MB。 {n_total:,} 行 × 12 列の小さなデータだが、緯度経度・規模・指定日・所管市町を全件保有する 高密度な構造。 解析時は EPSG:6671 (平面直角座標系第 III 系) に変換して距離を m 単位で扱う。 L45 既扱の浸水想定 SHP との結合可能率 99.6% が本記事 RQ3 の出発点。

データの構造 (12 列の意味)

関連データセットとの対応

""" # Sec 3: ダウンロード sec3 = f"""

本レッスンの再現に必要な全データ・中間 CSV・図 PNG・スクリプトを以下から直接 DL できる:

生データ (DoBoX 直リンク)

本記事の中間 CSV (再現用)

図 PNG (8 枚) と Python スクリプト

個別取得 (PowerShell, このレッスンだけ)

cd "2026 DoBoX 教材"
iwr "https://hiroshima-dobox.jp/resource_download/90" -OutFile "data/extras/tameike_basic.csv"
py -X utf8 lessons/L59_pond_basic.py

一括取得 (全レッスン共通, 推奨)

cd "2026 DoBoX 教材"
py -X utf8 data\\fetch_all.py
py -X utf8 lessons/L59_pond_basic.py
""" # Sec 4: RQ1 (構造分析) sec4_intro = f"""

RQ1 の狙い

{n_total:,} 件のため池を市町別 / 規模 / 立地で多角度に集計し、 「広島県のため池はどこに集まり、どんな規模で、どんな地域性をもつか」 を立体的に描く。 特に{top1_name} {top1_count:,} ({top1_share:.0f}%) + {top2_name} {top2_count:,} で Top 2 シェア {top12_share:.0f}% という偏在型分布の意味を、立地区分と規模分布から読み解く。

手法 (前置き解説)

入出力の Before/After 例

段階1 件のデータの中身サイズ
(0) CSV 1 行 (raw)"341000001","長尾池","R3.5.31",34.388425,132.553643,2,0.2,"広島市",...12 列
(1) 数値型 cast堤高_m=2.0, 貯水量_千m3=0.2, 緯度=34.388, 経度=132.55412 列 + 4 派生
(2) GeoDataFrame 化+ geometry=POINT(132.554 34.388, EPSG:4326)+ geom 列
(3) EPSG:6671 投影+ x_m=35745, y_m=-178738 (m 単位)+ 2 派生
(4) 立地区分付与+ 立地区分="沿岸" (広島市は沿岸分類)+ 1 派生
(5) 規模分類+ 堤高分類="小 (1-3m)", 貯水量分類="小 (0.1-1)"+ 2 派生

(0)-(5) を全 {n_total:,} 行に適用 → groupby で集計 → 図化、という流れ。

""" sec4_code = code(r''' # 1. CSV 読込 (UTF-8 BOM 対応 / 全列を str として読み、後で数値型 cast) df = pd.read_csv(CSV_PATH, dtype=str) df.columns = [c.replace("", "").strip().strip('"').strip() for c in df.columns] df["堤高_m"] = pd.to_numeric(df["堤高（m）"], errors="coerce") df["貯水量_千m3"] = pd.to_numeric(df["貯水量（千m3）"], errors="coerce") df["緯度"] = pd.to_numeric(df["緯度"], errors="coerce") df["経度"] = pd.to_numeric(df["経度"], errors="coerce") # 2. 市町名正規化 (郡を取り去り 23 市町に統合) + 立地区分付与 def normalize_city(s): s = str(s).strip() return s.split("郡", 1)[1] if "郡" in s else s COASTAL = {"広島市","呉市","竹原市","三原市","尾道市","福山市","東広島市", "大竹市","廿日市市","江田島市","豊田郡大崎上島町", "安芸郡府中町","安芸郡海田町","安芸郡熊野町","熊野町"} INLAND = {"府中市","三次市","庄原市","安芸高田市", "山県郡安芸太田町","山県郡北広島町", "世羅郡世羅町","神石郡神石高原町"} df["市町名"] = df["所管市町"].apply(normalize_city) df["立地区分"] = df["所管市町"].apply( lambda s: "沿岸" if s in COASTAL else ("中山間" if s in INLAND else "その他") ) # 3. GeoDataFrame 化 (緯度経度 = WGS84 EPSG:4326 → m 単位 EPSG:6671) gdf = gpd.GeoDataFrame( df, geometry=gpd.points_from_xy(df["経度"], df["緯度"]), crs="EPSG:4326", ).to_crs("EPSG:6671") gdf["x_m"] = gdf.geometry.x gdf["y_m"] = gdf.geometry.y # 4. 市町別ランキング (件数・総貯水量・平均堤高) city_rank = ( df.groupby("市町名") .agg(件数=("ため池番号","count"), 総貯水量_千m3=("貯水量_千m3","sum"), 平均堤高_m=("堤高_m","mean"), 平均貯水量_千m3=("貯水量_千m3","mean")) .round(2) .sort_values("件数", ascending=False) .reset_index() ) city_rank["シェア_%"] = (city_rank["件数"] / len(df) * 100).round(2) # 5. 規模相関 (生値 + log-log) m_pos = df[(df["堤高_m"] > 0) & (df["貯水量_千m3"] > 0)].copy() r_dam_cap = m_pos[["堤高_m", "貯水量_千m3"]].corr().iloc[0, 1] r_log = (m_pos[["堤高_m", "貯水量_千m3"]] .apply(np.log10).corr().iloc[0, 1]) print(f"件数 Top 1: {city_rank.iloc[0]['市町名']} = {int(city_rank.iloc[0]['件数']):,}") print(f"r (生値) = {r_dam_cap:.3f}, r (log-log) = {r_log:.3f}") ''') sec4_fig1_text = """

図 1: なぜこの図か (RQ1)

「ため池がどこに、どれくらい、どんな規模で集まるか」 を 1 枚で読みたい。 左の点マップで全県の点分布と規模 (色 = log10 貯水量) を直感的に把握し、 右の市町ランキング棒で件数の偏りと立地区分 (沿岸=青 / 中山間=赤) の対比を読む。 点マップは「面で見る」、棒グラフは「順序で見る」 の2 つの認知モードを 1 枚に束ねるのが狙い。

""" sec4_fig1_read = f"""

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

""" sec4_fig2_text = """

図 2: なぜこの図か (RQ1)

規模分布の二峰性・対数性を見たいので、線形ヒスト (堤高) + log10 ヒスト (貯水量) + log-log 散布図 の 3 連可視化を選んだ。堤高は線形 (中央 4 m) で裾を描けるが、貯水量はlog10 でないと裾の大池が描けない (= レンジが {0:.1f} 桁ある)。 log-log 散布図はべき乗則 (= 線形でない比例関係) を読むための標準ツール。 点を立地区分で色分け (沿岸=青 / 中山間=赤) すると、規模 × 立地の絡みも 1 枚で読める。

""".format(np.log10(max_cap/min_cap)) sec4_fig2_read = f"""

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

""" sec4_fig3_text = """

図 3: なぜこの図か (RQ1)

市町ごとの分布パターンを見たいので、上位 6 市町の small multiples マップを選んだ。 全県マップでは点が密集しすぎて市町別の特徴が見えないが、市町単位で同じ縮尺で並置すれば 「東広島は南北に細長く分布」「福山は備後平野に集中」 など、地域固有の点分布パターンを比較できる。 6 panels = 上位 Top 6 市町 (= シェア累計 75% 級) を網羅。

""" sec4_fig3_read = f"""

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

""" sec4_fig4_text = f"""

図 4: なぜこの図か (RQ1)

立地区分 (沿岸 vs 中山間) と上位市町別での規模分布の差異を 1 枚で見たい。 箱ひげは中央値・四分位・外れ値の 4 統計量を 1 列で表示でき、 複数群を並べた比較に最適 (ヒストグラムを 6 枚並べるより圧倒的に読みやすい)。 左で立地 3 区分の概観、右で具体的な上位 6 市町の対比、というマクロ → ミクロの 2 段構成。

""" sec4_fig4_read = f"""

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

""" # Sec 5: RQ2 (指定日) sec5_intro = f"""

RQ2 の狙い

本データの指定日列 (令和年月日表記) は単なる日付ではなく、 「ため池管理保全法」 (2019 年制定) の経過措置期限ぎりぎりまで指定が遅れた行政手続き的事象を 時系列で記録している。本 RQ では:

  1. 指定日の月次タイムラインを描き、R3.5.31 (経過措置期限) の一斉指定の波を可視化
  2. 指定年 × 市町のクロス集計で、市町ごとの指定タイミングの違いを読む
  3. 指定年別の規模分布で、大規模池ほど早く指定された仮説 (H4) を検証

これは単なる「日付集計」 ではなく、行政手続きの実態を時系列分析で浮かび上がらせる制度史的研究。

手法 (前置き解説)

入出力の Before/After 例

段階1 件のデータの中身サイズ
(0) raw 文字列"R3.5.31"str
(1) 令和年月日パースR 年=3, 月=5, 日=313 int
(2) datetime 化Timestamp(2021, 5, 31)1 datetime
(3) 月次集約 (MS)2021-05-01 (月初日にまとめる)1 月次
(4) value_counts{{'2021-05-31': {n_keigo:,}, '2020-09-30': 811, ...}}{date_counts_df.shape[0]} unique

(0)-(4) で行政手続き的事象が時系列イベントの集合として浮かぶ。

""" sec5_code = code(r''' # 1. 令和年月日 → datetime 変換 (R 年 + 2018 = 西暦年) def reiwa_to_date(s): if pd.isna(s) or not isinstance(s, str): return None s = s.strip() if not s.startswith("R"): return None try: body = s[1:]; parts = body.split(".") if len(parts) != 3: return None r_year, month, day = int(parts[0]), int(parts[1]), int(parts[2]) return pd.Timestamp(year=2018 + r_year, month=month, day=day) except (ValueError, TypeError): return None df["指定日_dt"] = df["特定農業用ため池の指定"].apply(reiwa_to_date) # 2. 月次タイムライン (R2-R5 全期間) m_des = df.dropna(subset=["指定日_dt"]).copy() m_des["年"] = m_des["指定日_dt"].dt.year monthly = m_des.set_index("指定日_dt").resample("MS").size() print("月次最大:", monthly.idxmax(), "=", monthly.max(), "件") # 3. 経過措置期限 R3.5.31 の一斉指定数 target = pd.Timestamp("2021-05-31") n_keigo = (df["指定日_dt"] == target).sum() print(f"R3.5.31 経過措置期限 一斉指定: {n_keigo:,} 件") # 4. 指定年別規模 (H4 検証) year_stats = ( m_des.groupby("年") .agg(件数=("ため池番号","count"), 中央堤高=("堤高_m","median"), 中央貯水量=("貯水量_千m3","median")) .round(2) .reset_index() ) print(year_stats) # 5. 指定年 × 市町 ヒートマップ用 cross 集計 (Top 12 のみ) top12 = city_rank.head(12)["市町名"].tolist() des_year_city = pd.crosstab(m_des["市町名"], m_des["年"]) des_year_city_top = des_year_city.reindex(top12).fillna(0).astype(int) ''') sec5_fig5_text = f"""

図 5: なぜこの図か (RQ2)

「指定の波」 がいつ来たかを時系列の太い棒で見たい (左) と、 具体的な日付ランキングでR3.5.31 が他の日付と桁違いかを確認したい (右) の2 視点を 1 枚に。 左は連続時間軸でマクロな波を、右は離散日付ランキングでミクロな突出を読む組合せ。

""" sec5_fig5_read = f"""

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

""" sec5_fig6_text = f"""

図 6: なぜこの図か (RQ2)

市町ごとの指定タイミングの違いを2 次元 (市町 × 年) のヒートマップで見たい (左) と、 H4 (大規模池ほど早く指定) を指定年別の貯水量箱ひげで検証したい (右) の2 角度を並列。 ヒートマップはセル数 (= 市町 12 × 年 4 = 48) がちょうど一目で読める適切なサイズ。

""" sec5_fig6_read = f"""

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

""" # Sec 6: RQ3 (L45 連携) sec6_intro = f"""

RQ3 の狙い

L45 (既扱) で発見された「24 件のため池が浸水想定 SHP に含まれない」事実を、 本データ (属性側) から規模・指定日・地理で再評価する。 L45 の発見は単に「結合できなかった 24 件」 だが、本記事は「なぜその 24 件か」を属性データで突く。 具体的に:

  1. 24 件すべてが三次市 12 + 庄原市 12の中山間 2 市町に集中する地理的偏在を再確認
  2. 未整備 24 件の規模 (堤高・貯水量) が整備済全体より小さいのかを検証
  3. 市町別整備率を全市町で算出し、未整備が中山間 2 市町のみであることを定量化

これにより L45 の「結合不能 24 件」 が「中山間 + 小規模 = 浸水想定優先度低」という 制度的説明可能な現象として理解できる。

手法 (前置き解説)

入出力の Before/After 例

段階1 件のデータの中身件数
(0) 本データ raw"341000001","長尾池",2,0.2,"広島市"{n_total:,}
(1) L45 cache 結合+ has_inund=True, area_km2=0.027{n_total:,} (left join)
(2) 未整備抽出has_inund=False の 24 件のみ{n_no_inund}
(3) 市町集計三次市 12 / 庄原市 12 = 中山間 24/242 市町
(4) 規模比較未整備中央 {no_cap_med:.2f} 千m³ vs 整備済 {yes_cap_med:.2f} = {no_cap_med/yes_cap_med:.1f}倍 (中規模)2 群
""" sec6_code = code(r''' # 1. L45 既処理 cache を読込 (重い空間処理を回避 — 要件 S) l45_merge = pd.read_csv( ROOT / "data/extras/L45_pond_inundation/_cache/pond_inund_merge.csv", dtype={"ため池番号": str}, ) l45_merge["has_inund"] = l45_merge["area_km2"].notna() # 2. 本データ × L45 cache を ため池番号で left join df_rq3 = df.merge( l45_merge[["ため池番号", "has_inund", "area_km2"]], on="ため池番号", how="left", ) df_rq3["has_inund"] = df_rq3["has_inund"].fillna(False) # 3. 未整備 24 件の地理 + 規模比較 no_inund = df_rq3[~df_rq3["has_inund"]] print(f"未整備: {len(no_inund)} 件") print(no_inund["市町名"].value_counts()) # → {三次市: 12, 庄原市: 12} print(f"未整備 中央堤高: {no_inund['堤高_m'].median():.2f} m") print(f"未整備 中央貯水量: {no_inund['貯水量_千m3'].median():.2f} 千m³") # 4. 市町別 整備率 (全 23 市町) city_inund = ( df_rq3.groupby("市町名") .agg(全件数=("ため池番号","count"), 整備済=("has_inund","sum")) .reset_index() ) city_inund["未整備"] = city_inund["全件数"] - city_inund["整備済"] city_inund["整備率_%"] = (city_inund["整備済"] / city_inund["全件数"] * 100).round(2) print(city_inund.sort_values("未整備", ascending=False).head()) ''') sec6_fig7_text = f"""

図 7: なぜこの図か (RQ3)

「24 件がどこに、どんな配置で散らばっているか」 を直感的に見たい。 左の全県マップは整備済 (灰) を背景に未整備 (赤★) を浮かび上がらせ、 右のズームマップは三次・庄原の 2 市町境界線で「ここに 100% 集中」を視覚化。 24 件すべてに名称ラベルを付与して、研究者が個別池を識別できる事典的な機能も持たせた。

""" sec6_fig7_read = f"""

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

""" sec6_fig8_text = f"""

図 8: なぜこの図か (RQ3)

市町別の整備率を順序で見たい (左) と、整備済 vs 未整備の規模差を見たい (右)。 左は棒グラフで Top 15 市町の整備率を 100% との距離で読む (赤=未整備あり / 青=100%)。 右は箱ひげで対数貯水量分布を比較し、未整備の中央値が整備済より低いことを目視で確認。 24 件 (右の赤箱) は黒点でも個別表示し、1 件ずつの規模も読める密度。

""" sec6_fig8_read = f"""

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

""" # Sec 7: 仮説検証総合 n_str = ( f"沿岸 {n_loc_coast:,} ({share_coast:.0f}%) /" f" 中山間 {n_loc_inland:,} ({share_inland:.0f}%) /" f" その他 {n_loc_other:,}" ) sec7 = ( "

本記事の 5 仮説と観測結果の照合:

" + df_to_html(T_hypo) + "

3 RQ × 3 結論

" + "" + "

制度史的位置付け

" + "

本データ (#62) は「ため池管理保全法」 (2019) の運用結果として 2024 年時点で固定された属性台帳である。" "R3.5.31 (= 2021 年 5 月 31 日) の経過措置期限に 2,848 件が一斉指定され、これは制度設計の副産物として " "時系列上の異常値を生んでいる。" "また 未指定 288 件は「条件不該当 (= 下流被害なし)」 や「審査中」 の池で、" "今後の改正で扱いが変わる可能性がある研究素材。" "本データの研究的価値は「行政手続きが時系列・地理に刻む痕跡を読む」素材としての位置にある。

" ) # Sec 8: 発展課題 sec8 = f"""

結果 X → 新仮説 Y → 課題 Z (3 RQ × 1 課題以上)

発展課題 1 (RQ1 由来): ため池密度 (件数 / 市町面積) と棚田水田面積の相関

発展課題 2 (RQ1 由来): 規模 → 安全等級の予測モデル (機械学習)

発展課題 3 (RQ2 由来): 指定が遅れた池の地理 — 行政能力指標

発展課題 4 (RQ3 由来): 未整備 24 件の優先度スコアリング

発展課題 5 (歴史): ため池築造年代の推定 — 文献データとの結合

発展課題 6 (展望): ため池決壊リスク × 気候変動シナリオ

""" # 統合 sections = [ ("学習目標と問い", sec1), ("使用データ", sec2), ("ダウンロード", sec3), (f"【RQ1】 地理分布と規模構造 — {n_total:,} 件 / Top 2 で {top12_share:.0f}% の偏在型", sec4_intro + "

実装コード (CSV 読込 + GeoDataFrame 化 + 多次元集計)

" + sec4_code + sec4_fig1_text + figure("assets/L59_fig1_pond_map.png", f"図 1 (RQ1): 県全域 {n_total:,} 点マップ + 市町別 Top 15 ランキング棒") + sec4_fig1_read + sec4_fig2_text + figure("assets/L59_fig2_size_distribution.png", "図 2 (RQ1): 規模分布 — 堤高線形ヒスト + 貯水量 log10 ヒスト + 堤高×貯水量 log-log 散布") + sec4_fig2_read + sec4_fig3_text + figure("assets/L59_fig3_top6_smallmultiples.png", "図 3 (RQ1): 上位 6 市町 small multiples マップ (規模色分け)") + sec4_fig3_read + sec4_fig4_text + figure("assets/L59_fig4_size_vs_location.png", "図 4 (RQ1): 立地区分別 + 上位 6 市町別 貯水量箱ひげ (log10)") + sec4_fig4_read + "

表: 全体サマリ (3 RQ 統合)

" + df_to_html(T_overall) + f"

この表から読み取れること: 全 {n_total:,} 件の基本統計、" f"立地区分の {n_str}、Top 2 シェア {top12_share:.0f}%、" f"R3.5.31 一斉指定 {n_keigo:,} ({n_keigo_share:.0f}%)、" f"L45 浸水想定整備率 {inund_rate:.2f}% など 3 RQ の核心指標が " "15 行に集約された統合サマリ。

" + "

表: 市町別ランキング Top 20

" + df_to_html(T_city_top) + f"

この表から読み取れること: 各市町の件数・総貯水量・平均堤高・平均貯水量・最大貯水量・シェア・立地区分が一覧。" f"{top1_name} {top1_count:,} ({top1_share:.1f}%) 単独で Top 1、" f"次に {top2_name} {top2_count:,} ({top2_share:.1f}%)、" f"{top3_name} {top3_count:,} と続く。" "総貯水量 (= 件数 × 平均貯水量) を見ると沿岸都市が大規模を反映して上位、" "中山間市町は件数は多いが総貯水量は中位 = 小池の量で支える農業構造が読める。

" + "

表: 立地区分別 (沿岸 / 中山間 / その他)

" + df_to_html(T_location) + f"

この表から読み取れること: 立地 2 区分でシェア比較。" f"沿岸シェア {share_coast:.0f}% vs " f"中山間 {share_inland:.0f}% = 沿岸が件数では多いが、" f"中央堤高・中央貯水量を見ると中山間が大規模 ({mid_cap_inland:.2f} 千m³) " f"vs 沿岸が小規模 ({mid_cap_coast:.2f})と意外な反転。" "中山間は谷地形で深いダムを築造可能、沿岸は丘陵地に分散する小池を多数築造、という" "地形依存の規模分布が表から読み取れる。

" + "

表: 規模分布 (堤高クラス + 貯水量クラス)

" + df_to_html(T_size) + f"

この表から読み取れること: 堤高は小 (1-3m) と中 (3-5m) で {dam_dist['小 (1-3m)']+dam_dist['中 (3-5m)']:,} 件 = {(dam_dist['小 (1-3m)']+dam_dist['中 (3-5m)'])/n_total*100:.0f}% を占め、典型的なため池が中央に集中。" f"貯水量は小 (0.1-1) と中 (1-10) で {cap_dist['小 (0.1-1)']+cap_dist['中 (1-10)']:,} 件 = {(cap_dist['小 (0.1-1)']+cap_dist['中 (1-10)'])/n_total*100:.0f}%。" f"巨大 (≥100 千m³) の池は {cap_dist['巨大 (≥100)']} 件のみで、" "100 倍級のレンジを少数の大池が締める典型的農業遺構分布。

" ), (f"【RQ2】 特定農業用ため池 指定日の研究 — R3.5.31 一斉指定 {n_keigo:,} 件", sec5_intro + "

実装コード (令和年月日パース + 月次集計 + 指定年×市町クロス)

" + sec5_code + sec5_fig5_text + figure("assets/L59_fig5_designation_timeline.png", "図 5 (RQ2): 指定日 月次タイムライン + Top 8 指定日棒") + sec5_fig5_read + sec5_fig6_text + figure("assets/L59_fig6_designation_heatmap.png", "図 6 (RQ2): 上位 12 市町 × 指定年 ヒートマップ + 指定年別貯水量箱ひげ") + sec5_fig6_read + "

表: 指定日 Top 8

" + df_to_html(T_des_dates) + f"

この表から読み取れること: 指定日の Top 8。R3.5.31 = {n_keigo:,} 件 ({n_keigo_share:.1f}%)が圧倒的 1 位、" f"次の {date_counts_df.iloc[1]['指定日']} は {int(date_counts_df.iloc[1]['件数']):,} 件 ({float(date_counts_df.iloc[1]['シェア_%']):.1f}%) で約 1/3 規模。" "Top 8 のうち年度末 (R2.3.31, R3.3.31) と半期末 (R2.9.30, R2.5.29) と R3.5.31 経過措置期限が並ぶ = " "「節目 + 期限の集中処理」という行政パターンが時系列に刻まれた。

" + "

表: 指定年別統計 (R2-R5)

" + df_to_html(T_year) + f"

この表から読み取れること: 各指定年の件数・規模統計。" f"R3 が突出 ({len(m_des[m_des['年'] == 2021]):,} 件 = {len(m_des[m_des['年'] == 2021])/n_total*100:.0f}%)、" f"R2 ({len(m_des[m_des['年'] == 2020]):,} 件) → R3 急増 → R4-R5 終息という3 段階構造。" f"中央規模を見るとR2 中央 {r2_size:.2f} 千m³ → R5 {r5_size:.2f}と単調減少 = " "大規模池ほど早く指定された制度運用が定量化された。

" + "

表: 上位 12 市町 × 指定年 (件数)

" + df_to_html(T_year_city) + f"

この表から読み取れること: 上位 12 市町の指定年別件数。R3 列がほぼ全市町で最大値を占める = " "経過措置期限の一斉処理が市町横断で発生。" "ただし東広島は R2 から段階的に処理、福山は R3 集中のように市町別パターン差がある。" "市町の指定能力 (担当者数・事務処理速度) の違いを反映。

" ), (f"【RQ3】 L45 浸水想定との対応 — 整備率 {inund_rate:.2f}% / 未整備 {n_no_inund} 件 中山間 100% 集中", sec6_intro + "

実装コード (L45 cache 読込 + ため池番号 join + 市町別整備率)

" + sec6_code + sec6_fig7_text + figure("assets/L59_fig7_inund_status_map.png", f"図 7 (RQ3): L45 浸水想定整備済 vs 未整備 {n_no_inund} 件全数マップ + 三次・庄原ズーム") + sec6_fig7_read + sec6_fig8_text + figure("assets/L59_fig8_inund_rate_size.png", "図 8 (RQ3): 市町別整備率 Top 15 + 整備済 vs 未整備 規模箱ひげ") + sec6_fig8_read + "

表: 市町別 浸水想定整備率 Top 15

" + df_to_html(T_inund_city) + f"

この表から読み取れること: 各市町の全件数・整備済・未整備・整備率。" f"三次市と庄原市のみ整備率 < 100% (各 -12 件)、他 21 市町は完全整備。" f"全県整備率 {inund_rate:.2f}% は実質「2 市町問題」 で、整備未了は制度的に説明可能な少数例外。" "未整備 24 件は下流被害が想定されない小規模・中山間池として浸水想定整備の対象外になっている。

" + f"

表: 浸水想定 SHP 未整備 {n_no_inund} 件 全数詳細

" + df_to_html(T_no_inund) + f"

この表から読み取れること: 未整備 24 件の全数詳細 — ため池番号・名称・市町・堤高・貯水量・指定日・緯度経度。" f"三次市 12 件 + 庄原市 12 件 = 24 件すべて中山間。" f"指定済 {n_no_inund_des} 件 / 未指定 {n_no_inund_no_des} 件で「指定はされているが浸水想定 SHP は未作成」 が大半 = " "行政の段階的整備の途上にある池群。" "緯度経度から個別池の位置を特定可能で、優先度スコアリング (発展課題 4) の出発点。

" ), ("仮説検証総合", sec7), ("発展課題", sec8), ] html = render_lesson( num=59, title="ため池基本情報 単独 3 研究例分析 — CSV 6,754 件から農業遺構の地理 × 制度史 × 整備状況を読む", tags=["L59", "ため池", "農業遺構", "ため池管理保全法", "RQ×3", "Format B", "CSV 解析", "規模分布", "log-log", "経過措置期限", "L45連携", "中山間集中", "瀬戸内型農業"], time="40 分", level="中級", data_label=f"DoBoX dataset 62 (CSV {n_total:,} 行 × 12 列, 1.1 MB) + L45 連携", sections=sections, script_filename="L59_pond_basic.py", ) OUT_HTML = LESSONS / "L59_pond_basic.html" OUT_HTML.write_text(html, encoding="utf-8") print(f" HTML: {OUT_HTML.name} ({len(html):,} chars)") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 終了 # ============================================================================= print(f"\n=== 完了: 全体 {time.time() - t_all:.1f}s ===") print(f"出力:") print(f" HTML: lessons/L59_pond_basic.html") print(f" Python: lessons/L59_pond_basic.py") print(f" 図: assets/L59_fig1-8_*.png (8 枚)") print(f" CSV: assets/L59_*.csv (11 表)") print(f" ZIP: assets/L59_pond_basic_raw.zip")