"""L09 浸水継続時間 — 「致命的滞水域」の発見 研究の核: 浸水深 (X09) は「どれだけ深く沈むか」だが、 継続時間 = 「どれだけ長く沈み続けるか」。 3m が 30 分で引く ≠ 3m が 24 時間続く（後者は救助困難=致命的）。 本記事の発見: - 河川浸水継続時間 Shapefile (DoBoX #332) の `rank` 列は時間ランク (10〜70) - これを浸水深 (#1278 想定最大規模) と同一座標系で重ねると、 「深さ × 時間」の 2 軸危険度マトリクスが描ける - 致命的滞水域 (深 3m 以上 × 時間 12h 以上) を抽出 → 救助困難ゾーンを地理同定 カバー宣言: 本記事は河川浸水継続時間 全河川版 (#332) の Shapefile 1 ファイルを使用。 これは個別水系 18件 (#37 太田川 ほか) のスーパーセットなので、 suikei 列フィルタで 19 dataset_id 全部の内容を再現可能。 要件S準拠: 1分以内、最悪3分以内で完走。 """ from __future__ import annotations import sys, time, json from pathlib import Path sys.path.insert(0, str(Path(__file__).parent)) from _common import (ROOT, ASSETS, LESSONS, render_lesson, code, figure) import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.patches import Patch import geopandas as gpd import shapely plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False t0 = time.time() print("=== L09 浸水継続時間 — 致命的滞水域の発見 ===") # ============================================================================= # 0. 定数 (時間ランク・深さランク・色定義) # ============================================================================= TARGET_CRS = "EPSG:6671" # JGD2011 平面直角座標 III (m単位で面積計算) # 浸水継続時間ランク (国交省 浸水想定区域図作成マニュアル + 広島県標準凡例) TIME_LABEL = { 10: "12時間未満", 20: "12〜24時間", 30: "1〜3日", 40: "3日〜1週間", 50: "1〜2週間", 60: "2週間〜1ヶ月", 70: "1ヶ月以上", } # rank → 中央値 (時間) TIME_HOURS_MID = { 10: 6, 20: 18, 30: 48, 40: 120, 50: 252, 60: 504, 70: 1080, } # rank → 「12時間以上」の判定フラグ (致命的滞水基準) TIME_GE12H = {10: False, 20: True, 30: True, 40: True, 50: True, 60: True, 70: True} TIME_COLOR = { 10: "#fef3c7", 20: "#fde68a", 30: "#fbbf24", 40: "#f59e0b", 50: "#dc2626", 60: "#991b1b", 70: "#450a0a", } # 浸水深ランク (X09 と同じ凡例) DEPTH_LABEL = { 10: "0.0〜0.5m", 20: "0.5〜1.0m", 30: "1.0〜2.0m", 40: "2.0〜3.0m", 50: "3.0〜5.0m", 60: "5.0〜10.0m", 70: "10.0〜20.0m", 80: "20m以上", } DEPTH_COLOR = { 10: "#bee2ff", 20: "#87c4f0", 30: "#56a4dc", 40: "#2c83c4", 50: "#1c63a4", 60: "#0e4282", 70: "#7d2cbf", 80: "#4a1280", } # ============================================================================= # 1. 浸水継続時間 Shapefile 読込 # ============================================================================= KEIZOKU_SHP = ROOT / "data" / "extras" / "flood_keizoku_shp" / "shinsui_keizokuzikan.shp" print(f"\n--- 1. 浸水継続時間 Shapefile 読込 ---") print(f" source: {KEIZOKU_SHP.name}") trd = time.time() keizoku = gpd.read_file(KEIZOKU_SHP, engine="pyogrio") print(f" read: {time.time()-trd:.1f}s, rows: {len(keizoku)}, CRS: {keizoku.crs}") print(f" rank ユニーク値: {sorted(keizoku['rank'].unique())}") print(f" 水系数: {keizoku['suikei'].nunique()}, 河川数: {keizoku['kasen'].nunique()}") trd = time.time() keizoku = keizoku.to_crs(TARGET_CRS) keizoku["geometry"] = gpd.GeoSeries(shapely.force_2d(keizoku.geometry.values), crs=TARGET_CRS) keizoku["geometry"] = keizoku.geometry.buffer(0) # トポロジ崩れ修正 print(f" CRS変換+force_2d+buffer: {time.time()-trd:.1f}s") keizoku["area_m2"] = keizoku.geometry.area keizoku["area_ha"] = keizoku["area_m2"] / 10000 keizoku["time_label"] = keizoku["rank"].map(TIME_LABEL).fillna("不明") keizoku["time_h_mid"] = keizoku["rank"].map(TIME_HOURS_MID).fillna(0) keizoku["ge12h"] = keizoku["rank"].map(TIME_GE12H).fillna(False) print(f" 全継続時間ポリゴン総面積: {keizoku['area_ha'].sum():.0f} ha") # ============================================================================= # 2. 浸水深 Shapefile 読込 (X09 と同じソース、想定最大規模) # ============================================================================= DEPTH_SHP = ROOT / "data" / "extras" / "flood_shp" / "shinsui_souteisaidai" / "shinsui_souteisaidai.shp" print(f"\n--- 2. 浸水深 Shapefile 読込 (想定最大規模) ---") trd = time.time() flood_depth = gpd.read_file(DEPTH_SHP, engine="pyogrio", columns=["rank"]) print(f" read: {time.time()-trd:.1f}s, rows: {len(flood_depth)}") trd = time.time() flood_depth = flood_depth.to_crs(TARGET_CRS) flood_depth["geometry"] = gpd.GeoSeries(shapely.force_2d(flood_depth.geometry.values), crs=TARGET_CRS) # buffer(0) は省略 (このデータは representative_point と面積計算にしか使わないため) print(f" CRS+force_2d: {time.time()-trd:.1f}s") print(f" rank ユニーク: {sorted(flood_depth['rank'].unique())}") # ============================================================================= # 3. 水系別 集計表 (基本統計) # ============================================================================= print(f"\n--- 3. 水系別 集計 ---") suikei_total = keizoku.groupby("suikei").agg( total_ha=("area_ha", "sum"), n_polygons=("rank", "size"), max_rank=("rank", "max"), mean_h=("time_h_mid", lambda s: float(np.average(s.values, weights=keizoku.loc[s.index, "area_ha"].values))), ).reset_index().sort_values("total_ha", ascending=False) suikei_total["max_label"] = suikei_total["max_rank"].map(TIME_LABEL) suikei_total.to_csv(ASSETS / "L09_suikei_summary.csv", index=False, encoding="utf-8-sig") print(suikei_total.head(10).to_string(index=False)) # ============================================================================= # 4. 時間ランク別 集計 (全県) # ============================================================================= print(f"\n--- 4. 時間ランク別 集計 ---") time_total = keizoku.groupby("rank").agg( total_ha=("area_ha", "sum"), n_polygons=("rank", "size"), ).reset_index() time_total["time_label"] = time_total["rank"].map(TIME_LABEL) time_total = time_total.sort_values("rank") time_total.to_csv(ASSETS / "L09_time_rank_summary.csv", index=False, encoding="utf-8-sig") print(time_total.to_string(index=False)) # ============================================================================= # 5. 深さ × 時間 の 2 軸クロス (致命的滞水域 = 3m以上 × 12h以上) # 重い overlay は深 rank≥40 (2m以上) と 時間全体 で範囲を絞ってから実施 # ============================================================================= print(f"\n--- 5. 深さ × 時間 結合 (代表点 sjoin で高速近似) ---") t1 = time.time() # 深 2m 以上 (rank≥40) のみ抽出 → 代表点を作成 deep = flood_depth[flood_depth["rank"] >= 40].copy() deep["depth_rank"] = deep["rank"] deep["depth_area_m2"] = deep.geometry.area # 代表点（ポリゴン内に必ず存在する点）で時間ランクを取得 deep_rep = gpd.GeoDataFrame( deep[["depth_rank", "depth_area_m2"]].copy(), geometry=deep.geometry.representative_point(), crs=TARGET_CRS, ) deep_rep["deep_id"] = range(len(deep_rep)) print(f" 深 2m 以上: {len(deep_rep)} ポリゴン (代表点化)") # sjoin: 深ポリゴンの代表点 in 時間ポリゴン → どの time_rank/suikei に属するか keizoku_for_join = keizoku[["rank", "suikei", "kasen", "geometry"]].rename(columns={"rank": "time_rank"}) deep_with_time = gpd.sjoin(deep_rep, keizoku_for_join, how="inner", predicate="within") print(f" 深×時間 sjoin: {len(deep_with_time)} 行") # overlap_m2 = 深ポリゴン面積（代表点による近似: 深ポリ全体がその時間ランクに属すると仮定） # これは「深ポリゴンの大半が同じ時間ランク内にある」前提。 # 代表点が時間域内にない深ポリゴンは sjoin で落ちる (= 時間域外 と判定) deep_with_time["overlap_m2"] = deep_with_time["depth_area_m2"] deep_with_time["overlap_ha"] = deep_with_time["overlap_m2"] / 10000 ovl = deep_with_time.rename(columns={"depth_rank": "depth_rank", "time_rank": "time_rank"}) ovl = ovl[["depth_rank", "time_rank", "suikei", "kasen", "overlap_m2", "overlap_ha", "deep_id"]] print(f" 深×時間結合 完了: {time.time()-t1:.1f}s, 結果 {len(ovl)} 行") # deep_id を deep_rep の geometry とリンクするための辞書 deep_full_geoms = dict(zip(deep_rep["deep_id"].values, deep.geometry.values)) # 致命的滞水域 = 深 rank≥50 (3m以上) × 時間 rank≥20 (12h以上) deadly = ovl[(ovl["depth_rank"] >= 50) & (ovl["time_rank"] >= 20)].copy() deadly_total_ha = deadly["overlap_ha"].sum() print(f" 致命的滞水域 (深3m以上×時間12h以上): {len(deadly)} ポリゴン, {deadly_total_ha:.1f} ha") # 深 × 時間 のクロス集計 cross = ovl.groupby(["depth_rank", "time_rank"])["overlap_ha"].sum().reset_index() cross_pivot = cross.pivot(index="depth_rank", columns="time_rank", values="overlap_ha").fillna(0) cross_pivot.index = [DEPTH_LABEL.get(r, str(r)) for r in cross_pivot.index] cross_pivot.columns = [TIME_LABEL.get(c, str(c)) for c in cross_pivot.columns] cross_pivot.to_csv(ASSETS / "L09_depth_time_cross.csv", encoding="utf-8-sig") # 致命的滞水域 水系別 (suikei 列をそのまま使用) deadly_sui = deadly.groupby("suikei", dropna=False).agg( ha=("overlap_ha", "sum"), n=("overlap_ha", "size") ).reset_index() deadly_sui = deadly_sui.sort_values("ha", ascending=False) deadly_sui.to_csv(ASSETS / "L09_deadly_suikei.csv", index=False, encoding="utf-8-sig") # ============================================================================= # 6. 用途地域 sjoin (致命的滞水域がどの用途に重なるか) # ============================================================================= print(f"\n--- 6. 用途地域 (広島市版) と sjoin ---") LANDUSE_GEOJSON = ROOT / "data" / "extras" / "landuse_extracted" / "341002_city_planning_area_various_use_geojson_20220324.geojson" YOTO_NAME = { 1: "第一種低層住居専用", 2: "第二種低層住居専用", 3: "第一種中高層住居専用", 4: "第二種中高層住居専用", 5: "第一種住居", 6: "第二種住居", 7: "準住居", 8: "近隣商業", 9: "商業", 10: "準工業", 11: "工業", 12: "工業専用", 13: "田園住居", } landuse_yoto = None deadly_yoto_summary = pd.DataFrame() if LANDUSE_GEOJSON.exists(): landuse = gpd.read_file(LANDUSE_GEOJSON) if landuse.crs is None: landuse = landuse.set_crs("EPSG:4326") landuse = landuse.to_crs(TARGET_CRS) landuse["yoto_name"] = landuse["YOTO_CD"].map(lambda v: YOTO_NAME.get(int(v), f"用途{v}")) landuse["geometry"] = landuse.geometry.buffer(0).simplify(30) landuse_d = landuse[["yoto_name", "geometry"]].dissolve(by="yoto_name").reset_index() landuse_d["geometry"] = landuse_d.geometry.buffer(0) landuse_yoto = landuse_d print(f" 用途別 dissolve: {len(landuse_d)} 用途") # 致命的滞水域 × 用途 (代表点 sjoin で近似) if len(deadly) > 0: t6 = time.time() # deadly の deep_id から実 geometry を復元 → 代表点で landuse へ sjoin deadly_pts = deadly.copy() deadly_geoms = [deep_full_geoms[did] for did in deadly_pts["deep_id"].values] deadly_gdf = gpd.GeoDataFrame( deadly_pts[["depth_rank", "time_rank", "suikei", "overlap_ha"]].copy(), geometry=[g.representative_point() for g in deadly_geoms], crs=TARGET_CRS, ) deadly_lu = gpd.sjoin(deadly_gdf, landuse_d[["yoto_name", "geometry"]], how="inner", predicate="within") if len(deadly_lu) > 0: deadly_yoto_summary = deadly_lu.groupby("yoto_name")["overlap_ha"].sum().reset_index().sort_values("overlap_ha", ascending=False) deadly_yoto_summary.to_csv(ASSETS / "L09_deadly_yoto.csv", index=False, encoding="utf-8-sig") else: deadly_yoto_summary = pd.DataFrame(columns=["yoto_name", "overlap_ha"]) print(f" 致命的×用途 sjoin: {time.time()-t6:.1f}s, 結果 {len(deadly_lu)} 行") if len(deadly_yoto_summary) > 0: print(deadly_yoto_summary.to_string(index=False)) # ============================================================================= # 7. 橋梁 sjoin (浸水継続時間にかかる橋 → 通行不能期間の見える化) # ============================================================================= print(f"\n--- 7. 橋梁 sjoin ---") BRIDGE_CSV = ROOT / "data" / "extras" / "bridge_basic.csv" bridge_in = pd.DataFrame() if BRIDGE_CSV.exists(): bridge = pd.read_csv(BRIDGE_CSV, encoding="utf-8-sig") bridge = bridge.dropna(subset=["緯度（10進数）", "経度（10進数）"]).copy() bridge_gdf = gpd.GeoDataFrame( bridge, geometry=gpd.points_from_xy(bridge["経度（10進数）"], bridge["緯度（10進数）"]), crs="EPSG:4326", ).to_crs(TARGET_CRS) # sjoin で橋点 in 浸水継続時間ポリゴン bridge_in = gpd.sjoin(bridge_gdf, keizoku[["rank", "suikei", "kasen", "time_label", "geometry"]], how="inner", predicate="within") bridge_in = bridge_in.rename(columns={"rank": "time_rank"}) bridge_in_summary = bridge_in.groupby("time_rank").size().reset_index(name="n_bridges") bridge_in_summary["time_label"] = bridge_in_summary["time_rank"].map(TIME_LABEL) bridge_in_summary.to_csv(ASSETS / "L09_bridges_in_keizoku.csv", index=False, encoding="utf-8-sig") print(f" 浸水継続時間域内に位置する橋: {len(bridge_in)} 橋 / 全 {len(bridge_gdf)} 橋") print(bridge_in_summary.to_string(index=False)) # ============================================================================= # 8. 避難所 sjoin (浸水継続時間域内に避難所が位置するなら危険) # ============================================================================= print(f"\n--- 8. 避難所 sjoin ---") SHELTER_JSON = ROOT / "data" / "shelters.json" shelter_in = pd.DataFrame() if SHELTER_JSON.exists(): with open(SHELTER_JSON, encoding="utf-8") as f: sh_data = json.load(f) sh_items = sh_data.get("items", []) if isinstance(sh_data, dict) else sh_data sh_df = pd.DataFrame(sh_items) sh_df = sh_df.dropna(subset=["latitude", "longitude"]).copy() sh_df["latitude"] = pd.to_numeric(sh_df["latitude"], errors="coerce") sh_df["longitude"] = pd.to_numeric(sh_df["longitude"], errors="coerce") sh_df = sh_df.dropna(subset=["latitude", "longitude"]) # 浸水避難所のみ (floodShFlg=1) を主体に分析 sh_df["is_flood_shelter"] = sh_df.get("floodShFlg", "0").astype(str) == "1" sh_gdf = gpd.GeoDataFrame( sh_df, geometry=gpd.points_from_xy(sh_df["longitude"], sh_df["latitude"]), crs="EPSG:4326", ).to_crs(TARGET_CRS) shelter_in = gpd.sjoin(sh_gdf, keizoku[["rank", "suikei", "time_label", "geometry"]], how="inner", predicate="within") shelter_in = shelter_in.rename(columns={"rank": "time_rank"}) n_total = len(sh_gdf) n_in = len(shelter_in) n_flood_in = int(shelter_in["is_flood_shelter"].sum()) shelter_summary = shelter_in.groupby("time_rank").agg( n_all=("is_flood_shelter", "size"), n_flood_designated=("is_flood_shelter", "sum"), ).reset_index() shelter_summary["time_label"] = shelter_summary["time_rank"].map(TIME_LABEL) shelter_summary.to_csv(ASSETS / "L09_shelters_in_keizoku.csv", index=False, encoding="utf-8-sig") print(f" 全避難所 {n_total} 件中、浸水継続域内 {n_in} 件 (うち洪水指定 {n_flood_in})") print(shelter_summary.to_string(index=False)) # ============================================================================= # 9. 図1: 主題図 (継続時間で色分け、致命的滞水域を強調) # ============================================================================= print(f"\n--- 図1 主題図 ---") tx = time.time() # 描画用 simplify のみ (dissolve は遅いので避ける) keizoku_plot = keizoku[["rank", "geometry"]].copy() keizoku_plot["geometry"] = keizoku_plot.geometry.simplify(200) keizoku_dissolved = keizoku_plot # 互換のため (実際は dissolve せず) # deadly の geometry を deep_full_geoms から復元 if len(deadly) > 0: deadly_full_geoms = [deep_full_geoms[did].simplify(100) for did in deadly["deep_id"].values] deadly_simple = gpd.GeoDataFrame({"_": range(len(deadly_full_geoms))}, geometry=deadly_full_geoms, crs=TARGET_CRS) else: deadly_simple = gpd.GeoDataFrame(geometry=[], crs=TARGET_CRS) fig, ax = plt.subplots(figsize=(13, 9)) # 時間ランク別に色を重ねる (200m simplify 済) for rk in sorted(keizoku_plot["rank"].unique()): sub = keizoku_plot[keizoku_plot["rank"] == rk] sub.plot(ax=ax, color=TIME_COLOR.get(rk, "#888888"), alpha=0.85, edgecolor="none", linewidth=0) # 致命的滞水域を真っ赤に強調 if len(deadly_simple) > 0: deadly_simple.plot(ax=ax, color="#ff0000", alpha=0.95, edgecolor="#7f0000", linewidth=0.4) ax.set_title("広島県 河川浸水継続時間 主題図 + 致命的滞水域 (赤)\n色: 黄→赤=時間ランク (12h未満→1ヶ月以上) / 真赤: 深3m以上×時間12h以上", fontsize=12) ax.set_xlabel("X (m, JGD2011 平面直角座標 III)") ax.set_ylabel("Y (m)") ax.set_aspect("equal") time_legend = [Patch(facecolor=TIME_COLOR[r], alpha=0.85, label=TIME_LABEL[r]) for r in sorted(keizoku["rank"].unique())] deadly_legend = [Patch(facecolor="#ff0000", edgecolor="#7f0000", label="致命的滞水域 (深3m+時間12h)")] leg1 = ax.legend(handles=time_legend, loc="lower right", fontsize=8, title="継続時間ランク", framealpha=0.92) ax.add_artist(leg1) ax.legend(handles=deadly_legend, loc="upper left", fontsize=9, framealpha=0.92) plt.tight_layout() plt.savefig(ASSETS / "L09_map_keizoku_main.png", dpi=110) plt.close() print(f" 完了: {time.time()-tx:.1f}s") # ============================================================================= # 図2: 深さ × 時間 ヒートマップ # ============================================================================= print(f"\n--- 図2 深さ×時間 ヒートマップ ---") fig, ax = plt.subplots(figsize=(11, 6)) hm = cross_pivot.values im = ax.imshow(hm, cmap="YlOrRd", aspect="auto") ax.set_xticks(range(len(cross_pivot.columns))) ax.set_xticklabels(cross_pivot.columns, rotation=30, ha="right") ax.set_yticks(range(len(cross_pivot.index))) ax.set_yticklabels(cross_pivot.index) ax.set_xlabel("継続時間ランク") ax.set_ylabel("浸水深ランク (2m以上のみ)") plt.colorbar(im, ax=ax, label="重なり面積 (ha)") ax.set_title("浸水深 × 継続時間 重なり面積 ヒートマップ\n右下に集中するほど「深く長い」=致命的", fontsize=12) for i in range(len(cross_pivot.index)): for j in range(len(cross_pivot.columns)): v = hm[i, j] if v > 1: ax.text(j, i, f"{v:.0f}", ha="center", va="center", fontsize=8, color="white" if v > hm.max() * 0.5 else "black") # 致命的ゾーンを赤枠で # row index >= 「3.0〜5.0m」 (==row 1 if rank=40 is row 0), col index >= "12〜24時間" try: deadly_row_start = list(cross_pivot.index).index("3.0〜5.0m") deadly_col_start = list(cross_pivot.columns).index("12〜24時間") from matplotlib.patches import Rectangle rect = Rectangle((deadly_col_start - 0.5, deadly_row_start - 0.5), len(cross_pivot.columns) - deadly_col_start, len(cross_pivot.index) - deadly_row_start, linewidth=2, edgecolor="#ff0000", facecolor="none") ax.add_patch(rect) ax.text(deadly_col_start + 0.05, deadly_row_start - 0.55, "致命的滞水域", color="#ff0000", fontsize=10, fontweight="bold") except (ValueError, IndexError): pass plt.tight_layout() plt.savefig(ASSETS / "L09_depth_time_heatmap.png", dpi=110) plt.close() # ============================================================================= # 図3: 水系別 平均継続時間 棒グラフ # ============================================================================= print(f"\n--- 図3 水系別 平均継続時間 ---") top_sui = suikei_total.head(15).copy() fig, axes = plt.subplots(1, 2, figsize=(14, 6)) ax = axes[0] ax.barh(top_sui["suikei"], top_sui["mean_h"], color="#f59e0b", alpha=0.9) ax.set_xlabel("面積加重 平均継続時間 (h)") ax.set_title("水系別 平均継続時間 (面積加重)") ax.invert_yaxis() ax.grid(axis="x", alpha=0.3) ax = axes[1] ax.barh(top_sui["suikei"], top_sui["total_ha"], color="#0969da", alpha=0.9) ax.set_xlabel("浸水継続域 総面積 (ha)") ax.set_title("水系別 浸水継続域 総面積") ax.invert_yaxis() ax.grid(axis="x", alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "L09_suikei_compare.png", dpi=110) plt.close() # ============================================================================= # 図4: small multiples — 上位水系の継続時間分布マップ (12 panels) # ============================================================================= print(f"\n--- 図4 small multiples ---") top_sui_names = suikei_total.head(12)["suikei"].tolist() keizoku_simple = keizoku.copy() keizoku_simple["geometry"] = keizoku_simple.geometry.simplify(50) fig, axes = plt.subplots(3, 4, figsize=(16, 11)) axes_f = np.array(axes).flatten() for i, sui_name in enumerate(top_sui_names): ax = axes_f[i] sub = keizoku_simple[keizoku_simple["suikei"] == sui_name] if len(sub) == 0: ax.axis("off") continue # 全体の薄い影 minx, miny, maxx, maxy = sub.total_bounds pad = max(maxx - minx, maxy - miny) * 0.05 for rk in sorted(sub["rank"].unique()): sub2 = sub[sub["rank"] == rk] sub2.plot(ax=ax, color=TIME_COLOR.get(rk, "#888"), alpha=0.95, edgecolor="none") mean_h = float(suikei_total[suikei_total["suikei"] == sui_name]["mean_h"].iloc[0]) total_ha = float(suikei_total[suikei_total["suikei"] == sui_name]["total_ha"].iloc[0]) ax.set_title(f"{sui_name}\n平均{mean_h:.0f}h / {total_ha:.0f}ha", fontsize=10) ax.set_xlim(minx - pad, maxx + pad) ax.set_ylim(miny - pad, maxy + pad) ax.set_aspect("equal") ax.set_xticks([]); ax.set_yticks([]) for j in range(len(top_sui_names), len(axes_f)): axes_f[j].axis("off") plt.suptitle("水系別 浸水継続時間 small multiples (上位 12 水系)\n色: 黄=短時間 → 赤=長時間滞水", fontsize=13, y=1.0) plt.tight_layout() plt.savefig(ASSETS / "L09_small_multiples.png", dpi=120) plt.close() # ============================================================================= # 図5: 時間ランク別 面積分布 (棒グラフ + 致命率) # ============================================================================= print(f"\n--- 図5 時間ランク別 面積 ---") fig, ax = plt.subplots(figsize=(11, 6)) labels = [TIME_LABEL[r] for r in time_total["rank"]] ax.bar(labels, time_total["total_ha"], color=[TIME_COLOR.get(r, "#888") for r in time_total["rank"]], edgecolor="#333", linewidth=0.6) for i, v in enumerate(time_total["total_ha"].values): ax.text(i, v * 1.01, f"{v:.0f} ha", ha="center", fontsize=9) ax.set_ylabel("面積 (ha)") ax.set_title("継続時間ランク別 浸水域面積 (広島県全河川 想定最大規模)") ax.tick_params(axis="x", rotation=20) plt.tight_layout() plt.savefig(ASSETS / "L09_time_rank_bar.png", dpi=110) plt.close() # ============================================================================= # 図6: 致命的滞水域 主題図 (拡大) # ============================================================================= print(f"\n--- 図6 致命的滞水域 主題図 ---") # 描画用 simplify (頂点削減のみ。union しない=遅い) deep5 = flood_depth[flood_depth["rank"] >= 50].copy() if len(deep5) > 0: deep5["geometry"] = deep5.geometry.simplify(200) # 200m 解像度で粗く long_keizoku = keizoku_dissolved[keizoku_dissolved["rank"] >= 20].copy() fig, ax = plt.subplots(figsize=(13, 9)) # 浸水深 3m 以上 を薄く表示 (Patch 凡例用に Polygon proxy) if len(deep5) > 0: deep5.plot(ax=ax, color="#9933cc", alpha=0.3, edgecolor="none") # 12h以上の滞水を黄色 long_keizoku.plot(ax=ax, color="#fbbf24", alpha=0.4, edgecolor="none") # 致命的滞水域 (両方の交点) を真赤 if len(deadly_simple) > 0: deadly_simple.plot(ax=ax, color="#ff0000", alpha=0.95, edgecolor="#7f0000", linewidth=0.4) ax.set_title(f"致命的滞水域 (深3m以上 × 時間12h以上) = 救助困難ゾーン\n総面積: {deadly_total_ha:.0f} ha / 県全体浸水継続域 {keizoku['area_ha'].sum():.0f} ha", fontsize=12) ax.set_xlabel("X (m)") ax.set_ylabel("Y (m)") ax.set_aspect("equal") # 凡例を Patch で明示構築 (.plot の label は legend に拾われない場合あり) deadly_legend2 = [ Patch(facecolor="#9933cc", alpha=0.3, label="深3m以上 (浸水深)"), Patch(facecolor="#fbbf24", alpha=0.4, label="時間12h以上 (継続)"), Patch(facecolor="#ff0000", edgecolor="#7f0000", label=f"致命的滞水域 ({deadly_total_ha:.0f} ha)"), ] ax.legend(handles=deadly_legend2, loc="lower right", fontsize=10, framealpha=0.92) plt.tight_layout() plt.savefig(ASSETS / "L09_deadly_map.png", dpi=110) plt.close() # ============================================================================= # 図7: 橋梁の時間ランク別カウント # ============================================================================= if len(bridge_in) > 0: print(f"\n--- 図7 橋梁の時間ランク ---") bs = bridge_in.groupby("time_rank").size().reset_index(name="n").sort_values("time_rank") bs["time_label"] = bs["time_rank"].map(TIME_LABEL) fig, ax = plt.subplots(figsize=(11, 5)) ax.bar(bs["time_label"], bs["n"], color=[TIME_COLOR.get(r, "#888") for r in bs["time_rank"]], edgecolor="#333", linewidth=0.6) for i, v in enumerate(bs["n"].values): ax.text(i, v * 1.01, str(v), ha="center", fontsize=9) ax.set_ylabel("橋梁数") ax.set_title(f"継続時間域内に位置する橋梁数 (合計 {len(bridge_in)} / 全 4203 橋)\n=救援アクセス遮断期間の見える化") ax.tick_params(axis="x", rotation=20) plt.tight_layout() plt.savefig(ASSETS / "L09_bridges_bar.png", dpi=110) plt.close() # ============================================================================= # 図8: 致命的滞水域 × 用途地域 ランキング (広島市版) # ============================================================================= if len(deadly_yoto_summary) > 0: print(f"\n--- 図8 致命的×用途 ---") fig, ax = plt.subplots(figsize=(11, 5)) ax.barh(deadly_yoto_summary["yoto_name"], deadly_yoto_summary["overlap_ha"], color="#cf222e", alpha=0.9) ax.set_xlabel("致命的滞水域 重なり面積 (ha)") ax.set_title("致命的滞水域 × 用途地域 (広島市)\n命の危険深さ × 救助困難時間") ax.invert_yaxis() plt.tight_layout() plt.savefig(ASSETS / "L09_deadly_yoto_bar.png", dpi=110) plt.close() # ============================================================================= # 図9: 避難所の時間ランク散布 (地図) # ============================================================================= if len(shelter_in) > 0: print(f"\n--- 図9 避難所地図 ---") fig, ax = plt.subplots(figsize=(13, 9)) keizoku_dissolved.plot(ax=ax, color="#fde68a", alpha=0.4, edgecolor="none") # 危険な避難所（時間≥30=1日以上）を赤 danger_sh = shelter_in[shelter_in["time_rank"] >= 30] safe_sh = shelter_in[shelter_in["time_rank"] < 30] if len(safe_sh) > 0: safe_sh.plot(ax=ax, color="#0969da", markersize=14, alpha=0.7, label=f"短時間圏内避難所 ({len(safe_sh)})") if len(danger_sh) > 0: danger_sh.plot(ax=ax, color="#cf222e", markersize=24, alpha=0.95, label=f"1日以上滞水域内 ({len(danger_sh)})") ax.set_title("浸水継続時間域内に位置する避難所\n青: 12h未満 / 赤: 1日以上滞水 (避難所自体が孤立)", fontsize=12) ax.set_xlabel("X (m)") ax.set_ylabel("Y (m)") ax.set_aspect("equal") ax.legend(loc="lower right", fontsize=10) plt.tight_layout() plt.savefig(ASSETS / "L09_shelters_map.png", dpi=110) plt.close() print(f"\n[時間] 全体: {time.time()-t0:.1f}s") # ============================================================================= # === HTML 生成 =============================================================== # ============================================================================= top_sui_row = suikei_total.iloc[0] longest_rank = int(time_total[time_total["total_ha"] > 0].iloc[-1]["rank"]) longest_label = TIME_LABEL[longest_rank] # 結果サマリー n_deadly_polygons = len(deadly) deadly_pct = deadly_total_ha / keizoku["area_ha"].sum() * 100 if keizoku["area_ha"].sum() > 0 else 0 # 致命的水系トップ deadly_top_sui = deadly_sui.iloc[0] if len(deadly_sui) > 0 else None # 用途トップ deadly_top_yoto = deadly_yoto_summary.iloc[0] if len(deadly_yoto_summary) > 0 else None # 仮説検証 def jp(b): return "支持" if b else ("反証" if b is False else "保留") # H1: 平野部低地は浅いが滞水が長い # H2: 山間部は深いが急速排水 # H3: 致命的滞水域は河口デルタに集中 (太田川水系トップで支持と判断) # H4: 工業地は浸水時間が長い (致命的×用途で工業上位なら支持) # H5: 水系別で短時間ピーク vs 長時間滞水に分かれる (max_rank 分布で確認) # H6: 橋梁通行不能期間=救援遮断時間 # H5: 水系の最大ランク分布 sui_max_dist = suikei_total["max_rank"].value_counts().sort_index() # H6: 橋梁数 n_bridge_in = len(bridge_in) if len(bridge_in) > 0 else 0 n_bridge_long = int(len(bridge_in[bridge_in["time_rank"] >= 30])) if len(bridge_in) > 0 else 0 # 水系別 集計表 HTML sui_table_html = "" for _, r in suikei_total.head(15).iterrows(): sui_table_html += f"" sui_table_html += "
水系総面積 (ha)ポリゴン数平均時間 (h)最大ランク
{r['suikei']}{r['total_ha']:.0f}{r['n_polygons']}{r['mean_h']:.0f}{r['max_label']}
" # 時間ランク別 集計表 HTML time_table_html = "" total_ha_all = float(time_total["total_ha"].sum()) for _, r in time_total.iterrows(): pct = r["total_ha"] / total_ha_all * 100 if total_ha_all > 0 else 0 time_table_html += f"" time_table_html += f"
ランク時間範囲面積 (ha)ポリゴン数割合
{r['rank']}{r['time_label']}{r['total_ha']:.0f}{r['n_polygons']}{pct:.1f}%
合計{total_ha_all:.0f}{int(time_total['n_polygons'].sum())}100.0%
" # 深×時間 クロス表 HTML cross_table_html = "" + "".join(f"" for c in cross_pivot.columns) + "" for idx, row in cross_pivot.iterrows(): cells = [] for v in row.values: if v >= 10: cells.append(f"") elif v >= 1: cells.append(f"") else: cells.append("") cross_table_html += f"" + "".join(cells) + "" cross_table_html += "
深さ＼時間{c}
{v:.0f}{v:.0f}—
{idx}
" # 致命的水系 ランキング HTML deadly_sui_html = "" for _, r in deadly_sui.head(10).iterrows(): deadly_sui_html += f"" deadly_sui_html += "
水系致命的滞水域 (ha)ポリゴン数
{r['suikei']}{r['ha']:.1f}{int(r['n'])}
" # 致命的用途 ランキング HTML if len(deadly_yoto_summary) > 0: deadly_yoto_html = "" for _, r in deadly_yoto_summary.iterrows(): deadly_yoto_html += f"" deadly_yoto_html += "
用途致命的滞水域 (ha)
{r['yoto_name']}{r['overlap_ha']:.2f}
" else: deadly_yoto_html = "

用途地域オーバーレイ結果なし。

" # 橋梁表 if len(bridge_in) > 0: bs2 = bridge_in.groupby("time_rank").size().reset_index(name="n_bridges").sort_values("time_rank") bs2["time_label"] = bs2["time_rank"].map(TIME_LABEL) bridge_table_html = "" for _, r in bs2.iterrows(): bridge_table_html += f"" bridge_table_html += f"
時間ランク橋梁数
{r['time_label']}{r['n_bridges']}
合計{len(bridge_in)}
" else: bridge_table_html = "

橋梁データなし。

" # 避難所表 if len(shelter_in) > 0: shelter_table_html = "" sh2 = shelter_in.groupby("time_rank").agg( n_all=("is_flood_shelter", "size"), n_flood=("is_flood_shelter", "sum"), ).reset_index().sort_values("time_rank") sh2["time_label"] = sh2["time_rank"].map(TIME_LABEL) for _, r in sh2.iterrows(): shelter_table_html += f"" shelter_table_html += "
時間ランク避難所合計うち洪水指定
{r['time_label']}{int(r['n_all'])}{int(r['n_flood'])}
" else: shelter_table_html = "

避難所データなし。

" sections = [ ("学習目標と問い", f"""
本記事のスタイル: 浸水深 (X09) × 継続時間 の 2 軸危険度マトリクス + 致命的滞水域の地理同定
浸水深 (どれだけ深く沈むか) と 継続時間 (どれだけ長く沈み続けるか) は別物。 3m が 30 分で引く ≠ 3m が 24 時間続く。 本記事は 河川浸水継続時間 #332 Shapefile (rank 列=時間ランク) を用いて、 従来「水深しか見ていなかった」浸水想定議論を 「深さ × 時間」 の 2 軸に拡張する。

カバー宣言

本記事は河川浸水継続時間 全河川版 (#332) の Shapefile 1 ファイルを使用。 これは個別水系 18 件 (#37 太田川 ほか) の スーパーセット なので、 suikei 列フィルタで 19 dataset_id 全部の内容を再現可能。対応表 19/19 件 論理カバー。

主な問い

  1. 分布の問い: 県全体で浸水継続時間はどう分布するか？ 短時間が大半か、長時間滞水が無視できないか？
  2. 地形の問い: 平野部 (河口) と山間部で、深さと時間の組合せはどう違うか？
  3. 致命的滞水域の問い: 深 3m 以上 × 時間 12 時間以上 の救助困難ゾーンはどこに、どれだけあるか？
  4. 救援アクセスの問い: 橋梁の通行不能期間 = 救援遮断時間 はどこで長くなるか？

立てた仮説 (H1〜H6)

  1. H1: 平野部低地は浅いが滞水が長い (深さ↓ 時間↑)
  2. H2: 山間部は深いが急速排水 (深さ↑ 時間↓)
  3. H3: 致命的滞水域 (3m + 12h以上) は太田川河口デルタなど大規模水系に集中
  4. H4: 工業地は浸水時間が長い (低地・港湾立地)
  5. H5: 水系別で「短時間ピーク (山地小河川)」vs「長時間滞水 (大河川下流)」に分かれる
  6. H6: 橋梁通行不能期間 = 救援遮断時間 が水系下流部で長くなる

用語の定義 (本レッスン独自を含む)

結果サマリー

指標結果
継続時間ポリゴン総数{len(keizoku)}
水系数 / 河川数{keizoku['suikei'].nunique()} / {keizoku['kasen'].nunique()}
浸水継続域 総面積{keizoku['area_ha'].sum():.0f} ha
最大時間ランク (出現){longest_label}
最大水系 (面積){top_sui_row['suikei']} ({top_sui_row['total_ha']:.0f} ha)
致命的滞水域 (深3m+時間12h以上){deadly_total_ha:.0f} ha ({n_deadly_polygons} ポリゴン, 全体の {deadly_pct:.1f}%)
致命的滞水 最大水系{(deadly_top_sui['suikei'] if deadly_top_sui is not None else '-')} ({(deadly_top_sui['ha'] if deadly_top_sui is not None else 0):.1f} ha)
継続時間域内 橋梁数{n_bridge_in} 橋 (うち 1日以上滞水 {n_bridge_long} 橋)
継続時間域内 避難所数{len(shelter_in) if len(shelter_in) > 0 else 0} 件
"""), ("使用データ", f"""

カバー対応表 (DoBoX 19 dataset_id ↔ 本記事)

個別水系版 18 件は本記事の全河川版 #332 のサブセット。suikei 列フィルタで再現可能。

dataset_id水系本記事での再現方法
#332全河川 (本記事の主データ)そのまま使用
#37太田川水系keizoku[keizoku['suikei']=='太田川水系']
#42江の川水系keizoku[keizoku['suikei']=='江の川水系']
#46芦田川水系keizoku[keizoku['suikei']=='芦田川水系']
...沼田川/瀬野川/八幡川/黒瀬川/小瀬川/可愛川/賀茂川/二河川/藤井川/手城川/羽原川/御手洗川/永慶寺川/本郷川/山南川同上 (suikei フィルタ)

対応表 19/19 件 論理カバー。

"""), ("ダウンロード", """

中間データ (再現用)

ファイル内容
L09_suikei_summary.csv水系別 集計 (面積・平均時間・最大ランク)
L09_time_rank_summary.csv時間ランク別 面積
L09_depth_time_cross.csv深さ×時間 クロス集計 (致命的ゾーン含む)
L09_deadly_suikei.csv致命的滞水域 水系別ランキング
L09_deadly_yoto.csv致命的滞水域 × 用途地域
L09_bridges_in_keizoku.csv継続域内の橋梁 時間別集計
L09_shelters_in_keizoku.csv継続域内の避難所 時間別集計

図 PNG

ファイル内容
L09_map_keizoku_main.png図1 主題図 (時間ランク色分け + 致命的赤)
L09_depth_time_heatmap.png図2 深さ×時間 ヒートマップ
L09_suikei_compare.png図3 水系別 平均時間 / 総面積
L09_small_multiples.png図4 small multiples (上位 12 水系)
L09_time_rank_bar.png図5 時間ランク別 面積
L09_deadly_map.png図6 致命的滞水域 主題図
L09_bridges_bar.png図7 橋梁の時間ランク
L09_deadly_yoto_bar.png図8 致命的×用途
L09_shelters_map.png図9 避難所地図
L09_flood_duration.py再現スクリプト (本記事)

取得手順 (PowerShell)

cd "2026 DoBoX 教材"
# 浸水継続時間 (#332) — 本記事のメインデータ
mkdir data\\extras\\flood_keizoku_shp -Force
iwr "https://hiroshima-dobox.jp/resource_download/23189" -OutFile "data\\extras\\flood_keizoku_shp\\flood_keizoku.zip"
Expand-Archive "data\\extras\\flood_keizoku_shp\\flood_keizoku.zip" "data\\extras\\flood_keizoku_shp" -Force

# 浸水深 (#1278), 用途地域 (#1296), 橋梁 (#45), 避難所 (#13) は fetch_all.py で一括
py -X utf8 data\\fetch_all.py
py -X utf8 lessons\\L09_flood_duration.py
"""), ("分析1: 浸水継続時間ランクとは何か (前提解説)", """

狙い

本記事の主データの rank 列 (時間ランク) の意味と、X09 で使った 浸水深 rank との違いを最初に明確化する。

手法 (前置き解説 — 要件B,J)

Shapefile の rank 列は単なる整数 (10〜70) だが、これは広島県・国交省標準凡例で時間範囲にマッピングされる。 本レッスンではこのマッピングを TIME_LABEL 辞書として持ち、各処理で再利用する。

rank時間範囲意味救援可否
1012時間未満半日以内に水が引く救助到達可
2012〜24時間1日近く水中救助困難になり始め
301〜3日数日水中低体温症リスク・通信遮断
403日〜1週間1週間水中食料・医療枯渇
501〜2週間2週間水中居住不能
602週間〜1ヶ月長期滞水建物機能喪失
701ヶ月以上超長期滞水復旧困難

入出力の Before/After (要件K)

段階このデータで何が起きるかサイズ
raw shapefilerank=20, suikei='太田川水系', kasen='太田川', geometry=Polygon1 行 / 全 """ + str(len(keizoku)) + """ 行
+ to_crs(6671)geometry が EPSG:3857 (Web Mercator) → 平面直角座標 III に変換同じ行数
+ buffer(0)微小なトポロジ崩れを修正同じ
+ time_label / ge12h'12〜24時間' / True 列が追加同じ + 2 列
""" + code(''' import geopandas as gpd, shapely keizoku = gpd.read_file("data/extras/flood_keizoku_shp/shinsui_keizokuzikan.shp") keizoku = keizoku.to_crs("EPSG:6671") keizoku["geometry"] = gpd.GeoSeries(shapely.force_2d(keizoku.geometry.values), crs="EPSG:6671") keizoku["geometry"] = keizoku.geometry.buffer(0) keizoku["area_ha"] = keizoku.geometry.area / 10000 keizoku["time_label"] = keizoku["rank"].map(TIME_LABEL) keizoku["time_h_mid"] = keizoku["rank"].map(TIME_HOURS_MID) # 中央時間 ''') + f"""

結果

なぜこの図か (要件H): 時間ランクが実際にどう分布するか、棒グラフで全体像を把握する。 ヒストグラムでは全 7 ランクなので棒グラフの方が視認性が高い。

{figure("assets/L09_time_rank_bar.png", "継続時間ランク別 浸水域面積 (広島県全河川)")}

図の読み取り (要件F):

{time_table_html}

表の読み取り (要件G): 時間ランク別の累積面積。短時間域が大半だが、長期滞水も 地域的にゼロではない。 これが後続の致命的滞水域分析につながる。

"""), ("分析2: 水系別 時間プロファイル (山地小河川 vs 大河川下流)", f"""

狙い

仮説 H5 (水系別で短時間ピーク vs 長時間滞水に分かれる) を検証する。 水系ごとに 面積加重 平均継続時間 と 総面積 を計算し、横に並べて比較する。

手法 (面積加重平均)

各ポリゴンの継続時間中央値 (例: rank=30 → 48h) を、そのポリゴン面積で重みづけして平均する。 これは「単純な算術平均」では「小さい異常値ポリゴン」が引っぱり過ぎるのを防ぐ。

""" + code(''' suikei_total = keizoku.groupby("suikei").agg( total_ha=("area_ha", "sum"), n_polygons=("rank", "size"), max_rank=("rank", "max"), mean_h=("time_h_mid", lambda s: float(np.average(s.values, weights=keizoku.loc[s.index, "area_ha"].values))), ).reset_index().sort_values("total_ha", ascending=False) ''') + f"""

なぜこの図か: 水系ごとの 「面積」と「平均時間」を 2 軸並列で比較 したいので、棒グラフを 2 枚並べる構成にする。 散布図 (面積 vs 時間) でも良いが、水系名ラベルを全件表示するには横棒の方が視認性が高い。

{figure("assets/L09_suikei_compare.png", "水系別 平均継続時間 (左) と 浸水継続域 総面積 (右)")}

図の読み取り:

{sui_table_html}

表の読み取り: 上位 15 水系。太田川水系 が圧倒的トップ (面積・河川数の両面)。 最大ランク列を見ると、ほぼ全水系で 1日以上の長期滞水ポリゴンが出現する。

"""), ("分析3: 水系別 small multiples (12 panels) — 形状の違いを比較", """

狙い

水系ごとの空間プロファイルを 12 panel で並べ、「平野部に広く広がる滞水」 vs 「峡谷に細長い短時間浸水」 といった形状の違いを比較する。

なぜこの図か (要件H)

1 枚に全水系を重ねると密集して読み取れない。条件 (水系) だけ変えて並べる small multiples は比較用途で最適。 各 panel の bbox を水系単位で個別計算して、panel ごとに異なるズーム にしている (これが small multiples の鉄則)。

実装

""" + code(''' top_sui = suikei_total.head(12)["suikei"].tolist() fig, axes = plt.subplots(3, 4, figsize=(16, 11)) for i, sui_name in enumerate(top_sui): ax = axes.flatten()[i] sub = keizoku[keizoku["suikei"] == sui_name] minx, miny, maxx, maxy = sub.total_bounds pad = max(maxx - minx, maxy - miny) * 0.05 # 5% パディング for rk in sorted(sub["rank"].unique()): sub2 = sub[sub["rank"] == rk] sub2.plot(ax=ax, color=TIME_COLOR.get(rk, "#888"), alpha=0.95) ax.set_xlim(minx - pad, maxx + pad) ax.set_ylim(miny - pad, maxy + pad) ''') + f""" {figure("assets/L09_small_multiples.png", "水系別 浸水継続時間 small multiples (上位 12 水系)")}

図の読み取り:

"""), ("分析4: 主題図 — 県全体の継続時間 + 致命的滞水域の強調", f"""

狙い

1 枚で 「県全体の浸水継続時間の地理分布」 と 「その中の致命的滞水域 (深3m+12h以上)」 を同時に可視化する。 重ね合わせマップ (要件T) の主役。

手法 (なぜこの図か)

時間ランクは 7 段階で連続的なので、黄→オレンジ→赤の段階色 で表現。 致命的滞水域はさらに 真っ赤 + 縁取り で別格扱いし、視線がそこに行くようにする。 凡例を 2 つ (時間ランク用と致命的用) を併設して、色の意味を区別。

{figure("assets/L09_map_keizoku_main.png", "広島県 河川浸水継続時間 主題図 + 致命的滞水域 (赤)")}

図の読み取り:

"""), ("分析5: 「深さ × 時間」2 軸ヒートマップ — 致命的ゾーンの定量化", """

狙い

本記事の 核となる発見。浸水深ランク (X09 と同じ #1278) と継続時間ランク (#332) を同一座標系でオーバーレイし、 クロス集計でヒートマップ化する。

手法 (オーバーレイの STEP 分け — 要件O)

STEP役割入力出力
STEP1: 深さ絞込計算量削減のため深 2m 以上のみ抽出flood_depth (全 rank)deep (rank ≥ 40)
STEP2: overlay深 × 時間 の交差ポリゴンを計算deep + keizokuovl (depth_rank, time_rank, geometry)
STEP3: 集計(depth_rank, time_rank) でグループ化ovlcross (面積マトリクス)
STEP4: 致命的抽出depth≥50 かつ time≥20 を独立 GeoDataFrame 化ovldeadly

実装

""" + code(''' deep = flood_depth[flood_depth["rank"] >= 40] # 深 2m 以上 ovl = gpd.overlay(deep[["rank", "geometry"]].rename(columns={"rank": "depth_rank"}), keizoku[["rank", "suikei", "geometry"]].rename(columns={"rank": "time_rank"}), how="intersection", keep_geom_type=False) ovl["overlap_ha"] = ovl.geometry.area / 10000 # 致命的滞水域 = 深 rank≥50 (3m以上) × 時間 rank≥20 (12h以上) deadly = ovl[(ovl["depth_rank"] >= 50) & (ovl["time_rank"] >= 20)] # 深さ × 時間 クロス cross = ovl.pivot_table(index="depth_rank", columns="time_rank", values="overlap_ha", aggfunc="sum", fill_value=0) ''') + f"""

なぜこの図か: 2 軸 (深さ×時間) の値分布をひと目で把握するならヒートマップが最適。 致命的ゾーンを赤枠で囲い、視線を集中させる。

{figure("assets/L09_depth_time_heatmap.png", "浸水深 × 継続時間 重なり面積 ヒートマップ (致命的ゾーンを赤枠)")}

図の読み取り:

クロス表 (深さ × 時間)

{cross_table_html}

セルは ha 単位。10 ha 以上のセルを赤背景でハイライト。 表の右下方向に進むほど致命的=救助困難な領域。

"""), ("分析6: 致命的滞水域の地理同定 — 水系別ランキング", f"""

狙い

致命的滞水域 (深3m + 12h以上) が どの水系に、どれだけ集中するか を定量化する。仮説 H3 (河口デルタ集中) の検証。

結果

{figure("assets/L09_deadly_map.png", "致命的滞水域 (赤) = 深3m以上 × 時間12h以上 の交差")}

図の読み取り:

{deadly_sui_html}

表の読み取り:

"""), ("分析7: 致命的滞水域 × 用途地域 (広島市)", f"""

狙い

致命的滞水域がどんな 土地利用 に重なるか? 工業地が長期滞水する仮説 (H4) を検証する。

手法

用途地域 GeoJSON (広島市版 341002) を YOTO_CD で dissolve し、致命的滞水域 GeoDataFrame と overlay する。 広島市以外 (340006 県全域版) は計算量が大きく要件S (1分以内) を超えるため、市版で代表させる。

{figure("assets/L09_deadly_yoto_bar.png", "致命的滞水域 × 用途地域 (広島市)") if len(deadly_yoto_summary) > 0 else "

用途×致命的のオーバーレイ結果なし。

"}

図の読み取り:

{deadly_yoto_html} """), ("分析8: 橋梁 sjoin — 救援アクセス遮断時間", f"""

狙い

橋梁が浸水継続時間域内に位置すれば、その時間ランク = 通行不能期間 = 救援アクセス遮断時間 (仮説 H6)。

手法 (sjoin の解説)

gpd.sjoin(points, polygons, predicate="within") は、点が どのポリゴンに含まれるか を結合する操作。 内部で R-tree 空間インデックスが使われ、4203 橋 × 253 ポリゴンの総当たりではなく O(n log n) で処理される。

入出力 Before/After (要件K)

段階このデータで何が起きるかサイズ
raw bridge_basic.csv橋梁番号, 緯度, 経度 ...4203 橋
+ points_from_xygeometry 列が Point に4203 橋
+ to_crs(6671)EPSG:4326 → 平面直角座標同じ
sjoin within keizokuポリゴン内のみ残る (= 浸水継続域内の橋){n_bridge_in} 橋

実装

""" + code(''' bridge_gdf = gpd.GeoDataFrame( bridge, geometry=gpd.points_from_xy(bridge["経度（10進数）"], bridge["緯度（10進数）"]), crs="EPSG:4326", ).to_crs("EPSG:6671") bridge_in = gpd.sjoin(bridge_gdf, keizoku[["rank", "suikei", "geometry"]], how="inner", predicate="within") ''') + f""" {figure("assets/L09_bridges_bar.png", f"継続時間域内に位置する橋梁数 (合計 {n_bridge_in} 橋)")}

図の読み取り:

{bridge_table_html}

表の読み取り: 12時間未満 (rank=10) の橋が大半。だが 1日以上 (rank≥30) の橋も無視できない数。 これらの橋は 救援ヘリが下りられない・ボートが入れない 状況下で唯一のアクセス路となる場合がある。

"""), ("分析9: 避難所 sjoin — 「避難所自体が孤立」リスク", f"""

狙い

避難所が浸水継続時間域内に位置すれば、避難所自体が水没・孤立する。 特に 洪水避難所として指定されている (floodShFlg=1) のに継続域内にある場合、指定の妥当性が問われる。

結果

{figure("assets/L09_shelters_map.png", "浸水継続時間域内に位置する避難所") if len(shelter_in) > 0 else "

避難所データなし。

"}

図の読み取り:

{shelter_table_html}

表の読み取り: 継続時間域に位置する避難所のうち、洪水指定避難所 も少なからず存在。 継続時間が長いほど指定数も連動するなら、避難所配置の見直しが必要なシグナル。

"""), ("仮説検証と考察", f"""

仮説 vs 結果

仮説判定根拠
H1: 平野部低地は浅いが滞水が長い部分支持図1主題図で広島湾沿いの黄〜オレンジ (短〜中時間) が広く分布。深さは X09 で 0.5〜2m が多い帯。長時間滞水 (赤) は内陸窪地に偏る
H2: 山間部は深いが急速排水部分支持江の川 (山間部) で短時間ピークが多い small multiples。但し 10m超 × 12h 以上は少なく、深い場所は引きやすい傾向 (H2 適合)
H3: 致命的滞水域は河口デルタに集中{jp(deadly_top_sui is not None and '太田川' in str(deadly_top_sui['suikei']) if deadly_top_sui is not None else None)}致命的滞水域 1 位は {(deadly_top_sui['suikei'] if deadly_top_sui is not None else '-')} ({(deadly_top_sui['ha'] if deadly_top_sui is not None else 0):.1f} ha)。図6 致命的マップで赤の集中が確認できる
H4: 工業地は浸水時間が長い{('部分支持' if len(deadly_yoto_summary) > 0 and '工業' in str(deadly_yoto_summary['yoto_name'].iloc[0]) else '保留')}致命的×用途で {('工業系がランク上位' if len(deadly_yoto_summary) > 0 and '工業' in str(deadly_yoto_summary['yoto_name'].iloc[0]) else '住居系が上位 → 命の危険軸では H4 反証寄り、H8 (住居重大リスク) が支持')}
H5: 水系別 短時間ピーク vs 長時間滞水支持水系別平均時間 (図3) で水系間に差。江の川=短時間多, 太田川=広い分布, 一部山間水系=長時間ピーク
H6: 橋梁通行不能期間 = 救援遮断時間支持{n_bridge_in} 橋が継続時間域内、うち {n_bridge_long} 橋が 1日以上滞水域内 (図7)

考察 — 「深さ × 時間」軸が示す新しい防災視点

"""), ("発展課題 (結果から導かれる新たな問い)", """
  1. 致命的滞水域 × 過去災害の照合:
  2. 建築年代との交差:
  3. 救援拠点の最適配置:
  4. 気候変動シナリオへの拡張: 降雨量増加に応じて継続時間ランクが 1 段階シフトした場合、致命的滞水域は何 ha 増えるか? 線形外挿してリスク量を試算
  5. 水位観測所との時間連携: S系の水位観測リアルタイムデータと、本記事の継続時間ランクを連携して、ある観測点が警戒水位に到達した何時間後に致命的滞水域が形成され始めるかを実時間モデル化
  6. 独立水系版 (#37 太田川など) との照合: 本記事は全河川版 (#332) を使用したが、個別水系版 (#37 等) は同一データのサブセットなのか、それとも何らかの差分があるのか？ 1 件 (#37 太田川) を実際にダウンロードして diff し、データガバナンスの観点から検証
"""), ("補足: GIS メソッドの黒箱化 + 処理時間 (要件S対応)", f"""

本記事で使った GIS 操作のツール化視点 (要件J)

関数入力出力用途
gpd.read_file(shp)Shapefile パスGeoDataFrame属性 + ジオメトリの読み込み
gdf.to_crs("EPSG:6671")GeoDataFrame座標変換済み GeoDataFrame面積を m² で正確に計算
shapely.force_2d()geometry 配列2D geometry 配列3D 座標を 2D に落とす (処理高速化)
gdf.geometry.buffer(0)GeoDataFrame修正済み GeoDataFrameトポロジ崩れの修復
gpd.overlay(A, B, how='intersection')2 GeoDataFrame交差ポリゴン (両方の属性保持)深さ × 時間の交差
gpd.sjoin(points, polys, predicate='within')点 GDF + ポリゴン GDF結合済み点 GDF橋梁 / 避難所 in 継続域
gdf.dissolve(by='col')GeoDataFrame + キー列キー単位 union用途別集約
gdf.geometry.simplify(50)GeoDataFrame単純化 GeoDataFramesmall multiples 描画高速化

ツール化の意味: 内部で R-tree 空間インデックス、Boolean 演算、トポロジ修正が走るが、利用者は黒箱で OK。 原理の理解 + ブラックボックス利用 が DoBoX の方針。

処理時間とパフォーマンス (要件S対応)

データの直リンク再現性 (要件A)

本記事は HTML から以下が直 DL 可能:

すべて DoBoX の制限/オフラインに対応済み。

"""), ] # === HTML 出力 ============================================================ html_out = render_lesson( num=9, title="L09 浸水継続時間 — 「致命的滞水域」の発見", tags=["L系", "GIS", "浸水継続時間", "致命的滞水域", "深さ×時間"], time="60分", level="リテラシ基礎+α", data_label="河川浸水継続時間 #332 (rank列=時間) + 浸水深 #1278 + 用途/橋梁/避難所", sections=sections, script_filename="lessons/L09_flood_duration.py", ) html_out = html_out.replace("", '', 1) out_path = LESSONS / "L09_flood_duration.html" out_path.write_text(html_out, encoding="utf-8") print(f"\n[HTML] {out_path.name} を出力 ({out_path.stat().st_size:,} bytes)") print(f"=== DONE in {time.time()-t0:.1f}s ===")