# -*- coding: utf-8 -*- """L72 緊急輸送道路 単独 3 研究例分析 — 広島県 緊急輸送道路 (dataset 1247) を 3 角度で解読 カバー宣言: 本記事は DoBoX のデータセット「緊急輸送道路」 (dataset_id = 1247) 1 件を 単独で取り上げ、 広島県の緊急輸送道路 (LineString, 4 階層 / 計 ~2,790 km) を 3 つの独立した研究角度 (RQ1 / RQ2 / RQ3) で並列に分析する。 「緊急輸送道路」 とは: 災害発生直後から救助・救急・医療・消火活動および緊急物資・人員等の 輸送のために緊急車両が通行する重要な道路。 都道府県地域防災計画で 指定され、 災害時の「生命線 (ライフライン)」として通行確保が 最優先される。 道路法上の道路 (高速・国道・県道・市町道) のうち、 防災拠点 (病院・警察・消防本部・自衛隊基地・港湾・空港・物資集積所) を相互に結ぶネットワークが指定対象。 広島県の緊急輸送道路は4 階層に区分されている: 第 1 次: 高速自動車国道・主要幹線国道 (国土幹線) 第 2 次: 主要地方道・一般国道 (県内幹線) 第 3 次: 一般県道 + 市町道 (補完路線) 補完線: 上記を補強する追加路線 総延長は約 2,790 km で、 県内道路網の約 5-10% が緊急輸送道路に 指定されている (推定)。 本記事は 橋梁 L66 / トンネル L67 / 道路法面 L71 / 道路規制 L50 / ハザード L11 と厳密に区別: L66 = 橋梁単独 4,203 件 (POINT, 中小河川クロス) L67 = トンネル単独 157 件 (POINT, 山岳貫通) L71 = 道路法面単独 12 区間 / 24 写真 (POINT, 沿道斜面保護) L50 = 道路規制単独 (今日 + 将来 2 dataset) L11 = トリプルハザード単独 (急傾斜 + 土石流 + 地すべり) L72 = 緊急輸送道路単独 (LineString, 災害時生命線, 4 階層) 本記事は他シリーズと合体しない。 RQ2/RQ3 で既扱データ (L11 警戒区域 + L66 橋梁 + L67 トンネル) を参照する形をとる。 研究の問い (3 RQ): RQ1 (主研究): 広島県の緊急輸送道路の構造 — 階層・延長・地理はどう描けるか? 4 階層 (1 次 / 2 次 / 3 次 / 補完) × 計 ~2,790 km の LineString ネットワークをセグメント数 × 延長 × 市町別カバー × 地理クラス (沿岸島嶼 / 平野都市 / 中山間山地)の 4 軸で集計し、 「県の 生命線ネットワーク」 の物理的形状を初めて系統的に記述する。 H1 = 第 2 次が県内総延長の過半数 (≥ 50%) を占める仮説。 RQ2 (副研究 1): 緊急輸送道路は災害脆弱性 — 浸水想定 + 土砂警戒区域 (L11 既知) にどれだけ晒されているか? 想定最大規模の浸水区域 (L8/L11 系の Shapefile) と土砂警戒区域 3 種 (急傾斜 / 土石流 / 地すべり) に対して、 緊急輸送道路の総延長のうち何 km / 何 % が重なるかを点サンプリング (50m 間隔) で計算し、 「災害時に 通行不能となる可能性のある区間」 を同定する。 H2 = 第 1 次は浸水・土砂双方で重なり率が低い (<10%) 仮説、 H3 = 全 4 階層を通じての浸水重なり率が土砂より高い仮説。 RQ3 (副研究 2): 緊急輸送道路上の橋梁・トンネル (L66/L67 既知) と 老朽橋の関係はどう描けるか? 緊急輸送道路 30m バッファ内に 入る橋梁数 + トンネル数を空間結合で同定し、 さらに老朽橋 (架設 1970 年以前)のうち緊急輸送道路上の件数を集計、 「BCP (事業継続計画) 脆弱箇所」 を実証的に同定する。 H4 = 緊急輸送道路上に橋梁の 30% 以上 + トンネルの 50% 以上 が乗る仮説、 H5 = 緊急輸送道路上の老朽橋密度が一般道より 同程度以上 (= BCP 脆弱箇所が定量化される) 仮説。 仮説 (5): H1 (RQ1, 階層延長): 緊急輸送道路 4 階層のうち第 2 次が総延長の 過半数 (≥ 50%)を占める。 これは第 1 次が高速道路・幹線国道に 限られ短距離、 第 2 次が県内主要地方道として最も長距離を占める という県管理道路の階層構造を反映する仮説。 H2 (RQ2, 第 1 次の防災優位): 第 1 次緊急輸送道路は (高速道路 中心のため高架・トンネルが多く) 浸水・土砂警戒区域との重なり率が 両方とも10% 未満。 これは第 1 次が「災害時最優先確保」 の 設計思想で整備された結果である仮説。 H3 (RQ2, 浸水 vs 土砂): 緊急輸送道路全体で、 浸水想定区域の重なり率は 土砂警戒区域 (3 種合計) の重なり率より高い。 これは平野部の 国道・県道が河川氾濫域を通過する一方、 土砂警戒区域は山地局所に 限られるためという地形的根拠の仮説。 H4 (RQ3, 橋梁・トンネル集中): 県内の橋梁 4,203 のうち30% 以上 (≥ 1,260 件)、 トンネル 157 のうち50% 以上 (≥ 78 件)が 緊急輸送道路 30m バッファ内に入る。 これは緊急輸送道路 = 主要 道路という性質上必然的だが、 「災害時にこれだけの構造物の機能を 維持する必要がある」 という BCP 規模の定量化となる仮説。 H5 (RQ3, 老朽橋の BCP リスク): 緊急輸送道路上の老朽橋 (架設 1970 年以前)は200 件以上存在し、 これは「災害時に最優先で 通行確保すべき道路上に、 耐震基準改定 (1980) 以前の橋梁が大量に 残存する」 という BCP 脆弱箇所の量的実証となる仮説。 要件 S 準拠 (1 分以内完走): - データ ZIP 1 件 (~390 KB) → 既存ローカル展開済 - 4 JSON 線データ → 630 LineString → ~2,790 km - 緊急輸送道路を 50m 間隔で点サンプル (~56,000 pts) → sjoin 高速判定 - L11 警戒区域 (既キャッシュ Shapefile) は polygon-in-point sjoin 1 秒 - L66/L67 既キャッシュ CSV 流用 (4,203 + 157 件) - 重い前処理 (浸水想定 Shapefile) は最初の 1 回だけ EPSG:6671 投影し GPKG キャッシュ → 2 回目以降は 1 秒で読込 - 本スクリプト全体で ~30 秒目標 要件 T 準拠 (位置情報あり = 地図必須): - RQ1: 階層別マップ + 階層別延長 棒 + 市町別 choropleth - RQ2: 浸水想定 + 土砂警戒区域 重ね合わせマップ + 階層別重なり率 - RQ3: 橋梁 + トンネル 重ね合わせマップ + 老朽橋集中エリア地図 要件 Q 準拠: 図 8+ / 表 12+ (3 RQ × 多角度: 構造 / 災害脆弱性 / 構造物リスク) データ仕様: - dataset 1247: 緊急輸送道路 (ZIP × 1 リソース) - resource 32491: 道路_緊急輸送道路_2022-09-26-T00:00:00 (~388 KB) - ZIP 内: 5 JSON - 05_kinkyu_route.json: 階層メタ (route_01〜06 + 00, color/weight/type) - 05_kinkyu_route_01.json: 第 1 次 (52 セグ / 470 km) - 05_kinkyu_route_02.json: 第 2 次 (379 セグ / 1,696 km) - 05_kinkyu_route_03.json: 第 3 次 (184 セグ / 599 km) - 05_kinkyu_route_04.json: 補完 (15 セグ / 25 km) - 各 LineString は { "e": 経度, "d": 緯度 } の点列 (NDJSON 風) - ライセンス: クリエイティブ・コモンズ表示 (CC-BY) メモリ対策: Figure ごとに plt.close('all') で確実に解放。 """ from __future__ import annotations import sys, time, zipfile, json from pathlib import Path sys.path.insert(0, str(Path(__file__).parent)) from _common import ROOT, ASSETS, LESSONS, render_lesson, code, figure, ensure_dataset import numpy as np import pandas as pd import geopandas as gpd from shapely.geometry import Point, LineString import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.lines import Line2D from matplotlib.patches import Patch plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False t_all = time.time() print("=== L72 緊急輸送道路 単独 3 研究例分析 ===", flush=True) # ============================================================================= # 0. 定数・パス # ============================================================================= TARGET_CRS = "EPSG:6671" # JGD2011 平面直角第 III 系 DATA_DIR = ROOT / "data" / "extras" / "L72_emergency_road" DATA_DIR.mkdir(parents=True, exist_ok=True) DATASET_ID = 1247 RESOURCE_ID = 32491 ZIP_NAME = "emergency_road.zip" EXTRACT_DIR_NAME = "340006_emergency_transport_road_20220908T000000" # 緊急輸送道路 4 階層 (本記事独自命名) RANK_DEF = { "01": {"label": "第1次 (高速・幹線国道)", "color": "#cf222e", "weight": 4, "short": "1次"}, "02": {"label": "第2次 (主要地方道・国道)", "color": "#0969da", "weight": 3, "short": "2次"}, "03": {"label": "第3次 (一般県道・市町道)", "color": "#1a7f37", "weight": 2, "short": "3次"}, "04": {"label": "補完線 (補強路線)", "color": "#cf6f00", "weight": 2, "short": "補完"}, } RANK_ORDER = ["01", "02", "03", "04"] # L11 警戒区域 (既キャッシュ) SED_DIR = ROOT / "data" / "extras" / "sediment_shp" SED_FILES = { "急傾斜": (SED_DIR / "kyukeisha" / "340006_sediment_disaster_hazard_area_steep_slope_20260427", "ky_"), "土石流": (SED_DIR / "doseki" / "340006_sediment_disaster_hazard_area_debris_flow_20260427", "dy_"), "地すべり": (SED_DIR / "jisuberi" / "340006_sediment_disaster_hazard_area_landslide_20260427", "jy_"), } # 想定最大規模 浸水想定 (L4-L11 系で扱われている) FLOOD_SHP = (ROOT / "data" / "extras" / "flood_shp" / "shinsui_souteisaidai" / "shinsui_souteisaidai.shp") FLOOD_CACHE = DATA_DIR / "_cache_flood_max.gpkg" # 行政界キャッシュ (L44 で生成済) ADMIN_GPKG = ROOT / "data" / "extras" / "L44_storm_surge" / "_cache" / "admin_diss.gpkg" # L66 橋梁 + L67 トンネル (既キャッシュ CSV) L66_CSV = ASSETS / "L66_all_bridges.csv" L67_CSV = ASSETS / "L67_all_tunnels.csv" # 緊急輸送道路 バッファ幅 (橋梁・トンネル所属判定 用) ROAD_BUFFER_M = 30.0 # 30m バッファ = 道路敷地 + 構造物近接帯 # 50m 間隔点サンプリング SAMPLE_M = 50.0 # CITY_CD → 市町名 CITY_NAME = { 101: "広島市中区", 102: "広島市東区", 103: "広島市南区", 104: "広島市西区", 105: "広島市安佐南区", 106: "広島市安佐北区", 107: "広島市安芸区", 108: "広島市佐伯区", 202: "呉市", 203: "竹原市", 204: "三原市", 205: "尾道市", 207: "福山市", 208: "府中市", 209: "三次市", 210: "庄原市", 211: "大竹市", 212: "東広島市", 213: "廿日市市", 214: "安芸高田市", 215: "江田島市", 302: "府中町", 304: "海田町", 307: "熊野町", 309: "坂町", 369: "安芸太田町", 462: "世羅町", } # 地理クラス (RQ1 用、 L67 と同分類 + 北広島町・大崎上島町補完) COASTAL_ISLAND = {"江田島市", "大崎上島町", "宮島町", "倉橋町"} CHUSANKAN_CITIES = { "庄原市", "三次市", "北広島町", "安芸太田町", "神石高原町", "世羅町", "府中市", "安芸高田市", "大崎上島町", } def geo_class(name): if name in COASTAL_ISLAND or "江田島" in name: return "沿岸島嶼" if name in CHUSANKAN_CITIES: return "中山間山地" return "平野・沿岸都市" # ============================================================================= # 1. データ取得 + 展開 # ============================================================================= print("\n[1] データ取得 + 展開", flush=True) t1 = time.time() zip_local = DATA_DIR / ZIP_NAME ensure_dataset(zip_local, resource_id=RESOURCE_ID, min_bytes=10000, label="L72 emergency_road.zip") extract_dir = DATA_DIR / EXTRACT_DIR_NAME if not extract_dir.exists() or not any(extract_dir.glob("*.json")): with zipfile.ZipFile(zip_local) as zf: zf.extractall(DATA_DIR) print(f" 展開: {zip_local.name}", flush=True) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 2. JSON → LineString GeoDataFrame # ============================================================================= print("\n[2] JSON 読込 → LineString", flush=True) t2 = time.time() # (1) 階層メタ JSON meta_path = extract_dir / "05_kinkyu_route.json" with open(meta_path, "r", encoding="utf-8") as f: meta_text = f.read() # 形式: 単一 dict が「,」 区切りで並ぶ NDJSON 風 → 配列にラップ meta_list = json.loads("[" + meta_text + "]") # meta_list 例: route_01..route_06 + route_00 (色/太さ/線種) # (2) 各階層 JSON (_01〜_04) を LineString に変換 def load_route_lines(idx): """route_idx.json (NDJSON 風) を LineString 配列に変換。""" p = extract_dir / f"05_kinkyu_route_{idx}.json" with open(p, "r", encoding="utf-8") as f: text = f.read() arr = json.loads("[" + text + "]") lines = [] for seg in arr: if len(seg) >= 2: coords = [(pt["e"], pt["d"]) for pt in seg] lines.append(LineString(coords)) return lines records = [] for idx in RANK_ORDER: lines = load_route_lines(idx) for ln in lines: records.append({"rank": idx, "rank_label": RANK_DEF[idx]["label"], "geometry": ln}) gdf_road = gpd.GeoDataFrame(records, crs="EPSG:4326").to_crs(TARGET_CRS) gdf_road["length_m"] = gdf_road.geometry.length n_seg = len(gdf_road) total_km = gdf_road["length_m"].sum() / 1000 print(f" セグメント: {n_seg} / 総延長: {total_km:.1f} km", flush=True) print(f" ({time.time()-t2:.1f}s)", flush=True) # ============================================================================= # 3. 50m 間隔の点サンプリング (RQ2 で使う) # ============================================================================= print("\n[3] 50m 間隔 点サンプリング", flush=True) t3 = time.time() points = [] for line_id, (_, row) in enumerate(gdf_road.iterrows()): L = row.geometry.length n = max(int(L / SAMPLE_M), 1) for i in range(n + 1): d = min(i * SAMPLE_M, L) pt = row.geometry.interpolate(d) points.append({"line_id": line_id, "rank": row["rank"], "d_m": d, "geometry": pt}) gdf_pts = gpd.GeoDataFrame(points, crs=TARGET_CRS) n_pts = len(gdf_pts) print(f" サンプル点: {n_pts}", flush=True) print(f" ({time.time()-t3:.1f}s)", flush=True) # ============================================================================= # 4. RQ1: 階層・延長・地理 — 構造分析 # ============================================================================= print("\n[4] RQ1 集計: 階層・延長・地理", flush=True) t4 = time.time() # (1) 階層別 延長 + セグメント数 T_rank = (gdf_road.groupby(["rank", "rank_label"]) .agg(セグメント数=("geometry", "count"), 延長_km=("length_m", lambda s: s.sum() / 1000)) .reset_index() .rename(columns={"rank": "階層"})) T_rank["シェア_%"] = (T_rank["延長_km"] / total_km * 100).round(1) T_rank["延長_km"] = T_rank["延長_km"].round(1) T_rank = T_rank.sort_values("階層").reset_index(drop=True) # (2) 市町別 延長 (sjoin で各セグメントが入る市町を特定) admin = gpd.read_file(ADMIN_GPKG).to_crs(TARGET_CRS) admin["市町名"] = admin["CITY_CD"].map(CITY_NAME).fillna( admin["CITY_CD"].astype(str)) # 点サンプルを admin に sjoin → 各点の市町 pts_city = gpd.sjoin(gdf_pts, admin[["CITY_CD", "市町名", "geometry"]], how="left", predicate="within") # index_right が NaN の点 (海上・市町外) は除外 pts_city = pts_city.dropna(subset=["市町名"]).copy() # 各点 = 50m → 市町別延長 (km) T_city = (pts_city.groupby("市町名") .size() .to_frame("点数") .reset_index()) T_city["延長_km"] = (T_city["点数"] * SAMPLE_M / 1000).round(1) T_city = T_city[["市町名", "延長_km", "点数"]].sort_values( "延長_km", ascending=False).reset_index(drop=True) T_city["地理クラス"] = T_city["市町名"].apply(geo_class) T_city["シェア_%"] = (T_city["延長_km"] / T_city["延長_km"].sum() * 100).round(1) # (3) 階層 × 市町 クロス (top-10 市町のみ) top_cities = T_city.head(10)["市町名"].tolist() T_rank_city = (pts_city[pts_city["市町名"].isin(top_cities)] .groupby(["市町名", "rank"]) .size().unstack(fill_value=0) * SAMPLE_M / 1000).round(1) T_rank_city.columns = [f"{RANK_DEF[c]['short']}_km" for c in T_rank_city.columns] T_rank_city = T_rank_city.reset_index() T_rank_city = T_rank_city.set_index("市町名").loc[top_cities].reset_index() # (4) 地理クラス別 延長 pts_city["地理クラス"] = pts_city["市町名"].apply(geo_class) T_geo = (pts_city.groupby("地理クラス") .size().to_frame("点数").reset_index()) T_geo["延長_km"] = (T_geo["点数"] * SAMPLE_M / 1000).round(1) T_geo["シェア_%"] = (T_geo["延長_km"] / T_geo["延長_km"].sum() * 100).round(1) T_geo = T_geo[["地理クラス", "延長_km", "点数", "シェア_%"]] # H1: 第 2 次が ≥ 50% か share_2 = float(T_rank.loc[T_rank["階層"] == "02", "シェア_%"].iloc[0]) h1_ok = share_2 >= 50.0 # 階層別 km km_rank = {r: float(T_rank.loc[T_rank["階層"] == r, "延長_km"].iloc[0]) for r in RANK_ORDER} print(f" 階層別 km: {km_rank}", flush=True) print(f" H1 (第2次 ≥ 50%): {h1_ok} (実 {share_2:.1f}%)", flush=True) print(f" ({time.time()-t4:.1f}s)", flush=True) # ============================================================================= # 5. RQ2: 災害脆弱性 — 浸水 + 土砂警戒区域との交差 # ============================================================================= print("\n[5] RQ2 集計: 災害脆弱性", flush=True) t5 = time.time() # (1) 浸水想定 (想定最大規模) — キャッシュ GPKG なら 1 秒、 無ければ最初の 1 回 5 秒 if FLOOD_CACHE.exists() and FLOOD_CACHE.stat().st_size > 1000: fl = gpd.read_file(FLOOD_CACHE) print(f" 浸水想定 (cache): {len(fl)} polys", flush=True) else: fl_raw = gpd.read_file(FLOOD_SHP, encoding="utf-8") fl = fl_raw[["geometry"]].to_crs(TARGET_CRS).copy() fl.to_file(FLOOD_CACHE, driver="GPKG") print(f" 浸水想定 (init cache): {len(fl)} polys", flush=True) # (2) sjoin: 各点が 浸水想定 / 急傾斜 / 土石流 / 地すべり に入るか def overlay_pts(pts, polys, label): """点が polys 内にあるかを sjoin で判定し、 重なり長さを集計。""" hits = gpd.sjoin(pts, polys[["geometry"]], how="left", predicate="within") flag = hits["index_right"].notna().values # 重複 (1 点が複数 poly に同時に入る) を line_id+d_m で uniq hit_df = hits[hits["index_right"].notna()].drop_duplicates( ["line_id", "d_m"]) # 階層別 by_rank = hit_df.groupby("rank").size().reindex(RANK_ORDER, fill_value=0) return flag[:len(pts)], hit_df, by_rank flood_flag, flood_hits, flood_by_rank = overlay_pts(gdf_pts, fl, "浸水") print(f" 浸水想定: {flood_flag.sum()} pts hit / {n_pts}", flush=True) sed_flags = {} sed_hits = {} sed_by_rank = {} sed_polys = {} # 図 4 で再利用するためキャッシュ for label, (d, key) in SED_FILES.items(): files = [p for p in d.rglob("*.shp") if key in p.stem] parts = [] for p in files: try: parts.append(gpd.read_file(p, encoding="utf-8")[["geometry"]]) except Exception as e: print(f" WARN {p.name}: {e}", flush=True) if not parts: sed_flags[label] = np.zeros(n_pts, dtype=bool) sed_hits[label] = pd.DataFrame() sed_by_rank[label] = pd.Series(0, index=RANK_ORDER) sed_polys[label] = None continue sed = gpd.GeoDataFrame(pd.concat(parts, ignore_index=True), crs=parts[0].crs).to_crs(TARGET_CRS) f, h, br = overlay_pts(gdf_pts, sed, label) sed_flags[label] = f sed_hits[label] = h sed_by_rank[label] = br sed_polys[label] = sed # キャッシュ print(f" {label}: {f.sum()} pts hit", flush=True) # (3) 土砂 3 種 「いずれか」 重なり sed_any = (sed_flags["急傾斜"] | sed_flags["土石流"] | sed_flags["地すべり"]) gdf_pts["浸水"] = flood_flag gdf_pts["急傾斜"] = sed_flags["急傾斜"] gdf_pts["土石流"] = sed_flags["土石流"] gdf_pts["地すべり"] = sed_flags["地すべり"] gdf_pts["土砂_いずれか"] = sed_any gdf_pts["災害_いずれか"] = flood_flag | sed_any # (4) 階層別 重なり率 def by_rank_summary(flag_col): res = (gdf_pts.groupby("rank")[flag_col].sum() .reindex(RANK_ORDER, fill_value=0)) return res rank_total = (gdf_pts.groupby("rank").size() .reindex(RANK_ORDER, fill_value=0)) flood_rank = by_rank_summary("浸水") sed_any_rank = by_rank_summary("土砂_いずれか") hazard_rank = by_rank_summary("災害_いずれか") ky_rank = by_rank_summary("急傾斜") ds_rank = by_rank_summary("土石流") js_rank = by_rank_summary("地すべり") T_rq2 = pd.DataFrame({ "階層": [RANK_DEF[r]["label"] for r in RANK_ORDER], "総延長_km": [round(rank_total[r] * SAMPLE_M / 1000, 1) for r in RANK_ORDER], "浸水_km": [round(flood_rank[r] * SAMPLE_M / 1000, 1) for r in RANK_ORDER], "浸水_%": [round(100 * flood_rank[r] / max(rank_total[r], 1), 1) for r in RANK_ORDER], "急傾斜_km": [round(ky_rank[r] * SAMPLE_M / 1000, 1) for r in RANK_ORDER], "土石流_km": [round(ds_rank[r] * SAMPLE_M / 1000, 1) for r in RANK_ORDER], "地すべり_km": [round(js_rank[r] * SAMPLE_M / 1000, 1) for r in RANK_ORDER], "土砂いずれか_km": [round(sed_any_rank[r] * SAMPLE_M / 1000, 1) for r in RANK_ORDER], "土砂いずれか_%": [round(100 * sed_any_rank[r] / max(rank_total[r], 1), 1) for r in RANK_ORDER], "災害いずれか_km": [round(hazard_rank[r] * SAMPLE_M / 1000, 1) for r in RANK_ORDER], "災害いずれか_%": [round(100 * hazard_rank[r] / max(rank_total[r], 1), 1) for r in RANK_ORDER], }) # 全体 (階層合算) all_total = T_rq2["総延長_km"].sum() all_flood_km = T_rq2["浸水_km"].sum() all_sed_km = T_rq2["土砂いずれか_km"].sum() all_hazard_km = T_rq2["災害いずれか_km"].sum() all_flood_pct = round(100 * all_flood_km / max(all_total, 1), 1) all_sed_pct = round(100 * all_sed_km / max(all_total, 1), 1) all_hazard_pct = round(100 * all_hazard_km / max(all_total, 1), 1) # H2: 第1次 浸水 < 10% AND 土砂 < 10% flood_1 = float(T_rq2.iloc[0]["浸水_%"]) sed_1 = float(T_rq2.iloc[0]["土砂いずれか_%"]) h2_ok = (flood_1 < 10.0) and (sed_1 < 10.0) # H3: 全体 浸水 % > 土砂 % h3_ok = all_flood_pct > all_sed_pct print(f" 全体 浸水 {all_flood_km:.1f} km ({all_flood_pct}%) / " f"土砂 {all_sed_km:.1f} km ({all_sed_pct}%) / " f"いずれか {all_hazard_km:.1f} km ({all_hazard_pct}%)", flush=True) print(f" H2 (第1次 浸水<10% AND 土砂<10%): {h2_ok} " f"(浸水 {flood_1}%, 土砂 {sed_1}%)", flush=True) print(f" H3 (全体 浸水 > 土砂): {h3_ok} " f"({all_flood_pct} vs {all_sed_pct})", flush=True) print(f" ({time.time()-t5:.1f}s)", flush=True) # ============================================================================= # 6. RQ3: 構造物リスク — 橋梁・トンネル・老朽橋 # ============================================================================= print("\n[6] RQ3 集計: 構造物リスク", flush=True) t6 = time.time() # (1) 橋梁・トンネル CSV → POINT df_b = pd.read_csv(L66_CSV, encoding="utf-8-sig") df_b = df_b.dropna(subset=["緯度（10進数）", "経度（10進数）"]) df_b = df_b[df_b["緯度（10進数）"] < 50].copy() gdf_b = gpd.GeoDataFrame( df_b, geometry=[Point(x, y) for x, y in zip(df_b["経度（10進数）"], df_b["緯度（10進数）"])], crs="EPSG:4326").to_crs(TARGET_CRS) df_t = pd.read_csv(L67_CSV, encoding="utf-8-sig") df_t = df_t.dropna(subset=["緯度（10進数）", "経度（10進数）"]) df_t = df_t[df_t["緯度（10進数）"] < 50].copy() gdf_t = gpd.GeoDataFrame( df_t, geometry=[Point(x, y) for x, y in zip(df_t["経度（10進数）"], df_t["緯度（10進数）"])], crs="EPSG:4326").to_crs(TARGET_CRS) n_bridge = len(gdf_b) n_tunnel = len(gdf_t) print(f" 橋梁: {n_bridge}, トンネル: {n_tunnel}", flush=True) # (2) 緊急輸送道路 30m バッファ (階層別 attr 保持) gdf_road["geometry_buf"] = gdf_road.geometry.buffer(ROAD_BUFFER_M) buf = gpd.GeoDataFrame( gdf_road[["rank", "rank_label", "geometry_buf"]].rename( columns={"geometry_buf": "geometry"}), geometry="geometry", crs=TARGET_CRS) # (3) 橋梁 within buffer (sjoin) — 階層情報も付与 b_on = gpd.sjoin(gdf_b, buf[["rank", "geometry"]], how="left", predicate="within") # 1 橋梁が複数階層 buffer に入る場合 → 最高位 (= 数字小さい) を採用 b_on_grp = b_on.groupby(b_on.index)["rank"].apply( lambda s: s.dropna().min() if s.notna().any() else np.nan) gdf_b["road_rank"] = b_on_grp.reindex(gdf_b.index) gdf_b["on_road"] = gdf_b["road_rank"].notna() t_on = gpd.sjoin(gdf_t, buf[["rank", "geometry"]], how="left", predicate="within") t_on_grp = t_on.groupby(t_on.index)["rank"].apply( lambda s: s.dropna().min() if s.notna().any() else np.nan) gdf_t["road_rank"] = t_on_grp.reindex(gdf_t.index) gdf_t["on_road"] = gdf_t["road_rank"].notna() n_b_on = int(gdf_b["on_road"].sum()) n_t_on = int(gdf_t["on_road"].sum()) b_on_pct = round(100 * n_b_on / n_bridge, 1) t_on_pct = round(100 * n_t_on / n_tunnel, 1) print(f" 橋梁 on 緊急輸送道路: {n_b_on}/{n_bridge} ({b_on_pct}%)", flush=True) print(f" トンネル on 緊急輸送道路: {n_t_on}/{n_tunnel} ({t_on_pct}%)", flush=True) # (4) 階層別 橋梁・トンネル 件数 def count_by_rank(gdf_x): s = gdf_x[gdf_x["on_road"]]["road_rank"].value_counts() return s.reindex(RANK_ORDER, fill_value=0).astype(int) b_by_rank = count_by_rank(gdf_b) t_by_rank = count_by_rank(gdf_t) # (5) 老朽橋 (架設 < 1970) on 緊急輸送道路 gdf_b["年代"] = gdf_b["架設年度"].astype(float, errors="ignore") gdf_b["is_old_pre1970"] = gdf_b["年代"] < 1970 gdf_b["is_old_pre1980"] = gdf_b["年代"] < 1980 # 耐震基準改定 1980 n_old70_total = int(gdf_b["is_old_pre1970"].sum()) n_old70_on = int((gdf_b["is_old_pre1970"] & gdf_b["on_road"]).sum()) n_old80_total = int(gdf_b["is_old_pre1980"].sum()) n_old80_on = int((gdf_b["is_old_pre1980"] & gdf_b["on_road"]).sum()) # (6) 階層別 老朽橋 (1970 以前) old_by_rank = (gdf_b[gdf_b["on_road"] & gdf_b["is_old_pre1970"]] ["road_rank"].value_counts() .reindex(RANK_ORDER, fill_value=0).astype(int)) # 階層別 BCP 表 T_rq3 = pd.DataFrame({ "階層": [RANK_DEF[r]["label"] for r in RANK_ORDER], "延長_km": [km_rank[r] for r in RANK_ORDER], "橋梁数": [int(b_by_rank[r]) for r in RANK_ORDER], "橋梁密度_件/km": [round(b_by_rank[r] / max(km_rank[r], 1), 2) for r in RANK_ORDER], "トンネル数": [int(t_by_rank[r]) for r in RANK_ORDER], "老朽橋_pre1970": [int(old_by_rank[r]) for r in RANK_ORDER], "老朽橋_%": [round(100 * old_by_rank[r] / max(b_by_rank[r], 1), 1) for r in RANK_ORDER], }) # H4: 橋梁 ≥ 30% AND トンネル ≥ 50% h4_ok = (b_on_pct >= 30.0) and (t_on_pct >= 50.0) # H5: 老朽橋 (pre1970) ≥ 200 件 h5_ok = n_old70_on >= 200 print(f" 老朽橋 (pre1970): on={n_old70_on}/{n_old70_total} (緊急輸送道路上)", flush=True) print(f" H4 (橋梁≥30% AND トンネル≥50%): {h4_ok}", flush=True) print(f" H5 (老朽橋 pre1970 ≥ 200): {h5_ok} ({n_old70_on})", flush=True) print(f" ({time.time()-t6:.1f}s)", flush=True) # ============================================================================= # 7. CSV 出力 # ============================================================================= print("\n[7] CSV 出力", flush=True) t7 = time.time() # (1) 緊急輸送道路 全セグメント (geometry除く) df_road_out = gdf_road.drop(columns=["geometry", "geometry_buf"]).copy() df_road_out["延長_km"] = (df_road_out["length_m"] / 1000).round(3) df_road_out.drop(columns=["length_m"]).to_csv( ASSETS / "L72_all_segments.csv", index=False, encoding="utf-8-sig") # (2) 階層別 / 市町別 / 地理クラス別 T_rank.to_csv(ASSETS / "L72_rank_summary.csv", index=False, encoding="utf-8-sig") T_city.to_csv(ASSETS / "L72_city_summary.csv", index=False, encoding="utf-8-sig") T_rank_city.to_csv(ASSETS / "L72_rank_x_city_top10.csv", index=False, encoding="utf-8-sig") T_geo.to_csv(ASSETS / "L72_geo_class.csv", index=False, encoding="utf-8-sig") # (3) RQ2 階層別災害重なり T_rq2.to_csv(ASSETS / "L72_hazard_overlap.csv", index=False, encoding="utf-8-sig") # (4) RQ3 階層別 BCP T_rq3.to_csv(ASSETS / "L72_bcp_summary.csv", index=False, encoding="utf-8-sig") # (5) 橋梁 + 老朽橋 on 緊急輸送道路 (個別) gdf_b_out = gdf_b[gdf_b["on_road"]][[ "橋梁番号", "施設名", "路線名", "道路種別", "架設年度", "延長(m)", "幅員(m)", "管理事務所名", "住所(市町)", "緯度（10進数）", "経度（10進数）", "判定区分", "is_old_pre1970", "is_old_pre1980", "road_rank"]].copy() gdf_b_out.rename(columns={"road_rank": "緊急輸送道路階層"}, inplace=True) gdf_b_out["緊急輸送道路階層"] = gdf_b_out["緊急輸送道路階層"].map( {r: RANK_DEF[r]["short"] for r in RANK_ORDER}) gdf_b_out.to_csv(ASSETS / "L72_bridges_on_emerg_road.csv", index=False, encoding="utf-8-sig") # (6) トンネル on 緊急輸送道路 gdf_t_out = gdf_t[gdf_t["on_road"]][[ "トンネル番号", "施設名", "路線名", "道路種別", "建設年度", "延長(m)", "幅員(m)", "管理事務所名", "住所(市町)", "緯度（10進数）", "経度（10進数）", "road_rank"]].copy() gdf_t_out.rename(columns={"road_rank": "緊急輸送道路階層"}, inplace=True) gdf_t_out["緊急輸送道路階層"] = gdf_t_out["緊急輸送道路階層"].map( {r: RANK_DEF[r]["short"] for r in RANK_ORDER}) gdf_t_out.to_csv(ASSETS / "L72_tunnels_on_emerg_road.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t7:.1f}s)", flush=True) # ============================================================================= # 8. 図の生成 (8 図) # ============================================================================= print("\n[8] 図の生成", flush=True) t8 = time.time() ASSETS.mkdir(parents=True, exist_ok=True) def save_fig(name, dpi=120): p = ASSETS / name plt.savefig(p, dpi=dpi, bbox_inches="tight", facecolor="white") plt.close('all') return p # ---- 図 1 (RQ1): 県全域 階層別 緊急輸送道路マップ ---- print(" fig1: 階層別 マップ", flush=True) fig, ax = plt.subplots(figsize=(13, 8)) admin.plot(ax=ax, color="#fff4e0", edgecolor="#888", linewidth=0.4, alpha=0.55) # 第3次/補完を先に薄く描画 → 第2次 → 第1次 が最前面 draw_order = ["04", "03", "02", "01"] for r in draw_order: sub = gdf_road[gdf_road["rank"] == r] sub.plot(ax=ax, color=RANK_DEF[r]["color"], linewidth=RANK_DEF[r]["weight"] * 0.7, alpha=0.95, zorder=2 + draw_order.index(r)) ax.set_xlim(-15000, 125000) ax.set_ylim(-220000, -90000) ax.set_aspect("equal") ax.set_title(f"図 1 (RQ1): 広島県 緊急輸送道路 4 階層マップ — " f"総 {total_km:.0f} km / " f"1次 {km_rank['01']:.0f} km / " f"2次 {km_rank['02']:.0f} km / " f"3次 {km_rank['03']:.0f} km / " f"補完 {km_rank['04']:.0f} km", fontsize=10.5) ax.set_xlabel("X (m, EPSG:6671)") ax.set_ylabel("Y (m, EPSG:6671)") patches = [ Line2D([0], [0], color=RANK_DEF[r]["color"], linewidth=RANK_DEF[r]["weight"] * 1.0, label=f"{RANK_DEF[r]['label']} ({km_rank[r]:.0f} km)") for r in RANK_ORDER ] ax.legend(handles=patches, loc="lower left", fontsize=10, title="階層 (色=識別)") plt.tight_layout() save_fig("L72_fig1_rank_map.png") # ---- 図 2 (RQ1): 階層別 セグメント数 + 延長 棒 + 地理クラス円 ---- print(" fig2: 階層・地理 構造", flush=True) fig, axes = plt.subplots(1, 3, figsize=(16, 5)) # 左: 階層別 延長 ax = axes[0] xs = np.arange(len(RANK_ORDER)) vals = [km_rank[r] for r in RANK_ORDER] cols = [RANK_DEF[r]["color"] for r in RANK_ORDER] ax.bar(xs, vals, color=cols, edgecolor="#333", linewidth=0.6, width=0.65) for x, v in zip(xs, vals): ax.text(x, v + total_km * 0.01, f"{v:.0f} km\n({100*v/total_km:.1f}%)", ha="center", fontsize=10, fontweight="bold") ax.set_xticks(xs) ax.set_xticklabels([RANK_DEF[r]["short"] for r in RANK_ORDER], fontsize=11) ax.set_ylabel("延長 (km)") ax.set_title(f"階層別 延長 — 総 {total_km:.0f} km\n" f"H1 第2次 {share_2:.1f}% (≥50%: {h1_ok})", fontsize=10.5) ax.grid(True, axis="y", alpha=0.3) # 中: 階層別 セグメント数 ax = axes[1] seg_vals = [int(T_rank.loc[T_rank['階層'] == r, 'セグメント数'].iloc[0]) for r in RANK_ORDER] ax.bar(xs, seg_vals, color=cols, edgecolor="#333", linewidth=0.6, width=0.65) for x, v in zip(xs, seg_vals): ax.text(x, v + max(seg_vals) * 0.02, f"{v}", ha="center", fontsize=10, fontweight="bold") ax.set_xticks(xs) ax.set_xticklabels([RANK_DEF[r]["short"] for r in RANK_ORDER], fontsize=11) ax.set_ylabel("セグメント数") ax.set_title(f"階層別 セグメント数 (計 {n_seg})", fontsize=10.5) ax.grid(True, axis="y", alpha=0.3) # 右: 地理クラス別 延長 (パイ) ax = axes[2] geo_labels = T_geo["地理クラス"].tolist() geo_vals = T_geo["延長_km"].tolist() geo_cols = ["#0969da", "#1a7f37", "#cf6f00"] wedges, texts, autotexts = ax.pie( geo_vals, labels=geo_labels, colors=geo_cols, autopct=lambda v: f"{v:.1f}%\n({v*sum(geo_vals)/100:.0f} km)", startangle=90, textprops={"fontsize": 10}) for t in autotexts: t.set_color("white") t.set_fontweight("bold") ax.set_title(f"地理クラス別 延長", fontsize=10.5) fig.suptitle("図 2 (RQ1): 階層別 延長 + セグメント数 + 地理クラス分布", fontsize=12.5, y=1.02) plt.tight_layout() save_fig("L72_fig2_rank_seg_geo.png") # ---- 図 3 (RQ1): 市町別 緊急輸送道路 延長 choropleth ---- print(" fig3: 市町別 choropleth", flush=True) fig, axes = plt.subplots(1, 2, figsize=(15, 7)) # 左: choropleth (km) ax = axes[0] admin_with = admin.merge(T_city[["市町名", "延長_km"]], on="市町名", how="left") admin_with["延長_km"] = admin_with["延長_km"].fillna(0) admin_with.plot(ax=ax, column="延長_km", cmap="YlOrRd", edgecolor="#444", linewidth=0.4, legend=True, legend_kwds={"label": "緊急輸送道路 延長 (km)", "orientation": "vertical", "shrink": 0.7}) # 第1次のみ重ねる (ハイライト) gdf_road[gdf_road["rank"] == "01"].plot( ax=ax, color=RANK_DEF["01"]["color"], linewidth=1.5, alpha=0.9, zorder=3) ax.set_xlim(-15000, 125000) ax.set_ylim(-220000, -90000) ax.set_aspect("equal") ax.set_title(f"市町別 緊急輸送道路 延長 (km, 全階層)\n" f"赤線 = 第1次 (高速・幹線国道) のみ重畳", fontsize=10.5) ax.set_xlabel("X (m)") ax.set_ylabel("Y (m)") # 右: 市町ランキング (top 15) ax = axes[1] top15 = T_city.head(15) ys = np.arange(len(top15)) ax.barh(ys, top15["延長_km"].values, color="#0969da", edgecolor="#333", linewidth=0.4) for y, v in zip(ys, top15["延長_km"].values): ax.text(v + 5, y, f"{v:.0f} km", va="center", fontsize=9.5) ax.set_yticks(ys) ax.set_yticklabels(top15["市町名"].values, fontsize=10) ax.invert_yaxis() ax.set_xlabel("延長 (km)") ax.set_title(f"市町別 ランキング (Top 15)", fontsize=10.5) ax.grid(True, axis="x", alpha=0.3) fig.suptitle("図 3 (RQ1): 市町別 緊急輸送道路カバー — 中山間 vs 沿岸都市", fontsize=12.5, y=1.02) plt.tight_layout() save_fig("L72_fig3_city_choropleth.png") # ---- 図 4 (RQ2): 緊急輸送道路 + 浸水想定 + 土砂警戒 重ね合わせマップ ---- # 大量ポリゴン (浸水 613 + 急傾斜 30k + 土石流 18k) は matplotlib が # 1 ポリゴン = 1 Patch を add_patch するため数十秒〜数分の描画時間がかかる。 # 解決策: 各ポリゴンの bbox を 1 km グリッドに丸めて 「グリッドセルの存否マスク」 # にした上で imshow で描画 (= 個別ポリゴンの輪郭は無視、 概観だけ示す)。 print(" fig4: 災害脆弱性 重ね合わせマップ", flush=True) GRID_M = 1000.0 xmin, ymin = -15000, -220000 xmax, ymax = 125000, -90000 nx = int((xmax - xmin) / GRID_M) + 1 ny = int((ymax - ymin) / GRID_M) + 1 def bbox_grid_mask(gdf_polys): """ポリゴン群の bbox を 1 km グリッドのブール mask に変換 (ベクトル化)。""" if gdf_polys is None or len(gdf_polys) == 0: return np.zeros((ny, nx), dtype=bool) b = gdf_polys.bounds i0 = np.maximum(((b["minx"].values - xmin) / GRID_M).astype(int), 0) i1 = np.minimum(((b["maxx"].values - xmin) / GRID_M).astype(int) + 1, nx) j0 = np.maximum(((b["miny"].values - ymin) / GRID_M).astype(int), 0) j1 = np.minimum(((b["maxy"].values - ymin) / GRID_M).astype(int) + 1, ny) m = np.zeros((ny, nx), dtype=bool) for k in range(len(b)): m[j0[k]:j1[k], i0[k]:i1[k]] = True return m fig, ax = plt.subplots(figsize=(13, 8)) admin.plot(ax=ax, color="#fff4e0", edgecolor="#888", linewidth=0.4, alpha=0.55) # 浸水 (青) mask_fl = bbox_grid_mask(fl) ax.imshow(mask_fl, extent=(xmin, xmax, ymin, ymax), origin="lower", aspect="equal", cmap=matplotlib.colors.ListedColormap(["#00000000", "#0969da"]), alpha=0.30, zorder=2, interpolation="nearest") # 土砂 (急傾斜 + 土石流) for label, c in [("急傾斜", "#cf222e"), ("土石流", "#7c3aed")]: mask = bbox_grid_mask(sed_polys.get(label)) if mask.any(): ax.imshow(mask, extent=(xmin, xmax, ymin, ymax), origin="lower", aspect="equal", cmap=matplotlib.colors.ListedColormap(["#00000000", c]), alpha=0.30, zorder=3, interpolation="nearest") # 緊急輸送道路 (全階層) — 630 セグだが線描画は速い for r in RANK_ORDER[::-1]: sub = gdf_road[gdf_road["rank"] == r] sub.plot(ax=ax, color=RANK_DEF[r]["color"], linewidth=RANK_DEF[r]["weight"] * 0.7, alpha=0.95, zorder=5 + RANK_ORDER.index(r)) ax.set_xlim(-15000, 125000) ax.set_ylim(-220000, -90000) ax.set_aspect("equal") ax.set_title(f"図 4 (RQ2): 緊急輸送道路 + 浸水想定 + 土砂警戒区域 重ね合わせ\n" f"(災害区域は 1 km グリッド近似表示; 重なり率は実形状で計算) — " f"浸水 {all_flood_pct}% / 土砂 {all_sed_pct}% / " f"いずれか {all_hazard_pct}%", fontsize=9.5) ax.set_xlabel("X (m, EPSG:6671)") ax.set_ylabel("Y (m, EPSG:6671)") patches = [ Patch(facecolor="#0969da", alpha=0.45, label="浸水想定 (想定最大規模)"), Patch(facecolor="#cf222e", alpha=0.45, label="急傾斜地崩壊警戒区域"), Patch(facecolor="#7c3aed", alpha=0.45, label="土石流警戒区域"), ] patches += [ Line2D([0], [0], color=RANK_DEF[r]["color"], linewidth=RANK_DEF[r]["weight"], label=RANK_DEF[r]["short"]) for r in RANK_ORDER ] ax.legend(handles=patches, loc="lower left", fontsize=9, title="災害区域 + 道路階層", ncol=2) plt.tight_layout() save_fig("L72_fig4_hazard_overlay.png") # ---- 図 5 (RQ2): 階層別 重なり率 グループ棒 ---- print(" fig5: 階層別 重なり率", flush=True) fig, axes = plt.subplots(1, 2, figsize=(15, 5.5)) # 左: 階層別 浸水 / 土砂 / いずれか % ax = axes[0] xs = np.arange(len(RANK_ORDER)) width = 0.27 vals_flood = [float(T_rq2.iloc[i]["浸水_%"]) for i in range(len(RANK_ORDER))] vals_sed = [float(T_rq2.iloc[i]["土砂いずれか_%"]) for i in range(len(RANK_ORDER))] vals_any = [float(T_rq2.iloc[i]["災害いずれか_%"]) for i in range(len(RANK_ORDER))] ax.bar(xs - width, vals_flood, width, color="#0969da", edgecolor="#333", linewidth=0.4, label="浸水想定") ax.bar(xs, vals_sed, width, color="#cf222e", edgecolor="#333", linewidth=0.4, label="土砂警戒 (3種)") ax.bar(xs + width, vals_any, width, color="#7c3aed", edgecolor="#333", linewidth=0.4, label="いずれか") for x, v in zip(xs - width, vals_flood): ax.text(x, v + 0.5, f"{v:.1f}", ha="center", fontsize=9) for x, v in zip(xs, vals_sed): ax.text(x, v + 0.5, f"{v:.1f}", ha="center", fontsize=9) for x, v in zip(xs + width, vals_any): ax.text(x, v + 0.5, f"{v:.1f}", ha="center", fontsize=9) ax.axhline(10, color="#666", linestyle="--", linewidth=1, label="10% 閾値 (H2)") ax.set_xticks(xs) ax.set_xticklabels([RANK_DEF[r]["short"] for r in RANK_ORDER], fontsize=11) ax.set_ylabel("延長重なり率 (%)") ax.set_title(f"階層別 災害区域 重なり率\n" f"H2 (1次 浸水<10% AND 土砂<10%): {h2_ok}", fontsize=10.5) ax.legend(fontsize=9, loc="upper left") ax.grid(True, axis="y", alpha=0.3) # 右: 全体 (合算) 浸水 vs 土砂 — H3 検証 ax = axes[1] labels = ["浸水想定\n(想定最大規模)", "土砂警戒\n(3種いずれか)", "災害いずれか\n(浸水 ∪ 土砂)"] vals3 = [all_flood_pct, all_sed_pct, all_hazard_pct] cols3 = ["#0969da", "#cf222e", "#7c3aed"] xs3 = np.arange(len(labels)) ax.bar(xs3, vals3, color=cols3, edgecolor="#333", linewidth=0.5, width=0.6) for x, v in zip(xs3, vals3): ax.text(x, v + 0.5, f"{v:.1f}%", ha="center", fontsize=11, fontweight="bold") ax.text(x, v / 2, f"{[all_flood_km, all_sed_km, all_hazard_km][int(x)]:.0f} km", ha="center", fontsize=10, color="white", fontweight="bold") ax.set_xticks(xs3) ax.set_xticklabels(labels, fontsize=10) ax.set_ylabel("緊急輸送道路 全延長に対する % ") ax.set_title(f"H3 検証: 全体 浸水 {all_flood_pct}% vs 土砂 {all_sed_pct}%\n" f"(H3 浸水>土砂: {h3_ok})", fontsize=10.5) ax.grid(True, axis="y", alpha=0.3) fig.suptitle("図 5 (RQ2): 階層別 + 全体 災害区域 重なり率", fontsize=12.5, y=1.02) plt.tight_layout() save_fig("L72_fig5_hazard_pct.png") # ---- 図 6 (RQ3): 緊急輸送道路 + 橋梁 + トンネル マップ ---- print(" fig6: 構造物 重ね合わせマップ", flush=True) fig, ax = plt.subplots(figsize=(13, 8)) admin.plot(ax=ax, color="#fff4e0", edgecolor="#888", linewidth=0.4, alpha=0.55) # 緊急輸送道路 (薄目に背景表示) for r in RANK_ORDER[::-1]: sub = gdf_road[gdf_road["rank"] == r] sub.plot(ax=ax, color=RANK_DEF[r]["color"], linewidth=RANK_DEF[r]["weight"] * 0.6, alpha=0.6, zorder=2 + RANK_ORDER.index(r)) # 橋梁 全件 (灰、 背景) gdf_b[~gdf_b["on_road"]].plot(ax=ax, color="#bbb", markersize=3, alpha=0.45, edgecolor="none", zorder=4) # 橋梁 on 緊急輸送 (赤、 前面) gdf_b[gdf_b["on_road"]].plot(ax=ax, color="#cf222e", markersize=6, alpha=0.85, edgecolor="#222", linewidth=0.2, zorder=5) # トンネル on 緊急輸送 (紫菱形、 最前面) gdf_t[gdf_t["on_road"]].plot(ax=ax, color="#7c3aed", markersize=70, marker="D", alpha=0.95, edgecolor="#000", linewidth=0.6, zorder=7) # トンネル off (薄) gdf_t[~gdf_t["on_road"]].plot(ax=ax, color="#aaa", markersize=22, marker="D", alpha=0.55, edgecolor="#444", linewidth=0.3, zorder=6) ax.set_xlim(-15000, 125000) ax.set_ylim(-220000, -90000) ax.set_aspect("equal") ax.set_title(f"図 6 (RQ3): 緊急輸送道路 + 橋梁・トンネル — " f"橋梁 on {n_b_on}/{n_bridge} ({b_on_pct}%) / " f"トンネル on {n_t_on}/{n_tunnel} ({t_on_pct}%)", fontsize=10) ax.set_xlabel("X (m, EPSG:6671)") ax.set_ylabel("Y (m, EPSG:6671)") patches = [ Line2D([0], [0], color=RANK_DEF["01"]["color"], linewidth=3, label="第1次 緊急輸送道路"), Line2D([0], [0], color=RANK_DEF["02"]["color"], linewidth=2, label="第2次"), Line2D([0], [0], color=RANK_DEF["03"]["color"], linewidth=1.5, label="第3次"), Line2D([0], [0], marker='o', color='w', markerfacecolor="#cf222e", markeredgecolor="#222", markersize=8, label=f"橋梁 on 緊急 ({n_b_on})"), Line2D([0], [0], marker='o', color='w', markerfacecolor="#bbb", markeredgecolor="none", markersize=6, label=f"橋梁 off 緊急 ({n_bridge - n_b_on})"), Line2D([0], [0], marker='D', color='w', markerfacecolor="#7c3aed", markeredgecolor="#000", markersize=10, label=f"トンネル on 緊急 ({n_t_on})"), Line2D([0], [0], marker='D', color='w', markerfacecolor="#aaa", markeredgecolor="#444", markersize=7, label=f"トンネル off 緊急 ({n_tunnel - n_t_on})"), ] ax.legend(handles=patches, loc="lower left", fontsize=8.5, ncol=2) plt.tight_layout() save_fig("L72_fig6_struct_overlay.png") # ---- 図 7 (RQ3): 階層別 橋梁数 + 老朽橋 (pre1970) 件数 棒 + 密度 ---- print(" fig7: 階層別 BCP 構造物", flush=True) fig, axes = plt.subplots(1, 3, figsize=(16, 5.5)) # 左: 階層別 橋梁・トンネル 件数 (積み) ax = axes[0] xs = np.arange(len(RANK_ORDER)) b_vals = [int(b_by_rank[r]) for r in RANK_ORDER] t_vals = [int(t_by_rank[r]) for r in RANK_ORDER] ax.bar(xs - 0.2, b_vals, 0.4, color="#cf222e", edgecolor="#333", linewidth=0.4, label="橋梁") ax.bar(xs + 0.2, t_vals, 0.4, color="#7c3aed", edgecolor="#333", linewidth=0.4, label="トンネル") for x, v in zip(xs - 0.2, b_vals): ax.text(x, v + max(b_vals) * 0.01, f"{v}", ha="center", fontsize=9.5, fontweight="bold") for x, v in zip(xs + 0.2, t_vals): ax.text(x, v + max(b_vals) * 0.01, f"{v}", ha="center", fontsize=9.5, fontweight="bold") ax.set_xticks(xs) ax.set_xticklabels([RANK_DEF[r]["short"] for r in RANK_ORDER], fontsize=11) ax.set_ylabel("件数") ax.set_title(f"階層別 橋梁 + トンネル (緊急輸送道路上)\n" f"橋梁 {n_b_on}/{n_bridge} / トンネル {n_t_on}/{n_tunnel}", fontsize=10) ax.legend(fontsize=10) ax.grid(True, axis="y", alpha=0.3) # 中: 階層別 橋梁密度 (件/km) ax = axes[1] b_dens = [round(b_by_rank[r] / max(km_rank[r], 1), 2) for r in RANK_ORDER] ax.bar(xs, b_dens, color="#0969da", edgecolor="#333", linewidth=0.4, width=0.65) for x, v in zip(xs, b_dens): ax.text(x, v + max(b_dens) * 0.01, f"{v}", ha="center", fontsize=10, fontweight="bold") ax.set_xticks(xs) ax.set_xticklabels([RANK_DEF[r]["short"] for r in RANK_ORDER], fontsize=11) ax.set_ylabel("橋梁密度 (件/km)") ax.set_title(f"階層別 橋梁密度 — 高密度 = 中小河川クロス多", fontsize=10) ax.grid(True, axis="y", alpha=0.3) # 右: 階層別 老朽橋 (pre1970) ax = axes[2] old_vals = [int(old_by_rank[r]) for r in RANK_ORDER] ax.bar(xs, old_vals, color="#cf6f00", edgecolor="#333", linewidth=0.4, width=0.65) for x, v in zip(xs, old_vals): ax.text(x, v + max(old_vals) * 0.02, f"{v}", ha="center", fontsize=10, fontweight="bold") ax.set_xticks(xs) ax.set_xticklabels([RANK_DEF[r]["short"] for r in RANK_ORDER], fontsize=11) ax.set_ylabel("老朽橋件数 (架設 < 1970)") ax.set_title(f"階層別 老朽橋 — BCP 脆弱箇所\n" f"H5: 緊急輸送道路上 老朽橋 {n_old70_on} 件 (≥200: {h5_ok})", fontsize=10) ax.grid(True, axis="y", alpha=0.3) fig.suptitle("図 7 (RQ3): 階層別 橋梁・トンネル + 老朽橋 BCP", fontsize=12.5, y=1.02) plt.tight_layout() save_fig("L72_fig7_bcp_struct.png") # ---- 図 8 (RQ3): 老朽橋 (pre1970) on 緊急輸送道路 マップ ---- print(" fig8: 老朽橋 BCP 脆弱箇所マップ", flush=True) fig, ax = plt.subplots(figsize=(13, 8)) admin.plot(ax=ax, color="#fff4e0", edgecolor="#888", linewidth=0.4, alpha=0.55) # 緊急輸送道路 (薄目) for r in RANK_ORDER[::-1]: sub = gdf_road[gdf_road["rank"] == r] sub.plot(ax=ax, color=RANK_DEF[r]["color"], linewidth=RANK_DEF[r]["weight"] * 0.5, alpha=0.4, zorder=2 + RANK_ORDER.index(r)) # 老朽橋 off (灰, 背景) old_off = gdf_b[gdf_b["is_old_pre1970"] & ~gdf_b["on_road"]] old_off.plot(ax=ax, color="#bbb", markersize=12, alpha=0.55, edgecolor="none", zorder=4) # 老朽橋 on 緊急輸送 (赤、 大きく前面) old_on = gdf_b[gdf_b["is_old_pre1970"] & gdf_b["on_road"]] old_on.plot(ax=ax, color="#cf222e", markersize=35, alpha=0.85, marker="o", edgecolor="#222", linewidth=0.5, zorder=6) ax.set_xlim(-15000, 125000) ax.set_ylim(-220000, -90000) ax.set_aspect("equal") ax.set_title(f"図 8 (RQ3): BCP 脆弱箇所マップ — 老朽橋 (架設 <1970) 全 " f"{n_old70_total} 件 / 緊急輸送道路上 {n_old70_on} 件 " f"({100*n_old70_on/max(n_old70_total,1):.1f}%)", fontsize=10) ax.set_xlabel("X (m, EPSG:6671)") ax.set_ylabel("Y (m, EPSG:6671)") patches = [ Line2D([0], [0], marker='o', color='w', markerfacecolor="#cf222e", markeredgecolor="#222", markersize=10, label=f"老朽橋 on 緊急輸送道路 ({n_old70_on})"), Line2D([0], [0], marker='o', color='w', markerfacecolor="#bbb", markeredgecolor="none", markersize=7, label=f"老朽橋 off 緊急輸送道路 ({len(old_off)})"), Line2D([0], [0], color=RANK_DEF["01"]["color"], linewidth=2.5, label="第1次 緊急輸送道路"), Line2D([0], [0], color=RANK_DEF["02"]["color"], linewidth=2, label="第2次"), Line2D([0], [0], color=RANK_DEF["03"]["color"], linewidth=1.5, label="第3次"), ] ax.legend(handles=patches, loc="lower left", fontsize=9) plt.tight_layout() save_fig("L72_fig8_old_bridges_bcp.png") print(f" ({time.time()-t8:.1f}s)", flush=True) # ============================================================================= # 9. 表データ・仮説検証 # ============================================================================= print("\n[9] 仮説検証", flush=True) t9 = time.time() def df_to_html(d): return d.to_html(index=False, classes="", border=0, escape=False, na_rep="-").replace(' style="text-align: right;"', "") def df_to_html_with_index(d): return d.to_html(classes="", border=0, escape=False, na_rep="-").replace(' style="text-align: right;"', "") # データセット仕様表 zip_size = zip_local.stat().st_size if zip_local.exists() else 0 T_dataset = pd.DataFrame([ ("dataset_id", str(DATASET_ID)), ("公式名", "緊急輸送道路"), ("公式説明", "広島県が管理する道路の緊急輸送道路情報"), ("リソース数", "1 (ZIP)"), ("リソース ID", str(RESOURCE_ID)), ("ZIP サイズ", f"{zip_size:,} byte (~{zip_size/1024:.0f} KB)"), ("ZIP 内 ファイル", "5 JSON (1 階層メタ + 4 階層 LineString)"), ("配信日", "2022-09-26 (公式 stamp)"), ("座標系 (元)", "WGS84 (EPSG:4326) → 本記事 EPSG:6671 で処理"), ("総延長", f"{total_km:.0f} km"), ("総セグメント数", f"{n_seg}"), ("階層数", "4 (本記事独自集約; 公式メタは 7 色)"), ("第1次 (高速・幹線国道)", f"{km_rank['01']:.0f} km / " f"{int(T_rank.loc[T_rank['階層']=='01', 'セグメント数'].iloc[0])} セグ"), ("第2次 (主要地方道・国道)", f"{km_rank['02']:.0f} km / " f"{int(T_rank.loc[T_rank['階層']=='02', 'セグメント数'].iloc[0])} セグ"), ("第3次 (一般県道・市町道)", f"{km_rank['03']:.0f} km / " f"{int(T_rank.loc[T_rank['階層']=='03', 'セグメント数'].iloc[0])} セグ"), ("補完線", f"{km_rank['04']:.0f} km / " f"{int(T_rank.loc[T_rank['階層']=='04', 'セグメント数'].iloc[0])} セグ"), ("ライセンス", "クリエイティブ・コモンズ表示 (CC-BY)"), ("URL", f"https://hiroshima-dobox.jp/datasets/{DATASET_ID}"), ("作成主体", "広島県 (土木建築局道路整備課・防災担当)"), ], columns=["項目", "値"]) # 全体サマリ表 T_overall = pd.DataFrame([ ("dataset", f"#{DATASET_ID} 緊急輸送道路"), ("総延長 (RQ1)", f"{total_km:.0f} km"), ("セグメント数", f"{n_seg}"), ("階層数", "4"), ("第1次 km (RQ1)", f"{km_rank['01']:.0f}"), ("第2次 km (RQ1)", f"{km_rank['02']:.0f}"), ("第3次 km (RQ1)", f"{km_rank['03']:.0f}"), ("補完 km (RQ1)", f"{km_rank['04']:.0f}"), ("第2次シェア % (RQ1)", f"{share_2:.1f}"), ("市町数 (RQ1)", f"{len(T_city)}"), ("最大カバー市町 (RQ1)", f"{T_city.iloc[0]['市町名']} ({T_city.iloc[0]['延長_km']:.0f} km)"), ("H1 (第2次≥50%) (RQ1)", "強支持" if h1_ok else "反証"), ("浸水重なり km (RQ2)", f"{all_flood_km:.0f} ({all_flood_pct}%)"), ("土砂重なり km (RQ2)", f"{all_sed_km:.0f} ({all_sed_pct}%)"), ("災害いずれか km (RQ2)", f"{all_hazard_km:.0f} ({all_hazard_pct}%)"), ("第1次 浸水 % (RQ2)", f"{flood_1}"), ("第1次 土砂 % (RQ2)", f"{sed_1}"), ("H2 (1次<10% AND 1次<10%) (RQ2)", "強支持" if h2_ok else "反証"), ("H3 (浸水>土砂) (RQ2)", "強支持" if h3_ok else "反証"), ("橋梁 on 緊急 (RQ3)", f"{n_b_on}/{n_bridge} ({b_on_pct}%)"), ("トンネル on 緊急 (RQ3)", f"{n_t_on}/{n_tunnel} ({t_on_pct}%)"), ("老朽橋 pre1970 on 緊急 (RQ3)", f"{n_old70_on}/{n_old70_total}"), ("老朽橋 pre1980 on 緊急 (RQ3)", f"{n_old80_on}/{n_old80_total}"), ("H4 (橋梁≥30% AND トンネル≥50%) (RQ3)", "強支持" if h4_ok else "反証"), ("H5 (老朽橋 pre1970 ≥200) (RQ3)", "強支持" if h5_ok else "反証"), ], columns=["指標", "値"]) T_overall.to_csv(ASSETS / "L72_overall.csv", index=False, encoding="utf-8-sig") # 仮説検証表 def jud(cond, ok="強支持", fail="反証"): return ok if cond else fail T_hypo = pd.DataFrame([ ("H1 第2次が総延長 ≥ 50% (RQ1)", f"観測 = 第2次 {km_rank['02']:.0f} km / " f"全 {total_km:.0f} km = {share_2:.1f}%", jud(h1_ok), f"H1 {jud(h1_ok)}: 第 2 次緊急輸送道路は{km_rank['02']:.0f} km " f"で総延長 {total_km:.0f} km の{share_2:.1f}%を占める。 " f"第 1 次 ({km_rank['01']:.0f} km, " f"{100*km_rank['01']/total_km:.1f}%) は高速道路 + 主要幹線国道に " f"限定される短距離ネットワークで、 県内主要都市 (広島・呉・福山・三原) を " f"結ぶ「骨格」 を形成。 第 2 次は主要地方道 + 一般国道 + 主要市町道で " f"県内全域を網羅する「肉付け」 ネットワークとして " f"第 1 次の {km_rank['02']/max(km_rank['01'], 1):.1f} 倍の延長を持つ。 " f"これは典型的な「太い骨格 + 細かい肉付け」 階層構造。"), ("H2 第1次は浸水<10% AND 土砂<10% (RQ2)", f"観測 = 第1次 浸水 {flood_1}% / 土砂 {sed_1}%", jud(h2_ok), f"H2 {jud(h2_ok)}: 第 1 次緊急輸送道路の重なり率は " f"浸水 {flood_1}% / 土砂 {sed_1}%。 " f"浸水は10% 未満 ({flood_1}%)に収まるが、 土砂は" f"{sed_1}%でわずかに 10% を超え、 H2 (両方とも < 10%) は" f"{('反証' if not h2_ok else '強支持')}となった。 " f"ただし第 2 次の土砂 {float(T_rq2.iloc[1]['土砂いずれか_%'])}% / " f"第 3 次の土砂 {float(T_rq2.iloc[2]['土砂いずれか_%'])}% と比較すると、 " f"第 1 次は依然として最低リスク階層である事は明確。 " f"高速道路中心で高架橋・トンネル・盛土区間が多く、 水位上昇の " f"影響を受けにくい構造であるが、 中国自動車道などの山岳貫通区間で " f"土砂警戒区域とのオーバーラップが発生する。 " f"GIS の polygon-on-line 判定では「重なり」 と判定されても、 " f"高架・トンネル構造のため実際の通行確保はもう少し優位 — " f"道路構造を加味した補正で実効重なり率はさらに下がる可能性がある。"), ("H3 全体 浸水重なり > 土砂重なり (RQ2)", f"観測 = 浸水 {all_flood_pct}% / 土砂 {all_sed_pct}%", jud(h3_ok), f"H3 {jud(h3_ok)}: 緊急輸送道路全体で浸水重なり率は " f"{all_flood_pct}% ({all_flood_km:.0f} km)、 " f"土砂警戒区域 (3 種いずれか) は{all_sed_pct}% ({all_sed_km:.0f} km)。 " f"差は{all_flood_pct - all_sed_pct:+.1f} ポイント。 " f"浸水は河川氾濫想定で沿岸都市・平野部の主要道路が広範囲に " f"指定区域を通過するため重なりが拡大する一方、 土砂は急傾斜地に " f"局所的に分布するためライン状道路との交差が物理的に短い。 " f"これは「浸水は面的脅威 / 土砂は線的脅威」 という性質の違いが " f"道路 × ハザード重なり率に直接反映された結果。"), ("H4 橋梁 ≥ 30% AND トンネル ≥ 50% on 緊急 (RQ3)", f"観測 = 橋梁 {b_on_pct}% / トンネル {t_on_pct}%", jud(h4_ok), f"H4 {jud(h4_ok)}: 県内橋梁 {n_bridge:,} のうち緊急輸送道路 " f"{ROAD_BUFFER_M:.0f}m バッファ内に入る橋梁は{n_b_on} 件 " f"({b_on_pct}%)、 トンネル {n_tunnel} 中" f"{n_t_on} 件 ({t_on_pct}%)。 " f"これは「災害時に通行確保しなければならない構造物」 の規模感を " f"定量化する。 とりわけトンネルの{t_on_pct}%が緊急輸送道路上に " f"位置するという事実は、 トンネル管理 (= 換気・照明・排水・耐震補強) が " f"県の BCP に直結することを意味する。 " f"橋梁 {n_b_on} 件 × トンネル {n_t_on} 件は「災害時優先点検対象 " f"のリスト」 として直接利用可能。"), ("H5 老朽橋 (pre1970) ≥ 200 on 緊急 (RQ3)", f"観測 = 老朽橋 (pre1970) on 緊急 {n_old70_on}", jud(h5_ok), f"H5 {jud(h5_ok)}: 緊急輸送道路上に架設 1970 年以前の老朽橋が " f"{n_old70_on} 件 (全老朽橋 {n_old70_total} 件の " f"{100*n_old70_on/max(n_old70_total, 1):.1f}%) 存在。 " f"これは耐震基準改定 (1980 年) 以前の構造物が県の生命線道路上に " f"大量に残存する事実を示し、 BCP 脆弱箇所の最重要リストとなる。 " f"参考: 1980 以前 (改定前耐震基準) で見ると " f"{n_old80_on} 件に拡大。 " f"これらの橋梁は災害時に倒壊・通行不能となれば「孤立集落 + 緊急車両 " f"進入不能」 の連鎖を引き起こすため、 計画的な耐震補強・架替えが " f"BCP の核心である。"), ], columns=["仮説", "観測値", "判定", "詳細解説"]) T_hypo.to_csv(ASSETS / "L72_hypothesis_check.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t9:.1f}s)", flush=True) # ============================================================================= # 10. HTML 生成 # ============================================================================= print("\n[10] HTML 生成", flush=True) t10 = time.time() # ----- セクション 1: 学習目標と問い ----- sec1 = f"""

本記事の対象 — 「緊急輸送道路」 1 件 単独分析

本記事は DoBoX のデータセット 「緊急輸送道路」 (dataset {DATASET_ID}) 1 件を 単独で取り上げ、 広島県の緊急輸送道路ネットワーク {n_seg} セグメント / 総延長 {total_km:.0f} km / 4 階層を 3 つの独立した研究角度で並列に分析する記事である。 他のシリーズ (橋梁 L66 / トンネル L67 / 道路法面 L71 / 道路規制 L50 / ハザード L11) と本記事は 合体しない。 RQ2 で警戒区域 (L11)、 RQ3 で橋梁・トンネル (L66/L67) を参照するが、 これは「緊急輸送道路の 災害脆弱性 + BCP 脆弱箇所」 を明らかにするための既扱データの従属的参照に 留め、 本記事の主軸はあくまで緊急輸送道路 1 dataset の分析である。

「緊急輸送道路」 とは:
災害発生直後から救助・救急・医療・消火活動および緊急物資・人員等の 輸送のために緊急車両が通行する重要な道路。 都道府県地域防災計画で 指定され、 災害時の「生命線 (ライフライン)」として通行確保が 最優先される。 道路法 (1952) 上の道路 (高速・国道・県道・市町道) のうち、 防災拠点 (病院・警察・消防本部・自衛隊基地・港湾・空港・物資集積所) を相互に結ぶネットワークが指定対象となる。

広島県の緊急輸送道路は4 階層で構成される: 本記事の主要発見 (3 RQ):

独自に定義する用語 (本記事限定)

研究の問い (3 RQ)

  1. RQ1 (主研究): 広島県の緊急輸送道路の構造 — 階層・延長・地理は どう描けるか? 4 階層 × 計 {total_km:.0f} km の LineString ネットワークを セグメント数 × 延長 × 市町別カバー × 地理クラスの 4 軸で集計し、 「県の生命線ネットワーク」 の物理的形状を初めて系統的に記述する。 H1 = 第 2 次が ≥ 50% を検証。
  2. RQ2 (副研究 1): 緊急輸送道路は災害脆弱性 — 浸水・土砂警戒に どれだけ晒されているか? 想定最大規模の浸水区域 + 土砂警戒区域 3 種に 対して点サンプリング sjoin で重なり延長を計算し、 「災害時通行不能 可能区間」 を同定する。 H2 = 第1次<10% / H3 = 浸水>土砂 を検証。
  3. RQ3 (副研究 2): 緊急輸送道路上の橋梁・トンネル・老朽橋は どう分布するか? 30m バッファで sjoin → 件数集計 → 老朽橋 (pre1970) を 抽出し、 BCP 脆弱箇所を同定。 H4 = 橋梁≥30% AND トンネル≥50% / H5 = 老朽橋≥200 を検証。

仮説 (5)

到達点

本記事を読み終えると、 (1) 県の緊急輸送道路 4 階層・総延長 {total_km:.0f} km・ {len(T_city)} 市町 のネットワーク構造を完全に俯瞰、 (2) 浸水想定 + 土砂警戒区域との重なり率 (浸水 {all_flood_pct}% / 土砂 {all_sed_pct}% / いずれか {all_hazard_pct}%) を定量把握、 (3) 緊急輸送道路上の橋梁 {n_b_on} 件・ トンネル {n_t_on} 件・老朽橋 {n_old70_on} 件 の BCP 脆弱箇所リストを 特定できる、 という 3 段階の知識が獲得できる。 これにより県の地域防災計画 (BCP) における道路インフラ管理の優先順位が研究者視点で見えるようになる。

""" # ----- セクション 2: 使用データ ----- sec2 = f"""

本研究で使う 1 つの dataset (1 ZIP リソース) を以下の表に示す。 本データは緊急輸送道路 LineString を JSON 配列形式 (NDJSON 風) で公開しており、 4 階層を別ファイルに分けて配信している点が大きな特徴。

データセット仕様

{df_to_html(T_dataset)}

ZIP 内 5 JSON の内訳

ファイル役割形式件数延長 (km)
05_kinkyu_route.json 階層メタ (色 / 太さ / 線種) dict 配列 (NDJSON 風) 7 階層メタ ({{name, color, weight, type}}) —
05_kinkyu_route_01.json 第 1 次 (高速・幹線国道) line array (NDJSON 風) {int(T_rank.loc[T_rank['階層']=='01', 'セグメント数'].iloc[0])} セグ {km_rank['01']:.1f}
05_kinkyu_route_02.json 第 2 次 (主要地方道・国道) line array (NDJSON 風) {int(T_rank.loc[T_rank['階層']=='02', 'セグメント数'].iloc[0])} セグ {km_rank['02']:.1f}
05_kinkyu_route_03.json 第 3 次 (一般県道・市町道) line array (NDJSON 風) {int(T_rank.loc[T_rank['階層']=='03', 'セグメント数'].iloc[0])} セグ {km_rank['03']:.1f}
05_kinkyu_route_04.json 補完線 line array (NDJSON 風) {int(T_rank.loc[T_rank['階層']=='04', 'セグメント数'].iloc[0])} セグ {km_rank['04']:.1f}

NDJSON 風の解読ポイント

各 JSON ファイルは厳密な JSON ではなく、 「dict / array」 が「,」 区切りで 複数並ぶ NDJSON 風の形式 (= JSON.parse でそのまま読めない)。 読み込みには「テキストを [ と ] でラップして配列化 → json.loads」 という工夫が必要。

# NDJSON 風を JSON 配列としてパース
with open(path, "r", encoding="utf-8") as f:
    text = f.read()
arr = json.loads("[" + text + "]")

各線ファイルの中身: arr は線セグメントの配列で、 各セグメントは [{{"e": 経度, "d": 緯度}}, ...] の点列。 e = easting (経度), d = degrees of latitude (緯度) と推定される。 これを shapely.LineString に変換する。

データの読み筋

""" # ----- セクション 3: ダウンロード ----- sec3 = f"""

本記事の再現に必要なすべてを直リンクで提供する。 HTML だけ読めば学習者が完全再現できることが目標 (要件 A)。

生データ (DoBoX 1 件)

このスクリプト本体

中間 CSV (本記事生成、再利用可)

図 (PNG, 直 DL 可)

""" # ----- セクション 4: RQ1 ----- sec4_code = ''' import json, zipfile from pathlib import Path import geopandas as gpd from shapely.geometry import LineString import pandas as pd DATA_DIR = Path("data/extras/L72_emergency_road") extract_dir = DATA_DIR / "340006_emergency_transport_road_20220908T000000" # (1) ZIP 展開 (既展開なら skip) zf = zipfile.ZipFile(DATA_DIR / "emergency_road.zip") zf.extractall(DATA_DIR) # (2) NDJSON 風を JSON 配列としてパースする helper def load_route_lines(idx): p = extract_dir / f"05_kinkyu_route_{idx}.json" with open(p, "r", encoding="utf-8") as f: text = f.read() arr = json.loads("[" + text + "]") # [...,] でラップ lines = [] for seg in arr: if len(seg) >= 2: # 1 点だけは LineString にできない coords = [(pt["e"], pt["d"]) for pt in seg] # e=経度, d=緯度 lines.append(LineString(coords)) return lines # (3) 4 階層を読み込み GeoDataFrame に統合 records = [] for idx in ["01", "02", "03", "04"]: for ln in load_route_lines(idx): records.append({"rank": idx, "geometry": ln}) gdf_road = gpd.GeoDataFrame(records, crs="EPSG:4326").to_crs("EPSG:6671") gdf_road["length_m"] = gdf_road.geometry.length # (4) 階層別集計 T_rank = (gdf_road.groupby("rank") .agg(セグメント数=("geometry", "count"), 延長_km=("length_m", lambda s: s.sum() / 1000)) .reset_index()) print(T_rank) ''' sec4 = f"""

狙い (RQ1)

RQ1 では「県の生命線ネットワーク」 の物理的構造を初めて系統的に記述する。 具体的には {n_seg} セグメント / 総延長 {total_km:.0f} km の緊急輸送道路を 4 階層 × 市町 × 地理クラスの 3 軸で集計し、 「どの階層がどこを通過するか」 を 1 枚で俯瞰できるようにする。 H1 (第 2 次 ≥ 50%) は階層構造の中心仮説を検証する。

手法 — 4 ステップ

  1. STEP 1: ZIP 展開 + NDJSON 風 JSON パース
    ZIP を zipfile.ZipFile.extractall() で展開。 5 JSON のうち 4 つの線データ (_01〜_04) はNDJSON 風(= dict / array が コンマ区切りで並ぶ非標準形式) なので、 "[" + text + "]" で ラップしてから json.loads() で読み込む。
  2. STEP 2: LineString 構築
    各セグメント = [{{"e": lon, "d": lat}}, ...] の点列。 shapely.LineString([(pt["e"], pt["d"]) for pt in seg]) で LineString に変換。 1 点だけのセグメント (= 始点 = 終点)は除外。
  3. STEP 3: 階層 attr 付与 + 投影
    ファイル名末尾 (_01〜_04) を rank 列として保持。 EPSG:4326 (WGS84) → EPSG:6671 (JGD2011 第 III 系) に to_crs() で投影変換 (距離計算が m 単位で正確になる)。
  4. STEP 4: 階層 / 市町 / 地理クラス 集計
    length_m をセグメント長として持ち、 groupby("rank") で階層別 km、 50m 間隔点サンプリング → 行政界 sjoin で市町別 km、 市町名から地理クラス (沿岸島嶼 / 平野・沿岸都市 / 中山間山地) に分類。

実装 (主要部のみ抜粋)

{code(sec4_code)}

結果 1: 階層別 緊急輸送道路マップ (図 1)

なぜこの図か: 4 階層の緊急輸送道路がどこを通過するかを県全域地図に 重ねて一目で確認したい。 階層を色 (赤=1次, 青=2次, 緑=3次, 橙=補完) で 区別し、 線の太さも階層に対応させて「骨格 → 肉付け → 補完」の階層構造を 直感できるようにする。

{figure("assets/L72_fig1_rank_map.png", f"図 1 (RQ1): 広島県 緊急輸送道路 4 階層マップ")}

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

結果 2: 階層別 延長 + セグメント + 地理クラス (図 2)

なぜこの図か: H1 (第 2 次 ≥ 50%) を直感検証するための階層別 km 棒、 セグメント数の比較、 地理クラスの分布 (沿岸都市 / 中山間 / 沿岸島嶼) を 1 枚で並べて、 県の緊急輸送道路の「形」 を 3 角度から把握する。

{figure("assets/L72_fig2_rank_seg_geo.png", "図 2 (RQ1): 階層別 延長 + セグメント + 地理クラス分布")}

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

結果 3: 市町別 緊急輸送道路カバー (図 3)

なぜこの図か: 「自分の市町にどれだけ緊急輸送道路が走っているか」 を 学習者が直感したい。 choropleth + ランキング棒で市町ごとのカバー量を 2 角度から見せる。 第 1 次のみハイライトすることで「骨格」 路線の通り抜け先も 分かる。

{figure("assets/L72_fig3_city_choropleth.png", "図 3 (RQ1): 市町別 緊急輸送道路カバー")}

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

結果 4: 階層別 + 市町別 サマリ表

階層別サマリ:

{df_to_html(T_rank)}

階層別 表から読み取れること: 第 2 次が {km_rank['02']:.0f} km / {share_2:.1f}%と最大、 第 1 次が {km_rank['01']:.0f} km / {100*km_rank['01']/total_km:.1f}% と次。 セグメント数も第 2 次が {int(T_rank.loc[T_rank['階層']=='02', 'セグメント数'].iloc[0])} と最多で、 県内主要地方道としての細分化が見える。

地理クラス別サマリ:

{df_to_html(T_geo)}

地理クラス 表から読み取れること: 沿岸都市 = "都市帯" として 緊急輸送道路の主要部を占める一方、 中山間 (人口疎) でも数百 km の延長を確保。 これは「人口集中地と中山間地の双方をカバーする」 防災設計の表れ。

市町別 ランキング (Top 15):

{df_to_html(T_city.head(15))}

市町別 表から読み取れること: 県北部の中山間市町 (庄原・三次・東広島・ 安芸太田町) が km 上位に並ぶのは面積の広さが反映された結果。 一方、 広島市 (8 区合計) は最も人口集中地で道路密度が高く、 都市部の緊急輸送道路総延長は 区合算で大きい。

階層 × 市町 (Top 10) クロス:

{df_to_html(T_rank_city)}

階層 × 市町 表から読み取れること: 中山間市町は第 2 次が支配的 (主要地方道 + 一般国道のみ)、 都市部は第 1 次 (高速・国道 2) が太く通過、 第 3 次は都市部の市町道として補完。 階層構造が地理に応じて使い分けられている ことが具体的に見える。

""" # ----- セクション 5: RQ2 ----- sec5_code = ''' # 50m 間隔の点サンプリング → sjoin で災害区域内判定 points = [] for line_id, (_, row) in enumerate(gdf_road.iterrows()): L = row.geometry.length n = max(int(L / 50.0), 1) for i in range(n + 1): d = min(i * 50.0, L) pt = row.geometry.interpolate(d) points.append({"line_id": line_id, "rank": row["rank"], "d_m": d, "geometry": pt}) gdf_pts = gpd.GeoDataFrame(points, crs="EPSG:6671") print(len(gdf_pts)) # 例: 56,000+ pts # 浸水想定 (想定最大規模) Shapefile を読込 → sjoin fl = gpd.read_file("data/extras/flood_shp/shinsui_souteisaidai/shinsui_souteisaidai.shp", encoding="utf-8")[["geometry"]].to_crs("EPSG:6671") hits = gpd.sjoin(gdf_pts, fl, how="left", predicate="within") hit_uniq = hits[hits["index_right"].notna()].drop_duplicates(["line_id", "d_m"]) flood_km = len(hit_uniq) * 50.0 / 1000 # 50m/点 × 点数 / 1000 = km print(f"浸水重なり: {flood_km:.0f} km") # 同様に土砂警戒区域 3 種 (急傾斜 / 土石流 / 地すべり) も sjoin # ...省略 (本記事スクリプト参照) # 階層別重なり率 = 各階層の hit 点数 / 総点数 × 100 rank_total = gdf_pts.groupby("rank").size() flood_by_rank = (hit_uniq.groupby("rank").size().reindex(rank_total.index, fill_value=0)) print((flood_by_rank / rank_total * 100).round(1)) ''' sec5 = f"""

狙い (RQ2)

RQ1 で「緊急輸送道路がどこを通過するか」 は分かったが、 これは位置情報のみ。 RQ2 では「緊急輸送道路がどれだけ災害区域に晒されているか」を空間関係で見る。 具体的には浸水想定 (想定最大規模) + 土砂警戒区域 3 種に対して点サンプリング sjoin で重なり延長を計算し、 「災害時通行不能可能区間」 を同定する。 H2 (第 1 次 < 10%) と H3 (浸水 > 土砂) の 2 仮説を検証する。

手法 — 50m 間隔 点サンプリング × sjoin

狙い: LineString と Polygon の直接交差 (intersection) は 計算が重い (~30 秒)。 代わりに緊急輸送道路を50m 間隔で点に 離散化 ({n_pts:,} 点) → sjoin (point-in-polygon) で災害区域内かを 判定 → 該当点数 × 50m で延長換算する。 これは1 秒未満で完了し、 誤差は ±50m 程度に留まる。

災害区分制度意味本研究との関連
浸水想定 (想定最大規模) 水防法 (改正 2015) 想定し得る最大規模降雨で起こる浸水 河川氾濫を通る平野部・沿岸部の道路面的脅威
急傾斜地崩壊警戒区域 土砂災害防止法 (2001) 30度以上の自然斜面の崩壊リスク区域 山岳道路の線的脅威。 法面 (L71) と直接対応
土石流警戒区域 土砂災害防止法 (2001) 渓流域での土石流リスク区域 渓流横断・並行路線で線的脅威
地すべり警戒区域 土砂災害防止法 (2001) 地すべりリスク区域 地質弱層を含む箇所の道路線的脅威

注意: 本記事は「重なり延長 (km, %)」を「災害時通行不能の可能性」 と 同定するが、 実際の通行可否は道路の構造 (高架 / 橋梁 / 路面標高) や降雨 シナリオ (実時間水位) に左右される。 GIS の polygon-in-polygon 判定は 「リスク存在の有無」 のスクリーニングとして使う。

実装 (主要部)

{code(sec5_code)}

結果 1: 緊急輸送道路 + 浸水想定 + 土砂警戒 重ね合わせマップ (図 4)

なぜこの図か: H2 と H3 を地図上で同時に直感検証するため、 4 階層の緊急輸送道路を線で、 浸水想定 + 土砂警戒 (急傾斜 + 土石流) を半透明 色塗りで重ねる。 「どの階層がどの災害区域に重なるか」 が空間的に分かる。

{figure("assets/L72_fig4_hazard_overlay.png", "図 4 (RQ2): 緊急輸送道路 + 浸水想定 + 土砂警戒区域 重ね合わせ")}

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

結果 2: 階層別 + 全体 重なり率 (図 5)

なぜこの図か: H2 (第 1 次 < 10%) と H3 (浸水 > 土砂) を 1 枚で 直接判定する。 左で階層 × 災害種別の % を棒で並べ、 右で全体合算 (浸水 / 土砂 / いずれか) を比較する。 10% 閾値ラインを入れて H2 が直感判定できる。

{figure("assets/L72_fig5_hazard_pct.png", "図 5 (RQ2): 階層別 + 全体 災害区域 重なり率")}

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

結果 3: 階層別 災害区域重なり 詳細表

階層別 重なり (RQ2 中核データ):

{df_to_html(T_rq2)}

階層別 重なり 表から読み取れること:

""" # ----- セクション 6: RQ3 ----- sec6_code = ''' # 緊急輸送道路 30m バッファを構築 gdf_road["geometry_buf"] = gdf_road.geometry.buffer(30) buf = gpd.GeoDataFrame(gdf_road[["rank", "geometry_buf"]].rename( columns={"geometry_buf": "geometry"}), geometry="geometry", crs="EPSG:6671") # L66 橋梁 + L67 トンネル CSV → POINT df_b = pd.read_csv("lessons/assets/L66_all_bridges.csv", encoding="utf-8-sig") df_b = df_b.dropna(subset=["緯度（10進数）", "経度（10進数）"]) df_b = df_b[df_b["緯度（10進数）"] < 50] gdf_b = gpd.GeoDataFrame(df_b, geometry=[Point(x, y) for x, y in zip(df_b["経度（10進数）"], df_b["緯度（10進数）"])], crs="EPSG:4326").to_crs("EPSG:6671") # sjoin: 橋梁 within 緊急輸送道路 30m バッファ b_on = gpd.sjoin(gdf_b, buf[["rank", "geometry"]], how="left", predicate="within") gdf_b["on_road"] = b_on.groupby(b_on.index)["rank"].apply( lambda s: s.notna().any()) print(f"橋梁 on 緊急輸送: {int(gdf_b['on_road'].sum())}/{len(gdf_b)}") # 老朽橋 (架設 < 1970) on 緊急輸送 gdf_b["is_old_pre1970"] = gdf_b["架設年度"].astype(float) < 1970 n_old_on = int((gdf_b["is_old_pre1970"] & gdf_b["on_road"]).sum()) print(f"老朽橋 on 緊急輸送: {n_old_on}") ''' sec6 = f"""

狙い (RQ3)

RQ2 で「災害区域との重なり」 という面的リスクを見たが、 BCP の もう一つの重要要素は「点的構造物リスク」 — 緊急輸送道路上にある橋梁・ トンネルが災害時に倒壊・通行不能となれば、 区域より広範な機能喪失を引き起こす。 RQ3 では緊急輸送道路 30m バッファ内に入る橋梁 + トンネル + 老朽橋を sjoin で同定し、 BCP 脆弱箇所のリストを定量化する。 H4 (≥30% 橋梁 + ≥50% トンネル) と H5 (老朽橋 ≥ 200) を検証する。

手法 — 30m バッファ + sjoin

狙い: 緊急輸送道路 LineString に対して30m バッファを作成 (= 道路敷地 + 構造物近接帯) し、 橋梁 POINT・トンネル POINT を gpd.sjoin(predicate="within") で重なり判定。 1 構造物が複数階層 バッファに重なる場合は最高階層 (= 数字最小)を採用。

項目値意味
バッファ幅 {ROAD_BUFFER_M:.0f} m 道路敷地 (典型的に 4-15m) + 構造物近接帯 (橋台・取付道路) の上限
橋梁 POINT L66 既扱 {n_bridge:,} 件 架設年度・延長・幅員・点検年度・判定区分付き
トンネル POINT L67 既扱 {n_tunnel} 件 建設年度・延長・幅員付き
老朽橋 (pre1970) 架設 < 1970 = {n_old70_total} 件 耐震基準改定 (1980) 以前 + 高度成長期前の構造物
老朽橋 (pre1980) 架設 < 1980 = {n_old80_total} 件 耐震基準改定 (1980) 以前の構造物 (= 旧耐震)

実装 (主要部)

{code(sec6_code)}

結果 1: 緊急輸送道路 + 橋梁 + トンネル マップ (図 6)

なぜこの図か: 「緊急輸送道路上にどれだけの橋梁・トンネルが乗っているか」 を 1 枚で把握したい。 橋梁 on (赤) / off (灰) を色分け、 トンネル on (紫菱形) / off (灰 菱形) を区別、 緊急輸送道路を背景線に置くことで「赤と紫が緊急輸送道路に乗る」 構図が直感できる。

{figure("assets/L72_fig6_struct_overlay.png", f"図 6 (RQ3): 緊急輸送道路 + 橋梁・トンネル")}

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

結果 2: 階層別 BCP 構造物 + 老朽橋 (図 7)

なぜこの図か: H4 と H5 を 3 角度から検証する: 階層別 橋梁・トンネル数 + 階層別 橋梁密度 (件/km) + 階層別 老朽橋。 これで「どの階層に構造物が集中するか」 「どの階層が老朽橋を多く抱えるか」 が同時に分かる。

{figure("assets/L72_fig7_bcp_struct.png", "図 7 (RQ3): 階層別 橋梁・トンネル + 老朽橋 BCP")}

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

結果 3: 老朽橋 BCP 脆弱箇所マップ (図 8)

なぜこの図か: 「老朽橋がどこに集中しているか」 を地理的に見せる。 緊急輸送道路 on (赤) / off (灰) を色分けすることで、 BCP 観点で優先補強すべき 箇所を地図上で直接見られる。

{figure("assets/L72_fig8_old_bridges_bcp.png", "図 8 (RQ3): 老朽橋 BCP 脆弱箇所マップ")}

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

結果 4: 階層別 BCP 構造物 詳細表

階層別 BCP 構造物 (RQ3 中核):

{df_to_html(T_rq3)}

階層別 BCP 表から読み取れること:

""" # ----- セクション 7: 仮説検証 ----- T_hypo_html = T_hypo.copy() T_hypo_html["判定"] = T_hypo_html["判定"].apply( lambda v: f'{v}') sec7 = f"""

仮説検証総合 (H1〜H5)

本記事冒頭で立てた 5 仮説の検証結果を以下にまとめる。 すべての仮説の検証根拠は本記事中の図表に明示されており、再現可能。

{df_to_html(T_hypo_html)}

主要発見の整理

本記事の独自貢献

  1. 「災害脆弱区間」 概念の定量化: 緊急輸送道路 LineString を 50m 間隔で 点サンプリング → sjoin で災害区域内判定 → 該当点数 × 50m で延長換算する 手法を導入。 LineString-Polygon 直接交差の30 倍高速 (1 秒 vs 30 秒)で 類似精度の集計が可能。
  2. 「BCP 脆弱箇所」 リストの作成: 緊急輸送道路 30m バッファ内の橋梁 ({n_b_on} 件) + トンネル ({n_t_on} 件) + 老朽橋 (pre1970, {n_old70_on} 件) を 個別 CSV で提供。 県の BCP 計画の耐震補強・架替えターゲットリストと して直接利用可能。
  3. 階層別 災害重なり率の体系化: 4 階層 × 4 災害区分 (浸水 + 急傾斜 + 土石流 + 地すべり) のクロス集計で、 「第 1 次 = 高架優位、 第 2 次・3 次 = 地形脆弱」 という設計差を空間データで実証。
  4. L66 + L67 + L11 との横断連携: 緊急輸送道路 (1 dataset) + 橋梁 (L66 既扱) + トンネル (L67 既扱) + 警戒区域 (L11 既扱) の 4 dataset を sjoin で組合わせ、 県の道路インフラ防災ネットを初めて統合的に 定量化した。
  5. NDJSON 風の独自パーサ: 公式 ZIP 内の JSON は非標準形式(= dict / array が「,」 区切り) で json.loads() でそのまま読めない。 "[" + text + "]" でラップして配列化するトリックを発見・公開。

本記事の限界

""" # ----- セクション 8: 発展課題 ----- sec8 = f"""

発展課題 — 結果 X → 新仮説 Y → 課題 Z 形式

発展課題 1 (RQ1 拡張): 緊急輸送道路 vs 全道路網のカバー率

発展課題 2 (RQ1 拡張): 防災拠点 (病院・消防本部・港湾)との接続性検証

発展課題 3 (RQ2 拡張): 河川氾濫 vs 高潮の同時被災シナリオ

発展課題 4 (RQ2 拡張): 降雨シナリオ別 通行可否シミュレーション

発展課題 5 (RQ3 拡張): L71 道路法面 + 第3次・補完路線の沿道斜面リスク統合

発展課題 6 (RQ3 拡張): 判定区分 (?) のマスク影響と老朽橋実態の補正

発展課題 7 (展望): 地震動 + 浸水 + 土砂 同時発生時の通行可能ネットワーク

""" # ----- 統合 ----- sections = [ ("学習目標と問い", sec1), ("使用データ", sec2), ("ダウンロード", sec3), (f"【RQ1】 緊急輸送道路の構造 — 4 階層 × {total_km:.0f} km / " f"第2次 {share_2:.0f}% / {len(T_city)} 市町", sec4), (f"【RQ2】 災害脆弱性 — 浸水 {all_flood_pct}% / 土砂 {all_sed_pct}% / " f"いずれか {all_hazard_pct}% ({all_hazard_km:.0f} km)", sec5), (f"【RQ3】 構造物リスク — 橋梁 {n_b_on}/{n_bridge} ({b_on_pct}%) / " f"トンネル {n_t_on}/{n_tunnel} ({t_on_pct}%) / " f"老朽橋 {n_old70_on}", sec6), ("仮説検証総合", sec7), ("発展課題", sec8), ] html = render_lesson( num=72, title=f"緊急輸送道路 単独 3 研究例分析 — " f"4 階層 / {total_km:.0f} km / " f"BCP 脆弱箇所 (橋梁 {n_b_on} + トンネル {n_t_on} + 老朽橋 {n_old70_on}) を読む", tags=["L72", "緊急輸送道路", "BCP", "事業継続計画", "災害脆弱区間", "バックアップルート", "第1次", "第2次", "第3次", "補完線", "RQ×3", "Format B", "geopandas", "LineString (JSON)", "NDJSON 風 パーサ", "L11連携 (警戒区域)", "L66連携 (橋梁)", "L67連携 (トンネル)", "L71連携 (法面)", "浸水想定", "急傾斜地崩壊", "土石流", "地すべり", "老朽橋", "耐震基準改定 (1980)"], time="50 分", level="中級", data_label=f"DoBoX dataset {DATASET_ID} (ZIP × 1, ~{zip_size/1024:.0f} KB)", sections=sections, script_filename="L72_emergency_road.py", ) OUT_HTML = LESSONS / "L72_emergency_road.html" OUT_HTML.write_text(html, encoding="utf-8") print(f" HTML: {OUT_HTML.name} ({len(html):,} chars)") print(f" ({time.time()-t10:.1f}s)", flush=True) # ============================================================================= # 終了 # ============================================================================= print(f"\n=== 完了 (合計 {time.time()-t_all:.1f} 秒) ===", flush=True)