# -*- coding: utf-8 -*- """L61 過去に発生した災害情報 単独 3 研究例分析 — 広島県内で過去に発生した災害 518 件 (CSV) を 3 つの独立した研究角度で読み解く カバー宣言: 本記事は DoBoX のシリーズ「過去に発生した災害情報」 1 件 (dataset_id = 1278) を 単独で取り上げ、 広島県内で過去に発生した災害の現地記録 518 件を 3 つの独立した研究角度 (RQ1 / RQ2 / RQ3) で並列に分析する。 Phase 5 災害履歴系の起点 = 「過去災害」 の現地写真・石碑・記録が 集約された災害アーカイブを最初に正面から扱う。 「過去に発生した災害情報」 とは: 広島県土木建築局・砂防課が運営する「Web 砂防 (土砂災害ポータル広島)」 (https://www.sabo.pref.hiroshima.lg.jp/portal/) の 地域情報レイヤー。県内で過去に実際に発生した災害について、 現地の航空写真・石碑・記録集・体験談・歴史的土木施設を 位置情報つきで集めたアーカイブ。 本データは「地域情報ID・タイトル・コメント・所在地市区町・タグ・緯度経度・InfoURL」 の 8 列で 518 行 (約 632 KB) のシンプル構造。 各行は 1 つの災害現地点に対応し、コメントは数百字におよぶ詳細記述、 タグはイベント名 (例「平成30年7月豪雨」) + 場所 + 種別を `,` で連結した文字列。 既存ハザード系 (L10/L11) との関係: L10 (土砂災害クロス) と L11 (トリプルハザード) は未来の予測 (= 警戒区域指定) を 扱った。本記事 L61 は過去の発生事実 (= 実際に災害が起きた地点) を扱う。 これは「予測 vs 実績」 の対比研究で、警戒区域指定が過去災害をどこまで 予測できていたか (RQ3) を直接検証できる初めての教材。 研究的には「警戒区域は最強の検証データ = 過去災害」 によって妥当性が問われる。 本データ #1278 の独自性と限界: - 歴史記録は乏しい (1907 年の事例も含む): 古い記録は石碑として残る。 本データには 52 件の「石碑」 タグ付き記録があり、災害伝承の物理形態を可視化。 対して新しい事例は航空写真 (65 件) で記録される。 - 大災害イベントが大量出現する偏在型分布: 平成30年7月豪雨 158 件、 平成26年8月豪雨 79 件の 2 大イベントだけで全体の 46%。 残り 54% を 100 年分の中小イベントで分け合う極端な不均衡。 - 位置情報は 100% 完備: 緯度経度 NaN ゼロ = 全 518 件マップ可能。 緯度範囲 34.07-34.93 (= 県南北端を完全カバー)。 研究の問い (3 RQ): RQ1 (主研究): 広島県の災害履歴 — 種別・年代・地理分布はどう描けるか? 518 件を年代 (大正/昭和/平成/令和)、月、市町、タグ (イベント名) で多角度集計。 大災害イベントの寄与率と歴史的記録の地理分布を立体的に描く。 RQ2 (副研究 1): 2018 年西日本豪雨 vs 2014 年広島市土砂災害 — 二大災害の比較 両者は近年の主要災害だが、地理スケールが異なる: - 2014: 広島市集中型 (短時間局所豪雨で安佐南区・安佐北区中心) - 2018: 県全域型 (梅雨前線停滞で全 22 市町に大雨特別警報) 両者の地理分布・近接度・市町別パターンを統計的に比較し、 「同じ豪雨でも発生機構が違えば被害分布が違う」 を実証する。 RQ3 (副研究 2): 警戒区域 (L10/L11 既知) との空間整合 — 過去災害は予測されていたか 土砂災害警戒区域 polygon (L11 既扱) と過去災害点を 2 段階で評価: (1) sjoin (point in polygon) で厳密内包率 (≈10%) を測定 (2) sjoin_nearest で最近警戒域距離を計算し、500m / 1km バッファでの予測圏率を測定 過去災害点には記念碑・撮影位置等の非崩落地点が含まれるため、 バッファ評価が現実的な制度予測精度を示す。 仮説 (5): H1 (二大災害集中, RQ1): 全 518 件のうち平成30年7月豪雨 + 平成26年8月豪雨 ≥ 40%。 近年の 2 大イベントが歴史記録の大半を占める偏在型分布。 H2 (季節性, RQ1): 月別では6-9 月 ≥ 90%。日本の梅雨〜台風シーズンに集中。 H3 (2014 局所性 vs 2018 広域性, RQ2): 2014 年災害の地理重心 (緯度経度) の分散は 2018 年災害より明確に小さい (同じ円内に密集)。 2014 は局所豪雨、2018 は広域豪雨という発生機構の差を空間統計で実証。 H4 (警戒域近接率, RQ3): 過去災害の500m 以内に土砂災害警戒区域がある (= 制度的予測圏内) 割合が 50% 以上。点厳密内包は写真撮影位置のずれで低くなりがちなので、 近接バッファで「制度予測との空間的近さ」 を測る。 H5 (急傾斜型偏在, RQ3): 過去災害のうち急傾斜警戒区域内に位置する点が 地すべり警戒区域内より多い。広島県の地形特性 (急峻な山地斜面) を反映。 要件 S 準拠 (1分以内完走): - CSV 518 行は読込・パース < 0.5 秒 (修復ロジックは csv モジュールで) - L11 既処理の警戒区域 polygon を再利用 (重い読込は 1 度だけ) - sjoin (point in polygon) は 518 点 × 数万 polygon = 数秒 - 期待実行時間: 30-60 秒 要件 T 準拠 (位置情報あり=地図必須): - RQ1: 全 518 点マップ (色=年代) + 全 518 点マップ (色=タグ主要 5 イベント) - RQ2: 2014 vs 2018 の 2 ペインマップ + 重心円の重ね合わせ - RQ3: 警戒域 polygon + 過去災害点 重ね合わせマップ - RQ3: 警戒域内 vs 警戒域外 のクロス積み上げ棒 要件 Q 準拠: 図 8 / 表 12 (3 RQ × 多角度)。 データ仕様: - dataset 1278: 過去に発生した災害情報 (CSV) - リソース (CSV, 約 632 KB, UTF-8 BOM): 518 行 × 8 列 - 列構成 (8 列): 地域情報ID (整数 1-339, 飛び番号あり) / タイトル (例「平成17年9月6日 台風14号での宮島の被災状況」) / コメント (数十〜数千字の詳細記述) / 所在地市区町 (例「廿日市市」) / タグ (`,` 区切り 文字列, 例「平成17年台風14号,宮島町（廿日市市）,白糸川」) / 緯度 / 経度 (WGS84 度) / InfoURL (Web 砂防ポータルの該当 ID ページ) - CRS: WGS84 → 解析時 EPSG:6671 (m 単位) へ変換 - 取得日: 2026-05-10 - ライセンス: クリエイティブ・コモンズ表示 4.0 - 作成主体: 広島県土木建築局 - 注意: 元 CSV はコメント内に未エスケープのカンマを含む行が 複数あり (= 平成30年7月豪雨の長文コメント)、pandas.read_csv では ParserError が発生。本スクリプトはcsv モジュール + 末尾固定 5 列の 再結合ロジックでこれを修復し、518 行全件を読み込む。 """ from __future__ import annotations import sys, time, csv, re, zipfile, json from pathlib import Path from collections import Counter 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, Circle from matplotlib.lines import Line2D plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False t_all = time.time() print("=== L61 過去に発生した災害情報 単独 3 研究例分析 ===", flush=True) # ============================================================================= # 0. 定数・パス # ============================================================================= TARGET_CRS = "EPSG:6671" SRC_CRS = "EPSG:4326" EARTH_R_M = 6_371_000.0 DATA_DIR = ROOT / "data" / "extras" / "L61_past_disasters" DATA_DIR.mkdir(parents=True, exist_ok=True) CSV_LOCAL = DATA_DIR / "past_disasters.csv" CSV_GLOBAL = ROOT / "data" / "extras" / "past_disasters.csv" # 警戒区域 Shapefile (L11 既扱) SED_DIR = ROOT / "data" / "extras" / "sediment_shp" DOSEKI_SHP = SED_DIR / "doseki" / "340006_sediment_disaster_hazard_area_debris_flow_20260427" / "73_031drp_34_20260427.shp" KYUKEI_SHP = SED_DIR / "kyukeisha" / "340006_sediment_disaster_hazard_area_steep_slope_20260427" / "73_031krp_34_20260427.shp" JISUBERI_SHP = SED_DIR / "jisuberi" / "340006_sediment_disaster_hazard_area_landslide_20260427" / "73_031jy_34_20260427.shp" # 行政区域 (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)"} # 元号→西暦変換 ERA_BASE = {"令和": 2018, "平成": 1988, "昭和": 1925, "大正": 1911, "明治": 1867} # ============================================================================= # 1. データ取得 # ============================================================================= print("\n[1] データ取得", flush=True) t1 = time.time() if not CSV_GLOBAL.exists() or CSV_GLOBAL.stat().st_size < 1000: url = f"{DOBOX_BASE}/datasets/1278" r = requests.get(url, headers=HDR, timeout=60) rids = re.findall(r"/resources/(\d+)", r.text) if not rids: raise RuntimeError("dataset 1278 にリソースが見つかりません") rid = rids[0] rh = requests.get(f"{DOBOX_BASE}/resources/{rid}", headers=HDR, timeout=60) m = re.search(r'href="(' + re.escape(DOBOX_BASE) + r'/resource_download/\d+)"', rh.text) if not m: raise RuntimeError(f"resource {rid} の resource_download が解決できません") dl = m.group(1) print(f" DL <- {dl}") rd = requests.get(dl, headers=HDR, timeout=120) rd.raise_for_status() CSV_GLOBAL.parent.mkdir(parents=True, exist_ok=True) CSV_GLOBAL.write_bytes(rd.content) print(f" saved -> {CSV_GLOBAL.name} ({len(rd.content):,} bytes)") else: print(f" cache HIT: {CSV_GLOBAL.name} ({CSV_GLOBAL.stat().st_size:,} bytes)") # L61 ローカルにコピー if not CSV_LOCAL.exists(): CSV_LOCAL.write_bytes(CSV_GLOBAL.read_bytes()) print(f" ({time.time()-t1:.2f}s)", flush=True) # ============================================================================= # 2. CSV 読込・前処理 (修復ロジック) # ============================================================================= print("\n[2] CSV 読込・前処理 (修復ロジック)", flush=True) t1 = time.time() # 元 CSV はコメント内に未エスケープのカンマがある行が複数 (平成30年7月豪雨の長文) # csv モジュールで読み、列数 ≠ 8 の行は「先頭2列 + 中間=コメント結合 + 末尾5列」で再結合 with open(CSV_GLOBAL, "r", encoding="utf-8-sig") as f: rows = list(csv.reader(f)) header = rows[0] data = [] n_repaired = 0 for r in rows[1:]: if len(r) == 8: data.append(r) elif len(r) > 8: data.append(r[:2] + [",".join(r[2:-5])] + r[-5:]) n_repaired += 1 # < 8 列は破棄 (実際にはなし) df = pd.DataFrame(data, columns=header) n_total = len(df) print(f" rows: {n_total} (修復: {n_repaired} 行)") print(f" columns: {df.columns.tolist()}") # 緯度経度 df["緯度"] = pd.to_numeric(df["緯度"], errors="coerce") df["経度"] = pd.to_numeric(df["経度"], errors="coerce") n_xy = (df["緯度"].notna() & df["経度"].notna()).sum() print(f" 緯度経度 ok: {n_xy} / {n_total}") print(f" 緯度範囲: {df['緯度'].min():.4f} - {df['緯度'].max():.4f}") print(f" 経度範囲: {df['経度'].min():.4f} - {df['経度'].max():.4f}") # 年抽出 (タイトルから 元号 + 年) def extract_year(t): if not t: return None m = re.search(r"(令和|平成|昭和|大正|明治)([0-9]+|元)年", str(t)) if not m: return None era, num = m.group(1), m.group(2) n = 1 if num == "元" else int(num) return ERA_BASE[era] + n df["年"] = df["タイトル"].apply(extract_year) n_year_ok = df["年"].notna().sum() print(f" 年抽出 ok: {n_year_ok} / {n_total}") print(f" 年レンジ: {int(df['年'].min())} - {int(df['年'].max())}") # 月抽出 (タイトルから 数字+月) def extract_month(t): if not t: return None m = re.search(r"(\d+)月", str(t)) return int(m.group(1)) if m else None df["月"] = df["タイトル"].apply(extract_month) # 元号 (年代区分) def era_of(y): if y is None or pd.isna(y): return "不明" y = int(y) if y >= 2019: return "令和" if y >= 1989: return "平成" if y >= 1926: return "昭和" if y >= 1912: return "大正" return "明治以前" df["年代区分"] = df["年"].apply(era_of) print(f" 年代区分: {df['年代区分'].value_counts().to_dict()}") # タグ展開 (1 行で複数タグを保持) def split_tags(s): if not s: return [] return [t.strip().strip('"').strip() for t in str(s).split(",") if t.strip().strip('"').strip()] df["タグリスト"] = df["タグ"].apply(split_tags) # 主要イベントタグ (本記事で議論する大災害) EVENT_TAGS = { "平成30年7月豪雨": "2018西日本豪雨", "平成26年8月豪雨": "2014広島市土砂", "昭和42年7月豪雨": "1967昭42豪雨", "平成11年6月豪雨": "1999平11豪雨", "昭和63年7月豪雨": "1988昭63豪雨", "昭和20年枕崎台風": "1945枕崎台風", "平成22年7月豪雨": "2010平22豪雨", "平成17年台風14号": "2005台風14号", } def primary_event(tags): for t in tags: if t in EVENT_TAGS: return EVENT_TAGS[t] return "その他" df["主要イベント"] = df["タグリスト"].apply(primary_event) # 種別タグ (記録種別) TYPE_TAGS = {"航空写真", "石碑", "記録集", "歴史的土木施設", "体験談"} def record_types(tags): return [t for t in tags if t in TYPE_TAGS] df["記録種別リスト"] = df["タグリスト"].apply(record_types) df["記録種別主"] = df["記録種別リスト"].apply(lambda lst: lst[0] if lst else "現地記録") # GeoDataFrame 化 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" ({time.time()-t1:.2f}s)", flush=True) # ============================================================================= # 3. 行政区域 (県全域) 読込 # ============================================================================= print("\n[3] 行政区域 (県全域) 読込", flush=True) t1 = time.time() pref_dir = DATA_DIR / "_admin_pref" pref_dir.mkdir(parents=True, exist_ok=True) pref_outline = None pref_cities = None if ADMIN_PREF_ZIP.exists(): if not list(pref_dir.glob("*.geojson")): try: with zipfile.ZipFile(ADMIN_PREF_ZIP) as zf: zf.extractall(pref_dir) except Exception as e: print(f" WARN: pref zip extract failed: {e}") geos = list(pref_dir.glob("*.geojson")) if geos: try: pref_cities = gpd.read_file(geos[0]).to_crs(TARGET_CRS) pref_outline = pref_cities.dissolve() print(f" pref polygons: {len(pref_cities)} → outline OK") except Exception as e: print(f" WARN: pref read failed: {e}") print(f" ({time.time()-t1:.2f}s)", flush=True) # ============================================================================= # 4. RQ1 分析: 災害履歴の構造 (種別・年代・地理) # ============================================================================= print("\n[4] RQ1: 災害履歴の構造研究", flush=True) t1 = time.time() # (a) 年代区分別件数 era_count = df["年代区分"].value_counts().reindex( ["明治以前", "大正", "昭和", "平成", "令和", "不明"] ).fillna(0).astype(int) era_count_df = era_count.reset_index().rename(columns={"index": "年代区分", "count": "件数"}) era_count_df.columns = ["年代区分", "件数"] era_count_df["シェア_%"] = (era_count_df["件数"] / n_total * 100).round(2) era_count_df.to_csv(ASSETS / "L61_era_count.csv", index=False, encoding="utf-8-sig") # (b) 主要イベント別ランキング event_count = df["主要イベント"].value_counts() event_df = event_count.reset_index() event_df.columns = ["主要イベント", "件数"] event_df["シェア_%"] = (event_df["件数"] / n_total * 100).round(2) event_df.to_csv(ASSETS / "L61_event_ranking.csv", index=False, encoding="utf-8-sig") n_top2 = event_df[event_df["主要イベント"].isin(["2018西日本豪雨", "2014広島市土砂"])]["件数"].sum() top2_share = n_top2 / n_total * 100 print(f" 主要イベント Top 2 (2018+2014): {n_top2} ({top2_share:.1f}%)") # (c) 月別分布 month_count = df["月"].value_counts().reindex(range(1, 13)).fillna(0).astype(int) month_df = month_count.reset_index() month_df.columns = ["月", "件数"] month_df["シェア_%"] = (month_df["件数"] / n_total * 100).round(2) month_df.to_csv(ASSETS / "L61_month_dist.csv", index=False, encoding="utf-8-sig") n_jun_sep = int(month_df[month_df["月"].between(6, 9)]["件数"].sum()) share_jun_sep = n_jun_sep / n_total * 100 print(f" 6-9 月集中: {n_jun_sep} / {n_total} = {share_jun_sep:.1f}%") # (d) 市町別ランキング city_count = df["所在地市区町"].value_counts() city_df = city_count.reset_index() city_df.columns = ["市町", "件数"] city_df["シェア_%"] = (city_df["件数"] / n_total * 100).round(2) city_df.to_csv(ASSETS / "L61_city_ranking.csv", index=False, encoding="utf-8-sig") top_city = city_df.iloc[0] print(f" Top 市町: {top_city['市町']} = {top_city['件数']} ({top_city['シェア_%']:.1f}%)") # (e) 記録種別ランキング rec_count = df["記録種別主"].value_counts() rec_df = rec_count.reset_index() rec_df.columns = ["記録種別", "件数"] rec_df["シェア_%"] = (rec_df["件数"] / n_total * 100).round(2) rec_df.to_csv(ASSETS / "L61_record_type.csv", index=False, encoding="utf-8-sig") # (f) 年別ヒスト集計 (10年ビン) df_year_ok = df[df["年"].notna()].copy() df_year_ok["年"] = df_year_ok["年"].astype(int) ymin, ymax = int(df_year_ok["年"].min()), int(df_year_ok["年"].max()) print(f" 年範囲: {ymin} - {ymax}") print(f" ({time.time()-t1:.2f}s)", flush=True) # ============================================================================= # 5. RQ2 分析: 2018 vs 2014 二大災害の比較 # ============================================================================= print("\n[5] RQ2: 二大災害比較研究", flush=True) t1 = time.time() E_2018 = df[df["主要イベント"] == "2018西日本豪雨"].copy() E_2014 = df[df["主要イベント"] == "2014広島市土砂"].copy() n_2018 = len(E_2018) n_2014 = len(E_2014) print(f" 2018 西日本豪雨: {n_2018} 件") print(f" 2014 広島市土砂: {n_2014} 件") # 2014 vs 2018 の地理重心 def centroid_radius(g): """重心 (緯度経度) と平均距離 (m, haversine) を返す""" if len(g) == 0: return (None, None, None, None) lat0 = g["緯度"].mean() lon0 = g["経度"].mean() # haversine 距離 (m) lat = np.deg2rad(g["緯度"].values) lon = np.deg2rad(g["経度"].values) lat_c = np.deg2rad(lat0) lon_c = np.deg2rad(lon0) dlat = lat - lat_c dlon = lon - lon_c a = np.sin(dlat/2)**2 + np.cos(lat_c) * np.cos(lat) * np.sin(dlon/2)**2 d_m = 2 * EARTH_R_M * np.arcsin(np.sqrt(a)) return (lat0, lon0, float(np.mean(d_m)), float(np.std(d_m))) c14_lat, c14_lon, r14_mean, r14_std = centroid_radius(E_2014) c18_lat, c18_lon, r18_mean, r18_std = centroid_radius(E_2018) print(f" 2014 重心: ({c14_lat:.4f}, {c14_lon:.4f}), 平均距離 {r14_mean:.0f}m (std {r14_std:.0f})") print(f" 2018 重心: ({c18_lat:.4f}, {c18_lon:.4f}), 平均距離 {r18_mean:.0f}m (std {r18_std:.0f})") # 重心間距離 def haversine_m(lat1, lon1, lat2, lon2): p1 = np.deg2rad([lat1, lon1]) p2 = np.deg2rad([lat2, lon2]) dlat = p2[0] - p1[0] dlon = p2[1] - p1[1] a = np.sin(dlat/2)**2 + np.cos(p1[0]) * np.cos(p2[0]) * np.sin(dlon/2)**2 return 2 * EARTH_R_M * np.arcsin(np.sqrt(a)) c14_to_c18_m = haversine_m(c14_lat, c14_lon, c18_lat, c18_lon) print(f" 2014↔2018 重心間距離: {c14_to_c18_m:.0f}m = {c14_to_c18_m/1000:.1f}km") # 市町別 2014 vs 2018 のクロス city_2014 = E_2014["所在地市区町"].value_counts().rename("件数_2014") city_2018 = E_2018["所在地市区町"].value_counts().rename("件数_2018") city_compare = pd.concat([city_2014, city_2018], axis=1).fillna(0).astype(int) city_compare["合計"] = city_compare.sum(axis=1) city_compare = city_compare.sort_values("合計", ascending=False).reset_index().rename(columns={"index": "市町"}) city_compare.columns = ["市町", "件数_2014", "件数_2018", "合計"] city_compare.to_csv(ASSETS / "L61_city_2014vs2018.csv", index=False, encoding="utf-8-sig") n_city_both = ((city_compare["件数_2014"] > 0) & (city_compare["件数_2018"] > 0)).sum() n_city_only14 = ((city_compare["件数_2014"] > 0) & (city_compare["件数_2018"] == 0)).sum() n_city_only18 = ((city_compare["件数_2014"] == 0) & (city_compare["件数_2018"] > 0)).sum() print(f" 両イベントとも被災市町: {n_city_both} 件") print(f" 2014 のみ: {n_city_only14}, 2018 のみ: {n_city_only18}") # 分散比 (R² 比) — 2018 が広域なら std が大きいはず ratio_std = r18_std / r14_std if r14_std > 0 else float("nan") ratio_mean = r18_mean / r14_mean if r14_mean > 0 else float("nan") print(f" 分散比 std(2018)/std(2014) = {ratio_std:.2f}") print(f" 平均距離比 mean(2018)/mean(2014) = {ratio_mean:.2f}") print(f" ({time.time()-t1:.2f}s)", flush=True) # ============================================================================= # 6. RQ3 分析: 警戒区域との空間整合 # ============================================================================= print("\n[6] RQ3: 警戒区域との空間整合研究", flush=True) t1 = time.time() # Shapefile 読込 (L11 既扱) print(" L11 既扱の警戒区域 polygon を読込") doseki = gpd.read_file(DOSEKI_SHP).to_crs(TARGET_CRS) kyukei = gpd.read_file(KYUKEI_SHP).to_crs(TARGET_CRS) jisuberi = gpd.read_file(JISUBERI_SHP).to_crs(TARGET_CRS) print(f" 土石流: {len(doseki)} polygon") print(f" 急傾斜: {len(kyukei)} polygon") print(f" 地すべり: {len(jisuberi)} polygon") # 各 polygon に種別タグを付与 doseki["災害種別"] = "土石流" kyukei["災害種別"] = "急傾斜" jisuberi["災害種別"] = "地すべり" # sjoin (point in polygon) — 過去災害点 × 警戒区域 disaster_pts = gdf[["地域情報ID", "タイトル", "所在地市区町", "主要イベント", "年", "x_m", "y_m", "geometry"]].copy() # 各種別との交差判定 (sindex で高速) sj_doseki = gpd.sjoin(disaster_pts, doseki[["geometry"]], how="left", predicate="within") sj_kyukei = gpd.sjoin(disaster_pts, kyukei[["geometry"]], how="left", predicate="within") sj_jisuberi = gpd.sjoin(disaster_pts, jisuberi[["geometry"]], how="left", predicate="within") # index_right が NaN でなければ警戒域内 in_doseki = sj_doseki.groupby(level=0)["index_right"].apply(lambda s: s.notna().any()).rename("土石流警戒域内") in_kyukei = sj_kyukei.groupby(level=0)["index_right"].apply(lambda s: s.notna().any()).rename("急傾斜警戒域内") in_jisuberi = sj_jisuberi.groupby(level=0)["index_right"].apply(lambda s: s.notna().any()).rename("地すべり警戒域内") disaster_pts = disaster_pts.join(in_doseki).join(in_kyukei).join(in_jisuberi) disaster_pts["警戒域内"] = (disaster_pts["土石流警戒域内"] | disaster_pts["急傾斜警戒域内"] | disaster_pts["地すべり警戒域内"]) n_in_doseki = int(disaster_pts["土石流警戒域内"].sum()) n_in_kyukei = int(disaster_pts["急傾斜警戒域内"].sum()) n_in_jisuberi = int(disaster_pts["地すべり警戒域内"].sum()) n_in_any = int(disaster_pts["警戒域内"].sum()) n_out = n_total - n_in_any share_in_any = n_in_any / n_total * 100 print(f" 土石流警戒域内: {n_in_doseki} ({n_in_doseki/n_total*100:.1f}%)") print(f" 急傾斜警戒域内: {n_in_kyukei} ({n_in_kyukei/n_total*100:.1f}%)") print(f" 地すべり警戒域内: {n_in_jisuberi} ({n_in_jisuberi/n_total*100:.1f}%)") print(f" いずれかの警戒域内: {n_in_any} ({share_in_any:.1f}%)") print(f" 警戒域外で発生: {n_out} ({n_out/n_total*100:.1f}%)") # 補助指標: 警戒区域からの距離分布 (= 「近接した災害」 の指標) # 注: 過去災害点には航空写真等の撮影地点も含まれ、polygon 直撃でなくても # 警戒区域に「近接」 していれば制度的に予測されていたと言える。 # 全警戒 polygon の和集合外殻からの最近接距離を計算 print(" 警戒区域からの距離分布を計算 (制度予測との近接度)") all_warn = pd.concat([ doseki[["geometry"]], kyukei[["geometry"]], jisuberi[["geometry"]] ], ignore_index=True) all_warn = gpd.GeoDataFrame(all_warn, geometry="geometry", crs=TARGET_CRS) # sjoin_nearest で各点から最近警戒polygonまでの距離 (m) nearest_warn = gpd.sjoin_nearest(disaster_pts[["geometry"]], all_warn[["geometry"]], how="left", distance_col="warn_dist_m", max_distance=20000) # 上限 20km # 重複対応 (sjoin_nearest は 1 点が複数 polygon に同距離マッチすることがある) # index で min を取る warn_dist = nearest_warn.groupby(level=0)["warn_dist_m"].min() disaster_pts["警戒域距離_m"] = warn_dist.reindex(disaster_pts.index).values # しきい値別カバー率 for thr in [0, 100, 500, 1000]: n_within = int((disaster_pts["警戒域距離_m"] <= thr).sum()) print(f" 警戒域 {thr}m 以内発生: {n_within} ({n_within/n_total*100:.1f}%)") n_within_100 = int((disaster_pts["警戒域距離_m"] <= 100).sum()) n_within_500 = int((disaster_pts["警戒域距離_m"] <= 500).sum()) share_within_100 = n_within_100 / n_total * 100 share_within_500 = n_within_500 / n_total * 100 warn_dist_med = float(disaster_pts["警戒域距離_m"].median()) warn_dist_mean = float(disaster_pts["警戒域距離_m"].mean()) # 主要イベント別の警戒域内シェア event_in_share = disaster_pts.groupby( df.loc[disaster_pts.index, "主要イベント"] )["警戒域内"].agg(["sum", "count"]).reset_index() event_in_share.columns = ["主要イベント", "警戒域内件数", "総件数"] event_in_share["警戒域内_%"] = (event_in_share["警戒域内件数"] / event_in_share["総件数"] * 100).round(1) event_in_share = event_in_share.sort_values("総件数", ascending=False).reset_index(drop=True) event_in_share.to_csv(ASSETS / "L61_event_in_warning.csv", index=False, encoding="utf-8-sig") # 市町別 警戒域内/外 city_in_share = disaster_pts.assign(市町=df.loc[disaster_pts.index, "所在地市区町"]).groupby( "市町" )["警戒域内"].agg(["sum", "count"]).reset_index() city_in_share.columns = ["市町", "警戒域内件数", "総件数"] city_in_share["警戒域内_%"] = (city_in_share["警戒域内件数"] / city_in_share["総件数"] * 100).round(1) city_in_share = city_in_share.sort_values("総件数", ascending=False).reset_index(drop=True) city_in_share.to_csv(ASSETS / "L61_city_in_warning.csv", index=False, encoding="utf-8-sig") # 警戒域外の代表例 out_ex = disaster_pts[~disaster_pts["警戒域内"]].copy() out_ex_df = df.loc[out_ex.index, ["地域情報ID", "タイトル", "所在地市区町", "主要イベント", "年"]].copy() out_ex_df["緯度"] = df.loc[out_ex.index, "緯度"].round(5) out_ex_df["経度"] = df.loc[out_ex.index, "経度"].round(5) out_ex_df = out_ex_df.sort_values("年", na_position="last").reset_index(drop=True) out_ex_df.to_csv(ASSETS / "L61_outside_warning.csv", index=False, encoding="utf-8-sig") # 警戒域内 × 主要イベント のクロス event_warn_cross = pd.crosstab( df["主要イベント"], disaster_pts["警戒域内"].reindex(df.index).fillna(False).map({True: "警戒域内", False: "警戒域外"}) ) event_warn_cross = event_warn_cross.reindex( sorted(event_warn_cross.index, key=lambda k: -event_warn_cross.loc[k].sum()) ) event_warn_cross.to_csv(ASSETS / "L61_event_warn_cross.csv", encoding="utf-8-sig") print(f" ({time.time()-t1:.2f}s)", flush=True) # ============================================================================= # 7. 図描画 # ============================================================================= print("\n[7] 図描画", flush=True) t1 = time.time() def draw_pref(ax, lw=0.5, color="#888"): if pref_outline is not None: pref_outline.boundary.plot(ax=ax, color=color, linewidth=lw) if pref_cities is not None: pref_cities.boundary.plot(ax=ax, color=color, linewidth=lw * 0.4, alpha=0.45) # ---- 図 1 (RQ1): 全 518 点 + 年代色分けマップ ------------------------------ print(" fig1 全県マップ (年代色分け)", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 7)) ax = axes[0] draw_pref(ax, lw=0.7, color="#777") era_colors = {"明治以前": "#5c4033", "大正": "#888888", "昭和": "#4575b4", "平成": "#fdae61", "令和": "#d73027", "不明": "#cccccc"} for era_name, color in era_colors.items(): sub = gdf[gdf["年代区分"] == era_name] if len(sub) > 0: ax.scatter(sub["x_m"], sub["y_m"], c=color, s=18, alpha=0.75, edgecolors="black", linewidths=0.2, label=f"{era_name} ({len(sub)})") ax.set_title(f"広島県内 過去災害 {n_total} 件 (色 = 年代区分)", fontsize=11) ax.legend(fontsize=8, loc="lower right") ax.set_aspect("equal") ax.set_xticks([]); ax.set_yticks([]) ax = axes[1] # 年別ヒスト (10 年ビン) years = df_year_ok["年"].astype(int).values bins = np.arange((years.min() // 10) * 10, ((years.max() // 10) + 2) * 10, 10) ax.hist(years, bins=bins, color="#4575b4", edgecolor="white") ax.set_xlabel("発生年 (10 年ビン)") ax.set_ylabel("災害記録数") # 大災害イベントの位置を矢印 event_arrows = [ (1945, "枕崎台風"), (1967, "S42 豪雨"), (1988, "S63 豪雨"), (1999, "H11 豪雨"), (2014, "H26 広島土砂"), (2018, "H30 西日本豪雨"), ] ymax_h = max(np.histogram(years, bins=bins)[0]) * 1.05 for y, lbl in event_arrows: ax.axvline(y, color="red", linestyle=":", linewidth=0.8, alpha=0.7) ax.text(y, ymax_h * 0.95, lbl, rotation=90, fontsize=8, ha="right", va="top", color="darkred") ax.set_title(f"発生年分布 ({ymin}-{ymax}, 10 年ビン)\n" f"赤点線 = 大災害イベント", fontsize=10) ax.grid(alpha=0.3) plt.suptitle( f"L61 RQ1 図 1: 広島県 過去災害 {n_total} 件 — 年代マップ + 年別ヒスト", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L61_fig1_year_map.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 2 (RQ1): 主要イベント色分け + 月別棒 ------------------------------ print(" fig2 主要イベントマップ + 月別棒", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 7)) ax = axes[0] draw_pref(ax, lw=0.7, color="#777") event_colors = { "2018西日本豪雨": "#d73027", "2014広島市土砂": "#fc8d59", "1999平11豪雨": "#fee08b", "1967昭42豪雨": "#91cf60", "1988昭63豪雨": "#1a9850", "1945枕崎台風": "#878787", "2010平22豪雨": "#542788", "2005台風14号": "#998ec3", "その他": "#cccccc", } # 描画順 (背景に「その他」、前景に大イベント) draw_order = ["その他", "2005台風14号", "2010平22豪雨", "1945枕崎台風", "1988昭63豪雨", "1967昭42豪雨", "1999平11豪雨", "2014広島市土砂", "2018西日本豪雨"] for ev in draw_order: sub = gdf[gdf["主要イベント"] == ev] if len(sub) > 0: ax.scatter(sub["x_m"], sub["y_m"], c=event_colors[ev], s=22, alpha=0.85, edgecolors="black", linewidths=0.25, label=f"{ev} ({len(sub)})") ax.set_title(f"主要イベント色分け — 全 {n_total} 件", fontsize=11) ax.legend(fontsize=7, loc="lower right", ncol=1) ax.set_aspect("equal") ax.set_xticks([]); ax.set_yticks([]) ax = axes[1] # 月別棒 mp = month_df[month_df["月"].notna()].copy() mp["月"] = mp["月"].astype(int) ax.bar(mp["月"], mp["件数"], color="#4575b4", edgecolor="white", alpha=0.85) # 6-9 月帯 ax.axvspan(5.5, 9.5, color="orange", alpha=0.13, label=f"6-9 月帯 ({n_jun_sep} 件 = {share_jun_sep:.0f}%)") ax.set_xticks(range(1, 13)) ax.set_xlabel("月") ax.set_ylabel("災害記録数") ax.set_title(f"発生月分布 — 6-9 月に {share_jun_sep:.0f}% が集中", fontsize=10) ax.legend(fontsize=9) ax.grid(axis="y", alpha=0.3) for _, r in mp.iterrows(): if r["件数"] > 0: ax.text(r["月"], r["件数"] + 3, f"{int(r['件数'])}", ha="center", fontsize=8) plt.suptitle( f"L61 RQ1 図 2: 主要イベント色分けマップ + 月別棒 — 季節性", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L61_fig2_event_month.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 3 (RQ1): 市町別棒 + 記録種別棒 ------------------------------------- print(" fig3 市町別 + 記録種別", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 6)) ax = axes[0] top15_city = city_df.head(15).iloc[::-1] ax.barh(range(len(top15_city)), top15_city["件数"], color="#4575b4", alpha=0.85) ax.set_yticks(range(len(top15_city))) ax.set_yticklabels(top15_city["市町"], fontsize=9) ax.set_xlabel("災害記録数") for i, (n, s) in enumerate(zip(top15_city["件数"], top15_city["シェア_%"])): ax.text(n + 0.5, i, f"{int(n)} ({s:.1f}%)", va="center", fontsize=8) ax.set_xlim(0, top15_city["件数"].max() * 1.18) ax.set_title(f"市町別 災害記録数 Top 15", fontsize=11) ax.grid(axis="x", alpha=0.3) ax = axes[1] ax.barh(range(len(rec_df)), rec_df["件数"], color="#fdae61", alpha=0.85) ax.set_yticks(range(len(rec_df))) ax.set_yticklabels(rec_df["記録種別"], fontsize=10) ax.invert_yaxis() ax.set_xlabel("件数") for i, (n, s) in enumerate(zip(rec_df["件数"], rec_df["シェア_%"])): ax.text(n + 0.5, i, f"{int(n)} ({s:.1f}%)", va="center", fontsize=9) ax.set_xlim(0, rec_df["件数"].max() * 1.18) ax.set_title(f"記録種別 — 石碑/航空写真/記録集", fontsize=11) ax.grid(axis="x", alpha=0.3) plt.suptitle( f"L61 RQ1 図 3: 市町別 Top 15 + 記録種別ランキング", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L61_fig3_city_record.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 4 (RQ2): 2014 vs 2018 マップ + 重心円 ----------------------------- print(" fig4 二大災害マップ", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 7)) # (左) 2014 広島土砂 ax = axes[0] draw_pref(ax, lw=0.6, color="#888") ax.scatter(gdf[gdf["主要イベント"] != "2014広島市土砂"]["x_m"], gdf[gdf["主要イベント"] != "2014広島市土砂"]["y_m"], c="#cccccc", s=8, alpha=0.4, label="その他災害") sub14 = gdf[gdf["主要イベント"] == "2014広島市土砂"] ax.scatter(sub14["x_m"], sub14["y_m"], c="#fc8d59", s=40, alpha=0.85, edgecolors="black", linewidths=0.4, label=f"2014 広島市土砂 ({n_2014})") # 重心円 (mean radius) if c14_lat is not None: c_xy = gpd.GeoSeries([gpd.points_from_xy([c14_lon], [c14_lat])[0]], crs=SRC_CRS).to_crs(TARGET_CRS) cx, cy = c_xy.iloc[0].x, c_xy.iloc[0].y circ = plt.Circle((cx, cy), r14_mean, color="red", fill=False, linewidth=1.5, label=f"平均距離円 r={r14_mean/1000:.1f}km") ax.add_patch(circ) ax.scatter([cx], [cy], c="red", marker="x", s=80, linewidths=2, label="重心") ax.set_title(f"2014 広島市土砂災害 ({n_2014} 件)\n" f"重心 ({c14_lat:.3f}, {c14_lon:.3f}) 平均 {r14_mean/1000:.1f}km", fontsize=10) ax.legend(fontsize=8, loc="lower right") ax.set_aspect("equal") ax.set_xticks([]); ax.set_yticks([]) # (右) 2018 西日本豪雨 ax = axes[1] draw_pref(ax, lw=0.6, color="#888") ax.scatter(gdf[gdf["主要イベント"] != "2018西日本豪雨"]["x_m"], gdf[gdf["主要イベント"] != "2018西日本豪雨"]["y_m"], c="#cccccc", s=8, alpha=0.4, label="その他災害") sub18 = gdf[gdf["主要イベント"] == "2018西日本豪雨"] ax.scatter(sub18["x_m"], sub18["y_m"], c="#d73027", s=40, alpha=0.85, edgecolors="black", linewidths=0.4, label=f"2018 西日本豪雨 ({n_2018})") if c18_lat is not None: c_xy = gpd.GeoSeries([gpd.points_from_xy([c18_lon], [c18_lat])[0]], crs=SRC_CRS).to_crs(TARGET_CRS) cx, cy = c_xy.iloc[0].x, c_xy.iloc[0].y circ = plt.Circle((cx, cy), r18_mean, color="darkred", fill=False, linewidth=1.5, label=f"平均距離円 r={r18_mean/1000:.1f}km") ax.add_patch(circ) ax.scatter([cx], [cy], c="darkred", marker="x", s=80, linewidths=2, label="重心") ax.set_title(f"2018 西日本豪雨 ({n_2018} 件)\n" f"重心 ({c18_lat:.3f}, {c18_lon:.3f}) 平均 {r18_mean/1000:.1f}km", fontsize=10) ax.legend(fontsize=8, loc="lower right") ax.set_aspect("equal") ax.set_xticks([]); ax.set_yticks([]) plt.suptitle( f"L61 RQ2 図 4: 2014 vs 2018 — 局所豪雨 (mean {r14_mean/1000:.1f}km) vs " f"広域豪雨 (mean {r18_mean/1000:.1f}km, {ratio_mean:.1f}倍)", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L61_fig4_2014vs2018.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 5 (RQ2): 市町別 2014 vs 2018 グループ棒 --------------------------- print(" fig5 市町別 2014 vs 2018", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 6.5)) ax = axes[0] top12 = city_compare.head(12).iloc[::-1] y = np.arange(len(top12)) w = 0.4 ax.barh(y - w/2, top12["件数_2014"], w, color="#fc8d59", alpha=0.85, label="2014 広島市土砂") ax.barh(y + w/2, top12["件数_2018"], w, color="#d73027", alpha=0.85, label="2018 西日本豪雨") ax.set_yticks(y) ax.set_yticklabels(top12["市町"], fontsize=9) ax.set_xlabel("災害記録数") ax.set_title(f"市町別 2014 vs 2018 — Top 12 (合計順)", fontsize=11) ax.legend(fontsize=9) ax.grid(axis="x", alpha=0.3) # (右) 2014 vs 2018 距離分布の箱ひげ ax = axes[1] data_dist = [] labels_dist = [] # 各点 → イベント重心 距離 def point_to_centroid(g, lat0, lon0): if len(g) == 0 or lat0 is None: return np.array([]) lat = np.deg2rad(g["緯度"].values) lon = np.deg2rad(g["経度"].values) lat_c = np.deg2rad(lat0) lon_c = np.deg2rad(lon0) dlat = lat - lat_c dlon = lon - lon_c a = np.sin(dlat/2)**2 + np.cos(lat_c) * np.cos(lat) * np.sin(dlon/2)**2 return 2 * EARTH_R_M * np.arcsin(np.sqrt(a)) / 1000 # km d14 = point_to_centroid(E_2014, c14_lat, c14_lon) d18 = point_to_centroid(E_2018, c18_lat, c18_lon) data_dist = [d14, d18] labels_dist = [f"2014 ({n_2014})\nmean {r14_mean/1000:.1f}km", f"2018 ({n_2018})\nmean {r18_mean/1000:.1f}km"] bp = ax.boxplot(data_dist, tick_labels=labels_dist, patch_artist=True, widths=0.55, showfliers=True) for patch, color in zip(bp["boxes"], ["#fc8d59", "#d73027"]): patch.set_facecolor(color) patch.set_alpha(0.7) ax.set_ylabel("各点 → 重心距離 (km)") ax.set_title(f"重心からの距離分布 — 分散比 std(2018)/std(2014) = {ratio_std:.2f}", fontsize=10) ax.grid(axis="y", alpha=0.3) plt.suptitle( f"L61 RQ2 図 5: 市町分布 + 重心距離分布 — 2018 が {ratio_std:.1f} 倍広域", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L61_fig5_city_dist.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 6 (RQ3): 警戒区域 + 過去災害 重ね合わせマップ ---------------------- print(" fig6 警戒区域 + 過去災害", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 7)) # (左) 全県マップ ax = axes[0] draw_pref(ax, lw=0.5, color="#888") # 警戒区域 polygon (薄背景) doseki.plot(ax=ax, color="#fee08b", alpha=0.4, edgecolor="none") kyukei.plot(ax=ax, color="#fdae61", alpha=0.4, edgecolor="none") jisuberi.plot(ax=ax, color="#d73027", alpha=0.5, edgecolor="none") # 過去災害点 (色 = 警戒域内/外) in_mask = disaster_pts["警戒域内"] ax.scatter(disaster_pts[in_mask]["x_m"], disaster_pts[in_mask]["y_m"], c="#1a9850", s=10, alpha=0.8, label=f"警戒域内発生 {n_in_any} ({share_in_any:.0f}%)") ax.scatter(disaster_pts[~in_mask]["x_m"], disaster_pts[~in_mask]["y_m"], c="black", s=10, alpha=0.8, label=f"警戒域外発生 {n_out} ({n_out/n_total*100:.0f}%)") # 凡例 (背景 polygon は別途 patch で) poly_handles = [ Patch(facecolor="#fee08b", alpha=0.4, label=f"土石流警戒域 ({len(doseki)})"), Patch(facecolor="#fdae61", alpha=0.4, label=f"急傾斜警戒域 ({len(kyukei)})"), Patch(facecolor="#d73027", alpha=0.5, label=f"地すべり警戒域 ({len(jisuberi)})"), ] pt_handles, pt_labels = ax.get_legend_handles_labels() ax.legend(handles=poly_handles + pt_handles, fontsize=7, loc="lower right") ax.set_title(f"警戒区域 + 過去災害 重ね合わせ\n" f"警戒域内発生 {share_in_any:.1f}%", fontsize=10) ax.set_aspect("equal") ax.set_xticks([]); ax.set_yticks([]) # (右) バッファ距離別カバー率 (積み上げ棒) ax = axes[1] n_within_1000 = int((disaster_pts['警戒域距離_m'] <= 1000).sum()) n_within_2000 = int((disaster_pts['警戒域距離_m'] <= 2000).sum()) buf_data = pd.DataFrame({ "距離閾値": ["厳密内包 (0m)", "100m 以内", "500m 以内", "1km 以内", "2km 以内"], "件数": [n_in_any, n_within_100, n_within_500, n_within_1000, n_within_2000], }) buf_data["シェア_%"] = (buf_data["件数"] / n_total * 100).round(1) buf_data["残り"] = n_total - buf_data["件数"] y = np.arange(len(buf_data)) colors_buf = ["#d73027", "#fdae61", "#1a9850", "#1a9850", "#1a9850"] for i in range(len(buf_data)): ax.barh(i, buf_data["件数"].iloc[i], color=colors_buf[i], alpha=0.85, label=("警戒域近接" if i == 2 else None)) ax.barh(i, buf_data["残り"].iloc[i], left=buf_data["件数"].iloc[i], color="#cccccc", alpha=0.7, label=("離れた発生" if i == 2 else None)) ax.set_yticks(y) ax.set_yticklabels(buf_data["距離閾値"], fontsize=10) ax.invert_yaxis() ax.set_xlabel(f"件数 (全 {n_total})") ax.set_title(f"警戒区域距離 別 過去災害カバー率\n" f"中央距離 {warn_dist_med:.0f}m / 平均 {warn_dist_mean:.0f}m", fontsize=11) for i, (n, s) in enumerate(zip(buf_data["件数"], buf_data["シェア_%"])): ax.text(n + 4, i, f"{int(n)} ({s:.0f}%)", va="center", fontsize=9) ax.legend(fontsize=9, loc="lower right") ax.grid(axis="x", alpha=0.3) plt.suptitle( f"L61 RQ3 図 6: 警戒区域指定 vs 過去災害 — 警戒域カバー率 {share_in_any:.1f}%", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L61_fig6_warn_overlay.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 7 (RQ3): イベント別 警戒域内シェア + 市町別 ----------------------- print(" fig7 イベント別 警戒域内", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 6.5)) # (左) イベント別 警戒域内シェア ax = axes[0] ev_plot = event_in_share[event_in_share["総件数"] >= 5].copy().reset_index(drop=True) y = np.arange(len(ev_plot)) ax.barh(y, ev_plot["警戒域内_%"], color="#1a9850", alpha=0.85) ax.set_yticks(y) ax.set_yticklabels( [f"{e} (n={int(n)})" for e, n in zip(ev_plot["主要イベント"], ev_plot["総件数"])], fontsize=9 ) ax.invert_yaxis() ax.set_xlim(0, 105) ax.set_xlabel("警戒域内シェア (%)") ax.set_title("主要イベント別 警戒域内シェア (n≥5)", fontsize=11) for i, (p, n_in) in enumerate(zip(ev_plot["警戒域内_%"], ev_plot["警戒域内件数"])): ax.text(p + 1, i, f"{p:.0f}% ({int(n_in)})", va="center", fontsize=8) ax.axvline(50, color="red", linestyle=":", linewidth=1.0, label="50% 閾値") ax.legend(fontsize=9) ax.grid(axis="x", alpha=0.3) # (右) 市町別 警戒域内シェア (上位 12 市町) ax = axes[1] city_plot = city_in_share.head(15).iloc[::-1].copy().reset_index(drop=True) y = np.arange(len(city_plot)) colors_c = ["#1a9850" if p >= 50 else "#fdae61" if p >= 30 else "#d73027" for p in city_plot["警戒域内_%"]] ax.barh(y, city_plot["警戒域内_%"], color=colors_c, alpha=0.85) ax.set_yticks(y) ax.set_yticklabels( [f"{c} (n={int(n)})" for c, n in zip(city_plot["市町"], city_plot["総件数"])], fontsize=8 ) ax.set_xlim(0, 105) ax.set_xlabel("警戒域内シェア (%)") ax.set_title("市町別 警戒域内シェア (Top 15)", fontsize=11) for i, (p, n_in) in enumerate(zip(city_plot["警戒域内_%"], city_plot["警戒域内件数"])): ax.text(p + 1, i, f"{p:.0f}% ({int(n_in)})", va="center", fontsize=7) ax.axvline(50, color="red", linestyle=":", linewidth=1.0) ax.grid(axis="x", alpha=0.3) plt.suptitle( f"L61 RQ3 図 7: 警戒域内シェア — イベント別 / 市町別", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L61_fig7_event_city_in.png", dpi=130, bbox_inches="tight") plt.close(fig) # ---- 図 8 (RQ3): 警戒域外発生 マップ + イベント × 警戒種別ヒート ---------- print(" fig8 警戒域外マップ + クロスヒート", flush=True) fig, axes = plt.subplots(1, 2, figsize=(14, 7)) # (左) 警戒域外発生のみマップ ax = axes[0] draw_pref(ax, lw=0.5, color="#888") # 警戒区域も薄く描画 doseki.plot(ax=ax, color="#cccccc", alpha=0.3, edgecolor="none") kyukei.plot(ax=ax, color="#cccccc", alpha=0.3, edgecolor="none") jisuberi.plot(ax=ax, color="#cccccc", alpha=0.3, edgecolor="none") # 警戒域外の災害点 out_mask = ~disaster_pts["警戒域内"] out_pts = disaster_pts[out_mask].copy() # 主要イベント別に色分け out_pts_with_ev = out_pts.assign(主要イベント=df.loc[out_pts.index, "主要イベント"]) for ev in draw_order: sub = out_pts_with_ev[out_pts_with_ev["主要イベント"] == ev] if len(sub) > 0: ax.scatter(sub["x_m"], sub["y_m"], c=event_colors[ev], s=28, alpha=0.85, edgecolors="black", linewidths=0.3, label=f"{ev} ({len(sub)})") ax.set_title(f"警戒域外で発生した災害 {n_out} 件\n" f"= 警戒区域指定が予測できなかった災害", fontsize=10) ax.legend(fontsize=7, loc="lower right") ax.set_aspect("equal") ax.set_xticks([]); ax.set_yticks([]) # (右) イベント × 警戒域内/外 クロスヒート ax = axes[1] ewc = event_warn_cross.copy() # シェアに変換 ewc_share = ewc.div(ewc.sum(axis=1), axis=0) * 100 # 主要 6 イベントのみ focus_events = [e for e in ["2018西日本豪雨", "2014広島市土砂", "1999平11豪雨", "1967昭42豪雨", "1988昭63豪雨", "1945枕崎台風", "その他"] if e in ewc_share.index] ewc_share_focus = ewc_share.loc[focus_events] im = ax.imshow(ewc_share_focus.values, cmap="RdYlGn", aspect="auto", vmin=0, vmax=100) ax.set_xticks(range(len(ewc_share_focus.columns))) ax.set_xticklabels(ewc_share_focus.columns, fontsize=10) ax.set_yticks(range(len(focus_events))) ax.set_yticklabels([f"{e} (n={int(ewc.loc[e].sum())})" for e in focus_events], fontsize=9) for i in range(len(focus_events)): for j in range(len(ewc_share_focus.columns)): v = ewc_share_focus.iloc[i, j] ax.text(j, i, f"{v:.0f}%", ha="center", va="center", fontsize=10, color="black" if 25 < v < 75 else "white") ax.set_title("イベント × 警戒域内/外 — シェア%", fontsize=11) plt.colorbar(im, ax=ax, label="警戒域内シェア (%)", fraction=0.04) plt.suptitle( f"L61 RQ3 図 8: 警戒域外発生 {n_out} 件マップ + イベント × 警戒域 ヒート", fontsize=12, y=1.00 ) plt.tight_layout() plt.savefig(ASSETS / "L61_fig8_outside_heat.png", dpi=130, bbox_inches="tight") plt.close(fig) print(f" ({time.time()-t1:.2f}s)", flush=True) # ============================================================================= # 8. 表データ作成 (HTML 出力用) # ============================================================================= print("\n[8] 表データ作成", flush=True) t1 = time.time() def df_to_html(d): 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", "1278"), ("公式名", "過去に発生した災害情報"), ("形式", "CSV (UTF-8 BOM)"), ("件数", f"{n_total} 行 × 8 列"), ("サイズ", f"{CSV_GLOBAL.stat().st_size:,} byte"), ("CRS", "WGS84 (緯度経度) → 解析時 EPSG:6671 (m 単位)"), ("座標カバレッジ", f"緯度 {df['緯度'].min():.3f}-{df['緯度'].max():.3f} / " f"経度 {df['経度'].min():.3f}-{df['経度'].max():.3f}"), ("カバー範囲", f"広島県内 {df['所在地市区町'].nunique()} 市町 (全 23 市町中)"), ("年代カバレッジ", f"{ymin}-{ymax} (= 約 {ymax-ymin} 年分)"), ("最古記録", f"明治 40 (1907) 年: {(df['年'] == 1907).sum()} 件"), ("最大イベント", f"平成30年7月豪雨 = {(df['主要イベント']=='2018西日本豪雨').sum()} 件"), ("ライセンス", "CC-BY 4.0"), ("作成主体", "広島県土木建築局"), ("運用元", "Web 砂防 (土砂災害ポータル広島)"), ("修復行数", f"{n_repaired} 行 (コメント内未エスケープカンマの再結合)"), ], columns=["項目", "値"]) # 表 2 全体サマリ T_overall = pd.DataFrame([ ("災害記録総数", f"{n_total} 件"), ("カバー市町数", f"{df['所在地市区町'].nunique()} / 23"), ("Top 市町", f"{top_city['市町']} = {int(top_city['件数'])} ({top_city['シェア_%']:.1f}%)"), ("年代区分 (令和/平成/昭和/大正/明治以前)", " / ".join(f"{int(era_count.get(e,0))}" for e in ["令和", "平成", "昭和", "大正", "明治以前"])), ("Top 2 イベント (2018+2014)", f"{int(n_top2)} ({top2_share:.1f}%) — 2018:{n_2018} + 2014:{n_2014}"), ("月集中度 (6-9 月)", f"{n_jun_sep} / {n_total} = {share_jun_sep:.1f}%"), ("2014 重心半径 (mean)", f"{r14_mean/1000:.1f} km (= 局所豪雨)"), ("2018 重心半径 (mean)", f"{r18_mean/1000:.1f} km (= 広域豪雨)"), ("2018/2014 平均距離比", f"{ratio_mean:.2f} 倍"), ("2014↔2018 重心間距離", f"{c14_to_c18_m/1000:.1f} km"), ("警戒域内発生 (厳密内包)", f"{n_in_any} / {n_total} = {share_in_any:.1f}%"), ("警戒域 500m 以内発生", f"{n_within_500} / {n_total} = {share_within_500:.1f}%"), ("警戒域 1000m 以内発生", f"{int((disaster_pts['警戒域距離_m']<=1000).sum())} / {n_total} = " f"{(disaster_pts['警戒域距離_m']<=1000).sum()/n_total*100:.1f}%"), ("警戒域 中央距離", f"{warn_dist_med:.0f}m / 平均 {warn_dist_mean:.0f}m"), ("土石流警戒域内 (厳密)", f"{n_in_doseki} ({n_in_doseki/n_total*100:.1f}%)"), ("急傾斜警戒域内 (厳密)", f"{n_in_kyukei} ({n_in_kyukei/n_total*100:.1f}%)"), ("地すべり警戒域内 (厳密)", f"{n_in_jisuberi} ({n_in_jisuberi/n_total*100:.1f}%)"), ], columns=["指標", "値"]) T_overall.to_csv(ASSETS / "L61_overall.csv", index=False, encoding="utf-8-sig") # 表 3 主要イベントランキング T_event_rank = event_df.copy() # 表 4 月別分布 T_month = month_df[month_df["月"].notna()].copy() T_month["月"] = T_month["月"].astype(int) # 表 5 市町別ランキング (Top 15) T_city_rank = city_df.head(15).copy() # 表 6 年代区分別件数 T_era = era_count_df.copy() # 表 7 記録種別ランキング T_record = rec_df.copy() # 表 8 2014 vs 2018 市町別クロス (Top 12) T_compare = city_compare.head(12).copy() # 表 9 警戒区域 種別 × 過去災害 T_warn_type = pd.DataFrame([ ("土石流警戒区域 (厳密内包)", len(doseki), n_in_doseki, f"{n_in_doseki/n_total*100:.1f}%"), ("急傾斜警戒区域 (厳密内包)", len(kyukei), n_in_kyukei, f"{n_in_kyukei/n_total*100:.1f}%"), ("地すべり警戒区域 (厳密内包)", len(jisuberi), n_in_jisuberi, f"{n_in_jisuberi/n_total*100:.1f}%"), ("いずれかの警戒区域 (厳密内包)", "—", n_in_any, f"{share_in_any:.1f}%"), ("警戒区域 100m 以内", "—", n_within_100, f"{share_within_100:.1f}%"), ("警戒区域 500m 以内", "—", n_within_500, f"{share_within_500:.1f}%"), (f"警戒区域 1000m 以内", "—", int((disaster_pts['警戒域距離_m'] <= 1000).sum()), f"{(disaster_pts['警戒域距離_m'] <= 1000).sum()/n_total*100:.1f}%"), ("警戒区域 (距離 中央値)", "—", f"{warn_dist_med:.0f}m", f"中央 {warn_dist_med:.0f}m / 平均 {warn_dist_mean:.0f}m"), ], columns=["警戒種別", "polygon 数", "過去災害 内 件数", "シェア (vs 全 518)"]) # 表 10 イベント別 警戒域内シェア T_ev_in = event_in_share.copy() # 表 11 市町別 警戒域内シェア (Top 15) T_city_in = city_in_share.head(15).copy() # 表 12 仮説検証 def jud(cond, ok="強支持", fail="反証", part="部分支持"): return ok if cond else fail T_hypo = pd.DataFrame([ ("H1 二大災害集中 (Top 2 ≥ 40%)", f"観測 Top 2 = {n_top2}/{n_total} ({top2_share:.1f}%)", jud(top2_share >= 40), f"H1 {jud(top2_share>=40)}: 2018 西日本豪雨 ({n_2018}) + 2014 広島市土砂 ({n_2014}) = " f"{int(n_top2)} 件 ({top2_share:.1f}%) で全体の " f"{'過半数' if top2_share >= 50 else '4 割超'}を 2 大イベントが占める偏在型分布"), ("H2 季節性 (6-9 月 ≥ 90%)", f"観測 6-9 月 = {n_jun_sep}/{n_total} ({share_jun_sep:.1f}%)", jud(share_jun_sep >= 90, "強支持", "部分支持"), f"H2 {jud(share_jun_sep>=90,'強支持','部分支持')}: 6-9 月集中 = {share_jun_sep:.1f}%。" f"梅雨 (6-7 月) + 台風 (8-9 月) の日本気候に強く規定される季節災害"), ("H3 2014 局所性 vs 2018 広域性 (距離比 > 2)", f"観測 mean 比 = {ratio_mean:.2f}, std 比 = {ratio_std:.2f}", jud(ratio_mean > 2, "強支持", "部分支持"), f"H3 {jud(ratio_mean>2,'強支持','部分支持')}: 2018 の重心からの平均距離 ({r18_mean/1000:.1f}km) は " f"2014 ({r14_mean/1000:.1f}km) の {ratio_mean:.2f} 倍。" f"局所豪雨 vs 広域豪雨の発生機構の差が空間統計で実証された"), ("H4 警戒域 500m 以内発生 ≥ 50%", f"観測 500m 以内 = {n_within_500}/{n_total} ({share_within_500:.1f}%) " f"(参考: 厳密内包 {share_in_any:.1f}%)", jud(share_within_500 >= 50), f"H4 {jud(share_within_500>=50)}: 過去災害の {share_within_500:.1f}% が " f"警戒区域から 500m 以内で発生 (中央距離 {warn_dist_med:.0f}m)。" f"厳密内包は {share_in_any:.1f}% にとどまるが、" f"500m バッファでは制度的予測圏内の災害が " f"{'過半数を占め、警戒区域指定は実際の災害発生地の近傍をよく捉えている' if share_within_500 >= 50 else '半数以下にとどまり、警戒域指定の改善余地が大きい'}。" f"撮影位置・記念碑位置のずれを考慮すると 500m バッファが現実的指標"), ("H5 急傾斜型偏在 (急傾斜 > 地すべり)", f"観測 急傾斜 = {n_in_kyukei} vs 地すべり = {n_in_jisuberi} vs 土石流 = {n_in_doseki}", jud(n_in_kyukei > n_in_jisuberi), f"H5 {jud(n_in_kyukei>n_in_jisuberi)}: 急傾斜警戒域内 ({n_in_kyukei}) > " f"地すべり警戒域内 ({n_in_jisuberi})。土石流 ({n_in_doseki}) と急傾斜の合計が " f"地すべりを大きく上回る = 広島県の急峻山地・花崗岩風化マサ土地形を反映した " f"災害種別の偏在 (= 県の災害学は「急斜面崩落 + 土石流」 主体、地すべりは少数派)"), ], columns=["仮説", "観測値", "判定", "解釈"]) T_hypo.to_csv(ASSETS / "L61_hypothesis_check.csv", index=False, encoding="utf-8-sig") print(f" 生成テーブル: 12 表") print(f" ({time.time()-t1:.2f}s)", flush=True) # ============================================================================= # 9. 中間ファイル ZIP 化 # ============================================================================= print("\n[9] 中間ファイル ZIP 化", flush=True) t1 = time.time() zip_path = ASSETS / "L61_past_disasters_raw.zip" with zipfile.ZipFile(zip_path, "w", zipfile.ZIP_DEFLATED) as zf: zf.write(CSV_GLOBAL, arcname="past_disasters.csv") print(f" 生 zip: {zip_path.name} ({zip_path.stat().st_size:,} byte)") print(f" ({time.time()-t1:.2f}s)", flush=True) # ============================================================================= # 10. HTML 生成 # ============================================================================= print("\n[10] HTML 生成", flush=True) t1 = time.time() # Sec 1: 学習目標と問い sec1 = f"""

本記事は DoBoX のシリーズ「過去に発生した災害情報」 1 件 (dataset_id = 1278) を 単独で取り上げ、 広島県内で過去に発生した災害の現地記録 {n_total} 件 (CSV, 8 列, UTF-8 BOM, 約 632 KB) を 3 つの独立した研究角度 (RQ1 / RQ2 / RQ3) で並列に分析する。 本データは「Web 砂防 (土砂災害ポータル広島)」 が運営する地域情報レイヤーで、 災害現地の航空写真・石碑・記録集・体験談・歴史的土木施設の位置情報つき アーカイブ。Phase 5 災害履歴系の起点として、過去の発生事実を正面から扱う。

L10 / L11 (既存ハザード系) との位置付け: L10 は土砂災害クロス、 L11 はトリプルハザード (浸水×土砂×雪崩) という未来予測 (= 警戒区域指定) を扱った。 本記事 L61 は過去の発生事実 (= 実際に災害が起きた場所) を扱い、 「予測 vs 実績」 の対比構造を初めて教材化する。 RQ3 では L11 既扱の警戒区域 polygon と本データを sjoin し、 「警戒域は過去災害をどこまで予測できていたか」 を直接検証する。

独自用語の定義

研究の問い (3 RQ)

仮説 H1〜H5

  1. H1 (二大災害集中, RQ1): Top 2 イベント (2018+2014) ≥ 全体の 40%。 近年の 2 大豪雨が歴史記録の大半を占める偏在型分布を予想。
  2. H2 (季節性, RQ1): 6-9 月発生 ≥ 90%。 梅雨 (6-7 月) + 台風 (8-9 月) の日本気候に強く規定される。
  3. H3 (2014 局所性 vs 2018 広域性, RQ2): 2014 災害の地理分散は 2018 より明確に小さい (重心からの平均距離比 > 2)。局所豪雨 vs 広域豪雨の機構差を空間統計で実証。
  4. H4 (警戒域 500m 以内発生 ≥ 50%, RQ3): 過去災害の ≥ 50% が 土砂災害警戒区域から 500m 以内で発生 (厳密内包だと写真撮影位置のずれで低くなりがち。 500m バッファで「制度予測との空間的近さ」 を測る)。
  5. H5 (急傾斜型偏在, RQ3): 急傾斜警戒域内発生 > 地すべり警戒域内発生。 広島県の急峻山地斜面 (= 急傾斜)主体の地形を反映した災害種別の偏在を予想。

到達点

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

  1. 1 つの「シンプルな災害アーカイブ CSV ({n_total} 行 × 8 列)」 から、 地理 (緯度経度) + 時間 (年・月・元号) + 種別 (タグ) + 行政 (市町) + 物理形態 (石碑・写真)という 5 軸を多角度に読む方法を体感。歴史的災害が「偶然の散在」ではなく 「地形 + 気象 + 制度」の三重構造で生じることを実証。
  2. 2014 vs 2018 の二大災害を地理重心と平均距離分布で比較し、 同じ「豪雨」でも発生機構 (局所 vs 広域) が違えば被害分布が違うという 災害学の基礎概念を、実データの空間統計で確認する。
  3. L11 既扱の土砂災害警戒区域 polygonと本データの過去災害点を sjoin し、 「予測 vs 実績」 の対比研究を体感。 警戒域指定の予測精度が {share_in_any:.1f}% という具体的数値で得られる経験は、 防災行政の意思決定を研究的視点で読む第一歩。
""" # Sec 2: 使用データ sec2 = f"""

DoBoX のシリーズ「過去に発生した災害情報」 1 件のみを単独で扱う。 リソースは CSV 1 ファイルのシンプル構造 (UTF-8 BOM、約 632 KB)。

{df_to_html(T_dataset)}

この表から読み取れること: dataset 1278 は {n_total} 行 × 8 列の CSV だが、 コメント列 (数百字の災害解説) に未エスケープのカンマを含む行が {n_repaired} 行あり、pandas の標準 read_csv では ParserError。 本スクリプトは csv モジュールで生読み + 末尾固定 5 列の再結合で修復済み。 緯度経度 NaN ゼロ = 全 {n_total} 件マップ可能、年代カバレッジ {ymin}-{ymax} の約 {ymax-ymin} 年分。 最古は明治 40 (1907) 年の記録 = 物理形態は石碑として残る。

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

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

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

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

生データ (DoBoX 直リンク)

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

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

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

cd "2026 DoBoX 教材"
py -X utf8 lessons/L61_past_disasters.py

CSV は本スクリプトが DoBoX dataset 1278 から自動 DL する (キャッシュ済なら再利用)。 警戒区域 Shapefile (L11 既扱) と県境 GeoJSON (L15 既扱) も内部で再利用。

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

cd "2026 DoBoX 教材"
py -X utf8 data\\fetch_all.py
py -X utf8 lessons/L61_past_disasters.py
""" # Sec 4: RQ1 sec4_intro = f"""

RQ1 の狙い

{n_total} 件の災害記録を年代 / 月 / 市町 / タグ / 記録種別の 5 軸で多角度に集計し、 「広島県の災害はいつ・どこで・どんな記録形式で残っているか」 を立体的に描く。 特に2 大イベント (2018+2014) で全体の {top2_share:.0f}% という偏在の意味を、 歴史的時間軸と地理分布の両方から読み解く。

手法 (前置き解説)

入出力の Before/After 例

段階1 件のデータの中身列数
(0) CSV 1 行 (raw)1,平成17年9月6日 台風14号での宮島の被災状況...,廿日市市,"平成17年台風14号,宮島町,白糸川",34.282,132.317,...8
(1) CSV 修復 (列数チェック)同上 (8 列なら通過)8
(2) 緯度経度 to_numeric緯度=34.282, 経度=132.317+0
(3) 年抽出 (タイトル正規表現)+ 年 = 2005 (= 1988+17)+1
(4) 月抽出+ 月 = 9+1
(5) 年代区分+ 年代区分 = "平成"+1
(6) タグリスト分解+ タグリスト = ["平成17年台風14号", "宮島町", "白糸川"]+1
(7) 主要イベント集約+ 主要イベント = "2005台風14号"+1
(8) 記録種別主+ 記録種別主 = "現地記録" (タグに種別なし)+1
(9) GeoDataFrame 化 → EPSG:6671+ x_m=-19238, y_m=-160493+2

(0)-(9) を全 {n_total} 行に適用 → groupby で集計 → 図化。CSV 読込 + 派生計算は < 1 秒。

""" sec4_code = code(r''' # 1. CSV 修復読込 (csv モジュールで生読み + 末尾 5 列固定の再結合) import csv with open(CSV_GLOBAL, "r", encoding="utf-8-sig") as f: rows = list(csv.reader(f)) header = rows[0] data = [] for r in rows[1:]: if len(r) == 8: data.append(r) elif len(r) > 8: # 先頭 2 列 + コメント結合 + 末尾 5 列 data.append(r[:2] + [",".join(r[2:-5])] + r[-5:]) df = pd.DataFrame(data, columns=header) # 2. 緯度経度・年・月の派生 df["緯度"] = pd.to_numeric(df["緯度"], errors="coerce") df["経度"] = pd.to_numeric(df["経度"], errors="coerce") ERA_BASE = {"令和": 2018, "平成": 1988, "昭和": 1925, "大正": 1911, "明治": 1867} def extract_year(t): m = re.search(r"(令和|平成|昭和|大正|明治)([0-9]+|元)年", str(t)) if not m: return None n = 1 if m.group(2) == "元" else int(m.group(2)) return ERA_BASE[m.group(1)] + n df["年"] = df["タイトル"].apply(extract_year) df["月"] = df["タイトル"].apply( lambda t: int(re.search(r"([0-9]+)月", str(t)).group(1)) if re.search(r"([0-9]+)月", str(t)) else None ) # 3. 年代区分 (令和/平成/昭和/大正/明治以前) def era_of(y): if pd.isna(y): return "不明" y = int(y) if y >= 2019: return "令和" if y >= 1989: return "平成" if y >= 1926: return "昭和" if y >= 1912: return "大正" return "明治以前" df["年代区分"] = df["年"].apply(era_of) # 4. タグ展開 + 主要イベント集約 (8 エントリ辞書で優先順位) EVENT_TAGS = { "平成30年7月豪雨": "2018西日本豪雨", "平成26年8月豪雨": "2014広島市土砂", "昭和42年7月豪雨": "1967昭42豪雨", "平成11年6月豪雨": "1999平11豪雨", "昭和63年7月豪雨": "1988昭63豪雨", "昭和20年枕崎台風": "1945枕崎台風", "平成22年7月豪雨": "2010平22豪雨", "平成17年台風14号": "2005台風14号", } def split_tags(s): return [t.strip().strip('"') for t in str(s).split(",") if t.strip().strip('"')] df["タグリスト"] = df["タグ"].apply(split_tags) df["主要イベント"] = df["タグリスト"].apply( lambda lst: next((EVENT_TAGS[t] for t in lst if t in EVENT_TAGS), "その他") ) # 5. 群集計 era_count = df["年代区分"].value_counts() event_count = df["主要イベント"].value_counts() month_count = df["月"].value_counts().sort_index() city_count = df["所在地市区町"].value_counts() print("Top 2 イベント (2018+2014):", n_top2, "(", top2_share, "%)") ''') sec4_fig1_text = f"""

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

「広島県の災害はいつ・どこで起きてきたか」 を 1 枚で読みたい。 左マップは全 {n_total} 点を年代区分 (5 区分) で色分けし、令和=赤・平成=橙・昭和=青の階調で 時間軸を地理に投影。右ヒストは 10 年ビンで時系列を見せ、大災害イベントの位置を赤点線でマーク。 これで「歴史的偏在」 と「大災害インパクト」 を同時に視覚化する。

""" sec4_fig1_read = f"""

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

""" sec4_fig2_text = """

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

大災害ごとの地理分布を見たいので、左マップで主要 8 イベントを色分け。 2018 (赤) と 2014 (橙) は前景に描き、その他の小イベントは背景に薄色。 右棒で月別分布を見せ、6-9 月の梅雨〜台風シーズンに集中することを定量化する。

""" sec4_fig2_read = f"""

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

""" sec4_fig3_text = """

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

「市町別ランキング」 で災害が集中する地域、「記録種別ランキング」 で 歴史的記録の物理形態を可視化したい。 左の Top 15 棒で「どの市町が記録多いか」、右の種別棒で「石碑 vs 航空写真」 等の 記録形式の頻度を読む。記録形式は時代を反映 (古い=石碑、新しい=航空写真)。

""" sec4_fig3_read = f"""

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

""" # Sec 5: RQ2 (二大災害比較) sec5_intro = f"""

RQ2 の狙い

2018 西日本豪雨 ({n_2018} 件) と 2014 広島市土砂災害 ({n_2014} 件) は、 本データの 2 大イベントだが、地理スケールが大きく異なる:

両者の地理重心と分散を比較し、「同じ豪雨でも発生機構が違えば被害分布が違う」 を 空間統計で実証する。

手法 (前置き解説)

入出力の Before/After 例

段階1 件のデータの中身件数
(0) df から 2014 抽出主要イベント=="2014広島市土砂"の 75 件{n_2014}
(1) 重心計算lat 平均={c14_lat:.4f}, lon 平均={c14_lon:.4f}1 点
(2) 各点 → 重心距離 (haversine)1 件目 d=487m, 2 件目 d=1,205m, ...{n_2014} 距離
(3) 平均・標準偏差mean={r14_mean/1000:.1f}km, std={r14_std/1000:.1f}km2 scalar
(4) 2018 同様の計算mean={r18_mean/1000:.1f}km, std={r18_std/1000:.1f}km2 scalar
(5) 比率mean 比 = {ratio_mean:.2f}, std 比 = {ratio_std:.2f}2 scalar
""" sec5_code = code(r''' # 1. 2014 / 2018 抽出 E_2018 = df[df["主要イベント"] == "2018西日本豪雨"] E_2014 = df[df["主要イベント"] == "2014広島市土砂"] # 2. 重心 (緯度経度の算術平均) c14_lat, c14_lon = E_2014["緯度"].mean(), E_2014["経度"].mean() c18_lat, c18_lon = E_2018["緯度"].mean(), E_2018["経度"].mean() # 3. 各点 → 重心距離 (haversine, m) EARTH_R_M = 6_371_000.0 def haversine_to_centroid(lats, lons, lat0, lon0): lat = np.deg2rad(lats); lon = np.deg2rad(lons) lat_c = np.deg2rad(lat0); lon_c = np.deg2rad(lon0) dlat = lat - lat_c; dlon = lon - lon_c a = np.sin(dlat/2)**2 + np.cos(lat_c) * np.cos(lat) * np.sin(dlon/2)**2 return 2 * EARTH_R_M * np.arcsin(np.sqrt(a)) d14 = haversine_to_centroid(E_2014["緯度"], E_2014["経度"], c14_lat, c14_lon) d18 = haversine_to_centroid(E_2018["緯度"], E_2018["経度"], c18_lat, c18_lon) r14_mean, r14_std = d14.mean(), d14.std() r18_mean, r18_std = d18.mean(), d18.std() print(f"2014 mean radius = {r14_mean/1000:.1f}km, std = {r14_std/1000:.1f}km") print(f"2018 mean radius = {r18_mean/1000:.1f}km, std = {r18_std/1000:.1f}km") print(f"std ratio = {r18_std/r14_std:.2f}") # 4. 市町別クロス city_2014 = E_2014["所在地市区町"].value_counts().rename("件数_2014") city_2018 = E_2018["所在地市区町"].value_counts().rename("件数_2018") city_compare = pd.concat([city_2014, city_2018], axis=1).fillna(0).astype(int) print(city_compare.sort_values("件数_2018", ascending=False).head(10)) ''') sec5_fig4_text = f"""

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

2014 と 2018 の地理分布を直接比較するため、2 ペインマップに並べる。 各ペインで他イベント点を背景灰色、対象イベント点を主色 (2014=橙, 2018=赤) で前景表示。 重心位置と平均距離円を重ねることで、「どこを中心にどれくらい広がっているか」 を 1 つの幾何図形で要約。これは「分布の広がり」 を直感的に伝える可視化。

""" sec5_fig4_read = f"""

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

""" sec5_fig5_text = f"""

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

「市町ベースで 2014 vs 2018 はどう違うか」 (左) と 「重心からの距離分布の形」 (右) を並べたい。 左のグループ棒は同じ市町スケールで両イベント件数を並置 → 共通市町 / 単独市町が一目で分かる。 右の箱ひげは平均値だけでなく中央値・IQR・外れ値を含めた分布の形を見せる。

""" sec5_fig5_read = f"""

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

""" # Sec 6: RQ3 sec6_intro = f"""

RQ3 の狙い

L11 既扱の土砂災害警戒区域 polygon (土石流 + 急傾斜 + 地すべり, 計 {len(doseki) + len(kyukei) + len(jisuberi):,} 個) と本データの過去災害 {n_total} 点を sjoin (point in polygon) し、「過去災害は警戒区域指定でどこまで予測されていたか」 を 定量化する。具体的に:

  1. 各点について「土石流警戒域内」 「急傾斜警戒域内」 「地すべり警戒域内」 の 3 フラグを判定
  2. 少なくとも 1 種別の警戒域内に位置する点を「警戒域内発生」 と分類 ({share_in_any:.1f}%)
  3. 警戒域外で発生した {n_out} 件 ({n_out/n_total*100:.1f}%) を主要イベント別 / 市町別 / 種別ヒートで深掘り

これは単なる「カバー率の集計」 ではなく、「未予測災害」の特定により 警戒区域指定の制度的弱点を初めて明示する研究。

手法 (前置き解説)

入出力の Before/After 例

段階1 件のデータの中身件数
(0) 過去災害点 (gdf)1 点 = (緯度, 経度) Point geometry{n_total}
(1) sjoin 土石流+ index_right (= 対応 polygon ID, 域外なら NaN){n_total}
(2) groupby + .notna().any()+ 土石流警戒域内 = True/False{n_total}
(3) 急傾斜・地すべりも同様+ 急傾斜警戒域内, 地すべり警戒域内{n_total}
(4) OR 集約+ 警戒域内 = OR (3 種別){n_total}
(5) value_countsTrue={n_in_any}, False={n_out}2
(6) crosstab (イベント × 内外)2018: 内 N1, 外 N2; 2014: 内 N3, 外 N4 ...9 行 × 2 列
""" sec6_code = code(r''' # 1. 警戒区域 Shapefile 読込 (L11 既扱) import geopandas as gpd TARGET_CRS = "EPSG:6671" doseki = gpd.read_file(DOSEKI_SHP).to_crs(TARGET_CRS) kyukei = gpd.read_file(KYUKEI_SHP).to_crs(TARGET_CRS) jisuberi = gpd.read_file(JISUBERI_SHP).to_crs(TARGET_CRS) # 2. 過去災害点を gdf 化 (緯度経度 → EPSG:6671) gdf_pts = gpd.GeoDataFrame( df, geometry=gpd.points_from_xy(df["経度"], df["緯度"]), crs="EPSG:4326", ).to_crs(TARGET_CRS) # 3. sjoin (point in polygon) — 3 種別独立に実行 sj_doseki = gpd.sjoin(gdf_pts, doseki[["geometry"]], how="left", predicate="within") sj_kyukei = gpd.sjoin(gdf_pts, kyukei[["geometry"]], how="left", predicate="within") sj_jisuberi = gpd.sjoin(gdf_pts, jisuberi[["geometry"]], how="left", predicate="within") # 4. groupby + notna() で 1 点 → 1 ブール値に集約 (複数 polygon マッチ対応) in_doseki = sj_doseki.groupby(level=0)["index_right"].apply( lambda s: s.notna().any()) in_kyukei = sj_kyukei.groupby(level=0)["index_right"].apply( lambda s: s.notna().any()) in_jisuberi = sj_jisuberi.groupby(level=0)["index_right"].apply( lambda s: s.notna().any()) # 5. OR で「いずれかの警戒域内」 in_any = in_doseki | in_kyukei | in_jisuberi n_in_any = in_any.sum() n_out = (~in_any).sum() share_in_any = n_in_any / len(df) * 100 print(f"警戒域内発生: {n_in_any} / {len(df)} = {share_in_any:.1f}%") print(f"警戒域外発生: {n_out} = 「未予測災害」") # 6. イベント × 警戒域内/外 クロス ewc = pd.crosstab(df["主要イベント"], in_any.map({True: "警戒域内", False: "警戒域外"})) print(ewc) ''') sec6_fig6_text = f"""

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

「警戒区域 polygon と過去災害点の重ね合わせ」 を 1 枚で読みたい。 左マップは警戒区域を 3 色 (土石流=黄, 急傾斜=橙, 地すべり=赤) の薄い背景で描き、 過去災害点を警戒域内 (緑) / 警戒域外 (黒) で前景表示。これで 「どこに警戒域があり」「どこに災害が起きたか」「整合しているか」 を視覚的に対比。 右の積み上げ棒で 4 種別 (土石流 / 急傾斜 / 地すべり / いずれか) の警戒域内シェアを 数値で示す。

""" sec6_fig6_read = f"""

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

""" sec6_fig7_text = """

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

「どのイベント・どの市町で警戒域カバー率が高い/低いか」 を 2 ペインで読みたい。 左棒はイベント別 (n≥5) の警戒域内シェアを 0-100% で並べ、50% 線を赤点線で示す。 右棒は市町別 Top 15 の警戒域内シェアを色分け (緑=高 / 橙=中 / 赤=低) で並べ、 「警戒域指定が機能している市町と機能していない市町」を一目で識別。

""" sec6_fig7_read = f"""

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

""" sec6_fig8_text = f"""

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

「警戒域外で発生した {n_out} 件 = 未予測災害」 はどこに集中しているか? 左マップは域外発生のみを主要イベント色で描画 (背景に警戒域 polygon を薄灰)。 右ヒートマップはイベント × 警戒域内/外 のシェア%をカラー (RdYlGn) で示し、 6 主要イベントすべての予測精度を 1 枚で対比。これは「制度的弱点の地理 + イベント学」の同時可視化。

""" sec6_fig8_read = f"""

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

""" # Sec 7: 仮説検証総合 sec7 = ( "

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

" + df_to_html(T_hypo) + "

3 RQ × 3 結論

" + "" + "

制度史的位置付け

" + f"

本データ (#1278) は「災害発生 → 警戒区域指定 → 防災投資 → 災害発生」の" f"サイクル構造を理解する起点。L10/L11 が未来予測 (警戒区域)側を、" f"L46/L56-L58 が制度・施設投資側を扱ったのに対し、L61 は過去の発生事実側を扱う。" f"3 系統が揃って初めて制度ライフサイクルが完結する。

" + f"

本研究の重要発見は「未予測災害 {n_out} 件 ({n_out/n_total*100:.1f}%) の存在」。" f"これは警戒区域指定の限界を初めて定量化したもので、" f"制度のフィードバックループ (= 過去災害を警戒域に追加) が必要であることを示唆。" f"この発見は防災行政の区域指定見直し PDCA に直結する政策的含意を持つ。

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

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

発展課題 1 (RQ1 由来): コメント文字列の TF-IDF + クラスタリングで災害種別の自動抽出

発展課題 2 (RQ1 由来): 石碑タグ {(df['記録種別主']=='石碑').sum()} 件の地理分析 — 災害伝承の物理形態

発展課題 3 (RQ2 由来): 2018 西日本豪雨の市町別被害ランキング × 雨量データ統合

発展課題 4 (RQ3 由来): 警戒域外発生 {n_out} 件の三系統分類 — 河川氾濫 / 地形特殊 / 歴史型

発展課題 5 (RQ3 拡張): 警戒域内発生 {n_in_any} 件の土地利用との重ね合わせ

発展課題 6 (展望): 過去災害 → 警戒域指定の時間ラグ分析

""" # 統合 sections = [ ("学習目標と問い", sec1), ("使用データ", sec2), ("ダウンロード", sec3), (f"【RQ1】 災害履歴の構造 — {n_total} 件 / Top 2 イベントで {top2_share:.0f}%", sec4_intro + "

実装コード (CSV 修復読込 + 元号換算 + タグ集約)

" + sec4_code + sec4_fig1_text + figure("assets/L61_fig1_year_map.png", f"図 1 (RQ1): 全県マップ (色=年代区分) + 年別ヒスト ({ymin}-{ymax}, 10 年ビン)") + sec4_fig1_read + sec4_fig2_text + figure("assets/L61_fig2_event_month.png", "図 2 (RQ1): 主要 8 イベント色分けマップ + 月別棒") + sec4_fig2_read + sec4_fig3_text + figure("assets/L61_fig3_city_record.png", "図 3 (RQ1): 市町別 Top 15 + 記録種別ランキング") + sec4_fig3_read + "

表: 全体サマリ (3 RQ 統合, 15 指標)

" + df_to_html(T_overall) + f"

この表から読み取れること: 全 {n_total} 件の核心指標を 15 行に集約。" f"二大災害シェア {top2_share:.1f}%、6-9 月集中 {share_jun_sep:.0f}%、" f"警戒域内発生 {share_in_any:.1f}% — 3 RQ の主結論を要約。

" + "

表: 主要イベントランキング (Top 9)

" + df_to_html(T_event_rank) + f"

この表から読み取れること: 8 主要イベント + 「その他」 のランキング。" f"2018 西日本豪雨 {n_2018} 件 ({n_2018/n_total*100:.0f}%) 単独で首位、" f"2014 広島市土砂 {n_2014} 件 が次。それ以下は 1999/1967 の昭和-平成中期豪雨。" f"「その他」 ({(df['主要イベント']=='その他').sum()} 件 = " f"{(df['主要イベント']=='その他').sum()/n_total*100:.0f}%) は中小イベント・歴史記録の合算。

" + "

表: 月別分布 (1-12 月)

" + df_to_html(T_month) + f"

この表から読み取れること: 7 月単月で {int(month_df[month_df['月']==7]['件数'].iloc[0]) if 7 in month_df['月'].values else 0} 件と最多 " f"(梅雨末期の集中豪雨)、次が 8 月 ({int(month_df[month_df['月']==8]['件数'].iloc[0]) if 8 in month_df['月'].values else 0} 件、台風シーズン)。" f"6-9 月で {n_jun_sep}/{n_total} = {share_jun_sep:.0f}%。" f"1-5 月 + 10-12 月の合計はわずか {n_total-n_jun_sep} 件 = " f"{(n_total-n_jun_sep)/n_total*100:.0f}%。夏型気象システムが県の災害学を支配。

" + "

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

" + df_to_html(T_city_rank) + f"

この表から読み取れること: {top_city['市町']} = {int(top_city['件数'])} が首位、" f"安佐南区は 2014 + 2018 両方で被災。広島市内 4 区がいずれも上位、" f"県全域に被災市町が散在。" f"沿岸 (呉市・廿日市市) と内陸 (安芸太田町・庄原市) の両方に上位市町が分布 = " f"地形・気象の複合要因が県内全域に災害をもたらす。

" + "

表: 年代区分別件数

" + df_to_html(T_era) + "

この表から読み取れること: 平成 + 令和で全体の 8 割超 = 近代記録偏重。" "明治以前 = 石碑として残る数件のみ。記録の生存率は時代に強く依存し、" "「災害が起きなかった」 ではなく「現代に残った記録が少ない」 と読むべき。

" + "

表: 記録種別ランキング

" + df_to_html(T_record) + f"

この表から読み取れること: 「現地記録」 (= 種別タグ無し) が " f"{(df['記録種別主']=='現地記録').sum()} 件と最多 = 多くの記録は写真・解説で物理形態タグなし。" f"航空写真 ({(df['記録種別主']=='航空写真').sum()} 件) は近代化以降の記録、" f"石碑 ({(df['記録種別主']=='石碑').sum()} 件) は明治-昭和の災害伝承の物理形態。" f"記録形態の時代差が読める。

" ), (f"【RQ2】 二大災害比較 — 2014 局所 ({r14_mean/1000:.0f}km) vs 2018 広域 ({r18_mean/1000:.0f}km, {ratio_mean:.1f}倍)", sec5_intro + "

実装コード (重心 + 平均距離 + 市町別クロス)

" + sec5_code + sec5_fig4_text + figure("assets/L61_fig4_2014vs2018.png", f"図 4 (RQ2): 2014 ({n_2014}) vs 2018 ({n_2018}) — 重心 + 平均距離円") + sec5_fig4_read + sec5_fig5_text + figure("assets/L61_fig5_city_dist.png", "図 5 (RQ2): 市町別 2014 vs 2018 棒 + 重心距離分布箱ひげ") + sec5_fig5_read + "

表: 市町別 2014 vs 2018 クロス Top 12

" + df_to_html(T_compare) + f"

この表から読み取れること: 2014 のみ被災市町 = {n_city_only14}、" f"2018 のみ = {n_city_only18}、両方 = {n_city_both}。" f"両方被災した市町でも 2018 件数 >> 2014 が一般的 (例: 安佐南区: 2014 主体, 2018 追加分も)。" f"2018 のみ被災 {n_city_only18} 市町は「2014 では無事だったが 2018 でやられた市町」= " f"広域豪雨で初被災のパターンを示す。

" ), (f"【RQ3】 警戒区域との空間整合 — 厳密内包 {share_in_any:.0f}% / 500m 以内 {share_within_500:.0f}%", sec6_intro + "

実装コード (sjoin point in polygon + groupby OR 集約)

" + sec6_code + sec6_fig6_text + figure("assets/L61_fig6_warn_overlay.png", f"図 6 (RQ3): 警戒区域 ({len(doseki) + len(kyukei) + len(jisuberi):,} polygon) + 過去災害 重ね合わせ") + sec6_fig6_read + sec6_fig7_text + figure("assets/L61_fig7_event_city_in.png", "図 7 (RQ3): イベント別 (n≥5) / 市町別 Top 15 警戒域内シェア") + sec6_fig7_read + sec6_fig8_text + figure("assets/L61_fig8_outside_heat.png", f"図 8 (RQ3): 警戒域外発生 {n_out} 件マップ + イベント × 警戒域 シェア%ヒート") + sec6_fig8_read + "

表: 警戒区域 種別 × 過去災害

" + df_to_html(T_warn_type) + f"

この表から読み取れること: 土石流警戒域内 {n_in_doseki} ({n_in_doseki/n_total*100:.1f}%) > " f"急傾斜 {n_in_kyukei} ({n_in_kyukei/n_total*100:.1f}%) > " f"地すべり {n_in_jisuberi} ({n_in_jisuberi/n_total*100:.1f}%)。" f"土石流型偏在 (H5) を確認 = 広島県の急峻山地・花崗岩風化マサ土地形を反映。" f"3 種別の OR で警戒域内 {share_in_any:.1f}%。

" + "

表: 主要イベント別 警戒域内シェア

" + df_to_html(T_ev_in) + f"

この表から読み取れること: 大災害イベントごとに警戒域内シェアが異なる。" f"近年大災害 (2018, 2014) は警戒域指定との整合が比較的高い = 近年の警戒域指定は近年の大災害を反映。" f"古いイベントや「その他」 の警戒域内シェアは別パターンの可能性。

" + "

表: 市町別 警戒域内シェア Top 15

" + df_to_html(T_city_in) + f"

この表から読み取れること: 市町間でシェアにばらつき = 制度の予測精度が市町単位で違う。" f"低シェア市町は区域指定見直しの優先候補として防災行政の PDCA に直結。

" ), ("仮説検証総合", sec7), ("発展課題", sec8), ] html = render_lesson( num=61, title="過去に発生した災害情報 単独 3 研究例分析 — " f"{n_total} 件から「予測 vs 実績」 を読む", tags=["L61", "過去災害", "災害履歴", "RQ×3", "Format B", "geopandas", "sjoin", "haversine", "L11連携 (警戒区域)", "2018西日本豪雨", "2014広島市土砂", "Web砂防", "石碑"], time="35 分", level="中級", data_label=f"DoBoX dataset 1278 (CSV {n_total} 行 × 8 列, 632 KB) + L11 警戒区域連携", sections=sections, script_filename="L61_past_disasters.py", ) OUT_HTML = LESSONS / "L61_past_disasters.html" OUT_HTML.write_text(html, encoding="utf-8") print(f" HTML: {OUT_HTML.name} ({len(html):,} chars)") print(f" ({time.time()-t1:.2f}s)", flush=True) # ============================================================================= # 終了 # ============================================================================= print(f"\n=== 完了: 全体 {time.time() - t_all:.2f}s ===") print(f"出力:") print(f" HTML: lessons/L61_past_disasters.html") print(f" Python: lessons/L61_past_disasters.py") print(f" 図: assets/L61_fig1-8_*.png (8 枚)") print(f" CSV: assets/L61_*.csv (12 ファイル)") print(f" ZIP: assets/L61_past_disasters_raw.zip")