# -*- coding: utf-8 -*- """ X08: 河川浸水想定区域 × 避難所立地 — 「避難所自体が水没するか?」 点 in ポリゴン判定研究 ================================================================================ DoBoX の以下 2 系統を空間的に重ね合わせ、4,065 件の避難所点が 「想定最大規模」「計画規模」の浸水想定区域ポリゴンに **入っているかどうか** を 1 件ごとに Ray casting で判定する **点 in ポリゴン判定研究** (X 水準, リテラシレベル)。 【本記事のスタイル: 点 in ポリゴン判定研究】 4,065 避難所の **どれが** 浸水想定域内に立地するかを 1 件ごとに判定し、 **危険避難所をランキング化** する。X07 の「面積率研究」(ポリゴン × ポリゴン)と 対比的に、本記事は「点 × ポリゴン」という別軸の空間結合を扱う。 散布図行列・PCA・ロジ回帰・k-means は使わない。 **危険避難所ランキング棒・市町別棒・ボーダー判定・フラグ×立地クロス・標高散布** の 5 図を主役に据える。 データ: - DoBoX #42 (避難所情報, JSON, 4,065 件; 既取得 data/shelters.json) - DoBoX #295 (河川浸水想定区域 計画規模 全河川, Shapefile, Web Mercator) - DoBoX #313 (河川浸水想定区域 想定最大規模 全河川, Shapefile, Web Mercator) - 市町別 平均標高 (シンボリック値, 既知の地理データから手動定義) 仮説 H1〜H5: H1: 全 4,065 避難所のうち 25〜45% が想定最大規模浸水想定区域内に立地する H2: 「危険避難所」(浸水域内 × 想定浸水深推定 > 標高差) は沿岸都市に集中する H3: 「ボーダー避難所」(計画規模で安全 / 想定最大規模で危険) が一定数 (50〜500件) 存在する H4: 災害種別フラグ「洪水対応(floodShFlg=1)」と「浸水想定区域内立地」は逆相関する (フラグが立っている避難所ほど安全な場所に置かれているはず) H5: 避難所が密集する市町ほど、危険避難所の絶対数も多い (= 避難所数 vs 危険避難所数の正相関) 主要可視化 (5 図): 図1 (主役): 危険避難所 上位 30 ランキング水平棒 (市町ラベル付) 図2: 市町別 危険避難所件数ランキング棒 図3: ボーダー避難所判定 (計画 OK / 最大 NG) 4 区分件数比較棒 図4: 洪水対応フラグ × 浸水域内立地 クロス集計 (4 セルヒートマップ) 図5: 避難所立地標高 (市町平均 m) vs 想定浸水域内フラグ 散布 + 箱ひげ 実装: - 点 in ポリゴン判定 = Ray casting (純 Python + numpy ベクトル化) - bbox プリフィルタで 4,065 点 × ~600 ポリゴン × 2 規模 = 488 万判定を秒単位に - 浸水 Shapefile を 1 度ロード後に逐次処理 (メモリ < 500 MB) 実行: cd "2026 DoBoX 教材" py -X utf8 lessons/X08_flood_shelter_in_polygon.py 制約: - 並列処理なし、メモリ < 500 MB """ from __future__ import annotations from pathlib import Path import json import math import warnings import zipfile import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import shapefile # pyshp from _common import (ROOT, ASSETS, LESSONS, render_lesson, code, figure, data_recipe, ensure_dataset) plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False warnings.filterwarnings("ignore", category=UserWarning) # === パス定義 ================================================================= DATA_DIR = ROOT / "data" / "extras" FLOOD_DIR = DATA_DIR / "flood_shp" FLOOD_DIR.mkdir(parents=True, exist_ok=True) SHELTERS_JSON = ROOT / "data" / "shelters.json" OUT_JUDGE = ASSETS / "X08_shelter_judge.csv" # 4,065 行 × 判定結果列 OUT_DANGER = ASSETS / "X08_danger_top.csv" # 危険避難所 上位 (浸水域内 × 低標高) OUT_CITY = ASSETS / "X08_city_summary.csv" # 市町別 集計 OUT_BORDER = ASSETS / "X08_border_breakdown.csv" # 計画 vs 最大 4区分内訳 OUT_FLAG = ASSETS / "X08_flag_cross.csv" # 洪水対応フラグ × 立地 クロス表 OUT_TRACK = ASSETS / "X08_track_one_shelter.csv" # 1 件追跡 (要件K) # 追加分析6-10 用 中間 CSV (3 本追加 = 計 9 本) OUT_MISMATCH = ASSETS / "X08_mismatch.csv" # フラグ無 × 浸水域内 の潜在リスク OUT_CAPSUM = ASSETS / "X08_capacity_summary.csv" # 市町別 収容力 × 危険度 OUT_CITYRISK = ASSETS / "X08_city_risk_rate.csv" # 市町別 危険率ヒートマップ用 PNG_RANK = ASSETS / "X08_danger_rank.png" # 主役: 危険避難所上位30 マルチ属性ランキング PNG_RANK_PANELS = ASSETS / "X08_danger_rank_panels.png" # 補助: 上位30 を 3 パネル並列 (capacity / 標高 / フラグ) PNG_CITY = ASSETS / "X08_city_danger_bar.png" # 市町別 危険件数 PNG_BORDER = ASSETS / "X08_border_breakdown.png" # ボーダー4区分 PNG_FLAG = ASSETS / "X08_flag_cross.png" # フラグ×立地クロス PNG_ELEV = ASSETS / "X08_elev_scatter.png" # 標高 vs 立地 # 追加分析6-10 用 PNG (5 枚追加 = 計 10 枚) PNG_MAP_DANGER = ASSETS / "X08_map_danger.png" # 危険度の地理マッピング PNG_MAP_BORDER = ASSETS / "X08_map_border.png" # border_class 別 地理マッピング PNG_MISMATCH = ASSETS / "X08_mismatch.png" # 5 災害種別フラグ × 浸水立地ミスマッチ PNG_CAP_DANGER = ASSETS / "X08_capacity_danger.png" # 収容力 × 危険度 バブル PNG_CITY_HEAT = ASSETS / "X08_city_heatmap.png" # 市町別 危険率ヒートマップ (地理) # === 0. 投影変換ユーティリティ =============================================== WGS_R = 6378137.0 # WGS84 球面近似 半径 (m) def merc_to_ll(x: float, y: float) -> tuple[float, float]: """Web Mercator (EPSG:3857) → WGS84 経緯度 (厳密式)。""" lon = math.degrees(x / WGS_R) lat = math.degrees(2 * math.atan(math.exp(y / WGS_R)) - math.pi / 2) return lon, lat # === 1. データ取得 (DoBoX 直叩き) ============================================= print("=== 1. データ取得 ===") FLOOD_RID = { "計画規模": (295, 23061, "shinsui_keikaku.zip", "shinsui_keikaku"), "想定最大規模": (313, 23118, "shinsui_souteisaidai.zip", "shinsui_souteisaidai"), } for label, (_, rid, fn, _stem) in FLOOD_RID.items(): p = FLOOD_DIR / fn ensure_dataset(p, resource_id=rid, min_bytes=1_000_000, label=f"浸水/{label}") # 避難所 JSON は前回までの取得物を再利用 (X03/L03 でも使用) if not SHELTERS_JSON.exists(): raise FileNotFoundError( f"data/shelters.json が見当たりません。L03 のスクリプトで先に取得してください。" ) # === 2. 避難所 JSON を読み込み 4,065 件をテーブル化 =========================== print("\n=== 2. 避難所 JSON 読込 ===") with open(SHELTERS_JSON, encoding="utf-8") as f: shelter_obj = json.load(f) shelters_raw = shelter_obj["items"] print(f" 総避難所数: {len(shelters_raw)} 件") shel_df = pd.DataFrame([{ "facility_id": s.get("facilityId"), "name": s.get("name"), "address": s.get("address01", ""), "lat": float(s["latitude"]), "lon": float(s["longitude"]), "muni": s.get("municipalityName", ""), "muni_cd": s.get("municipalityCd", ""), "capacity": s.get("capacity"), "flood_flg": int(s.get("floodShFlg") or 0), "sediment_flg": int(s.get("sedimentDisasterShFlg") or 0), "storm_flg": int(s.get("stormSurgeShFlg") or 0), "eq_flg": int(s.get("earthquakeShFlg") or 0), "tsunami_flg": int(s.get("tsunamiShFlg") or 0), } for s in shelters_raw]) print(f" shel_df: {len(shel_df)} 行 × {len(shel_df.columns)} 列") print(f" 洪水対応フラグ立ち: {shel_df['flood_flg'].sum()}/{len(shel_df)}" f" ({shel_df['flood_flg'].mean()*100:.1f}%)") print(f" 市町数: {shel_df['muni'].nunique()}") # === 2.5 市町 → 市 (大くくり) 集約 + シンボリック平均標高テーブル =========== # 広島市の各区を「広島市」に集約しつつも、市町別ランキングは元粒度で見る。 # 平均標高はシンボリック (公開地形データから既知の概算値を手動で定義)。 # ここでは「役所所在地付近の標高」を採用 (m, 概数)。 CITY_ELEV_M = { # 広島市 8 区 (デルタ平野) "広島市中区": 3, "広島市東区": 6, "広島市南区": 4, "広島市西区": 4, "広島市安佐南区": 30, "広島市安佐北区": 110, "広島市安芸区": 25, "広島市佐伯区": 18, # 周辺都市 "呉市": 8, "竹原市": 5, "三原市": 6, "尾道市": 5, "福山市": 7, "府中市": 80, "三次市": 165, "庄原市": 230, "大竹市": 6, "東広島市": 220, "廿日市市": 8, "安芸高田市": 175, "江田島市": 5, "海田町": 4, "府中町": 8, "坂町": 4, "熊野町": 220, "安芸太田町": 350, "北広島町": 470, "大崎上島町": 30, "世羅町": 380, "神石高原町": 480, # 念の為のデフォルト (未定義市町用) } DEFAULT_ELEV = 50 # 未収録時のフォールバック shel_df["elev_m"] = shel_df["muni"].map(CITY_ELEV_M).fillna(DEFAULT_ELEV).astype(float) print(f" 標高未定義 (DEFAULT_ELEV={DEFAULT_ELEV} 適用): " f"{(shel_df['elev_m'] == DEFAULT_ELEV).sum()} 件") # === 3. 浸水 Shapefile を読み込み、ポリゴン bbox を経緯度で保持 =============== print("\n=== 3. 浸水 Shapefile 読込 (Web Mercator → 経緯度) ===") flood_data: dict[str, list[dict]] = {} for label, (_, _, fn, stem) in FLOOD_RID.items(): zp = FLOOD_DIR / fn ext = FLOOD_DIR / stem ext.mkdir(exist_ok=True) with zipfile.ZipFile(zp) as z: for n in z.namelist(): if not (ext / Path(n).name).exists(): z.extract(n, ext) shp_path = next(ext.glob("*.shp")) sf = shapefile.Reader(str(shp_path)) fields = [f[0] for f in sf.fields[1:]] polys = [] for shp_idx, (rec, shp) in enumerate(zip(sf.records(), sf.shapes())): if not shp.points: continue bx0, by0, bx1, by1 = shp.bbox lon0, lat0 = merc_to_ll(bx0, by0) lon1, lat1 = merc_to_ll(bx1, by1) parts = list(shp.parts) + [len(shp.points)] rings = [] for k in range(len(shp.parts)): seg = shp.points[parts[k]:parts[k + 1]] lons = [merc_to_ll(p[0], p[1])[0] for p in seg] lats = [merc_to_ll(p[0], p[1])[1] for p in seg] rings.append((np.asarray(lons, dtype=np.float64), np.asarray(lats, dtype=np.float64))) polys.append({ "shp_idx": shp_idx, "kasen": rec[fields.index("kasen")] if "kasen" in fields else "?", "suikei": rec[fields.index("suikei")] if "suikei" in fields else "?", "lon_min": lon0, "lon_max": lon1, "lat_min": lat0, "lat_max": lat1, "rings": rings, }) flood_data[label] = polys print(f" {label}: {len(polys)} ポリゴン (河川数 = " f"{len(set(p['kasen'] for p in polys))})") sf.close() # === 4. Ray casting (ベクトル化) ============================================= def points_in_ring_vec(xs_pt: np.ndarray, ys_pt: np.ndarray, ring_xs: np.ndarray, ring_ys: np.ndarray) -> np.ndarray: """N 点が 1 リング (M 頂点) の内側か。N が小さい場合の安定実装。""" n_pts = len(xs_pt) inside = np.zeros(n_pts, dtype=bool) m = len(ring_xs) j_idx = (np.arange(m) - 1) % m for i in range(m): yi = ring_ys[i]; yj = ring_ys[j_idx[i]] xi = ring_xs[i]; xj = ring_xs[j_idx[i]] cond_y = (yi > ys_pt) != (yj > ys_pt) if not cond_y.any(): continue x_int = (xj - xi) * (ys_pt - yi) / (yj - yi + 1e-15) + xi cond = cond_y & (xs_pt < x_int) inside ^= cond return inside def points_in_any_polygon_set(xs_pt: np.ndarray, ys_pt: np.ndarray, polys: list) -> np.ndarray: """N 点 vs P 個のポリゴン集合 → bool 配列 (N,)。 bbox プリフィルタ + 該当点だけ Ray casting。""" n_pts = len(xs_pt) inside = np.zeros(n_pts, dtype=bool) for p in polys: # bbox プリフィルタ (ポリゴン bbox 外の点は確定で外) cand = (~inside) & \ (xs_pt >= p["lon_min"]) & (xs_pt <= p["lon_max"]) & \ (ys_pt >= p["lat_min"]) & (ys_pt <= p["lat_max"]) if not cand.any(): continue idx = np.where(cand)[0] sub_x = xs_pt[idx] sub_y = ys_pt[idx] sub_in = np.zeros(len(idx), dtype=bool) for ring_x, ring_y in p["rings"]: if len(ring_x) < 3: continue still = ~sub_in if not still.any(): break r = points_in_ring_vec(sub_x[still], sub_y[still], ring_x, ring_y) # 偶奇規則: 1 つでも入れば inside (穴は無視, 近似) sub_in[still] |= r inside[idx] |= sub_in return inside # === 5. 4,065 避難所 × 2 規模 を一括判定 ===================================== print("\n=== 5. 点 in ポリゴン 判定 (4,065 × 2 規模) ===") xs = shel_df["lon"].values.astype(np.float64) ys = shel_df["lat"].values.astype(np.float64) import time as _t for label in ["計画規模", "想定最大規模"]: t0 = _t.time() in_flood = points_in_any_polygon_set(xs, ys, flood_data[label]) col = "in_keikaku" if label == "計画規模" else "in_max" shel_df[col] = in_flood.astype(int) n = int(in_flood.sum()) print(f" {label}: {n}/{len(shel_df)} 件 浸水域内 " f"({n/len(shel_df)*100:.1f}%, {_t.time()-t0:.1f}s)") # 派生フラグ shel_df["only_max_in"] = ((shel_df["in_max"] == 1) & (shel_df["in_keikaku"] == 0)).astype(int) shel_df["both_in"] = ((shel_df["in_max"] == 1) & (shel_df["in_keikaku"] == 1)).astype(int) shel_df["both_safe"] = ((shel_df["in_max"] == 0) & (shel_df["in_keikaku"] == 0)).astype(int) # OUT_JUDGE 保存は danger_score / border_class 追加後 (セクション 8 末尾) に行う # === 6. 危険避難所スコア定義 ================================================ # 「危険」= 浸水域内 (想定最大規模) AND 標高が低い # danger_score = (浸水域内なら 1, それ以外 0) × (50 - elev_m を 0〜1 に正規化) # ここでは elev_m が小さいほど危険、計画+最大両方で水没 (both_in=1) なら強い加点。 # シンプルな線形合成スコア (主観だが透明) DANGER_ELEV_REF = 30.0 # この高度を超えれば標高による危険減衰が始まる shel_df["danger_score"] = ( shel_df["in_max"] * (1.0 + 0.5 * shel_df["both_in"]) # 計画規模も内なら 1.5 倍 * np.clip(1.0 - shel_df["elev_m"] / DANGER_ELEV_REF, 0.05, 1.0) ).round(4) danger = shel_df[shel_df["in_max"] == 1].copy() danger = danger.sort_values("danger_score", ascending=False).reset_index(drop=True) danger_top = danger.head(30).copy() danger_top.to_csv(OUT_DANGER, index=False, encoding="utf-8-sig") print(f"\n=== 6. 危険避難所スコア (想定最大規模) ===") print(f" 浸水域内: {len(danger)} 件 (= H1 検証分母)") print(f" 上位 30 (danger_score 降順):") print(danger_top[["name", "muni", "elev_m", "in_keikaku", "in_max", "flood_flg", "danger_score"]].head(10).to_string(index=False)) # === 7. 市町別 集計 ========================================================== print("\n=== 7. 市町別 集計 ===") city_summary = shel_df.groupby("muni").agg( n_total=("facility_id", "count"), n_in_max=("in_max", "sum"), n_in_keikaku=("in_keikaku", "sum"), n_only_max=("only_max_in", "sum"), n_flood_flg=("flood_flg", "sum"), elev_m=("elev_m", "first"), ).reset_index() city_summary["pct_in_max"] = (city_summary["n_in_max"] / city_summary["n_total"] * 100).round(2) city_summary["pct_only_max"] = (city_summary["n_only_max"] / city_summary["n_total"] * 100).round(2) city_summary = city_summary.sort_values("n_in_max", ascending=False).reset_index(drop=True) city_summary.to_csv(OUT_CITY, index=False, encoding="utf-8-sig") print(city_summary.head(10).to_string(index=False)) print(f" → {OUT_CITY.name}") # === 8. ボーダー避難所 4 区分 ================================================ # A: 両方 OUT (安全) ; B: 計画 IN / 最大 IN (両 NG) ; # C: 計画 OUT / 最大 IN (= ボーダー) ; D: 計画 IN / 最大 OUT (異常: 想定最大規模が縮むのは稀) print("\n=== 8. ボーダー判定 4 区分 ===") def classify(row): k, m = row["in_keikaku"], row["in_max"] if k == 0 and m == 0: return "A 両規模で安全" if k == 1 and m == 1: return "B 両規模で水没" if k == 0 and m == 1: return "C ボーダー (最大規模のみ水没)" return "D 計画のみ水没 (異常)" shel_df["border_class"] = shel_df.apply(classify, axis=1) border_breakdown = shel_df["border_class"].value_counts().reset_index() border_breakdown.columns = ["区分", "件数"] border_breakdown["割合(%)"] = (border_breakdown["件数"] / len(shel_df) * 100).round(2) border_breakdown = border_breakdown.sort_values("区分").reset_index(drop=True) border_breakdown.to_csv(OUT_BORDER, index=False, encoding="utf-8-sig") print(border_breakdown.to_string(index=False)) n_border = int((shel_df["border_class"] == "C ボーダー (最大規模のみ水没)").sum()) print(f" ボーダー避難所件数 (H3): {n_border} 件") # OUT_JUDGE を 19+2 列で書き出し (border_class, danger_score を追加後) shel_df.to_csv(OUT_JUDGE, index=False, encoding="utf-8-sig") print(f" → {OUT_JUDGE.name} ({len(shel_df)} 行 × {len(shel_df.columns)} 列)") # === 9. 洪水対応フラグ × 立地 クロス集計 (2x2) ============================== print("\n=== 9. 洪水対応フラグ × 想定最大規模立地 クロス ===") ct = pd.crosstab( shel_df["flood_flg"].map({0: "洪水対応無", 1: "洪水対応有"}), shel_df["in_max"].map({0: "浸水域外", 1: "浸水域内"}), margins=True, margins_name="合計", ) ct.to_csv(OUT_FLAG, encoding="utf-8-sig") print(ct) # χ² 風の独立性指標 (簡易, scipy なし) def expected_chi2(crosstab: pd.DataFrame) -> tuple[float, pd.DataFrame]: """2x2 の期待度数と Pearson χ² 値を計算 (合計行/列込み crosstab を渡す)。""" obs = crosstab.iloc[:-1, :-1].values.astype(float) row_sum = obs.sum(axis=1) col_sum = obs.sum(axis=0) total = obs.sum() exp = np.outer(row_sum, col_sum) / total chi2 = float(((obs - exp) ** 2 / exp).sum()) df_exp = pd.DataFrame(exp, index=crosstab.index[:-1], columns=crosstab.columns[:-1]) return chi2, df_exp chi2_val, exp_df = expected_chi2(ct) print(f"\n Pearson χ² = {chi2_val:.2f} (df=1, p<0.001 の閾値 ≈ 10.83)") # 「フラグありかつ域内」の率 vs 「フラグなしかつ域内」の率 flag_in_rate = shel_df[shel_df["flood_flg"] == 1]["in_max"].mean() noflag_in_rate = shel_df[shel_df["flood_flg"] == 0]["in_max"].mean() print(f" フラグあり: 浸水域内率 {flag_in_rate*100:.1f}%") print(f" フラグなし: 浸水域内率 {noflag_in_rate*100:.1f}%") print(f" 逆相関?: {flag_in_rate < noflag_in_rate} (H4 検証)") # === 10. 1件追跡 (要件K) ===================================================== print("\n=== 10. 1件追跡 (海田町の代表避難所) ===") # 海田町の代表的な公民館を 1 つ拾う (浸水域内に立地している可能性が高い) target_muni = "海田町" cand = shel_df[shel_df["muni"] == target_muni].copy() # できれば「公民館 OR 集会所 OR センター」を含む key_match = cand[cand["name"].str.contains("公民館|センター|集会|小学校", na=False)] if len(key_match): cand = key_match # in_max == 1 を優先 in_max_cand = cand[cand["in_max"] == 1] if len(in_max_cand): target = in_max_cand.iloc[0] else: target = cand.iloc[0] track = pd.DataFrame([ {"段階": "1. JSON 取得", "値": "DoBoX #42 避難所情報 (4,065 件)", "内訳": "data/shelters.json から1件抽出"}, {"段階": "2. 1件選択", "値": str(target["name"]), "内訳": f"市町={target['muni']} / 住所={target['address']}"}, {"段階": "3. 経緯度抽出", "値": f"({target['lon']:.5f}, {target['lat']:.5f})", "內訳" if False else "内訳": "JSON の latitude / longitude 文字列を float 化"}, {"段階": "4. 計画規模 判定", "値": "IN" if target["in_keikaku"] == 1 else "OUT", "內訳" if False else "内訳": f"Ray casting × {len(flood_data['計画規模'])} ポリゴン"}, {"段階": "5. 想定最大規模 判定", "値": "IN" if target["in_max"] == 1 else "OUT", "内訳": f"Ray casting × {len(flood_data['想定最大規模'])} ポリゴン"}, {"段階": "6. 標高 lookup", "値": f"{target['elev_m']:.0f} m", "内訳": f"市町別代表標高 ({target['muni']} = {target['elev_m']:.0f} m)"}, {"段階": "7. ボーダー区分", "値": str(target["border_class"]), "内訳": "(計画, 最大) の 4 区分パターン"}, {"段階": "8. 洪水対応フラグ", "値": "有" if target["flood_flg"] == 1 else "無", "内訳": "JSON の floodShFlg (DoBoX 側の災害種別フラグ)"}, {"段階": "9. 危険スコア", "値": f"{target['danger_score']:.3f}", "内訳": "in_max × (1 + 0.5 × both_in) × clip(1 - elev/30, 0.05, 1.0)"}, ]) track = track.rename(columns={"內訳": "内訳"}).fillna("") # Drop accidental 內訳 col if exists if "內訳" in track.columns: track = track.drop(columns=["內訳"]) track.to_csv(OUT_TRACK, index=False, encoding="utf-8-sig") print(track.to_string(index=False)) # === 11. 仮説判定 =========================================================== print("\n=== 11. 仮説判定 ===") n_total = len(shel_df) n_in_max = int(shel_df["in_max"].sum()) pct_in_max = n_in_max / n_total * 100 h1_supported = 25 <= pct_in_max <= 45 # H2: 沿岸都市の危険避難所占有率 (近似: 海沿い 5 市町に集中するか) COASTAL_CITIES = ["広島市中区", "広島市南区", "広島市西区", "呉市", "尾道市", "福山市", "三原市", "竹原市", "大竹市", "廿日市市", "海田町", "府中町", "坂町", "江田島市", "大崎上島町"] n_danger = len(danger) n_coastal_danger = int(danger["muni"].isin(COASTAL_CITIES).sum()) coastal_share = n_coastal_danger / max(n_danger, 1) * 100 h2_supported = coastal_share > 50 h3_supported = 50 <= n_border <= 500 # H4: フラグあり群とフラグなし群で「浸水域内率」を比較 h4_supported = flag_in_rate < noflag_in_rate # H5: 市町別 (避難所数 vs 危険避難所数) の相関 corr_h5 = float(np.corrcoef(city_summary["n_total"], city_summary["n_in_max"])[0, 1]) h5_supported = corr_h5 > 0.6 verdict = { "H1_浸水域内率25-45%": (h1_supported, f"{pct_in_max:.1f}%"), "H2_沿岸集中>50%": (h2_supported, f"沿岸シェア {coastal_share:.1f}% / 危険 {n_danger}件中"), "H3_ボーダー50-500件": (h3_supported, f"{n_border}件"), "H4_フラグ逆相関": (h4_supported, f"フラグあり {flag_in_rate*100:.1f}% vs なし {noflag_in_rate*100:.1f}%"), "H5_全件×危険正の相関": (h5_supported, f"r={corr_h5:.3f}"), } for k, (ok, msg) in verdict.items(): print(f" {k}: {'支持' if ok else '反証'} ({msg})") findings_text = "\n".join( f" {k}: {'支持' if ok else '反証'} ({msg})" for k, (ok, msg) in verdict.items() ) # === 12. 図1 (主役): 危険避難所 上位30 マルチ属性ランキング =================== # 旧版は danger_score で赤濃度を変えるだけ → ほぼ全バーが似た色で情報密度が低かった。 # 新版は 1 枚に「収容力 (バー長) × 市町 (色) × 標高 (右数値) × 危険要因 (記号)」を # 圧縮表示し、見ただけで「どこが大きい・低地・複合危険」を判別できる。 print("\n=== 12. 図1 (主役): 危険避難所 上位30 マルチ属性ランキング ===") top = danger.head(30).iloc[::-1].reset_index(drop=True) # 下から上へ昇順 ypos = np.arange(len(top)) # capacity 欠損は中央値で補完 (バー長 0 だと「危険なのに見えない」になるため) cap_med = float(shel_df["capacity"].dropna().median()) top["cap_disp"] = top["capacity"].fillna(cap_med).astype(float) # 市町 → カテゴリカル色 unique_munis = list(top["muni"].unique()) cmap_muni = plt.get_cmap("tab20") muni_color = {m: cmap_muni(i % 20) for i, m in enumerate(unique_munis)} bar_colors_muni = [muni_color[m] for m in top["muni"]] # 5 災害フラグ FLAG_COLS = ["flood_flg", "sediment_flg", "storm_flg", "eq_flg", "tsunami_flg"] top["n_hazard"] = top[FLAG_COLS].sum(axis=1) fig, ax = plt.subplots(figsize=(13, 10)) ax.barh(ypos, top["cap_disp"], color=bar_colors_muni, edgecolor="#222", linewidth=0.5) # y 軸ラベル: 順位 + 名 + 市町 labels = [f"#{len(top)-i:>2} {r['name'][:16]} ({r['muni']})" for i, (_, r) in enumerate(top.iterrows())] ax.set_yticks(ypos) ax.set_yticklabels(labels, fontsize=9) # 注: 順位は danger 降順なので、グラフ最上部 (i=len-1) が「1 位 = 最も危険」 # (上で iloc[::-1] 反転済 → 下から #30, ..., 上が #1) # バー右端に補助マーカー + 標高表示 xmax = top["cap_disp"].max() mark_off_star = xmax * 0.015 # ★ オフセット mark_off_tri = xmax * 0.05 # ▲ mark_off_dot = xmax * 0.085 # ● text_off = xmax * 0.13 # 数値表示 for i, (_, r) in enumerate(top.iterrows()): barlen = r["cap_disp"] # ★: 洪水フラグ無し (制度上は洪水避難所ではない = 特に運用矛盾) if r["flood_flg"] == 0: ax.text(barlen + mark_off_star, ypos[i], "★", va="center", ha="center", fontsize=14, color="#cf222e", fontweight="bold") # ▲: 標高 ≦ 5m (低地) if r["elev_m"] <= 5: ax.text(barlen + mark_off_tri, ypos[i], "▲", va="center", ha="center", fontsize=12, color="#bf5700") # ●: 5 災害フラグ全部 1 (= 複合避難所) if r["n_hazard"] == 5: ax.text(barlen + mark_off_dot, ypos[i], "●", va="center", ha="center", fontsize=12, color="#1a7f37") # 標高 + capacity 数値 cap_int = int(r["cap_disp"]) cap_lbl = f"{cap_int:,}人" + ("*" if pd.isna(r["capacity"]) else "") ax.text(barlen + text_off, ypos[i], f"標高{r['elev_m']:.0f}m {cap_lbl} score={r['danger_score']:.2f}", va="center", fontsize=8.2, color="#333") ax.set_xlabel("収容力 capacity (人) ※ * は欠損中央値補完") ax.set_xlim(0, xmax * 1.55) ax.set_title(f"図1 (主役): 危険避難所 上位 30 マルチ属性ランキング (浸水域内 {len(danger)}件中)\n" "バー長=収容力 / バー色=市町 / ★=洪水フラグ無 ▲=標高≦5m ●=5災害フラグ全立 右数値=標高/収容/score") ax.grid(axis="x", alpha=0.3) # 市町凡例 (上位 8 のみ; 残りは "他") muni_count = top["muni"].value_counts() top_munis = list(muni_count.head(8).index) handles = [plt.Rectangle((0, 0), 1, 1, color=muni_color[m]) for m in top_munis] labels_leg = [f"{m} ({muni_count[m]})" for m in top_munis] if len(unique_munis) > 8: labels_leg.append(f"他 {len(unique_munis)-8} 市町") handles.append(plt.Rectangle((0, 0), 1, 1, color="#ccc")) ax.legend(handles, labels_leg, loc="lower right", fontsize=8.5, ncol=2, title="市町 (上位30 中の出現件数)", title_fontsize=9) plt.tight_layout() plt.savefig(PNG_RANK, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_RANK.name}") print(f" マーカー集計: ★(洪水フラグ無)={int((top['flood_flg']==0).sum())} / " f"▲(標高≦5m)={int((top['elev_m']<=5).sum())} / " f"●(5フラグ全立)={int((top['n_hazard']==5).sum())}") # === 12B. 図1B (補助): 上位30 を 3 パネル並列 (capacity / 標高 / フラグ heatmap) # 1 図に 4 軸を詰めると密になりすぎる場合の「展開図」。同じ上位30 を行で揃え、 # 左=収容力 (市町色) / 中=標高 (低いほど赤) / 右=5災害フラグの 1/0 ヒートマップ。 print("\n=== 12B. 図1B (補助): 上位30 3 パネル並列 ===") fig, axes = plt.subplots(1, 3, figsize=(15, 9), gridspec_kw={"width_ratios": [1.4, 1.0, 1.0]}, sharey=True) # --- 左パネル: capacity (市町色) --- ax = axes[0] ax.barh(ypos, top["cap_disp"], color=bar_colors_muni, edgecolor="#222", linewidth=0.4) labels_p = [f"#{len(top)-i:>2} {r['name'][:14]}" for i, (_, r) in enumerate(top.iterrows())] ax.set_yticks(ypos) ax.set_yticklabels(labels_p, fontsize=8.5) for i, (_, r) in enumerate(top.iterrows()): cap_int = int(r["cap_disp"]) ax.text(r["cap_disp"] + top["cap_disp"].max()*0.01, ypos[i], f"{cap_int:,}", va="center", fontsize=7.5, color="#333") ax.set_xlabel("収容力 (人)") ax.set_title("(左) capacity / 色=市町") ax.grid(axis="x", alpha=0.3) # --- 中央パネル: 標高 (低いほど赤) --- ax = axes[1] elev_vals = top["elev_m"].values.astype(float) elev_max_ref = 30.0 # これ以上は青寄り (=危険減衰) elev_norm = np.clip(1.0 - elev_vals / elev_max_ref, 0.0, 1.0) # 低標高ほど 1 elev_colors = plt.get_cmap("RdYlBu")(1.0 - elev_norm) # 低い=赤, 高い=青 ax.barh(ypos, elev_vals, color=elev_colors, edgecolor="#222", linewidth=0.4) for i, v in enumerate(elev_vals): ax.text(v + 1, ypos[i], f"{v:.0f}m", va="center", fontsize=8, color="#333") ax.axvline(5, color="#cf222e", lw=0.8, ls="--", alpha=0.7) ax.text(5, len(top)-0.3, " 5m", color="#cf222e", fontsize=8, va="top") ax.set_xlabel("市町平均 標高 (m)") ax.set_title("(中) 標高 / 赤=低地 (≦5m に破線)") ax.grid(axis="x", alpha=0.3) # --- 右パネル: 5 災害フラグ heatmap --- ax = axes[2] flag_mat = top[FLAG_COLS].values.astype(float) # (30, 5) # 0 はほぼ白、1 は濃緑 ax.imshow(flag_mat, aspect="auto", cmap="Greens", vmin=0, vmax=1.5, extent=[-0.5, 4.5, -0.5, len(top)-0.5], origin="lower") ax.set_xticks(range(5)) ax.set_xticklabels(["洪水", "土砂", "高潮", "地震", "津波"], fontsize=9) ax.set_xlim(-0.5, 4.5) # 数値 (1) と、洪水フラグ無は赤×で強調 for i in range(len(top)): for j in range(5): v = flag_mat[i, j] if v == 1: ax.text(j, i, "1", va="center", ha="center", fontsize=8.5, color="#0a3d0a", fontweight="bold") else: ax.text(j, i, "0", va="center", ha="center", fontsize=8, color="#bbb") # 洪水フラグ無 (= 第0列が0) の行に右端 ★ for i in range(len(top)): if flag_mat[i, 0] == 0: ax.text(4.7, i, "★", va="center", ha="center", fontsize=11, color="#cf222e", fontweight="bold") ax.set_xlim(-0.5, 5.2) # y 軸グリッドを行間に for i in range(len(top)+1): ax.axhline(i-0.5, color="#eee", lw=0.4, zorder=0) ax.set_title("(右) 災害種別フラグ 5列 / ★=洪水無") fig.suptitle(f"図1B (補助): 危険避難所 上位30 を 3 軸並列で展開 (左=capacity / 中=標高 / 右=5災害フラグ)\n" f"行は図1 と同じ並び (上=#1 最も危険)", y=1.005, fontsize=12) plt.tight_layout() plt.savefig(PNG_RANK_PANELS, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_RANK_PANELS.name}") # === 13. 図2: 市町別 危険避難所件数ランキング棒 ============================== print("\n=== 13. 図2: 市町別 危険避難所件数ランキング ===") top_city = city_summary.head(20).iloc[::-1].reset_index(drop=True) fig, ax = plt.subplots(figsize=(11, 8)) ypos = np.arange(len(top_city)) ax.barh(ypos, top_city["n_in_max"], color="#cf222e", edgecolor="#222", linewidth=0.4, label="想定最大規模 浸水域内") ax.barh(ypos, top_city["n_in_keikaku"], color="#0969da", edgecolor="#222", linewidth=0.4, alpha=0.85, label="計画規模 浸水域内") for i, (_, r) in enumerate(top_city.iterrows()): ax.text(r["n_in_max"] + 1, ypos[i], f"{r['n_in_max']}/{r['n_total']} ({r['pct_in_max']:.1f}%)", va="center", fontsize=8.5) ax.set_yticks(ypos) ax.set_yticklabels(top_city["muni"], fontsize=10) ax.set_xlabel("浸水域内立地避難所数 (件)") ax.set_title(f"図2: 市町別 危険避難所件数 上位 20 (全 {city_summary['muni'].nunique()} 市町中)\n" f"想定最大規模(赤) と 計画規模(青) を重ね描き / 数値は最大/全件 (%)") ax.legend(loc="lower right") ax.grid(axis="x", alpha=0.3) plt.tight_layout() plt.savefig(PNG_CITY, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_CITY.name}") # === 14. 図3: ボーダー避難所 4区分件数 比較棒 ================================ print("\n=== 14. 図3: ボーダー判定 4 区分 ===") fig, axes = plt.subplots(1, 2, figsize=(13, 5.5), gridspec_kw={"width_ratios": [1.2, 1.0]}) ax = axes[0] b = border_breakdown.copy() # 行順を固定 (D 異常を最後に) order = ["A 両規模で安全", "B 両規模で水没", "C ボーダー (最大規模のみ水没)", "D 計画のみ水没 (異常)"] b = b.set_index("区分").reindex(order, fill_value=0).reset_index() xpos = np.arange(len(b)) bar_colors = ["#1a7f37", "#cf222e", "#bf8700", "#888"] ax.bar(xpos, b["件数"], color=bar_colors, edgecolor="#222", linewidth=0.4) for i, (_, r) in enumerate(b.iterrows()): ax.text(xpos[i], r["件数"] + max(b["件数"]) * 0.015, f"{r['件数']:,}\n({r['割合(%)']:.2f}%)", ha="center", fontsize=10) ax.set_xticks(xpos) ax.set_xticklabels([s.split(" ", 1)[0] for s in b["区分"]], fontsize=11) ax.set_ylabel("避難所件数") ax.set_title(f"図3-(a): 計画 × 最大 4 区分内訳 (n={len(shel_df):,})\n" "C = 計画OK・最大NG = 「ボーダー避難所」") ax.grid(axis="y", alpha=0.3) # 右パネル: ボーダー / 全危険 の内訳円風 横棒 ax = axes[1] total_max = b.loc[b["区分"].isin(["B 両規模で水没", "C ボーダー (最大規模のみ水没)"]), "件数"].sum() border_n = b.loc[b["区分"] == "C ボーダー (最大規模のみ水没)", "件数"].iloc[0] both_n = b.loc[b["区分"] == "B 両規模で水没", "件数"].iloc[0] ax.barh([0], [both_n], color="#cf222e", edgecolor="#222", label="B 両規模で水没") ax.barh([0], [border_n], left=[both_n], color="#bf8700", edgecolor="#222", label=f"C ボーダー (最大のみ): {border_n} 件") ax.set_yticks([0]) ax.set_yticklabels(["想定最大規模で水没\n(B + C)"]) ax.set_xlabel("避難所件数") ax.set_title(f"図3-(b): 想定最大規模 浸水域内 {total_max:,} 件 の内訳\n" f"ボーダー比率 = {border_n}/{total_max} " f"= {border_n/max(total_max,1)*100:.1f}%") ax.legend(loc="lower right") ax.grid(axis="x", alpha=0.3) plt.tight_layout() plt.savefig(PNG_BORDER, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_BORDER.name}") # === 15. 図4: フラグ × 立地 クロス ヒートマップ ============================== print("\n=== 15. 図4: 洪水対応フラグ × 想定最大規模立地 クロス ===") fig, ax = plt.subplots(figsize=(8, 5.5)) mat = ct.iloc[:-1, :-1].values # 2×2 row_labels = ct.index[:-1] col_labels = ct.columns[:-1] im = ax.imshow(mat, cmap="YlOrRd", aspect="auto") for i in range(mat.shape[0]): for j in range(mat.shape[1]): share = mat[i, j] / mat.sum() * 100 ax.text(j, i, f"{mat[i, j]:,}\n({share:.1f}%)", ha="center", va="center", fontsize=12, color="#fff" if mat[i, j] > mat.max() * 0.5 else "#222", weight="bold") ax.set_xticks(np.arange(len(col_labels))) ax.set_xticklabels(col_labels, fontsize=11) ax.set_yticks(np.arange(len(row_labels))) ax.set_yticklabels(row_labels, fontsize=11) ax.set_title(f"図4: 洪水対応フラグ × 想定最大規模 浸水域内立地 (2×2 クロス)\n" f"χ² = {chi2_val:.2f} (p<0.001 閾値 ≈ 10.83) / " f"フラグ有内率 {flag_in_rate*100:.1f}% vs フラグ無内率 {noflag_in_rate*100:.1f}%") fig.colorbar(im, ax=ax, shrink=0.7, label="件数") plt.tight_layout() plt.savefig(PNG_FLAG, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_FLAG.name}") # === 16. 図5: 標高 vs 浸水立地 散布 + 箱ひげ ================================ print("\n=== 16. 図5: 標高 vs 立地 ===") fig, axes = plt.subplots(1, 2, figsize=(14, 6.0), gridspec_kw={"width_ratios": [1.4, 1.0]}) # 散布 (jitter) ax = axes[0] rng_j = np.random.default_rng(20260502) jitter_x = rng_j.uniform(-0.4, 0.4, len(shel_df)) y_elev = shel_df["elev_m"].values x_in = shel_df["in_max"].values + jitter_x colors_elev = np.where(shel_df["in_max"] == 1, "#cf222e", "#1a7f37") ax.scatter(x_in, y_elev, s=8, c=colors_elev, alpha=0.4, edgecolor="none") ax.set_xticks([0, 1]) ax.set_xticklabels(["浸水域外", "浸水域内"]) ax.set_ylabel("市町別 平均標高 (m, シンボリック値)") ax.set_title("図5-(a): 各避難所の立地標高 vs 想定最大規模 域内/域外 (ジッタ散布, n=4,065)") ax.grid(alpha=0.3) ax.set_ylim(0, 520) # 箱ひげ比較 (浸水域内 vs 域外) ax = axes[1] elev_in = shel_df.loc[shel_df["in_max"] == 1, "elev_m"].values elev_out = shel_df.loc[shel_df["in_max"] == 0, "elev_m"].values bp = ax.boxplot([elev_in, elev_out], vert=True, patch_artist=True, widths=0.5, labels=[f"浸水域内\n(n={len(elev_in)})", f"浸水域外\n(n={len(elev_out)})"]) for patch, color in zip(bp["boxes"], ["#fdd", "#dfd"]): patch.set_facecolor(color) patch.set_edgecolor("#222") ax.axhline(np.median(shel_df["elev_m"]), color="#0969da", ls="--", lw=1.0, label=f"全体中央値 {np.median(shel_df['elev_m']):.0f} m") ax.set_ylabel("市町平均標高 (m)") ax.set_title(f"図5-(b): 立地標高分布 比較\n" f"中央値: 域内 {np.median(elev_in):.0f}m vs 域外 {np.median(elev_out):.0f}m") ax.legend(loc="upper right") ax.grid(axis="y", alpha=0.3) plt.tight_layout() plt.savefig(PNG_ELEV, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_ELEV.name}") # === 16-A. 追加分析6: 危険度の地理マッピング ================================ print("\n=== 16-A. 図6: 危険度の地理マッピング ===") fig, ax = plt.subplots(figsize=(11.5, 8.5)) # 浸水域外 (in_max=0) を薄く背景に置く mask_out = shel_df["in_max"] == 0 mask_in = shel_df["in_max"] == 1 ax.scatter(shel_df.loc[mask_out, "lon"], shel_df.loc[mask_out, "lat"], s=6, c="#cccccc", alpha=0.35, edgecolor="none", label=f"浸水域外 (n={int(mask_out.sum()):,})") # 浸水域内 (in_max=1) のみ濃く danger_score で色塗り sc = ax.scatter(shel_df.loc[mask_in, "lon"], shel_df.loc[mask_in, "lat"], s=22, c=shel_df.loc[mask_in, "danger_score"], cmap="RdYlGn_r", vmin=0.0, vmax=1.5, alpha=0.85, edgecolor="#222", linewidth=0.25) cbar = plt.colorbar(sc, ax=ax, shrink=0.78, pad=0.02) cbar.set_label("danger_score (0=安全 → 1.5=最大危険)") ax.set_xlabel("経度 (°E)") ax.set_ylabel("緯度 (°N)") ax.set_title(f"図6: 危険度の地理マッピング (4,065 避難所, 浸水域内 {int(mask_in.sum()):,} 件のみ濃く描画)\n" "色 = danger_score (RdYlGn 反転: 緑→黄→赤)") ax.grid(alpha=0.3) ax.legend(loc="lower left", fontsize=9) ax.set_aspect(1 / np.cos(np.deg2rad(34.4))) # 広島県中緯度の縦横比 plt.tight_layout() plt.savefig(PNG_MAP_DANGER, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_MAP_DANGER.name}") # === 16-B. 追加分析7: border_class 別 地理マッピング ========================= print("\n=== 16-B. 図7: border_class 別 地理マッピング ===") fig, ax = plt.subplots(figsize=(11.5, 8.5)) bc_styles = { "A 両規模で安全": {"color": "#1a7f37", "alpha": 0.30, "size": 7, "label_short": "A 安全"}, "B 両規模で水没": {"color": "#cf222e", "alpha": 0.85, "size": 18, "label_short": "B 両危険"}, "C ボーダー (最大規模のみ水没)": {"color": "#bf8700", "alpha": 0.85, "size": 18, "label_short": "C ボーダー"}, "D 計画のみ水没 (異常)": {"color": "#888888", "alpha": 0.85, "size": 22, "label_short": "D 異常"}, } # 描画順: A を背景, B/C/D を前面 for cls in ["A 両規模で安全", "B 両規模で水没", "C ボーダー (最大規模のみ水没)", "D 計画のみ水没 (異常)"]: style = bc_styles[cls] sub = shel_df[shel_df["border_class"] == cls] if len(sub) == 0: continue ax.scatter(sub["lon"], sub["lat"], s=style["size"], c=style["color"], alpha=style["alpha"], edgecolor="#222" if cls != "A 両規模で安全" else "none", linewidth=0.25, label=f"{style['label_short']} (n={len(sub):,})") ax.set_xlabel("経度 (°E)") ax.set_ylabel("緯度 (°N)") ax.set_title("図7: border_class 別 地理マッピング (4,065 避難所)\n" "A=両規模で安全 / B=両規模で水没 / C=ボーダー(最大のみ) / D=異常") ax.grid(alpha=0.3) ax.legend(loc="lower left", fontsize=10, framealpha=0.92) ax.set_aspect(1 / np.cos(np.deg2rad(34.4))) plt.tight_layout() plt.savefig(PNG_MAP_BORDER, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_MAP_BORDER.name}") # === 16-C. 追加分析8: 5 災害種別フラグ × 浸水域立地のミスマッチ ============== print("\n=== 16-C. 図8: 災害種別フラグ × 浸水域立地ミスマッチ ===") # 「洪水フラグなし」かつ「想定最大規模浸水域内」= 制度上洪水対応外なのに実は浸水域に立地 mismatch = shel_df[(shel_df["flood_flg"] == 0) & (shel_df["in_max"] == 1)].copy() mismatch = mismatch.sort_values("danger_score", ascending=False).reset_index(drop=True) # CSV 保存 (主要列のみ) mismatch_save = mismatch[ ["facility_id", "name", "muni", "lat", "lon", "elev_m", "in_keikaku", "in_max", "border_class", "flood_flg", "sediment_flg", "storm_flg", "eq_flg", "tsunami_flg", "capacity", "danger_score"] ].copy() mismatch_save.to_csv(OUT_MISMATCH, index=False, encoding="utf-8-sig") n_mismatch = len(mismatch) print(f" ミスマッチ件数 (洪水フラグ無 × 浸水域内): {n_mismatch:,}/{len(shel_df):,}" f" ({n_mismatch/len(shel_df)*100:.2f}%)") print(f" → {OUT_MISMATCH.name}") # 5 災害フラグごとに「フラグ無 × 浸水域内」ミスマッチ件数を集計 mismatch_by_flag = pd.DataFrame([ {"災害種別": "洪水", "フラグ列": "flood_flg", "ミスマッチ(フラグ無×浸水内)": int(((shel_df["flood_flg"]==0) & (shel_df["in_max"]==1)).sum()), "対応有×浸水内": int(((shel_df["flood_flg"]==1) & (shel_df["in_max"]==1)).sum())}, {"災害種別": "土砂", "フラグ列": "sediment_flg", "ミスマッチ(フラグ無×浸水内)": int(((shel_df["sediment_flg"]==0) & (shel_df["in_max"]==1)).sum()), "対応有×浸水内": int(((shel_df["sediment_flg"]==1) & (shel_df["in_max"]==1)).sum())}, {"災害種別": "高潮", "フラグ列": "storm_flg", "ミスマッチ(フラグ無×浸水内)": int(((shel_df["storm_flg"]==0) & (shel_df["in_max"]==1)).sum()), "対応有×浸水内": int(((shel_df["storm_flg"]==1) & (shel_df["in_max"]==1)).sum())}, {"災害種別": "地震", "フラグ列": "eq_flg", "ミスマッチ(フラグ無×浸水内)": int(((shel_df["eq_flg"]==0) & (shel_df["in_max"]==1)).sum()), "対応有×浸水内": int(((shel_df["eq_flg"]==1) & (shel_df["in_max"]==1)).sum())}, {"災害種別": "津波", "フラグ列": "tsunami_flg", "ミスマッチ(フラグ無×浸水内)": int(((shel_df["tsunami_flg"]==0) & (shel_df["in_max"]==1)).sum()), "対応有×浸水内": int(((shel_df["tsunami_flg"]==1) & (shel_df["in_max"]==1)).sum())}, ]) print(mismatch_by_flag.to_string(index=False)) # 図: 左パネル = 5災害ミスマッチ件数バー / 右パネル = 洪水ミスマッチを地図にプロット (赤) fig, axes = plt.subplots(1, 2, figsize=(14.5, 7.0), gridspec_kw={"width_ratios": [1.0, 1.4]}) ax = axes[0] xpos = np.arange(len(mismatch_by_flag)) ax.bar(xpos - 0.2, mismatch_by_flag["ミスマッチ(フラグ無×浸水内)"], width=0.4, color="#cf222e", edgecolor="#222", linewidth=0.4, label="フラグ無 × 浸水域内 (潜在リスク)") ax.bar(xpos + 0.2, mismatch_by_flag["対応有×浸水内"], width=0.4, color="#bf8700", edgecolor="#222", linewidth=0.4, label="フラグ有 × 浸水域内 (運用上の盲点)") for i, r in mismatch_by_flag.iterrows(): ax.text(i - 0.2, r["ミスマッチ(フラグ無×浸水内)"] + 8, f"{r['ミスマッチ(フラグ無×浸水内)']:,}", ha="center", fontsize=9, color="#cf222e") ax.text(i + 0.2, r["対応有×浸水内"] + 8, f"{r['対応有×浸水内']:,}", ha="center", fontsize=9, color="#bf8700") ax.set_xticks(xpos) ax.set_xticklabels(mismatch_by_flag["災害種別"], fontsize=10) ax.set_ylabel("浸水域内 (in_max=1) 件数") ax.set_title("図8-(a): 5 災害種別フラグ × 浸水域立地 ミスマッチ\n" "赤=制度外なのに浸水内 / 黄=指定有なのに浸水内") ax.legend(loc="upper right", fontsize=9) ax.grid(axis="y", alpha=0.3) ax = axes[1] # 全避難所を背景に ax.scatter(shel_df["lon"], shel_df["lat"], s=4, c="#cccccc", alpha=0.30, edgecolor="none", label=f"全避難所 (n={len(shel_df):,})") # 洪水フラグ無 × 浸水域内 を赤で強調 ax.scatter(mismatch["lon"], mismatch["lat"], s=18, c="#cf222e", alpha=0.85, edgecolor="#222", linewidth=0.25, label=f"洪水フラグ無 × 浸水域内 (n={n_mismatch:,})") ax.set_xlabel("経度 (°E)") ax.set_ylabel("緯度 (°N)") ax.set_title(f"図8-(b): 「制度上は洪水避難所ではないが浸水域に立地」 {n_mismatch:,} 件の地理分布\n" "= 洪水時に開設指定が無いのに地形的には水没する潜在リスク") ax.grid(alpha=0.3) ax.legend(loc="lower left", fontsize=9) ax.set_aspect(1 / np.cos(np.deg2rad(34.4))) plt.tight_layout() plt.savefig(PNG_MISMATCH, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_MISMATCH.name}") # === 16-D. 追加分析9: 収容力 × 危険度 バブル散布 (市町別) ===================== print("\n=== 16-D. 図9: 収容力 × 危険度 バブル散布 (市町別集約) ===") # capacity 欠損は中央値補完 cap_med = shel_df["capacity"].median(skipna=True) shel_df["capacity_filled"] = shel_df["capacity"].fillna(cap_med) cap_summary = shel_df.groupby("muni").agg( n_shelters=("facility_id", "count"), capacity_total=("capacity_filled", "sum"), capacity_mean=("capacity_filled", "mean"), danger_score_mean=("danger_score", "mean"), danger_score_max=("danger_score", "max"), n_in_max=("in_max", "sum"), elev_m=("elev_m", "first"), ).reset_index() cap_summary["pct_in_max"] = (cap_summary["n_in_max"] / cap_summary["n_shelters"] * 100).round(2) cap_summary["risk_capacity"] = ( cap_summary["capacity_total"] * cap_summary["danger_score_mean"] ).round(1) cap_summary = cap_summary.sort_values("risk_capacity", ascending=False).reset_index(drop=True) cap_summary.to_csv(OUT_CAPSUM, index=False, encoding="utf-8-sig") print(cap_summary.head(10).to_string(index=False)) print(f" → {OUT_CAPSUM.name}") # 単回帰直線 (市町別 capacity_mean vs danger_score_mean) x_cap = cap_summary["capacity_mean"].values.astype(float) y_dg = cap_summary["danger_score_mean"].values.astype(float) slope, intercept = np.polyfit(x_cap, y_dg, 1) corr_cd = float(np.corrcoef(x_cap, y_dg)[0, 1]) print(f" 単回帰: danger_score_mean = {slope:.5f} × capacity_mean + {intercept:.4f}" f" (r={corr_cd:.3f})") fig, ax = plt.subplots(figsize=(11.5, 8.0)) # バブル半径 = sqrt(避難所件数) × 18 radii = np.sqrt(cap_summary["n_shelters"].values) * 18 # 全市町を散布。色 = pct_in_max sc = ax.scatter(x_cap, y_dg, s=radii, c=cap_summary["pct_in_max"].values, cmap="Reds", vmin=0, vmax=max(cap_summary["pct_in_max"].max(), 30), alpha=0.78, edgecolor="#222", linewidth=0.5) cbar = plt.colorbar(sc, ax=ax, shrink=0.78, pad=0.02) cbar.set_label("市町別 浸水域内率 pct_in_max (%)") # 上位 8 市町に名前ラベル top_label = cap_summary.head(8) for _, r in top_label.iterrows(): ax.annotate(r["muni"], (r["capacity_mean"], r["danger_score_mean"]), fontsize=9, alpha=0.95, xytext=(6, 6), textcoords="offset points") # 単回帰直線 xs_line = np.linspace(x_cap.min(), x_cap.max(), 100) ax.plot(xs_line, slope * xs_line + intercept, color="#0969da", lw=1.6, ls="--", label=f"単回帰: y = {slope:.5f}x + {intercept:.3f} (r={corr_cd:.3f})") # 「右上象限」(平均より大 & 平均より高危険) を点線で示す mean_x = x_cap.mean() mean_y = y_dg.mean() ax.axvline(mean_x, color="#888", lw=0.8, ls=":") ax.axhline(mean_y, color="#888", lw=0.8, ls=":") ax.text(mean_x * 1.02, ax.get_ylim()[1] * 0.95, "右上 = 「大きいけど危ない」象限", fontsize=9, color="#444") ax.set_xlabel("市町別 平均 capacity (人, 中央値補完済)") ax.set_ylabel("市町別 平均 danger_score") ax.set_title(f"図9: 収容力 × 危険度 バブル散布 (市町別, n={len(cap_summary)} 市町)\n" "バブル半径 = √避難所件数 × 18 / 色 = 浸水域内率 (%)") ax.legend(loc="upper right", fontsize=9) ax.grid(alpha=0.3) plt.tight_layout() plt.savefig(PNG_CAP_DANGER, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_CAP_DANGER.name}") # === 16-E. 追加分析10: 市町別 危険率ヒートマップ (地理) ======================= print("\n=== 16-E. 図10: 市町別 危険率ヒートマップ (地理) ===") # 市町ごとに重心 (lat 平均, lon 平均) を算出 city_geo = shel_df.groupby("muni").agg( lat_mean=("lat", "mean"), lon_mean=("lon", "mean"), n_shelters=("facility_id", "count"), n_in_max=("in_max", "sum"), ).reset_index() city_geo["risk_rate"] = (city_geo["n_in_max"] / city_geo["n_shelters"] * 100).round(2) city_geo = city_geo.sort_values("risk_rate", ascending=False).reset_index(drop=True) city_geo.to_csv(OUT_CITYRISK, index=False, encoding="utf-8-sig") print(city_geo.head(10).to_string(index=False)) print(f" → {OUT_CITYRISK.name}") fig, ax = plt.subplots(figsize=(11.5, 8.5)) # 背景: 全避難所を薄く ax.scatter(shel_df["lon"], shel_df["lat"], s=2.5, c="#dddddd", alpha=0.30, edgecolor="none") # 市町中心円: 半径 = sqrt(避難所件数) × 5, 色 = 危険率 radii_city = np.sqrt(city_geo["n_shelters"].values) * 30 sc = ax.scatter(city_geo["lon_mean"], city_geo["lat_mean"], s=radii_city, c=city_geo["risk_rate"].values, cmap="Reds", vmin=0, vmax=max(city_geo["risk_rate"].max(), 30), alpha=0.80, edgecolor="#222", linewidth=0.6) cbar = plt.colorbar(sc, ax=ax, shrink=0.78, pad=0.02) cbar.set_label("市町別 危険率 risk_rate = n_in_max / n_shelters (%)") # 上位 10 市町に名前ラベル for _, r in city_geo.head(10).iterrows(): ax.annotate(f"{r['muni']}\n{r['risk_rate']:.0f}%", (r["lon_mean"], r["lat_mean"]), fontsize=8.5, alpha=0.95, xytext=(7, 7), textcoords="offset points", bbox=dict(boxstyle="round,pad=0.15", fc="white", ec="#888", alpha=0.85, lw=0.4)) ax.set_xlabel("経度 (°E)") ax.set_ylabel("緯度 (°N)") ax.set_title(f"図10: 市町別 危険率ヒートマップ (地理) — {len(city_geo)} 市町\n" "円中心 = 市町重心 (避難所平均) / 半径 = √避難所件数 × 30 / 色 = 危険率") ax.grid(alpha=0.3) ax.set_aspect(1 / np.cos(np.deg2rad(34.4))) plt.tight_layout() plt.savefig(PNG_CITY_HEAT, dpi=140, bbox_inches="tight") plt.close() print(f" → {PNG_CITY_HEAT.name}") # === 17. HTML 構築 ========================================================== print("\n=== 17. HTML 組み立て ===") # 表示用テーブル city_disp = city_summary.head(20)[ ["muni", "n_total", "n_in_max", "n_in_keikaku", "n_only_max", "pct_in_max", "n_flood_flg", "elev_m"] ].copy() city_disp.columns = ["市町", "全避難所", "想定最大域内", "計画域内", "ボーダー(C)", "最大域内率(%)", "洪水対応有", "市町平均標高(m)"] city_html = city_disp.to_html(index=False, border=0) danger_disp = danger.head(15)[ ["name", "muni", "elev_m", "in_keikaku", "in_max", "flood_flg", "danger_score"] ].copy() danger_disp.columns = ["避難所名", "市町", "標高(m)", "計画", "最大", "洪水対応", "危険スコア"] danger_disp["計画"] = danger_disp["計画"].map({0: "OUT", 1: "IN"}) danger_disp["最大"] = danger_disp["最大"].map({0: "OUT", 1: "IN"}) danger_disp["洪水対応"] = danger_disp["洪水対応"].map({0: "無", 1: "有"}) danger_html = danger_disp.to_html(index=False, border=0) border_disp = border_breakdown.copy() border_html = border_disp.to_html(index=False, border=0) flag_html = ct.to_html(border=0) track_html = track.to_html(index=False, border=0) verdict_rows = [] for k, (ok, msg) in verdict.items(): verdict_rows.append({"仮説": k, "判定": "支持" if ok else "反証", "根拠": msg}) verdict_html = pd.DataFrame(verdict_rows).to_html(index=False, border=0) # 追加分析 6-10 用 表データ mismatch_disp = mismatch.head(15)[ ["name", "muni", "elev_m", "border_class", "sediment_flg", "storm_flg", "eq_flg", "tsunami_flg", "danger_score"] ].copy() mismatch_disp.columns = ["避難所名", "市町", "標高(m)", "border_class", "土砂", "高潮", "地震", "津波", "危険スコア"] mismatch_html = mismatch_disp.to_html(index=False, border=0) mismatch_flag_html = mismatch_by_flag.to_html(index=False, border=0) cap_disp = cap_summary.head(15)[ ["muni", "n_shelters", "capacity_total", "capacity_mean", "danger_score_mean", "pct_in_max", "risk_capacity"] ].copy() cap_disp.columns = ["市町", "避難所数", "総capacity", "平均capacity", "平均危険スコア", "浸水域内率(%)", "リスク収容人数"] cap_html = cap_disp.to_html(index=False, border=0) cityrisk_disp = city_geo.head(15)[ ["muni", "n_shelters", "n_in_max", "risk_rate", "lat_mean", "lon_mean"] ].copy() cityrisk_disp.columns = ["市町", "避難所数", "浸水域内件数", "危険率(%)", "重心緯度", "重心経度"] cityrisk_html = cityrisk_disp.to_html(index=False, border=0) # 追加分析6-10 用 集計値 (HTML 内で使うため事前計算) n_in_max_total = int(shel_df["in_max"].sum()) mismatch_pct = n_mismatch / max(n_in_max_total, 1) * 100 # 浸水内のうち何 % が「フラグ無」か border_class_n = shel_df["border_class"].value_counts().to_dict() border_A_n = int(border_class_n.get("A 両規模で安全", 0)) border_B_n = int(border_class_n.get("B 両規模で水没", 0)) border_C_n = int(border_class_n.get("C ボーダー (最大規模のみ水没)", 0)) border_D_n = int(border_class_n.get("D 計画のみ水没 (異常)", 0)) risk_rate_top_muni = city_geo.iloc[0]["muni"] risk_rate_top_val = float(city_geo.iloc[0]["risk_rate"]) risk_rate_top_n = int(city_geo.iloc[0]["n_shelters"]) risk_capacity_top = cap_summary.iloc[0] n_city_above50 = int((city_geo["risk_rate"] >= 50).sum()) n_city_above30 = int((city_geo["risk_rate"] >= 30).sum()) sections = [] # --- 1. 学習目標と問い ---------------------------------------------------- sections.append(("学習目標と問い", f"""
【本記事のスタイル: 点 in ポリゴン判定研究】 4,065 避難所のうち どれが 浸水想定域内に立地しているかを 1 件ごとに Ray casting 法で判定し、危険避難所をランキング化する。 X07 が「ポリゴン × ポリゴンの面積率研究」だったのに対し、 本記事は「点 × ポリゴン」という別軸の空間結合。 PCA・散布図行列・ロジ回帰・k-means は使わず、 危険避難所ランキング棒・市町別棒・ボーダー判定棒・フラグ×立地クロス・標高散布 の 5 図を主役に据える。

本レッスンは 避難所点 4,065 件 (DoBoX #42) と 河川浸水想定区域 (計画規模・想定最大規模 各 1 件, DoBoX #295/#313) を 同一の経緯度系 (WGS84) に揃え、 4,065 × 2 規模 = 8,130 件の点判定を行う。 得られた判定は次の 5 つの仮説と照合する。

このレッスンで答えたい問い (1 文)

「広島県内の避難所 4,065 件のうち何件が河川浸水想定区域内に立地しており、 そのうちで 本当に避難してはいけない 危険避難所はどこか？」

立てた仮説 H1〜H5

独自定義の用語 (このレッスン専用 — 要件M)

到達点

  1. Ray casting で 1 点が 1 多角形の内側か外側かを判定できる (= 1 関数 10 行)。
  2. 4,065 点 × 1,029 ポリゴンの総当たり判定を bbox プリフィルタで 秒単位に高速化できる。
  3. 異なる 座標系 (CRS) の点とポリゴンを、共通の経緯度系に揃えることで判定を統一できる。
  4. 「避難所自体が水没する」という反直観的な事実を、4,065 件の数字で確認できる。
  5. 計画規模 vs 想定最大規模の差を 「ボーダー避難所」として可視化し、 通常運用では見えない過去最大級リスクを発見できる。

なぜ「危険スコア」を主指標に選んだか (要件 H/I)

候補指標長所限界判定
(A) in_max のみ (二値)解釈最簡水深と標高の差を区別できない補助指標
(B) 危険スコア (in_max × 標高ペナルティ) ★連続値で順位付け可能標高はシンボリック近似採用 (主指標)
(C) 浸水深 × 標高差物理的に正確浸水深データが本シェープにない不採用
(D) 容量加重スコア政策的に重要capacity 欠損 535/4,065派生指標として採用

(B) を主指標、(A)(D) を補助指標とする。 浸水深の代替として 「両規模で水没 (both_in)」を 1.5 倍重み付けで取り込み、 「常時危険」と「最大規模時のみ危険」を順位上で区別している。

""")) # --- 2. 使用データ -------------------------------------------------------- sections.append(("使用データ", f"""

本レッスンは DoBoX 3 系統を使う。 避難所は 4,065 件すべて、浸水は 39 件のうち 「全河川」集約版 2 件のみを使う (個別水系 37 件は X07 と同じ理由でメタ集計のみ)。

原データ — 避難所情報 (1 件, 4,065 行)

ID名称形式件数主要列
#42避難所情報JSON{len(shel_df):,} name, latitude, longitude, municipalityName, floodShFlg, capacity

原データ — 河川浸水想定区域情報 (39 件のうち 2 件を使用)

ID名称形式規模ポリゴン数
#295計画規模 全河川Shapefile 計画規模 (1/100〜1/30){len(flood_data['計画規模']):,}
#313想定最大規模 全河川Shapefile 想定最大規模 (1/1000){len(flood_data['想定最大規模']):,}

派生データ — 市町別代表標高 (シンボリック)

本記事では浸水深の 絶対値データが Shapefile に含まれないため、 標高との比較は 市町ごとの代表標高 (役所所在地付近, m, 整数概算) を ローカル定数として定義する手法を採る。 {len(CITY_ELEV_M)} 市町を網羅、未定義は {DEFAULT_ELEV} m とフォールバック。 これは シンボリック値であり、避難所 1 件ごとの実標高ではない点に注意 (= 限界)。

サイズ・次元の整理 (要件 L)

段行/列/サイズ役割
原 避難所 JSON{len(shel_df):,} 行 × 30 列 1 行 = 1 避難所 (ID, 名称, lat, lon, 各種フラグ)
原 浸水 Shapefile (計画+最大) {len(flood_data['計画規模'])} + {len(flood_data['想定最大規模'])} = {len(flood_data['計画規模'])+len(flood_data['想定最大規模'])} ポリゴン 1 行 = 1 浸水ポリゴン
判定後 shel_df {len(shel_df):,} 行 × {len(shel_df.columns)} 列 各避難所に in_keikaku, in_max, only_max_in, border_class, danger_score を追加
市町別集計 city_summary {len(city_summary)} 行 × {len(city_summary.columns)} 列 市町ごとに件数集計 (主集計テーブル)
判定総回数 (粗計算) {len(shel_df):,} 点 × ({len(flood_data['計画規模'])+len(flood_data['想定最大規模'])}) ポリゴン = {len(shel_df) * (len(flood_data['計画規模']) + len(flood_data['想定最大規模'])):,} 回 bbox プリフィルタなしの上限

※ 表示の都合で「上位 30」「上位 20 市町」と書くが、 全件は 常に {len(shel_df):,} 避難所 / {len(city_summary)} 市町で集計済み (要件L)。

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

原データ (DoBoX, 直リンク・直 DL)

{data_recipe([ {"label": "避難所情報", "dataset_id": 42, "resource_id": None, "file": "data/shelters.json", "format": "JSON", "size": "~3 MB"}, {"label": "計画規模 全河川", "dataset_id": 295, "resource_id": 23061, "file": "data/extras/flood_shp/shinsui_keikaku.zip", "format": "Shapefile (zip)", "size": "22 MB"}, {"label": "想定最大規模 全河川", "dataset_id": 313, "resource_id": 23118, "file": "data/extras/flood_shp/shinsui_souteisaidai.zip", "format": "Shapefile (zip)", "size": "38 MB"}, ])}

本レッスン生成の中間データ (HTML から直 DL)

図 PNG (HTML から直 DL)

再現スクリプト

X08_flood_shelter_in_polygon.py を以下で実行:

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

初回実行時に data/extras/flood_shp/ へ 浸水 2 件を自動取得する。所要 30〜60 秒 (DL は逐次)。

""")) # --- 4. 分析1: 1件追跡 + Ray casting 解説 (要件K) ----------------------- sections.append(("分析1: 4,065 避難所と 2 規模浸水を経緯度系で揃える (要件K 1件追跡)", f"""

狙い

避難所点 (lat, lon) は WGS84 経緯度で配信されているが、 浸水 Shapefile は EPSG:3857 (Web Mercator) 単位 m で配信されている。 両者を WGS84 経緯度に揃えて重ね合わせる。 本記事では geopandas/pyproj を使わず、球面近似式 1 本だけで変換する。

手法 (要件 B/J: ツール化視点で簡潔に)

STEP 分け (要件 O)

段役割入力出力
STEP1: Ray casting1 点が 1 多角形の内側か(x, y, ring)True/False
STEP2: bbox プリフィルタ外接矩形外を即外側判定点 + ポリゴン bbox候補のみ抽出
STEP3: ベクトル化判定4,065 点を一括処理numpy 配列bool 配列
STEP4: 規模ごとに繰り返し計画/最大 で 2 周flood_data dictin_keikaku, in_max

実装

""" + code(r''' import math import numpy as np WGS_R = 6378137.0 def merc_to_ll(x, y): """Web Mercator (EPSG:3857) → WGS84 経緯度 (厳密式)""" lon = math.degrees(x / WGS_R) lat = math.degrees(2 * math.atan(math.exp(y / WGS_R)) - math.pi/2) return lon, lat def points_in_ring_vec(xs_pt, ys_pt, ring_xs, ring_ys): """N 点が 1 リング (M 頂点) の内側かをベクトル化 Ray casting で判定。""" n_pts = len(xs_pt) inside = np.zeros(n_pts, dtype=bool) m = len(ring_xs) j_idx = (np.arange(m) - 1) % m for i in range(m): yi, yj = ring_ys[i], ring_ys[j_idx[i]] xi, xj = ring_xs[i], ring_xs[j_idx[i]] cond_y = (yi > ys_pt) != (yj > ys_pt) # y 方向の交差検査 if not cond_y.any(): continue x_int = (xj-xi) * (ys_pt-yi) / (yj-yi+1e-15) + xi inside ^= cond_y & (xs_pt < x_int) # XOR で内外を反転 return inside ''') + f"""

結果 (表と読み取り) — {target_muni} 1 件追跡 (要件K)

「JSON 取得 → 1 件選択 → 経緯度抽出 → 計画判定 → 最大判定 → 標高 → ボーダー区分 → フラグ確認 → 危険スコア」の 9 段階 を、{target['name']} 1 件で具体値を追って示す。

表1: {target['muni']} {target['name']} 段階表 (Before → After)

{track_html}

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

""")) # --- 5. 分析2: 点 in ポリゴン判定 = Ray casting + bbox 高速化 ----------- sections.append(("分析2: 点 in ポリゴン判定 (Ray casting + bbox プリフィルタ) (要件J)", f"""

狙い

避難所点 4,065 件が 計画規模 {len(flood_data['計画規模'])} ポリゴンと 想定最大規模 {len(flood_data['想定最大規模'])} ポリゴンのどれかに入っているかを、 ライブラリ (geopandas/shapely) を使わずに 純 Python + numpy だけで判定する。

手法 (要件B/J)

直感的説明 — Ray casting

判定したい点 P から右方向に半直線 (光線, ray) を伸ばし、ポリゴン外周との 交差回数を数える。 奇数なら P は内側、偶数 (0 含む) なら外側。 これは「外周を 1 回横切るたびに内→外、または外→内」という単純な原理に基づく。 凹多角形でも穴がない限り正しく動作する。

アルゴリズム (4 ステップ)

  1. 避難所点を numpy 配列化: xs = np.array([s['lon'] for s in shelters]) 等。
  2. 規模ごとにポリゴンをループ: 計画 → 最大 の 2 周。
  3. 各ポリゴンで bbox プリフィルタ: ポリゴン bbox 外の点は確定で外側 → Ray casting スキップ。
  4. 残った候補に対してリング単位で Ray casting: 1 ポリゴンに複数のリング (parts) がある 場合は 奇偶 XOR で内外を集約。

入出力の形 (=このツールで何ができるか)

計算量と高速化

限界 (=このツールで何ができないか)

実装

""" + code(r''' def points_in_any_polygon_set(xs_pt, ys_pt, polys): """N 点 vs P ポリゴン集合 → bool 配列。bbox プリフィルタ付き。""" n_pts = len(xs_pt) inside = np.zeros(n_pts, dtype=bool) for p in polys: # bbox プリフィルタ (このポリゴンの bbox 外なら確定で外) cand = (~inside) & \ (xs_pt >= p["lon_min"]) & (xs_pt <= p["lon_max"]) & \ (ys_pt >= p["lat_min"]) & (ys_pt <= p["lat_max"]) if not cand.any(): continue idx = np.where(cand)[0] sub_x, sub_y = xs_pt[idx], ys_pt[idx] sub_in = np.zeros(len(idx), dtype=bool) for ring_x, ring_y in p["rings"]: still = ~sub_in if not still.any(): break r = points_in_ring_vec(sub_x[still], sub_y[still], ring_x, ring_y) sub_in[still] |= r inside[idx] |= sub_in return inside # 4,065 避難所 × 2 規模を一括判定 xs = shel_df["lon"].values.astype(np.float64) ys = shel_df["lat"].values.astype(np.float64) shel_df["in_keikaku"] = points_in_any_polygon_set(xs, ys, flood_data["計画規模"]).astype(int) shel_df["in_max"] = points_in_any_polygon_set(xs, ys, flood_data["想定最大規模"]).astype(int) ''') + f"""

結果 (表と読み取り) — 4,065 × 2 規模 判定サマリ

規模浸水域内件数占有率処理時間 (参考)
計画規模{int(shel_df['in_keikaku'].sum()):,} {shel_df['in_keikaku'].mean()*100:.1f}%~5 秒
想定最大規模{int(shel_df['in_max'].sum()):,} {shel_df['in_max'].mean()*100:.1f}%~10 秒

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

""")) # --- 6. 分析3: 主役図 危険避難所マルチ属性ランキング ------------------------ _top30 = danger.head(30).copy() _top30_n_floodless = int((_top30["flood_flg"] == 0).sum()) _top30_n_lowland = int((_top30["elev_m"] <= 5).sum()) _top30_n_allflag = int((_top30[["flood_flg", "sediment_flg", "storm_flg", "eq_flg", "tsunami_flg"]].sum(axis=1) == 5).sum()) _top30_cap_total = int(_top30["capacity"].fillna( shel_df["capacity"].dropna().median()).sum()) _top30_n_munis = int(_top30["muni"].nunique()) sections.append(("分析3: 主役図 — 危険避難所 上位 30 マルチ属性ランキング (要件H)", f"""

狙い

4,065 件のうち想定最大規模浸水域内に立地する {int(shel_df['in_max'].sum())} 件をすべて並べると 横棒が長くなりすぎる。「危険スコア」で 上位 30を抽出し、1 枚に複数属性を圧縮表示する。 本記事の 主役図。学習者は他の図を見る前に、この図で 「具体的にどの避難所が一番危なく、なぜそれが危ないのか (低地? 大規模? 洪水避難所未指定?)」を地名レベルで把握する。

旧版 (赤濃淡だけ) の問題点

旧 図1 は danger_score をバー長と赤濃度で表現したが、 上位 30 はほぼ全件が score > 0.85 の「高スコア帯」に集中するため バーがほぼ同じ長さ・同じ赤色になり、視覚的に区別がつきにくかった。 1 枚に 1 指標では情報密度が低い。

新版: マルチ属性化 (1 枚に 5 軸)

視覚チャネルマッピング先狙い
バー長 (x 軸)capacity (収容人数)「同じ危険でも何人巻き込まれるか」を即座に把握
バー色 (カテゴリ)市町 (tab20 カラー)同じ市町の連続出現 = 「危険集中エリア」が一目で分かる
右マーカー ★flood_flg = 0 (洪水避難所未指定)「制度上は洪水避難所ではないのに浸水域内」= 運用矛盾の典型
右マーカー ▲elev_m ≦ 5m (低地)河口・デルタ平野立地の即時識別
右マーカー ●5 災害フラグ全立 (複合避難所)「全災害対応」を看板に掲げた避難所が河川氾濫で水没する皮肉
右数値テキスト標高 / capacity / score具体値の確認

なぜマルチ属性が必要か (4 案比較)

方式長所限界判定
(A) 旧版: 赤濃淡 score バーシンプル, 順位は明快上位30 がほぼ同色 → 識別不能不採用 (情報密度不足)
(B) capacity バー × 市町色 × 補助マーカー ★5 軸を 1 枚で表示, 市町集中も即見えるマーカー記号の凡例が必要採用 (新主役)
(C) 散布図 (標高×score)2 変数同時名前ラベルが重なる, 順位感がない不採用
(D) 並列ボード 3 パネル ★各属性を分離して比較1 枚に集約しないので「全体感」がない図1B として採用 (補助)

実装

""" + code(r''' import matplotlib.pyplot as plt import numpy as np danger = shel_df[shel_df["in_max"] == 1].sort_values("danger_score", ascending=False) top = danger.head(30).iloc[::-1].reset_index(drop=True) ypos = np.arange(len(top)) # capacity 欠損は中央値補完 (バー長 0 を回避) cap_med = float(shel_df["capacity"].dropna().median()) top["cap_disp"] = top["capacity"].fillna(cap_med).astype(float) # 市町 → カテゴリカル色 unique_munis = list(top["muni"].unique()) cmap_muni = plt.get_cmap("tab20") muni_color = {m: cmap_muni(i % 20) for i, m in enumerate(unique_munis)} bar_colors = [muni_color[m] for m in top["muni"]] FLAG_COLS = ["flood_flg", "sediment_flg", "storm_flg", "eq_flg", "tsunami_flg"] top["n_hazard"] = top[FLAG_COLS].sum(axis=1) fig, ax = plt.subplots(figsize=(13, 10)) ax.barh(ypos, top["cap_disp"], color=bar_colors, edgecolor="#222", linewidth=0.5) xmax = top["cap_disp"].max() for i, (_, r) in enumerate(top.iterrows()): barlen = r["cap_disp"] if r["flood_flg"] == 0: ax.text(barlen + xmax*0.015, ypos[i], "★", va="center", ha="center", fontsize=14, color="#cf222e", fontweight="bold") if r["elev_m"] <= 5: ax.text(barlen + xmax*0.05, ypos[i], "▲", va="center", ha="center", fontsize=12, color="#bf5700") if r["n_hazard"] == 5: ax.text(barlen + xmax*0.085, ypos[i], "●", va="center", ha="center", fontsize=12, color="#1a7f37") ax.text(barlen + xmax*0.13, ypos[i], f"標高{r['elev_m']:.0f}m {int(r['cap_disp']):,}人 " f"score={r['danger_score']:.2f}", va="center", fontsize=8.2, color="#333") ''') + f"""

結果 (図と読み取り) — 主役図

{figure("assets/X08_danger_rank.png", "図1 (主役): 危険避難所 上位 30 マルチ属性ランキング — バー長=収容力 / バー色=市町 / ★=洪水フラグ無 / ▲=標高≦5m / ●=5災害フラグ全立")}

表2: 危険避難所 上位 15 (危険スコア降順)

{danger_html}

この図と表から読み取れること (5 軸を統合した読み):

図1B (補助): 上位30 を 3 軸並列展開

図1 は 1 枚に 5 軸を圧縮しているため、属性ごとの比較には別の見方が要る。 同じ上位 30 の 行並びを揃えて 3 パネルに分けると、 「順位は高いが capacity は小さい」「順位は中位だが標高 0m」「災害フラグが洪水のみ」 といった 属性間の不一致が一目で分かる。

{figure("assets/X08_danger_rank_panels.png", "図1B (補助): 上位30 を 3 パネル並列 (左=capacity[市町色] / 中=標高[低いほど赤] / 右=5災害フラグ heatmap)")}

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

""")) # --- 7. 分析4: 市町別ランキング棒 ----------------------------------------- sections.append(("分析4: 市町別 危険避難所件数ランキング (H5 検証)", f"""

狙い

図1 は個別避難所の順位だが、政策単位は市町。 市町ごとに「危険避難所が何件あるか」を集計し、上位 20 を棒グラフで並べる。 これにより 市町行政の備えるべき件数が明示される。

手法

結果 (表と読み取り)

表3: 市町別 集計 上位 20

{city_html}

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

結果 (図と読み取り)

{figure("assets/X08_city_danger_bar.png", "図2: 市町別 危険避難所件数 上位 20 (赤=想定最大, 青=計画)")}

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

""")) # --- 8. 分析5: ボーダー判定 4 区分 --------------------------------------- sections.append(("分析5: ボーダー避難所 4 区分内訳 (H3 検証)", f"""

狙い

計画規模 vs 想定最大規模の 判定差を 4 区分に分解する: A (両方安全), B (両方水没), C (ボーダー = 計画OK / 最大NG), D (異常 = 計画NG / 最大OK)。 本セクションの主役は 区分 C = 「過去最大級降雨で初めて使えなくなる」避難所。

4 区分の意味

区分計画規模想定最大規模意味政策含意
AOUTOUT両規模で安全常時開設可能
BININ両規模で水没洪水対応指定取り消し検討
COUTINボーダー過去最大級時のみ代替必要
DINOUT異常 (理論上稀)判定エラーまたは縮小済み堤防

手法

結果 (表と読み取り)

表4: 4 区分内訳

{border_html}

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

結果 (図と読み取り)

{figure("assets/X08_border_breakdown.png", "図3: 4 区分内訳 (左) と 想定最大規模水没のボーダー比率 (右)")}

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

""")) # --- 9. 分析6: フラグ × 立地 クロス + 標高散布 --------------------------- sections.append(("分析6: 洪水対応フラグ × 立地 クロス、標高 vs 立地 (H4 検証)", f"""

狙い

避難所 JSON には floodShFlg という「洪水時に開設対象」を示す 0/1 フラグがある。 このフラグと 点 in ポリゴン判定の整合性を検証する。 仮説 H4: フラグ有 → 立地が安全 (= 浸水域内率が低い) はず。

手法

図4 (フラグ × 立地 クロス)

図5 (標高 vs 立地)

実装

""" + code(r''' ct = pd.crosstab( shel_df["flood_flg"].map({0: "洪水対応無", 1: "洪水対応有"}), shel_df["in_max"].map({0: "浸水域外", 1: "浸水域内"}), margins=True, margins_name="合計", ) def expected_chi2(ct): obs = ct.iloc[:-1, :-1].values.astype(float) row, col, tot = obs.sum(1), obs.sum(0), obs.sum() exp = np.outer(row, col) / tot return float(((obs - exp)**2 / exp).sum()), exp chi2, _ = expected_chi2(ct) print(f"χ² = {chi2:.2f} (df=1, p<0.001 閾値 ≈ 10.83)") ''') + f"""

結果 (表と読み取り)

表5: 洪水対応フラグ × 想定最大規模立地 (2×2 クロス)

{flag_html}

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

結果 (図と読み取り) — 図4 フラグ × 立地

{figure("assets/X08_flag_cross.png", "図4: 洪水対応フラグ × 想定最大規模立地 2×2 クロスヒートマップ")}

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

結果 (図と読み取り) — 図5 標高 vs 立地

{figure("assets/X08_elev_scatter.png", "図5: 立地標高 vs 想定最大規模 域内/域外 (左: ジッタ散布, 右: 箱ひげ)")}

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

""")) # --- 9-A. 追加分析7: 危険度の地理マッピング ------------------------------ sections.append(("分析7: 危険度の地理マッピング (4,065 点を danger_score で色塗り)", f"""

狙い

分析3〜6 までは「ランキング棒」「集計棒」「クロス表」「ジッタ散布」と 非地理的な可視化だった。しかし「避難所がどこに立地するか」は本質的に地理問題。 4,065 件の (lat, lon) を直接散布し、危険度の連続値 (danger_score) で色塗りすれば、 沿岸/内陸/山間の分布パターンが一目で読める。

なぜこの図か (要件H, 4 案比較)

方式長所限界判定
(A) 全 {len(shel_df):,} 点 単色散布地理分布のみ把握danger 強弱が出ない不採用
(B) 浸水域内のみ danger_score で色塗り ★「どこが危険か」が連続色で読める背景に全体感が要る採用
(C) ヒートマップ (KDE)密度を平滑化個別避難所が消える不採用
(D) ベース地図 + マーカー市町境界が見える外部地図ライブラリ依存不採用 (純matplotlib縛り)

手法

実装

""" + code(r''' import numpy as np import matplotlib.pyplot as plt mask_in = shel_df["in_max"] == 1 fig, ax = plt.subplots(figsize=(11.5, 8.5)) # 背景: 浸水域外を薄く ax.scatter(shel_df.loc[~mask_in, "lon"], shel_df.loc[~mask_in, "lat"], s=6, c="#cccccc", alpha=0.35, edgecolor="none") # 前景: 浸水域内のみ danger_score で色塗り sc = ax.scatter(shel_df.loc[mask_in, "lon"], shel_df.loc[mask_in, "lat"], s=22, c=shel_df.loc[mask_in, "danger_score"], cmap="RdYlGn_r", vmin=0.0, vmax=1.5, alpha=0.85, edgecolor="#222", linewidth=0.25) plt.colorbar(sc, ax=ax, label="danger_score (0=安全 → 1.5=最大危険)") ax.set_aspect(1 / np.cos(np.deg2rad(34.4))) ''') + f"""

結果 (図と読み取り)

{figure("assets/X08_map_danger.png", "図6: 危険度の地理マッピング — 浸水域内のみ danger_score で色塗り (RdYlGn 反転)")}

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

""")) # --- 9-B. 追加分析8: border_class 別 地理マッピング ------------------------ sections.append(("分析8: border_class 別 地理マッピング (A/B/C 沿岸-内陸-山間パターン)", f"""

狙い

分析5 (図3) で 4 区分 A/B/C/D の件数を棒グラフで見たが、地理的にどう散らばるかは分からなかった。 本図では同じ 4,065 点を border_class で色分けして 1 枚の地図に重ね、 「B (常時危険) は河口部に / C (ボーダー) は中流部に / A (安全) は山間に」 というクラスタが 目視で読み取れるかを検証する。

なぜこの図か (要件H)

図6 (前節) は連続値 danger_score を色塗りしたため 「程度」は分かるが、 「種類」(両規模水没 vs ボーダー) は区別できない。一方、図3 は 4 区分件数だけで 地理パターンが消える。両者の中間として、離散カテゴリ × 地理座標を重ねる本図を選んだ。

手法

結果 (図と読み取り)

{figure("assets/X08_map_border.png", "図7: border_class 別 地理マッピング — A=緑(背景) / B=赤(両規模水没) / C=黄(ボーダー) / D=灰(異常)")}

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

""")) # --- 9-C. 追加分析9: 5 災害種別フラグ × 浸水域立地ミスマッチ ------------- sections.append(("分析9: 5 災害種別フラグ × 浸水域立地のミスマッチ (制度 vs 地形)", f"""

狙い

避難所 JSON には 5 種類の災害対応フラグがある (洪水/土砂/高潮/地震/津波)。 分析6 では 洪水フラグ × 浸水域内のみ見たが、「洪水フラグなし」かつ「浸水域内」の避難所は 制度上は洪水避難所ではないが、実は地形的に浸水域に立地する潜在リスク群。 さらに 5 災害分を一括比較すれば、フラグ運用の盲点が見える。

なぜこの図か (要件H)

図4 (フラグ × 立地 2×2 クロス) は洪水フラグだけの分析。 本節では 5 災害種別を横並びで比較し、さらに 制度外なのに浸水内の {n_mismatch:,} 件を 地図にプロットすることで「制度上の指定」と「地形上のリスク」の乖離を可視化する。

手法

結果 (表と読み取り) — 5 災害種別フラグ集計

表6: 5 災害種別フラグ × 想定最大規模浸水域内立地

{mismatch_flag_html}

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

表7: 洪水フラグ無 × 浸水域内 ミスマッチ 上位 15 件 (危険スコア降順)

{mismatch_html}

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

結果 (図と読み取り)

{figure("assets/X08_mismatch.png", "図8: (a) 5 災害種別フラグ × 浸水内 件数 / (b) 洪水フラグ無 × 浸水内の地理分布")}

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

""")) # --- 9-D. 追加分析10: 収容力 × 危険度 バブル散布 ------------------------- sections.append(("分析10: 収容力 × 危険度 バブル散布 (市町別 — 大きいけど危ない避難所の特定)", f"""

狙い

これまでの分析は「件数」中心。しかし政策的には 「何人が避難できなくなるか」が本質。 本節では capacity (収容人数) を取り込み、市町別に capacity × danger_score を集計。 2 軸散布図で 「右上象限」 = 大きいけど危ない 市町を特定する。

なぜこの図か (要件H)

方式長所限界判定
(A) 棒グラフ (件数のみ)シンプル容量を反映しない図2 で既出
(B) 散布 (capacity × danger_score) ★2 変数同時 + 回帰直線市町数 = {len(cap_summary)} で点が混雑採用
(C) 散布 (避難所 1 件 × danger)4,065 点全部市町判別不能不採用
(D) 棒グラフ (リスク収容人数)1 変数で順位付けcapacity と danger の関係が消える表で代替

手法

実装

""" + code(r''' import numpy as np import matplotlib.pyplot as plt cap_med = shel_df["capacity"].median(skipna=True) shel_df["capacity_filled"] = shel_df["capacity"].fillna(cap_med) cap_summary = shel_df.groupby("muni").agg( n_shelters=("facility_id", "count"), capacity_total=("capacity_filled", "sum"), capacity_mean=("capacity_filled", "mean"), danger_score_mean=("danger_score", "mean"), n_in_max=("in_max", "sum"), ).reset_index() cap_summary["pct_in_max"] = cap_summary["n_in_max"] / cap_summary["n_shelters"] * 100 cap_summary["risk_capacity"] = ( cap_summary["capacity_total"] * cap_summary["danger_score_mean"] ) # 単回帰 x = cap_summary["capacity_mean"].values y = cap_summary["danger_score_mean"].values slope, intercept = np.polyfit(x, y, 1) r = np.corrcoef(x, y)[0, 1] ''') + f"""

結果 (表と読み取り)

表8: リスク収容人数 ランキング 上位 15 (capacity_total × 平均危険スコア)

{cap_html}

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

結果 (図と読み取り)

{figure("assets/X08_capacity_danger.png", "図9: 収容力 × 危険度 バブル散布 (市町別) — 右上象限 = 大きいけど危ない")}

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

""")) # --- 9-E. 追加分析11: 市町別 危険率ヒートマップ (地理) ------------------- sections.append(("分析11: 市町別 危険率ヒートマップ (地理) (危険率 高い市町を一目で)", f"""

狙い

分析4 (図2) の市町別棒グラフは 件数 での順位だが、件数大 = 母集団大 のことが多い。 占有率 (危険率 = n_in_max / n_shelters) の方が 地形リスクの本質を映す。 本節では 占有率を地図上の円の色、避難所件数を円の半径として、 両指標を 1 図に圧縮する。

なぜこの図か (要件H)

方式長所限界判定
(A) 市町別 危険率 棒グラフ順位明快地理が消える図2 で件数版あり
(B) 地理ヒートマップ ★色=率, 半径=件数, 1 図で 2 指標市町境界が出ない採用
(C) コロプレス (市町ポリゴン色塗)正確な市町境界市町境界 Shapefile が必要不採用 (本記事範囲外)
(D) 散布図 (件数 × 危険率)2 軸明示地理が消える図9 で類似版あり

手法

結果 (表と読み取り)

表9: 市町別 危険率 ランキング 上位 15

{cityrisk_html}

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

結果 (図と読み取り)

{figure("assets/X08_city_heatmap.png", "図10: 市町別 危険率ヒートマップ (地理) — 円の色=危険率 / 円の半径=√避難所件数")}

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

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

仮説判定表

{verdict_html}

主要発見の物語化

  1. 「避難所自体が水没する」は実証された (H1 {('支持' if h1_supported else '反証')}, {pct_in_max:.1f}%)。 広島県の避難所 4,065 件のうち {n_in_max:,} 件 ({pct_in_max:.1f}%) が想定最大規模浸水想定区域内に立地。 これは {int(round(100/max(pct_in_max,0.01)))} 件に 1 件の頻度で「避難所自体が水没する」ことを意味する。 学校や公民館は住宅地中心に置かれるため、住宅地が浸水する地理を反映した結果。
  2. 沿岸都市集中 (H2 {'支持' if h2_supported else '反証'}, 沿岸シェア {coastal_share:.1f}%)。 {('上位 30 危険避難所のうち過半数が沿岸 5 市町に集中。河川氾濫はデルタ平野・河口部で最も広範囲になるため自然な結果。' if h2_supported else '沿岸/内陸が混在しており、沿岸集中の単純化は不正確。河川中流の三次市・庄原市にも危険避難所が分散している。')}
  3. ボーダー避難所 = {n_border} 件 (H3 {('支持' if h3_supported else '反証')})。 計画規模では安全だが想定最大規模で水没する避難所が {n_border} 件存在。 これらは 「常時閉鎖」ではなく「警報時切替」という運用設計が必要。 2018 年西日本豪雨級の事象で初めて顕在化するリスク。
  4. 洪水対応フラグの整合性 (H4 {('支持' if h4_supported else '反証')})。 {('フラグ有の方が浸水域内率が低く、行政の指定方針は機能している。ただし「フラグ有 × 浸水域内」'+str(int(ct.loc['洪水対応有', '浸水域内']))+' 件は運用上の盲点。' if h4_supported else 'フラグ有の方が浸水域内率が逆に高い ('+f'{flag_in_rate*100:.1f}% > {noflag_in_rate*100:.1f}%'+')。これは「洪水対応指定 = 洪水時に開設」が必ずしも「浸水しない安全な場所」を意味しないことを示す。指定の意味の再定義が必要。')}
  5. 避難所数と危険件数の相関 (H5 r={corr_h5:.3f}, {('支持' if h5_supported else '反証')})。 {('避難所が多い市町ほど危険件数も多い = 母集団効果が支配。占有率では市町ランキングが大きく変わるため、件数と占有率を併記する必要がある。' if h5_supported else '避難所数と危険件数の単純な比例は成立しない。市町ごとに地形条件が大きく違い、件数では測れない地形リスクが存在。')}

方法論的限界 (この記事を「鵜呑みにしない」ためのチェック)

主要発見 (1 段落要約)

広島県の避難所 4,065 件のうち {pct_in_max:.1f}% ({n_in_max:,} 件) が想定最大規模浸水想定区域内に立地し、 そのうち {n_border} 件は計画規模では安全な「ボーダー避難所」である。 危険避難所の上位 30 件は標高 30m 以下に集中し、{', '.join(danger.head(30)['muni'].value_counts().head(3).index.tolist())} の沿岸都市群に偏在する。 洪水対応フラグ有の浸水域内率は {flag_in_rate*100:.1f}% で、フラグ無 ({noflag_in_rate*100:.1f}%) より{('低く、指定方針は機能している' if h4_supported else '高く、指定の意味の再定義が必要')}。 χ² = {chi2_val:.1f} で {'統計的に有意' if chi2_val > 10.83 else '独立に近い'}。 「避難所自体が水没するか」という反直観的問いに、4,065 件の点 in ポリゴン判定で具体地名レベルの答えが出た。

""")) # --- 11. 発展課題 ----------------------------------------------------------- sections.append(("発展課題 (結果X → 新仮説Y → 課題Z)", f"""

本レッスンの結果から、次の 5 つの発展課題が論理的に導かれる:

課題1: 避難所 1 件ごとの実標高で危険スコアを精緻化

課題2: 浸水深カテゴリを取り込んだリスクスコア

課題3: 容量加重リスク (= 何人が避難できなくなるか)

課題4: 高潮浸水・津波浸水との合成リスク

課題5: 危険避難所からの「最近傍安全避難所」距離

""")) # === HTML 書き出し =========================================================== data_label = ( f"避難所 {len(shel_df):,} 件 (DoBoX #42 JSON) + 浸水想定 計画規模 " f"{len(flood_data['計画規模'])} ポリゴン + 想定最大規模 " f"{len(flood_data['想定最大規模'])} ポリゴン (DoBoX 直 DL)" ) html = render_lesson( num=80, title="X08: 河川浸水想定区域 × 避難所立地 — 「避難所自体が水没するか?」 点 in ポリゴン判定研究", tags=["X-tier", "点 in ポリゴン判定研究", "Ray casting", "bbox プリフィルタ", "避難所", "浸水想定区域", "危険避難所", "L01"], time="60-90分", level="リテラシレベル (L01)", data_label=data_label, sections=sections, script_filename="X08_flood_shelter_in_polygon.py", ) # data-draft="1" data-stier="X" を に挿入 html = html.replace( '
', '
', ) out_html = LESSONS / "X08_flood_shelter_in_polygon.html" out_html.write_text(html, encoding="utf-8") print(f"\n→ {out_html}") print("\n=== 主要発見 ===") print(findings_text) print("\n=== 完了 ===")