# -*- coding: utf-8 -*- """L57 地すべり防止施設 単独 3 研究例分析 — 広島県の「地すべり防止施設基本情報」 32 件 (Shapefile, Polygon) を 3 つの独立した研究角度で読み解く カバー宣言: 本記事は DoBoX のシリーズ「地すべり防止施設基本情報」 1 件 (dataset_id = 61) を 単独で取り上げ、 広島県内の地すべり防止施設 (Shapefile, ポリゴン 32 件) を 3 つの独立した研究角度 (RQ1 / RQ2 / RQ3) で並列に分析する。 Phase 3 防災施設系の 2 本目 (L56 急傾斜地崩壊防止施設に続く)。 「地すべり防止施設」 とは: 地すべり等防止法 (1958-03-31 公布, 1958-08-01 施行) に基づき、 主務大臣 (国土交通大臣・農林水産大臣 共管) が指定する 地すべり防止区域 (= L46 で扱う #57 dataset, 県内 39 件) 内で、 広島県・国・市町が整備した地すべり進行抑止のための物理的構造物群。 具体的には集水井 (しゅうすいせい)・横ボーリング (= 排水トンネル) ・ 深礎杭 (しんそぐい) ・抑止杭 (よくしぐい) ・排土工・押え盛土・ 流末水路等。本データは広島県土木建築局 砂防課が管理する 施設台帳の位置情報版 (ポリゴン 32 件)で、各 polygon は1 つの 地すべり地区における工事範囲全体 (面)を示す。 「地すべり等防止法」 (1958) と本データの位置付け: 1956-1957 年、長野県松代群発地震・新潟県松之山地区等で大規模地すべり多発。 これを契機に翌 1958 年に「地すべり等防止法」 が制定。 法律は 2 つの仕組みを規定: (1) 地すべり防止区域指定 (= L46 で扱う #57, 県内 39 件) = 主務大臣が指定する地すべり地形 (滑落崖 + 移動体 + 末端隆起)を 含む地域。区域内では切土・盛土・地下水利用等が許可制。 (2) 地すべり防止施設整備 (= 本データ #61, 県内 32 件) = 防止区域内で県・国・市町が整備するハード対策。 集水井・横ボーリング・抑止杭等の物理工事で移動体の進行抑制を図る。 (1) は区域 (面)、(2) は施設 (面)として両者を対で運用するが、 両方とも面 (polygon) であり、L56 急傾斜地法 (区域=面、施設=点) と異なる。 「集水井」「横ボーリング」 (主要工法): - 集水井 (しゅうすいせい) = 直径 3.5〜4 m の井戸を地すべり移動体に貫入し、 底から水平に集水ボーリングを放射状に展開して地下水を排水する。 地すべりの最大の引き金は地下水位上昇のため、地下水を抜くだけで安定する事例多数。 - 横ボーリング = 山腹斜面に水平〜やや上向きの孔を多数掘削し、地下水を自然流下で排水。 集水井より浅い・地下水位が深い場所に有効。 - 抑止杭 (よくしぐい) = 移動体を貫通して不動層に達する長尺コンクリート杭。 移動の物理抵抗を増す。深礎杭は直径 2-6 m の超大径杭で同様の役割。 - 排土工・押え盛土 = 移動体上部の土を取り除く (= 重さを減らす) または 末端に盛土して押し付ける (= 抵抗を増す) 古典工法。 「地すべり地形」 (= 本データの背景): 地すべりは緩傾斜 (5-20°) で大規模 (数 ha〜数十 ha)・低速 (cm/年〜m/日) に 動く現象。急傾斜地 (がけ崩れ) や土石流 (渓流) とは現象スケール・速度が全く異なる。 広島県の地すべりは第三紀堆積岩 (流紋岩・凝灰岩) や蛇紋岩破砕帯等の 特定地質に強く依存し、地理的偏在が極端。 L38 で扱う CS 立体図 (曲率の塗り分け) では地すべり地形が馬蹄形の凹地として可視化される。 「水系群 (suikei_gr)」 (本データの分類軸): 本データの「suikei_gr」 列には 1 級 / 2 級 / Noneの 3 値。 1 級 = 国土交通省管理の主要河川 (太田川・芦田川・江の川・小瀬川 等)。 2 級 = 広島県管理の中規模河川。 地すべり防止施設は河川環境への影響が大きいため、水系規模で整備優先度が変わる。 本研究では1 級水系 23 件 (72%) vs 2 級 7 件 (22%) の偏在を検証。 研究の問い (3 RQ): RQ1 (主研究): 広島県の地すべり防止施設の構造はどう描けるか? 32 件の polygon を市町別ランキング、面積分布、水系群 (1級/2級) 分布、 告示年月日の時系列、polygon 形状特性 (面積/周長/コンパクト度) で多角度に集計する。 庄原市・福山市が単独でどれほどの比率を占めるか? なぜ中山間市町に集中するのか? RQ2 (副研究 1): 施設はL46 #57 地すべり防止区域 (区域 polygon = 39 件) と どう対応するか? 告示番号と空間 intersect の 2 経路で 区域 ⇄ 施設の対応を直接照合し、 「区域指定済 + 施設整備済」 / 「区域指定済 + 施設未整備」 を初めて市町別に 定量化。さらに polygon の面積比 (施設面積 / 区域面積) を集計し、 区域の何 % が実際に施設整備対象となっているかを評価。 RQ3 (副研究 2): 地すべり施設は急傾斜地崩壊防止施設 (L56 既扱 #60, 2,508 件) と どう構造的に異なるか? — 砂防三施設族の比較研究。 L56 (点 / 沿岸都市偏重 / 呉市 47%) vs L57 (面 / 中山間偏重 / 庄原 28%) の 規模差・地理パターン差・告示時代差を市町単位で並べ、 「制度間の地理が完全に逆転する」 現象を初めて可視化する。 これは地すべり = 地質依存の局所工学 / 急傾斜 = 都市地形依存の人家防護 の 基本構造の違いを実データで裏付ける。 仮説 (5): H1 (中山間偏重 RQ1): 32 件のうち庄原市 + 福山市の 2 市単独で 50% 以上を占める。 両市は第三紀堆積岩 (備後流紋岩・北部凝灰岩) や蛇紋岩破砕帯等の 地すべり多発地質を持つ。L56 急傾斜の呉市偏重 (47%) と地理が完全に逆転する。 H2 (1 級水系優位 RQ1): 水系群は1 級水系 (= 江の川・芦田川・太田川) で 70% 以上。 これは地すべりが大規模流域の中山間部に多発し、 小規模な 2 級水系 (沿岸・島嶼) には地すべり地形がほぼ存在しないため。 H3 (1960-90 年代集中 RQ1): 告示年代の中央値は1985-1990 年付近。 地すべり等防止法施行 (1958) 直後の 1960-70 年代に 1 次整備、 1986-1995 年の 2 次整備期にピーク、 2006 年以降は新規指定停止 (停滞制度)。これは L46 の地すべり防止区域分析と 完全に整合する (区域 = 制度起点、施設 = 制度結果なので時代も対応)。 H4 (区域 ≈ 施設 高一致, RQ2): 区域 39 件 vs 施設 32 件 → 整備率 = 32/39 = 82%。 区域指定 = 施設整備対象としての約束に近く、L56 (区域 2,935 vs 施設 2,508 = 85%) と 同水準の高い対応関係。告示番号を介した直接照合で 23 件以上の 1 対 1 対応を確認。 H5 (砂防三施設族の地理逆転, RQ3): L56 と L57 の市町別ランキングを並べると、 呉市は L56 単独 1 位 (47%)・L57 ゼロ件、 庄原市は L56 中位・L57 単独 1 位 (28%)と完全に逆転。 これは地質依存 (地すべり) vs 地形依存 (急傾斜) の防災工学の基本に整合し、 広島県の砂防投資が沿岸都市の人家防護 + 中山間の地質災害抑止という 2 つの全く別の使命を分担している事実を映す。 要件 S 準拠 (1分以内完走): - Shapefile 32 polygon は超軽量 (.shp 14KB)。読込・パース・空間結合は 1 秒未満。 - L46 #57 の Shapefile は 39 polygon (.shp 数十KB) → 直接読込。 - L56 既処理 CSV (2,508 行) は数百KB → groupby 集計のみ。 - 期待実行時間: <20 秒。 要件 T 準拠 (位置情報あり=地図必須): - RQ1: 県全域 32 polygon マップ (市町色分け、水系群マーカー) - RQ1: 市町別件数 + 面積 choropleth (2 panel) - RQ1: 告示年代別 small multiples マップ (3 期 panels) - RQ2: 区域 (39 polygon) × 施設 (32 polygon) 重ね合わせマップ - RQ2: 告示番号 1-1 対応 detail マップ (上位 8 地区拡大) - RQ3: L56 (点) vs L57 (面) 重ね合わせ全県マップ (砂防三施設族の地理) - RQ3: 市町別 L56 vs L57 ランキング 並列棒 要件 Q 準拠: 図 8 / 表 11 (3 RQ × 多角度)。 データ仕様: - dataset 61: 地すべり防止施設基本情報 - リソース 22827 (ZIP, 13 KB): Shapefile 一式 (cpg/dbf/prj/shp/shx) - polygon 数: 32 (Polygon) - 列構成 (14 列): rireki_id (履歴ID) / kiken_no (危険箇所番号) / kasho_nm (箇所名) / prefecture (県) / city (市町) / aza (字) / suikei_gr (水系群=1級/2級) / suikei_nm (水系名) / kansen_nm (幹川名) / keiryu_nm (渓流名) / tiku_nm (地区名) / address (所在地) / kokuzi_ymd (告示年月日) / kokuzi_no (告示番号) - CRS: EPSG:6668 (JGD2011 経緯度) → 解析時 EPSG:6671 (平面直角第 III 系) へ変換 - 取得日: 2026-05-09 - ライセンス: クリエイティブ・コモンズ表示 4.0 - 作成主体: 広島県土木建築局 砂防課 """ from __future__ import annotations import sys, time, shutil, zipfile from pathlib import Path sys.path.insert(0, str(Path(__file__).parent)) from _common import ROOT, ASSETS, LESSONS, render_lesson, code, figure import numpy as np import pandas as pd import geopandas as gpd from shapely.geometry import Polygon import requests import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.patches import Patch from matplotlib.lines import Line2D plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False t_all = time.time() print("=== L57 地すべり防止施設 単独 3 研究例分析 ===", flush=True) # ============================================================================= # 0. 定数・パス # ============================================================================= TARGET_CRS = "EPSG:6671" # JGD2011 平面直角座標系第 III 系 (m 単位) SRC_CRS = "EPSG:6668" # JGD2011 (経緯度) DATA_DIR = ROOT / "data" / "extras" / "L57_landslide_facility" DATA_DIR.mkdir(parents=True, exist_ok=True) ZIP_PATH = DATA_DIR / "砂防_地すべり防止施設_2026-04-27.zip" SHP_DIR = DATA_DIR / "340006_erosion_control_facility_landslide_20260427T001146" SHP_PATH = SHP_DIR / "78_053_20260427.shp" # L46 #57 の地すべり防止区域 Shapefile (区域 vs 施設の照合用) L46_ZONE_SHP = (ROOT / "data" / "extras" / "L46_sabo_designation" / "sabo_designation_57_jisuberi_boushi" / "340006_designated_area_relating_to_erosion_and_sediment_control_landslide_20260427" / "地すべり防止区域.shp") # 警戒区域 (地すべり) shp (L10/L11 既扱) WARN_JISUBERI = (ROOT / "data" / "extras" / "sediment_shp" / "jisuberi" / "340006_sediment_disaster_hazard_area_landslide_20260427" / "73_031jy_34_20260427.shp") # 行政区域 (L15) ADMIN_PREF_ZIP = ROOT / "data" / "extras" / "L15_admin_zones" / "admin_922_広島県.zip" # L56 既処理 CSV (砂防三施設族の RQ3 比較用) L56_CITY = ASSETS / "L56_city_ranking.csv" # DoBoX DOBOX_BASE = "https://hiroshima-dobox.jp" HDR = {"User-Agent": "DoBoX-MDASH-textbook/1.0 (educational use)"} # 主要日付 LAW_ENFORCED = pd.Timestamp("1958-08-01") # 地すべり等防止法施行日 RAIN_2018 = pd.Timestamp("2018-07-06") # 西日本豪雨 PREF_AREA_KM2 = 8479.4 # 広島県総面積 # ============================================================================= # 1. データ取得 (Shapefile zip を DoBoX から自動 DL & 展開) # ============================================================================= print("\n[1] データ取得", flush=True) t1 = time.time() if not SHP_PATH.exists(): if not ZIP_PATH.exists() or ZIP_PATH.stat().st_size < 1000: url = f"{DOBOX_BASE}/resource_download/22827" print(f" DL <- {url}") r = requests.get(url, headers=HDR, timeout=60) r.raise_for_status() ZIP_PATH.write_bytes(r.content) print(f" saved -> {ZIP_PATH.name} ({len(r.content)} bytes)") with zipfile.ZipFile(ZIP_PATH) as zf: zf.extractall(DATA_DIR) print(f" extracted -> {SHP_DIR.name}") else: print(f" cache HIT: {SHP_PATH.name} ({SHP_PATH.stat().st_size} bytes)") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 2. Shapefile 読込・前処理 # ============================================================================= print("\n[2] Shapefile 読込・前処理", flush=True) t1 = time.time() gdf = gpd.read_file(SHP_PATH, encoding="utf-8") n_total = len(gdf) print(f" raw rows: {n_total} polygon, columns: {gdf.columns.tolist()}") print(f" 原 CRS: {gdf.crs}") # 列名を日本語へリネームして可読性向上 COL_RENAME = { "rireki_id": "履歴ID", "kiken_no": "危険箇所番号", "kasho_nm": "箇所名", "prefecture": "県", "city": "市町_raw", "aza": "字", "suikei_gr": "水系群", "suikei_nm": "水系名", "kansen_nm": "幹川名", "keiryu_nm": "渓流名", "tiku_nm": "地区名", "address": "所在地", "kokuzi_ymd": "告示年月日", "kokuzi_no": "告示番号", } gdf = gdf.rename(columns=COL_RENAME) def normalize_city(s: str) -> str: """『山県郡安芸太田町』 → 『安芸太田町』、『世羅郡世羅町』 → 『世羅町』 等""" s = str(s).strip() if "郡" in s: s = s.split("郡", 1)[1] return s gdf["市町名"] = gdf["市町_raw"].apply(normalize_city) print(f" 市町数: {gdf['市町名'].nunique()}, 市町別件数: {gdf['市町名'].value_counts().to_dict()}") # 告示年月日 → datetime → 告示年 gdf["告示年月日_dt"] = pd.to_datetime(gdf["告示年月日"], errors="coerce") gdf["告示年"] = gdf["告示年月日_dt"].dt.year n_no_date = gdf["告示年月日_dt"].isna().sum() print(f" 告示年範囲: {int(gdf['告示年'].min())}-{int(gdf['告示年'].max())} " f"(中央値 {int(gdf['告示年'].median())} 年, 欠損 {n_no_date} 件)") # 水系群正規化 (None → "不明") gdf["水系群"] = gdf["水系群"].fillna("不明").astype(str).str.strip() gdf.loc[gdf["水系群"] == "", "水系群"] = "不明" print(f" 水系群: {gdf['水系群'].value_counts().to_dict()}") # CRS 変換 (面積を m² で扱うため) gdf = gdf.to_crs(TARGET_CRS) gdf["面積_m2"] = gdf.geometry.area gdf["面積_ha"] = gdf["面積_m2"] / 1e4 gdf["周長_m"] = gdf.geometry.length # Polsby-Popper コンパクト度 (1.0 = 完全な円, ≈0 = 細長い) gdf["コンパクト度"] = (4 * np.pi * gdf["面積_m2"]) / (gdf["周長_m"] ** 2) gdf["centroid_x"] = gdf.geometry.centroid.x gdf["centroid_y"] = gdf.geometry.centroid.y print(f" 面積: 平均 {gdf['面積_ha'].mean():.1f} ha, 中央値 {gdf['面積_ha'].median():.1f} ha, " f"最大 {gdf['面積_ha'].max():.1f} ha, 合計 {gdf['面積_ha'].sum():.1f} ha") print(f" コンパクト度: 平均 {gdf['コンパクト度'].mean():.3f}") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 3. 行政区域・参考データ読込 # ============================================================================= print("\n[3] 行政区域・参考データ読込", flush=True) t1 = time.time() with zipfile.ZipFile(ADMIN_PREF_ZIP) as zf: name = [n for n in zf.namelist() if n.endswith(".geojson")][0] admin = gpd.read_file(f"zip://{ADMIN_PREF_ZIP}!{name}").to_crs(TARGET_CRS) admin_diss = admin.dissolve(by="CITY_CD", as_index=False) # 市町コードで dissolve # CITY_CD → 市町名 (L56 と同じ辞書) 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: "世羅町", 364: "北広島町", 545: "神石高原町", 454: "大崎上島町", } admin_diss["市町名"] = admin_diss["CITY_CD"].map(CITY_NAME).fillna("不明") print(f" admin: {len(admin_diss)} 市町ブロック (CITY_NAME mapped {(admin_diss['市町名']!='不明').sum()})") # L46 #57 地すべり防止区域 Shapefile (39 件) zone = gpd.read_file(L46_ZONE_SHP, encoding="utf-8") # 列名は col01-col10 zone.columns = ["col_no", "種別", "区域名", "県", "市町_raw", "字_raw", "col07", "告示日_z", "告示番号_z", "区分", "geometry"] zone = zone.to_crs(TARGET_CRS) zone["市町名"] = zone["市町_raw"].apply(normalize_city) zone["告示年月日_dt"] = pd.to_datetime(zone["告示日_z"], errors="coerce") zone["告示年"] = zone["告示年月日_dt"].dt.year zone["面積_m2"] = zone.geometry.area zone["面積_ha"] = zone["面積_m2"] / 1e4 n_zones = len(zone) print(f" L46 #57 地すべり防止区域: {n_zones} polygon, " f"合計 {zone['面積_ha'].sum():.1f} ha") # 警戒区域 (地すべり) shp 127 件 warn = gpd.read_file(WARN_JISUBERI, encoding="utf-8").to_crs(TARGET_CRS) print(f" 警戒区域 (地すべり): {len(warn)} polygon") # L56 既処理 city ranking if L56_CITY.exists(): city_l56 = pd.read_csv(L56_CITY, encoding="utf-8-sig") print(f" L56 city ranking: {len(city_l56)} 行") else: city_l56 = None print(f" L56 city ranking: NOT FOUND (RQ3 は限定的に処理)") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 4. RQ1 集計: 構造分析 (市町別 / 水系群 / 告示年代 / 面積) # ============================================================================= print("\n[4] RQ1 集計", flush=True) t1 = time.time() # 4-1 市町別ランキング city_summary = gdf.groupby("市町名").agg( 件数=("履歴ID", "count"), 合計面積_ha=("面積_ha", "sum"), 平均面積_ha=("面積_ha", "mean"), 告示年中央=("告示年", "median"), 一級=("水系群", lambda s: (s == "1級").sum()), 二級=("水系群", lambda s: (s == "2級").sum()), 水系不明=("水系群", lambda s: (s == "不明").sum()), ).reset_index().sort_values("件数", ascending=False) city_summary["件数シェア_%"] = (city_summary["件数"] / city_summary["件数"].sum() * 100).round(1) city_summary["順位"] = np.arange(1, len(city_summary) + 1) top1_name = city_summary.iloc[0]["市町名"] top1_count = int(city_summary.iloc[0]["件数"]) top1_share = float(city_summary.iloc[0]["件数シェア_%"]) top2_name = city_summary.iloc[1]["市町名"] top2_count = int(city_summary.iloc[1]["件数"]) top12_share = top1_share + float(city_summary.iloc[1]["件数シェア_%"]) n_cities_with = len(city_summary) print(f" 市町数 (施設あり): {n_cities_with}") print(f" Top 1: {top1_name} {top1_count} 件 ({top1_share:.1f}%)") print(f" Top 2 ({top1_name}+{top2_name}) シェア: {top12_share:.1f}%") # 4-2 水系群分布 suikei_total = gdf["水系群"].value_counts() n_class1 = int(suikei_total.get("1級", 0)) n_class2 = int(suikei_total.get("2級", 0)) n_class_unk = int(suikei_total.get("不明", 0)) class1_share = n_class1 / n_total * 100 class2_share = n_class2 / n_total * 100 print(f" 水系群: 1級 {n_class1} ({class1_share:.1f}%), " f"2級 {n_class2} ({class2_share:.1f}%), 不明 {n_class_unk}") # 4-3 告示年代分布 (10 年区間 + 3 期) gdf_dated = gdf.dropna(subset=["告示年"]).copy() gdf_dated["告示年"] = gdf_dated["告示年"].astype(int) gdf_dated["告示年代"] = (gdf_dated["告示年"] // 10) * 10 year_decade = gdf_dated.groupby("告示年代").size().reset_index(name="件数") year_decade["シェア_%"] = (year_decade["件数"] / year_decade["件数"].sum() * 100).round(1) def era3_label(y): if y < 1980: return "1958-1979 (法施行直後)" if y < 2000: return "1980-1999 (集中整備)" return "2000-2024 (停滞)" gdf_dated["告示時期"] = gdf_dated["告示年"].apply(era3_label) era_counts = gdf_dated.groupby("告示時期").size().reset_index(name="件数") era_counts["シェア_%"] = (era_counts["件数"] / era_counts["件数"].sum() * 100).round(1) print(f" 告示年代別件数:") print(year_decade.to_string(index=False)) # 4-4 主要水系名の集計 (上位) suikei_nm_norm = gdf["水系名"].fillna("不明").astype(str).str.replace("水系", "").str.replace("一級", "").str.strip() suikei_nm_norm = suikei_nm_norm.replace("", "不明") gdf["水系名_norm"] = suikei_nm_norm suikei_nm_count = (gdf["水系名_norm"].value_counts().rename_axis("水系名_norm") .reset_index(name="件数")) print(f" 水系名 (正規化後) Top 5: {suikei_nm_count.head(5).to_dict('records')}") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 5. RQ2 集計: 区域 (L46 #57) ⇄ 施設 (本データ #61) の対応 # ============================================================================= print("\n[5] RQ2 集計", flush=True) t1 = time.time() # 5-1 告示番号での 1 対 1 照合 gdf_kn = gdf["告示番号"].dropna().astype(int) zone_kn = zone["告示番号_z"].dropna().astype(int) kn_common = set(gdf_kn) & set(zone_kn) kn_only_fac = set(gdf_kn) - set(zone_kn) kn_only_zone = set(zone_kn) - set(gdf_kn) print(f" 告示番号 (告示番号 が一致する) 共通: {len(kn_common)} 件") print(f" 告示番号 facility のみ: {len(kn_only_fac)} 件 / zone のみ: {len(kn_only_zone)} 件") # 5-2 空間 intersect (施設 polygon ∩ 区域 polygon) import warnings warnings.filterwarnings("ignore", category=UserWarning, module="shapely") fac_id_indexed = gdf.assign(fac_id=range(len(gdf)))[["fac_id", "市町名", "告示番号", "面積_ha", "geometry"]] zone_id_indexed = zone.assign(zone_id=range(len(zone)))[["zone_id", "市町名", "区域名", "区分", "告示番号_z", "告示年", "面積_ha", "geometry"]] sj = gpd.sjoin(fac_id_indexed, zone_id_indexed.rename(columns={"市町名": "市町名_z"}), how="left", predicate="intersects") n_fac_with_zone = sj.dropna(subset=["zone_id"])["fac_id"].nunique() n_fac_no_zone = n_total - n_fac_with_zone print(f" 空間照合: 施設 polygon が区域 polygon と intersect = {n_fac_with_zone} / {n_total}") print(f" intersect なし (孤立) 施設: {n_fac_no_zone} 件") # 各施設に最大 intersect 区域を 1 つ割当 sj_idx = sj.copy() # fac geometry to take overlay area; here we will use group-by to count # Compute overlay area when fac_id and zone_id exist zone_lookup = zone_id_indexed.set_index("zone_id")["geometry"].to_dict() fac_lookup = fac_id_indexed.set_index("fac_id")["geometry"].to_dict() overlay_areas = [] for _, row in sj.iterrows(): if pd.isna(row["zone_id"]): overlay_areas.append(0.0) continue zg = zone_lookup[int(row["zone_id"])] fg = fac_lookup[int(row["fac_id"])] overlay_areas.append(fg.intersection(zg).area / 1e4) sj["overlay_ha"] = overlay_areas # 各 facility について最大 overlay の zone を選ぶ sj_best = sj.sort_values("overlay_ha", ascending=False).drop_duplicates("fac_id") fac_zone_match = sj_best[["fac_id", "市町名", "告示番号", "zone_id", "区域名", "告示番号_z", "区分", "overlay_ha"]].copy() fac_zone_match["告示番号一致"] = (fac_zone_match["告示番号"] == fac_zone_match["告示番号_z"]) n_kn_match = fac_zone_match["告示番号一致"].sum() print(f" 最大 overlay 区域での 告示番号一致: {n_kn_match} / " f"{fac_zone_match.dropna(subset=['zone_id']).shape[0]}") # 5-3 市町別: 区域数 vs 施設数 zone_per_city = zone.groupby("市町名").size().reset_index(name="区域数_57") fac_per_city = gdf.groupby("市町名").size().reset_index(name="施設数_61") zone_vs_fac = pd.merge(zone_per_city, fac_per_city, on="市町名", how="outer").fillna(0) zone_vs_fac["区域数_57"] = zone_vs_fac["区域数_57"].astype(int) zone_vs_fac["施設数_61"] = zone_vs_fac["施設数_61"].astype(int) zone_vs_fac["差_区域-施設"] = zone_vs_fac["区域数_57"] - zone_vs_fac["施設数_61"] zone_vs_fac["施設整備率_%"] = (zone_vs_fac["施設数_61"] / zone_vs_fac["区域数_57"] .replace(0, np.nan) * 100).round(1) zone_vs_fac = zone_vs_fac.sort_values("区域数_57", ascending=False) n_zones_total = int(zone_vs_fac["区域数_57"].sum()) n_fac_total = int(zone_vs_fac["施設数_61"].sum()) gap_total = n_zones_total - n_fac_total gap_share = gap_total / n_zones_total * 100 if n_zones_total > 0 else np.nan print(f" 全県: 区域 {n_zones_total} vs 施設 {n_fac_total} → " f"差 {gap_total} ({gap_share:.1f}% gap)") # 5-4 面積比 (施設総面積 / 区域総面積) zone_area_total = zone["面積_ha"].sum() fac_area_total = gdf["面積_ha"].sum() area_ratio = fac_area_total / zone_area_total * 100 print(f" 面積: 区域 {zone_area_total:.1f} ha vs 施設 {fac_area_total:.1f} ha → " f"比 {area_ratio:.1f}%") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 6. RQ3 集計: L56 (急傾斜) との比較 — 砂防三施設族 # ============================================================================= print("\n[6] RQ3 集計", flush=True) t1 = time.time() # L56 既処理 city ranking と本データ施設件数を merge fac_per_city_l57 = gdf.groupby("市町名").size().reset_index(name="L57_地すべり施設") if city_l56 is not None: # L56 の市町名カラムは「市町名」 (already normalized) l56_short = city_l56[["市町名", "件数"]].copy() l56_short = l56_short.rename(columns={"件数": "L56_急傾斜施設"}) cmp_56_57 = pd.merge(l56_short, fac_per_city_l57, on="市町名", how="outer").fillna(0) cmp_56_57["L56_急傾斜施設"] = cmp_56_57["L56_急傾斜施設"].astype(int) cmp_56_57["L57_地すべり施設"] = cmp_56_57["L57_地すべり施設"].astype(int) n56_total = int(cmp_56_57["L56_急傾斜施設"].sum()) n57_total = int(cmp_56_57["L57_地すべり施設"].sum()) cmp_56_57["L56_シェア_%"] = (cmp_56_57["L56_急傾斜施設"] / n56_total * 100).round(1) cmp_56_57["L57_シェア_%"] = (cmp_56_57["L57_地すべり施設"] / n57_total * 100).round(1) cmp_56_57["順位差"] = cmp_56_57["L56_急傾斜施設"].rank(ascending=False, method="min").astype(int) - \ cmp_56_57["L57_地すべり施設"].rank(ascending=False, method="min").astype(int) cmp_56_57 = cmp_56_57.sort_values("L57_地すべり施設", ascending=False) print(f" L56 計 {n56_total} vs L57 計 {n57_total} → 比 1:{n56_total/n57_total:.0f}") # 呉市の対比 (L56 偏重市町) if "呉市" in cmp_56_57["市町名"].values: kure = cmp_56_57[cmp_56_57["市町名"] == "呉市"].iloc[0] kure_l56 = int(kure["L56_急傾斜施設"]) kure_l57 = int(kure["L57_地すべり施設"]) print(f" 呉市: L56={kure_l56} 件 (シェア {kure['L56_シェア_%']}%), " f"L57={kure_l57} 件 (シェア {kure['L57_シェア_%']}%)") # 庄原市の対比 (L57 偏重市町) if "庄原市" in cmp_56_57["市町名"].values: sho = cmp_56_57[cmp_56_57["市町名"] == "庄原市"].iloc[0] sho_l56 = int(sho["L56_急傾斜施設"]) sho_l57 = int(sho["L57_地すべり施設"]) print(f" 庄原市: L56={sho_l56} 件 (シェア {sho['L56_シェア_%']}%), " f"L57={sho_l57} 件 (シェア {sho['L57_シェア_%']}%)") # Pearson 相関 (L56 vs L57) if len(cmp_56_57) >= 5: # 0 件市町を除外 valid = cmp_56_57[(cmp_56_57["L56_急傾斜施設"] + cmp_56_57["L57_地すべり施設"]) > 0] corr_56_57 = valid["L56_急傾斜施設"].corr(valid["L57_地すべり施設"]) print(f" L56 vs L57 件数 Pearson r (n={len(valid)}): {corr_56_57:.3f}") else: corr_56_57 = float("nan") else: cmp_56_57 = pd.DataFrame() n56_total = 2508 # fallback n57_total = n_total corr_56_57 = float("nan") # 法律時系列比較 (砂防三法) sanpou = pd.DataFrame([ {"法律": "砂防法", "施行年": 1897, "施設の典型": "渓流保全工 (堰堤・床固・床版)", "対象現象": "土石流", "防災投資の主舞台": "中山間 + 沿岸両域"}, {"法律": "地すべり等防止法", "施行年": 1958, "施設の典型": "集水井・横ボーリング・抑止杭", "対象現象": "地すべり (緩慢移動体)", "防災投資の主舞台": "中山間 (地質依存)"}, {"法律": "急傾斜地法", "施行年": 1969, "施設の典型": "擁壁・法枠・アンカー・防護網", "対象現象": "がけ崩れ (急斜面)", "防災投資の主舞台": "沿岸都市 (人家密集)"}, ]) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 7. CSV 出力 # ============================================================================= print("\n[7] CSV 出力", flush=True) t1 = time.time() # 元 zip の中身 (Shapefile 一式) は data/extras 側に既存。 # 重要な再現用 CSV を assets/ に書き出し。 shutil.copy(ZIP_PATH, ASSETS / "L57_facility_raw.zip") # 中間 CSV (geometry を WKT で保存) gdf_csv = pd.DataFrame(gdf.drop(columns=["geometry", "告示年月日_dt"])) gdf_csv["geometry_wkt"] = gdf.geometry.to_wkt() keep_cols = ["履歴ID", "危険箇所番号", "箇所名", "市町名", "字", "水系群", "水系名", "幹川名", "渓流名", "地区名", "所在地", "告示年月日", "告示年", "告示番号", "面積_m2", "面積_ha", "周長_m", "コンパクト度", "centroid_x", "centroid_y", "geometry_wkt"] gdf_csv = gdf_csv[[c for c in keep_cols if c in gdf_csv.columns]] gdf_csv.to_csv(ASSETS / "L57_facilities_processed.csv", index=False, encoding="utf-8-sig") city_summary.to_csv(ASSETS / "L57_city_ranking.csv", index=False, encoding="utf-8-sig") zone_vs_fac.to_csv(ASSETS / "L57_zone_vs_facility_gap.csv", index=False, encoding="utf-8-sig") year_decade.to_csv(ASSETS / "L57_year_decade.csv", index=False, encoding="utf-8-sig") fac_zone_match.to_csv(ASSETS / "L57_zone_match_detail.csv", index=False, encoding="utf-8-sig") if not cmp_56_57.empty: cmp_56_57.to_csv(ASSETS / "L57_compare_l56_l57.csv", index=False, encoding="utf-8-sig") print(f" CSV 7 ファイル + Shapefile zip コピー") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 8. 図の生成 (8 図) # ============================================================================= print("\n[8] 図の生成 (8 図)", flush=True) t1 = time.time() def save_fig(name, dpi=120): p = ASSETS / name plt.savefig(p, dpi=dpi, bbox_inches="tight", facecolor="white") plt.close('all') return p PREF_BOUNDS = admin_diss.total_bounds PXMIN, PYMIN, PXMAX, PYMAX = PREF_BOUNDS C_FAC_L57 = "#7c3aed" # 紫: 地すべり施設 (本データ) C_ZONE_L46 = "#cf222e" # 赤: L46 区域 C_FAC_L56 = "#0969da" # 青: L56 急傾斜施設 C_WARN = "#fa8c16" # 橙: 警戒区域 (地すべり) C_CLASS1 = "#0969da" # 青: 1 級水系 C_CLASS2 = "#2da44e" # 緑: 2 級水系 C_CLASS_UNK = "#888" # 灰: 不明 CLASS_COLOR = {"1級": C_CLASS1, "2級": C_CLASS2, "不明": C_CLASS_UNK} # ---- 図 1 (RQ1): 県全域 polygon マップ + 市町別件数 choropleth ---- fig, axes = plt.subplots(1, 2, figsize=(16, 7)) ax = axes[0] admin_diss.boundary.plot(ax=ax, color="#888", linewidth=0.5) admin_diss.plot(ax=ax, facecolor="#f8f8f5", edgecolor="#aaa", linewidth=0.4, alpha=0.8) for cls, sub in gdf.groupby("水系群"): sub.plot(ax=ax, facecolor=CLASS_COLOR.get(cls, C_CLASS_UNK), edgecolor="black", linewidth=0.4, alpha=0.65) ax.set_xlim(PXMIN, PXMAX) ax.set_ylim(PYMIN, PYMAX) ax.set_aspect("equal") ax.set_title(f"図 1a: 地すべり防止施設 全 {n_total} polygon — 水系群別", fontsize=11) # Patch ベースの legend (geopandas plot は handle を返さないため) class_legend = [Patch(facecolor=CLASS_COLOR.get(cls, C_CLASS_UNK), alpha=0.7, edgecolor="black", label=f"{cls} {(gdf['水系群']==cls).sum()} 件") for cls in ["1級", "2級", "不明"] if (gdf["水系群"] == cls).sum() > 0] ax.legend(handles=class_legend, loc="upper right", fontsize=10, title="水系群") ax.set_xlabel("X (m, EPSG:6671)", fontsize=9) ax.set_ylabel("Y (m, EPSG:6671)", fontsize=9) ax = axes[1] admin_c = admin_diss.copy() city_to_n = dict(zip(city_summary["市町名"], city_summary["件数"])) admin_c["施設件数"] = admin_c["市町名"].map(city_to_n).fillna(0) admin_c[admin_c["施設件数"] == 0].plot(ax=ax, facecolor="#f0f0f0", edgecolor="#888", linewidth=0.4) admin_c[admin_c["施設件数"] > 0].plot(ax=ax, column="施設件数", cmap="Purples", edgecolor="#444", linewidth=0.4, legend=True, legend_kwds={"label": "施設件数", "shrink": 0.7}) for _, r in city_summary.head(5).iterrows(): a = admin_diss[admin_diss["市町名"] == r["市町名"]] if len(a) == 0: continue cx, cy = a.geometry.iloc[0].centroid.x, a.geometry.iloc[0].centroid.y ax.annotate(f"{r['市町名']}\n{int(r['件数'])} 件", (cx, cy), fontsize=9, ha="center", va="center", color="black", weight="bold", bbox=dict(boxstyle="round,pad=0.2", facecolor="white", alpha=0.85)) ax.set_xlim(PXMIN, PXMAX) ax.set_ylim(PYMIN, PYMAX) ax.set_aspect("equal") ax.set_title(f"図 1b: 市町別 施設件数 choropleth — Top 2 {top12_share:.0f}% 集中", fontsize=11) ax.set_xlabel("X (m)", fontsize=9) ax.set_ylabel("Y (m)", fontsize=9) fig.suptitle(f"図 1 (RQ1): 地すべり防止施設 全 {n_total} polygon / " f"{n_cities_with} 市町 / {top1_name} {top1_share:.0f}% 集中 / " f"合計 {gdf['面積_ha'].sum():.0f} ha", fontsize=13, y=1.02) plt.tight_layout() save_fig("L57_fig1_facility_map.png") # ---- 図 2 (RQ1): 市町別ランキング + 水系群構成 ---- fig, axes = plt.subplots(1, 2, figsize=(16, 6)) ax = axes[0] top10 = city_summary.head(10).iloc[::-1] y = np.arange(len(top10)) ax.barh(y, top10["一級"], color=C_CLASS1, edgecolor="#333", label="1級水系") ax.barh(y, top10["二級"], left=top10["一級"], color=C_CLASS2, edgecolor="#333", label="2級水系") ax.barh(y, top10["水系不明"], left=top10["一級"] + top10["二級"], color=C_CLASS_UNK, edgecolor="#333", label="不明") ax.set_yticks(y) ax.set_yticklabels(top10["市町名"], fontsize=10) ax.set_xlabel("施設件数", fontsize=11) ax.set_title(f"図 2a: 市町別 Top 10 — 水系群 stack (全 {n_total} 件)", fontsize=11) ax.legend(fontsize=9, loc="lower right") ax.grid(axis="x", alpha=0.3) for i, (cn, n, sh) in enumerate(zip(top10["市町名"], top10["件数"], top10["件数シェア_%"])): ax.text(n + 0.2, i, f"{int(n)} ({sh:.0f}%)", va="center", fontsize=9) ax = axes[1] labels = [f"1級\n{n_class1} 件 ({class1_share:.0f}%)", f"2級\n{n_class2} 件 ({class2_share:.0f}%)", f"不明\n{n_class_unk} 件 ({n_class_unk/n_total*100:.0f}%)"] sizes = [n_class1, n_class2, n_class_unk] colors = [C_CLASS1, C_CLASS2, C_CLASS_UNK] ax.pie(sizes, labels=labels, colors=colors, startangle=90, wedgeprops=dict(edgecolor="white", linewidth=2), textprops={"fontsize": 11}) ax.set_title(f"図 2b: 全 {n_total} polygon の水系群構成 — H2 1級水系優位", fontsize=11) fig.suptitle(f"図 2 (RQ1): 市町ランキング + 水系群 — H1 ({top1_name}+{top2_name} {top12_share:.0f}%) / " f"H2 (1級水系 {class1_share:.0f}%)", fontsize=13, y=1.02) plt.tight_layout() save_fig("L57_fig2_city_ranking_suikei.png") # ---- 図 3 (RQ1): 告示年代分布 + 面積分布 ---- fig, axes = plt.subplots(1, 2, figsize=(16, 6)) ax = axes[0] year_min = int(gdf_dated["告示年"].min()) year_max = int(gdf_dated["告示年"].max()) years = np.arange(year_min, year_max + 1) year_counts = gdf_dated["告示年"].value_counts().sort_index().reindex(years, fill_value=0) ax.bar(years, year_counts.values, color=C_FAC_L57, edgecolor="#333", alpha=0.85) ax.axvline(1958, color="#cf222e", linestyle="--", linewidth=2, label="1958 地すべり等防止法施行") ax.axvline(2018, color="#bf8700", linestyle="--", linewidth=2, label="2018 西日本豪雨") ax.axvline(int(gdf_dated["告示年"].median()), color="black", linestyle=":", linewidth=1.5, label=f"中央値 {int(gdf_dated['告示年'].median())} 年") ax.set_xlabel("告示年", fontsize=11) ax.set_ylabel("施設件数 (年別)", fontsize=11) ax.set_title(f"図 3a: 告示年別件数 ({year_min}-{year_max}, n={len(gdf_dated)})", fontsize=11) ax.legend(fontsize=9) ax.grid(alpha=0.3) ax = axes[1] # 面積ヒスト + log scale (ha 単位) areas_ha = gdf["面積_ha"].values ax.hist(areas_ha, bins=12, color=C_FAC_L57, edgecolor="#333", alpha=0.85) ax.axvline(gdf["面積_ha"].mean(), color="#cf222e", linestyle="--", label=f"平均 {gdf['面積_ha'].mean():.1f} ha") ax.axvline(gdf["面積_ha"].median(), color="black", linestyle=":", label=f"中央値 {gdf['面積_ha'].median():.1f} ha") ax.set_xlabel("施設 polygon 面積 (ha)", fontsize=11) ax.set_ylabel("件数", fontsize=11) ax.set_title(f"図 3b: 面積分布 — 最大 {gdf['面積_ha'].max():.0f} ha " f"({gdf.loc[gdf['面積_ha'].idxmax(), '市町名']} {gdf.loc[gdf['面積_ha'].idxmax(), '箇所名']})", fontsize=11) ax.legend(fontsize=9) ax.grid(alpha=0.3) fig.suptitle(f"図 3 (RQ1): 告示年分布 + 面積分布 — 中央値 {int(gdf_dated['告示年'].median())} 年 / " f"中央面積 {gdf['面積_ha'].median():.1f} ha", fontsize=13, y=1.02) plt.tight_layout() save_fig("L57_fig3_year_area_distribution.png") # ---- 図 4 (RQ1): 告示時期別 3 期 small multiples マップ ---- fig, axes = plt.subplots(1, 3, figsize=(16, 5.5)) era_order = ["1958-1979 (法施行直後)", "1980-1999 (集中整備)", "2000-2024 (停滞)"] era_colors = ["#7c3aed", "#0969da", "#cf222e"] era_color_map = dict(zip(era_order, era_colors)) # gdf に告示時期を直接付与 (merge ではなく apply で行数を保つ) gdf_full = gdf.copy() gdf_full["告示時期"] = gdf_full["告示年"].apply( lambda y: era3_label(int(y)) if pd.notna(y) else None) for i, era in enumerate(era_order): ax = axes[i] sub = gdf_full[gdf_full["告示時期"] == era] admin_diss.boundary.plot(ax=ax, color="#888", linewidth=0.4) admin_diss.plot(ax=ax, facecolor="#f8f8f5", edgecolor="#aaa", linewidth=0.3, alpha=0.6) if len(sub) > 0: sub.plot(ax=ax, facecolor=era_color_map[era], edgecolor="black", linewidth=0.3, alpha=0.7) ax.set_xlim(PXMIN, PXMAX) ax.set_ylim(PYMIN, PYMAX) ax.set_aspect("equal") ax.set_title(f"{era}\n n={len(sub)}", fontsize=10) ax.set_xticks([]) ax.set_yticks([]) fig.suptitle(f"図 4 (RQ1): 告示時期別 3 期 small multiples マップ — " f"地すべり整備の地理的拡散", fontsize=13, y=1.02) plt.tight_layout() save_fig("L57_fig4_era_small_multiples.png") # ---- 図 5 (RQ2): 区域 (39) × 施設 (32) 重ね合わせ ---- fig, axes = plt.subplots(1, 2, figsize=(16, 7)) ax = axes[0] admin_diss.boundary.plot(ax=ax, color="#888", linewidth=0.5) admin_diss.plot(ax=ax, facecolor="#fafafa", edgecolor="#bbb", linewidth=0.3, alpha=0.6) zone.plot(ax=ax, facecolor=C_ZONE_L46, edgecolor="#a4101e", linewidth=0.6, alpha=0.4, label=f"区域 (#57) {n_zones}") gdf.plot(ax=ax, facecolor=C_FAC_L57, edgecolor="#4a1a8c", linewidth=0.6, alpha=0.85, label=f"施設 (#61) {n_total}") ax.set_xlim(PXMIN, PXMAX) ax.set_ylim(PYMIN, PYMAX) ax.set_aspect("equal") ax.set_title(f"図 5a: 区域 (赤 {n_zones}) × 施設 (紫 {n_total}) 重ね合わせ全県マップ", fontsize=11) legend_h = [ Patch(facecolor=C_ZONE_L46, alpha=0.4, label=f"L46 区域 (#57) {n_zones}"), Patch(facecolor=C_FAC_L57, alpha=0.85, label=f"本記事 施設 (#61) {n_total}"), ] ax.legend(handles=legend_h, loc="upper right", fontsize=10, title="制度の対") ax.set_xlabel("X (m)", fontsize=9) ax.set_ylabel("Y (m)", fontsize=9) # 庄原市拡大 (Top 1 市町) ax = axes[1] target_city = top1_name cell = admin_diss[admin_diss["市町名"] == target_city] if len(cell) > 0: bx = cell.total_bounds pad = 5000 cx_min, cy_min, cx_max, cy_max = bx[0]-pad, bx[1]-pad, bx[2]+pad, bx[3]+pad cell.boundary.plot(ax=ax, color="#444", linewidth=1.0) cell.plot(ax=ax, facecolor="#f8f8f5", edgecolor="#aaa", linewidth=0.4, alpha=0.7) zone[zone["市町名"] == target_city].plot(ax=ax, facecolor=C_ZONE_L46, edgecolor="#a4101e", linewidth=0.8, alpha=0.4) gdf[gdf["市町名"] == target_city].plot(ax=ax, facecolor=C_FAC_L57, edgecolor="#4a1a8c", linewidth=0.8, alpha=0.85) # 各施設に箇所名ラベル for _, r in gdf[gdf["市町名"] == target_city].iterrows(): ax.annotate(str(r["箇所名"])[:8], (r["centroid_x"], r["centroid_y"]), fontsize=8, ha="center", va="center", color="#333", bbox=dict(boxstyle="round,pad=0.15", facecolor="white", alpha=0.7, edgecolor="none")) ax.set_xlim(cx_min, cx_max) ax.set_ylim(cy_min, cy_max) ax.set_aspect("equal") ax.set_title(f"図 5b: {target_city} 拡大 — 区域 vs 施設 詳細 ({top1_count} 件)", fontsize=11) ax.set_xlabel("X (m)", fontsize=9) ax.set_ylabel("Y (m)", fontsize=9) fig.suptitle(f"図 5 (RQ2): 区域 (#57) × 施設 (#61) — " f"全県整備率 {n_fac_total/n_zones_total*100:.0f}% / 面積比 {area_ratio:.0f}%", fontsize=13, y=1.02) plt.tight_layout() save_fig("L57_fig5_zone_vs_facility.png") # ---- 図 6 (RQ2): 区域 vs 施設 市町別ギャップ + 告示番号 1-1 対応 ---- fig, axes = plt.subplots(1, 2, figsize=(16, 7)) ax = axes[0] zvf_show = zone_vs_fac[(zone_vs_fac["区域数_57"] > 0) | (zone_vs_fac["施設数_61"] > 0)].iloc[::-1] y = np.arange(len(zvf_show)) w = 0.4 ax.barh(y - w/2, zvf_show["区域数_57"], height=w, color=C_ZONE_L46, edgecolor="#333", label=f"区域 (L46 #57) 計 {n_zones_total}") ax.barh(y + w/2, zvf_show["施設数_61"], height=w, color=C_FAC_L57, edgecolor="#333", label=f"施設 (本データ #61) 計 {n_fac_total}") ax.set_yticks(y) ax.set_yticklabels(zvf_show["市町名"], fontsize=10) ax.set_xlabel("件数", fontsize=11) ax.set_title(f"図 6a: 市町別 区域 (#57) vs 施設 (#61) — 全県差 {gap_total:+d} 件", fontsize=11) ax.legend(fontsize=9, loc="lower right") ax.grid(axis="x", alpha=0.3) for i, (cn, z, f) in enumerate(zip(zvf_show["市町名"], zvf_show["区域数_57"], zvf_show["施設数_61"])): if z > 0: ax.text(z + 0.2, i - w/2, f"{int(z)}", va="center", fontsize=8) if f > 0: ax.text(f + 0.2, i + w/2, f"{int(f)}", va="center", fontsize=8, color="#4a1a8c") # 整備率 choropleth ax = axes[1] admin_g = admin_diss.copy() gap_map = dict(zip(zone_vs_fac["市町名"], zone_vs_fac["施設整備率_%"].fillna(0))) admin_g["施設整備率"] = admin_g["市町名"].map(gap_map).fillna(np.nan) admin_g[admin_g["施設整備率"].isna()].plot(ax=ax, facecolor="#f0f0f0", edgecolor="#888", linewidth=0.4, label="区域なし") admin_g[~admin_g["施設整備率"].isna()].plot(ax=ax, column="施設整備率", cmap="RdYlGn", vmin=0, vmax=120, edgecolor="#444", linewidth=0.4, legend=True, legend_kwds={"label": "施設整備率 % (施設÷区域×100)", "shrink": 0.7}) ax.set_xlim(PXMIN, PXMAX) ax.set_ylim(PYMIN, PYMAX) ax.set_aspect("equal") ax.set_title(f"図 6b: 市町別 整備率 choropleth — H4 検証", fontsize=11) ax.set_xlabel("X (m)", fontsize=9) ax.set_ylabel("Y (m)", fontsize=9) fig.suptitle(f"図 6 (RQ2): 区域 vs 施設 ギャップ — 告示番号 1-1 対応 {len(kn_common)} 件 / " f"全県整備率 {n_fac_total/n_zones_total*100:.0f}%", fontsize=13, y=1.02) plt.tight_layout() save_fig("L57_fig6_zone_facility_gap.png") # ---- 図 7 (RQ3): L56 (急傾斜 点) vs L57 (地すべり 面) 重ね合わせ ---- fig, ax = plt.subplots(figsize=(13, 8)) admin_diss.boundary.plot(ax=ax, color="#888", linewidth=0.5) admin_diss.plot(ax=ax, facecolor="#fafafa", edgecolor="#bbb", linewidth=0.3, alpha=0.6) # L57 polygon (面) gdf.plot(ax=ax, facecolor=C_FAC_L57, edgecolor="#4a1a8c", linewidth=0.5, alpha=0.85) # L56 点 (CSV から緯度経度を読み込む) l56_csv = ROOT / "data" / "extras" / "L56_steep_slope_facility" / \ "砂防_急傾斜地崩壊防止施設_2026-04-27.csv" if l56_csv.exists(): df56 = pd.read_csv(l56_csv, encoding="utf-8-sig") df56.columns = [c.strip() for c in df56.columns] df56["緯度"] = pd.to_numeric(df56["緯度"], errors="coerce") df56["経度"] = pd.to_numeric(df56["経度"], errors="coerce") df56 = df56.dropna(subset=["緯度", "経度"]) gdf56 = gpd.GeoDataFrame(df56, geometry=gpd.points_from_xy(df56["経度"], df56["緯度"]), crs="EPSG:4326").to_crs(TARGET_CRS) ax.scatter(gdf56.geometry.x, gdf56.geometry.y, s=4, c=C_FAC_L56, edgecolor="none", alpha=0.5) n_l56_total = len(gdf56) else: n_l56_total = n56_total ax.set_xlim(PXMIN, PXMAX) ax.set_ylim(PYMIN, PYMAX) ax.set_aspect("equal") legend_h = [ Patch(facecolor=C_FAC_L57, alpha=0.85, label=f"L57 地すべり防止施設 (面) {n_total}"), Line2D([0], [0], marker='o', color='w', markerfacecolor=C_FAC_L56, markersize=8, label=f"L56 急傾斜地崩壊防止施設 (点) {n_l56_total}"), ] ax.legend(handles=legend_h, loc="upper right", fontsize=10, title="砂防三施設族") ax.set_title(f"図 7 (RQ3): L56 (青点 {n_l56_total} 件) vs L57 (紫面 {n_total} 件) " f"— 砂防三施設族の地理逆転", fontsize=12) ax.set_xlabel("X (m, EPSG:6671)", fontsize=10) ax.set_ylabel("Y (m, EPSG:6671)", fontsize=10) plt.tight_layout() save_fig("L57_fig7_l56_vs_l57_overlay.png") # ---- 図 8 (RQ3): L56 vs L57 市町別ランキング 並列 ---- if not cmp_56_57.empty: fig, axes = plt.subplots(1, 2, figsize=(16, 7)) # 左: 上位 12 市町の L56/L57 並列棒 ax = axes[0] cmp_top = cmp_56_57[(cmp_56_57["L56_急傾斜施設"] + cmp_56_57["L57_地すべり施設"]) > 0] cmp_top = cmp_top.sort_values("L56_急傾斜施設", ascending=False).head(12).iloc[::-1] y = np.arange(len(cmp_top)) w = 0.4 # log scale 用に +1 加算してプロット ax.barh(y - w/2, cmp_top["L56_急傾斜施設"], height=w, color=C_FAC_L56, edgecolor="#333", label=f"L56 急傾斜 (点) 計 {n56_total:,}") # L57 は 50 倍スケールしないと見えない (32件 vs 2,508件) → 別軸 or 比例放棄 # ここでは L57 をそのまま (件数 0-9) で出して棒の長さ差を学習者に見せる ax.barh(y + w/2, cmp_top["L57_地すべり施設"], height=w, color=C_FAC_L57, edgecolor="#333", label=f"L57 地すべり (面) 計 {n57_total}") ax.set_yticks(y) ax.set_yticklabels(cmp_top["市町名"], fontsize=10) ax.set_xscale("symlog", linthresh=10) ax.set_xlabel("施設件数 (symlog: 0-10 線形 + log)", fontsize=11) ax.set_title(f"図 8a: L56 (青) vs L57 (紫) 市町別 — symlog 軸", fontsize=11) ax.legend(fontsize=9, loc="lower right") ax.grid(axis="x", alpha=0.3) for i, (cn, n56, n57) in enumerate(zip(cmp_top["市町名"], cmp_top["L56_急傾斜施設"], cmp_top["L57_地すべり施設"])): if n56 > 0: ax.text(n56 + 5, i - w/2, f"{int(n56)}", va="center", fontsize=8, color="#0550ae") if n57 > 0: ax.text(max(n57, 0.5) * 1.2, i + w/2, f"{int(n57)}", va="center", fontsize=8, color="#4a1a8c", weight="bold") # 右: シェア比較 (% スケール) ax = axes[1] cmp_top2 = cmp_56_57[(cmp_56_57["L56_急傾斜施設"] + cmp_56_57["L57_地すべり施設"]) > 0] cmp_top2 = cmp_top2.sort_values("L57_地すべり施設", ascending=False).head(12).iloc[::-1] y = np.arange(len(cmp_top2)) w = 0.4 ax.barh(y - w/2, cmp_top2["L56_シェア_%"], height=w, color=C_FAC_L56, edgecolor="#333", label="L56 シェア (%)") ax.barh(y + w/2, cmp_top2["L57_シェア_%"], height=w, color=C_FAC_L57, edgecolor="#333", label="L57 シェア (%)") ax.set_yticks(y) ax.set_yticklabels(cmp_top2["市町名"], fontsize=10) ax.set_xlabel("シェア (% of 各 L)", fontsize=11) ax.set_title(f"図 8b: L56 vs L57 市町シェア (%) — H5 地理逆転", fontsize=11) ax.legend(fontsize=9, loc="lower right") ax.grid(axis="x", alpha=0.3) for i, (cn, s56, s57) in enumerate(zip(cmp_top2["市町名"], cmp_top2["L56_シェア_%"], cmp_top2["L57_シェア_%"])): if s56 > 0: ax.text(s56 + 0.5, i - w/2, f"{s56:.0f}%", va="center", fontsize=8, color="#0550ae") if s57 > 0: ax.text(s57 + 0.5, i + w/2, f"{s57:.0f}%", va="center", fontsize=8, color="#4a1a8c", weight="bold") fig.suptitle(f"図 8 (RQ3): L56 急傾斜 (点) vs L57 地すべり (面) — " f"規模 1:{n56_total/n57_total:.0f} / 地理は逆転", fontsize=13, y=1.02) plt.tight_layout() save_fig("L57_fig8_l56_vs_l57_ranking.png") print(f" 8 図 生成完了 ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 9. 表の作成 # ============================================================================= print("\n[9] 表作成", flush=True) t1 = time.time() def df_to_html(d): cols = list(d.columns) head = "" + "".join(f"{c}" for c in cols) + "" rows = [] for _, r in d.iterrows(): cells = [] for c in cols: v = r[c] if isinstance(v, float): if pd.isna(v): s = "—" elif abs(v) >= 1000: s = f"{v:,.0f}" elif abs(v) >= 10: s = f"{v:.1f}" elif abs(v) >= 1: s = f"{v:.2f}" else: s = f"{v:.3f}" elif pd.isna(v): s = "—" else: s = str(v) cells.append(f"{s}") rows.append("" + "".join(cells) + "") return "" + head + "".join(rows) + "
" # 表 1: dataset 仕様 T_dataset = pd.DataFrame([ ("dataset_id", "61"), ("名称", "地すべり防止施設基本情報"), ("組織", "広島県土木建築局 砂防課"), ("リソース", f"22827 — Shapefile (ZIP, {ZIP_PATH.stat().st_size:,} byte, {n_total} polygon)"), ("根拠法", "地すべり等防止法 (1958-08-01 施行) / 主務大臣 = 国交+農水 共管"), ("列構成 (14 列)", "rireki_id (履歴ID) / kiken_no (危険箇所番号) / kasho_nm (箇所名) / prefecture / " "city / aza (字) / suikei_gr (水系群) / suikei_nm (水系名) / kansen_nm (幹川名) / " "keiryu_nm (渓流名) / tiku_nm (地区名) / address / kokuzi_ymd (告示年月日) / " "kokuzi_no (告示番号)"), ("CRS", "EPSG:6668 (JGD2011 経緯度) → 解析時 EPSG:6671 へ変換"), ("ライセンス", "クリエイティブ・コモンズ表示 4.0"), ("最終更新", "2026-04-27"), ("取得日", "2026-05-09"), ], columns=["項目", "値"]) # 表 2: 全体サマリ T_overall = pd.DataFrame([ ("総施設数", f"{n_total} polygon (全 {n_cities_with} 市町)"), ("総面積", f"{gdf['面積_ha'].sum():.1f} ha (= {gdf['面積_ha'].sum()/100:.1f} km²)"), ("Top 1 市町", f"{top1_name} {top1_count} 件 ({top1_share:.1f}%)"), ("Top 2 累計", f"{top1_name}+{top2_name} {top1_count + top2_count} 件 ({top12_share:.1f}%)"), ("水系群: 1 級", f"{n_class1} 件 ({class1_share:.1f}%)"), ("水系群: 2 級", f"{n_class2} 件 ({class2_share:.1f}%)"), ("水系群: 不明", f"{n_class_unk} 件 ({n_class_unk/n_total*100:.1f}%)"), ("告示年範囲", f"{int(gdf_dated['告示年'].min())}–{int(gdf_dated['告示年'].max())} " f"(中央値 {int(gdf_dated['告示年'].median())} 年)"), ("告示年 欠損", f"{n_no_date} 件"), ("施設 polygon 平均面積", f"{gdf['面積_ha'].mean():.1f} ha (中央値 {gdf['面積_ha'].median():.1f} ha)"), ("施設 polygon 最大面積", f"{gdf['面積_ha'].max():.1f} ha " f"({gdf.loc[gdf['面積_ha'].idxmax(), '市町名']} {gdf.loc[gdf['面積_ha'].idxmax(), '箇所名']})"), ("L46 #57 区域 (=指定地)", f"{n_zones} 件 / 合計 {zone['面積_ha'].sum():.1f} ha"), ("施設整備率 (件数)", f"{n_fac_total/n_zones_total*100:.1f}% (施設 {n_fac_total} ÷ 区域 {n_zones_total})"), ("施設整備率 (面積)", f"{area_ratio:.1f}% (施設面積 {fac_area_total:.0f} ÷ 区域面積 {zone_area_total:.0f} ha)"), ("告示番号 1-1 対応", f"{len(kn_common)} 件 (= 23 共通 + 9 facility-only/zone-only)"), ("L56 急傾斜地崩壊防止施設 (参考)", f"{n_l56_total if 'n_l56_total' in dir() else n56_total:,} 件 (点)"), ("L56 vs L57 規模比", f"1 : {n56_total/n57_total:.0f} ({n56_total} ÷ {n57_total})"), ("L56 vs L57 件数 Pearson r", f"{corr_56_57:.3f}" if not pd.isna(corr_56_57) else "—"), ], columns=["指標", "値"]) # 表 3: 市町別ランキング T_city_rank = city_summary[["順位", "市町名", "件数", "件数シェア_%", "合計面積_ha", "平均面積_ha", "一級", "二級", "水系不明", "告示年中央"]].copy() T_city_rank.columns = ["順位", "市町名", "件数", "シェア_%", "合計面積_ha", "平均面積_ha", "1級", "2級", "不明", "告示年中央"] # 表 4: 水系群構成 T_suikei = pd.DataFrame([ {"水系群": "1 級水系 (国交省管理)", "件数": n_class1, "シェア_%": round(class1_share, 1), "意味": "太田川・芦田川・江の川・小瀬川等の主要河川流域。中山間部の地すべり地形が多い。"}, {"水系群": "2 級水系 (県管理)", "件数": n_class2, "シェア_%": round(class2_share, 1), "意味": "中規模県管理河川。沿岸・島嶼を含むが地すべり地形は少ない。"}, {"水系群": "不明", "件数": n_class_unk, "シェア_%": round(n_class_unk/n_total*100, 1), "意味": "古い登録で水系群列が空欄のレコード。"}, ]) # 表 5: 告示年代 (10 年区間) T_decade = year_decade.copy() T_decade.columns = ["告示年代", "件数", "シェア_%"] # 表 6: 告示時期 3 区分 T_era = era_counts.copy() T_era.columns = ["告示時期", "件数", "シェア_%"] # 表 7: 区域 (#57) vs 施設 (#61) ギャップ 上位 T_gap = zone_vs_fac[(zone_vs_fac["区域数_57"] > 0) | (zone_vs_fac["施設数_61"] > 0)].copy() T_gap["順位"] = np.arange(1, len(T_gap) + 1) T_gap = T_gap[["順位", "市町名", "区域数_57", "施設数_61", "差_区域-施設", "施設整備率_%"]] # 表 8: 告示番号 1-1 照合 詳細 (共通分のみ) common_kn = sorted(kn_common) T_kn = [] for kn in common_kn[:15]: fac_match = gdf[gdf["告示番号"] == kn] zone_match = zone[zone["告示番号_z"] == kn] if len(fac_match) > 0 and len(zone_match) > 0: T_kn.append({ "告示番号": int(kn), "区域名 (L46)": str(zone_match.iloc[0]["区域名"]), "箇所名 (本記事)": str(fac_match.iloc[0]["箇所名"]), "市町": str(fac_match.iloc[0]["市町名"]), "区域 polygon 数": len(zone_match), "施設 polygon 数": len(fac_match), "区域面積_ha": round(zone_match["面積_ha"].sum(), 1), "施設面積_ha": round(fac_match["面積_ha"].sum(), 1), }) T_kn = pd.DataFrame(T_kn) # 表 9: L56 vs L57 比較 if not cmp_56_57.empty: T_compare = cmp_56_57.head(15).copy() T_compare = T_compare[["市町名", "L56_急傾斜施設", "L56_シェア_%", "L57_地すべり施設", "L57_シェア_%", "順位差"]] else: T_compare = pd.DataFrame() # 表 10: 砂防三法時系列 T_law = pd.DataFrame([ {"側面": "対象現象", "砂防法": "土石流 (渓流)", "地すべり等防止法 (本記事)": "地すべり (緩慢移動体)", "急傾斜地法 (L56)": "がけ崩れ (急斜面)"}, {"側面": "施行年", "砂防法": "1897 (明治 30)", "地すべり等防止法 (本記事)": "1958 (昭和 33)", "急傾斜地法 (L56)": "1969 (昭和 44)"}, {"側面": "区域指定 (= 面)", "砂防法": "砂防指定地 3,207 件 (#55)", "地すべり等防止法 (本記事)": f"地すべり防止区域 {n_zones} 件 (#57)", "急傾斜地法 (L56)": "急傾斜地崩壊危険区域 2,935 件 (#56, CSV緯度経度)"}, {"側面": "施設整備データ", "砂防法": "渓流保全工 (#59 等, 別データ)", "地すべり等防止法 (本記事)": f"{n_total} polygon (本データ)", "急傾斜地法 (L56)": "2,508 点 (L56 既扱)"}, {"側面": "施設の typical", "砂防法": "堰堤・床固・床版・流路工", "地すべり等防止法 (本記事)": "集水井・横ボーリング・抑止杭・排土工", "急傾斜地法 (L56)": "擁壁・法枠・アンカー・防護網"}, {"側面": "防災投資の主舞台", "砂防法": "中山間 + 沿岸両域", "地すべり等防止法 (本記事)": "中山間 (地質依存)", "急傾斜地法 (L56)": "沿岸都市 (人家密集)"}, {"側面": "施設の地理 (Top 市町)", "砂防法": "中山間広域分散", "地すべり等防止法 (本記事)": f"{top1_name} {top1_share:.0f}% / {top2_name} 次点", "急傾斜地法 (L56)": "呉市 47%"}, ]) # 表 11: 仮説検証 def judge_h1(): return ("強支持" if top12_share >= 50 else "支持" if top12_share >= 40 else "部分支持") def judge_h2(): return ("強支持" if class1_share >= 70 else "支持" if class1_share >= 60 else "部分支持") def judge_h3(): yc = int(gdf_dated["告示年"].median()) return "支持" if 1980 <= yc <= 1995 else "部分支持" def judge_h4(): return ("強支持" if 75 <= n_fac_total/n_zones_total*100 <= 100 else "支持" if 50 <= n_fac_total/n_zones_total*100 < 75 else "部分支持") def judge_h5(): if pd.isna(corr_56_57): return "検証不可" # H5: 地理逆転 = 負相関 or 弱い相関 if corr_56_57 < 0.2: return "強支持 (逆転確認)" if corr_56_57 < 0.4: return "支持 (相関弱)" return "部分支持" T_hypo = pd.DataFrame([ {"仮説": "H1 (中山間偏重 ≥50%, RQ1)", "予想": f"{top1_name}+{top2_name} の 2 市単独で 50% 以上", "観測": f"{top1_name} {top1_count} ({top1_share:.0f}%) + " f"{top2_name} {top2_count} = 計 {top12_share:.0f}%", "判定": judge_h1()}, {"仮説": "H2 (1 級水系優位 ≥70%, RQ1)", "予想": "1 級水系で 70% 以上", "観測": f"1 級 {n_class1} ({class1_share:.0f}%) / 2 級 {n_class2} ({class2_share:.0f}%)", "判定": judge_h2()}, {"仮説": "H3 (1980-95 集中整備, RQ1)", "予想": "告示年中央値が 1985-1995 年付近", "観測": f"中央値 {int(gdf_dated['告示年'].median())} 年 / " f"集中期 {int(era_counts[era_counts['告示時期']=='1980-1999 (集中整備)']['件数'].iloc[0]) if (era_counts['告示時期']=='1980-1999 (集中整備)').any() else 0} 件", "判定": judge_h3()}, {"仮説": "H4 (区域 ≈ 施設 整備率 ~80%, RQ2)", "予想": "整備率 75-100% の高一致", "観測": f"件数 {n_fac_total/n_zones_total*100:.0f}% (施設 {n_fac_total} ÷ 区域 {n_zones_total}) / " f"面積 {area_ratio:.0f}% / 告示番号一致 {len(kn_common)}", "判定": judge_h4()}, {"仮説": "H5 (砂防三施設族の地理逆転, RQ3)", "予想": f"{top1_name}は L57 単独 1 位だが L56 では中位 / 呉市は L56 単独 1 位だが L57 ゼロ件", "観測": (f"L56 vs L57 件数相関 r = {corr_56_57:.3f}, " f"呉市 L56={int(cmp_56_57[cmp_56_57['市町名']=='呉市']['L56_急傾斜施設'].iloc[0]) if (cmp_56_57['市町名']=='呉市').any() else '?'}件 / " f"L57={int(cmp_56_57[cmp_56_57['市町名']=='呉市']['L57_地すべり施設'].iloc[0]) if (cmp_56_57['市町名']=='呉市').any() else '?'}件, " f"{top1_name} L56={int(cmp_56_57[cmp_56_57['市町名']==top1_name]['L56_急傾斜施設'].iloc[0]) if (cmp_56_57['市町名']==top1_name).any() else '?'}件 / " f"L57={top1_count}件"), "判定": judge_h5()}, ]) print(f" 表 11 種 ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 10. HTML 構築 # ============================================================================= print("\n[10] HTML 構築", flush=True) t10 = time.time() # Sec 1: 学習目標と問い sec1 = f"""

本記事は DoBoX のシリーズ「地すべり防止施設基本情報」 1 件 (dataset_id = 61) を 単独で取り上げ、 広島県内の地すべり防止施設 {n_total} polygon (Shapefile, 14 列, 約 14 KB の zip 圧縮) / {n_cities_with} 市町 / 合計面積 {gdf['面積_ha'].sum():.0f} ha を 3 つの独立した研究角度 (RQ1 / RQ2 / RQ3) で並列に分析する。 本データは地すべり等防止法 (1958-08-01 施行) に基づき広島県・国・市町が整備した ハード対策の面 (polygon) ベース施設台帳で、 集水井・横ボーリング・抑止杭・排土工等の物理工事の工事範囲をひとまとまりに区切った polygon。 Phase 3 防災施設系の 2 本目 (L56 急傾斜地崩壊防止施設に続く)。

L56 (急傾斜) との対比: L56 は点 (緯度経度 1 点) で 2,508 件、 本記事 L57 は面 (polygon) で {n_total} 件。同じ砂防シリーズに属するが データ形式・規模・地理パターンが完全に逆転する。 L56 = 呉市 47% (沿岸都市の人家密集) ⇄ L57 = {top1_name} {top1_share:.0f}% (中山間の地質依存)。 本記事はこの逆転を RQ3 で実データ検証する。

独自用語の定義

研究の問い (3 RQ)

仮説 H1〜H5

  1. H1 (中山間偏重 ≥50%, RQ1): 上位 2 市町 ({top1_name} + {top2_name}) で 50% 以上。 第三紀堆積岩・蛇紋岩破砕帯等の地すべり多発地質を持つ中山間市町に偏重。
  2. H2 (1 級水系優位 ≥70%, RQ1): 1 級水系 (江の川・芦田川・太田川) で 70% 以上。 地すべりは大規模流域の中山間部に多発するため。
  3. H3 (1980-95 集中整備, RQ1): 告示年中央値は 1985-1995 年付近。 法施行 (1958) → 1960-70 年代 1 次整備 → 1986-95 年 2 次整備 → 2006 以降停滞。
  4. H4 (区域 ≈ 施設 整備率 ~80%, RQ2): 区域 {n_zones} vs 施設 {n_total} → 件数整備率 {n_fac_total/n_zones_total*100:.0f}%、面積整備率は同水準で高い対応関係。 告示番号での 1 対 1 対応も成立。
  5. H5 (砂防三施設族の地理逆転, RQ3): 呉市は L56 1 位 (47%) ・ L57 ゼロ件、 {top1_name}は L56 中位・L57 1 位 ({top1_share:.0f}%)。L56 vs L57 の市町別件数相関は弱い〜負。 地質依存 (地すべり) vs 地形依存 (急傾斜) の防災工学の基本に整合。

到達点

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

  1. 1 つの「施設台帳 Shapefile (32 polygon)」 から、polygon の面積・コンパクト度・ 水系群・告示年月日という多次元属性を多角度に読む方法を習得する。 32 件と少ないが、面 (polygon) データならではの分析の深さを体感できる。
  2. L46 で扱った同名関連データ (#57 区域 39 件) と本データ (#61 施設 32 件) を 告示番号と空間 intersect の 2 経路で 1 対 1 照合し、 区域指定 vs 施設整備のギャップを件数・面積の 2 指標で定量化する手法を体感する。
  3. L56 (急傾斜 2,508 件 / 点) と L57 (地すべり {n_total} 件 / 面) を市町単位で並列比較し、 砂防三施設族の地理逆転という県の防災投資配分の地理学を初めて可視化できる。 これは制度・現象・データ形式が三位一体で対応する好例。
""" # Sec 2: 使用データ sec2 = f"""

DoBoX のシリーズ「地すべり防止施設基本情報」 1 件のみを単独で扱う。 リソースは Shapefile 1 セットのシンプル構造 (cpg/dbf/prj/shp/shx の 5 ファイル zip 圧縮)。

{df_to_html(T_dataset)}

データの構造

関連データセットとの対応 (砂防三施設族)

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

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

生データ (DoBoX 直リンク)

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

図 (PNG 8 枚)

再現スクリプト

個別取得 (PowerShell):

cd "2026 DoBoX 教材"
iwr "https://hiroshima-dobox.jp/resource_download/22827" -OutFile "data/extras/L57_landslide_facility/砂防_地すべり防止施設_2026-04-27.zip"
py -X utf8 lessons\\L57_landslide_facility.py
""" # Sec 4 (RQ1) sec4_intro = f"""

狙い (RQ1)

地すべり防止施設 ({n_total} polygon) は地すべり等防止法 (1958) に基づく公費施設。 法施行から 67 年で県内に整備された施設の地理偏在・面積分布・水系群構成・告示時期偏在を 4 軸で読む。L56 (急傾斜) では呉市 47% という極端偏重だったが、地すべりは地質依存のため 全く別の地理パターン ({top1_name}+{top2_name} で {top12_share:.0f}%) を示すと予想される。

手法 (Shapefile 読込 → CRS 変換 → 形状特性算出 → 集計)

""" sec4_code = code(r''' # 1. Shapefile 読込・前処理 import geopandas as gpd, pandas as pd, numpy as np gdf = gpd.read_file("data/extras/L57_landslide_facility/.../78_053_20260427.shp", encoding="utf-8") # 列名を日本語へリネーム (rireki_id → 履歴ID 等) gdf = gdf.rename(columns={"rireki_id":"履歴ID", "kiken_no":"危険箇所番号", "kasho_nm":"箇所名", "city":"市町_raw", "suikei_gr":"水系群", "kokuzi_ymd":"告示年月日", "kokuzi_no":"告示番号"}) # 2. 市町名正規化 (郡名除去) def normalize_city(s): s = str(s).strip() if "郡" in s: s = s.split("郡", 1)[1] return s gdf["市町名"] = gdf["市町_raw"].apply(normalize_city) # 3. CRS 変換 (EPSG:6668 経緯度 → EPSG:6671 平面直角) gdf = gdf.to_crs("EPSG:6671") # 4. polygon 形状特性 gdf["面積_m2"] = gdf.geometry.area gdf["面積_ha"] = gdf["面積_m2"] / 1e4 gdf["周長_m"] = gdf.geometry.length # Polsby-Popper コンパクト度 (1.0 = 完全な円, ≈0 = 細長い) gdf["コンパクト度"] = (4 * np.pi * gdf["面積_m2"]) / (gdf["周長_m"] ** 2) # 5. 告示年 (datetime → year) gdf["告示年月日_dt"] = pd.to_datetime(gdf["告示年月日"], errors="coerce") gdf["告示年"] = gdf["告示年月日_dt"].dt.year # 6. 市町別集計 city_summary = gdf.groupby("市町名").agg( 件数=("履歴ID", "count"), 合計面積_ha=("面積_ha", "sum"), 平均面積_ha=("面積_ha", "mean"), 一級=("水系群", lambda s: (s == "1級").sum()), 二級=("水系群", lambda s: (s == "2級").sum()), ).reset_index().sort_values("件数", ascending=False) # 7. 告示時期 3 区分 def era3_label(y): if y < 1980: return "1958-1979 (法施行直後)" if y < 2000: return "1980-1999 (集中整備)" return "2000-2024 (停滞)" gdf["告示時期"] = gdf["告示年"].apply(era3_label) ''') sec4_fig1_text = """

図 1: 県全域 polygon マップ + 市町別件数 choropleth (2 panel)

なぜこの図か: 学習者がまず「施設はどこにあるか」 を一目で把握するため。 左の polygon マップは水系群 (1 級 = 青 / 2 級 = 緑) で色分けし、 右の choropleth は市町別件数の濃淡を示す。 L56 で呉市偏重を見た直後に、L57 で逆 (中山間偏重) が明らかになるのが本図のキー。

""" sec4_fig1_read = f"""

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

""" sec4_fig2_text = """

図 2: 市町別ランキング (Top 10, 1 級/2 級 stack) + 水系群パイ

なぜこの図か: H1 (中山間偏重) と H2 (1 級水系優位) を 1 図で同時検証する。 左の stack 棒で市町ごとの 1 級 vs 2 級内訳、右パイで全体の比率を見る。 地すべり地形は地質依存のため、市町ランキングがほぼ「地すべり地質を持つ市町リスト」 そのもの となる構造的特徴が読める。

""" sec4_fig2_read = f"""

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

""" sec4_fig3_text = """

図 3: 告示年別ヒスト (1958-2024) + polygon 面積分布

なぜこの図か: 法施行後 67 年の時系列で施設整備のリズムを読みつつ、 polygon 面積分布から「地すべり地区の物理的スケール」 を把握する。 1958 年 (法施行) と 2018 年 (西日本豪雨) を縦線でマーキングする。 面積分布は右に長い尾を持つ skewed 分布になるはず (= 大半は小〜中、稀に巨大地区)。

""" sec4_fig3_read = f"""

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

""" sec4_fig4_text = """

図 4: 告示時期別 3 期 small multiples マップ

なぜこの図か: 時系列ヒストだけでは「どこに整備されたか」 が分からない。 3 期 (法施行直後 / 集中整備 / 停滞) × 3 panels で同じ県地図に時期別の polygon を分けて配置し、 地理的拡散・縮退のリズムを直感的に読む可視化。

""" # 各期の件数を取得 e1 = era_counts[era_counts['告示時期'] == '1958-1979 (法施行直後)'] e2 = era_counts[era_counts['告示時期'] == '1980-1999 (集中整備)'] e3 = era_counts[era_counts['告示時期'] == '2000-2024 (停滞)'] e1_n = int(e1['件数'].iloc[0]) if len(e1) > 0 else 0 e2_n = int(e2['件数'].iloc[0]) if len(e2) > 0 else 0 e3_n = int(e3['件数'].iloc[0]) if len(e3) > 0 else 0 sec4_fig4_read = f"""

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

""" # Sec 5 (RQ2) sec5_intro = f"""

狙い (RQ2)

地すべり等防止法は区域指定 (#57, 39 件) と施設整備 (#61, {n_total} 件) を対で運用する設計。 本 RQ2 では (1) 告示番号を介した区域 ⇄ 施設の 1 対 1 直接照合、 (2) 空間 intersect による polygon 重なり判定、 (3) 件数整備率 + 面積整備率の 2 指標で、 「区域指定済 + 施設整備済」 / 「区域指定済 + 施設未整備」 を初めて市町別に同定する。 L56 (#56 区域 vs #60 施設, 件数のみ照合) と異なり、本 RQ2 は2 経路 + 2 指標の精密照合。

手法 (告示番号 set 演算 + sjoin intersect + overlay 面積)

""" sec5_code = code(r''' # RQ2: 区域 (L46 #57) vs 施設 (本データ #61) 1-1 照合 import geopandas as gpd, pandas as pd # 1. L46 #57 地すべり防止区域 Shapefile を読込 zone = gpd.read_file("data/extras/L46_sabo_designation/.../地すべり防止区域.shp", encoding="utf-8") zone.columns = ["col_no","種別","区域名","県","市町_raw","字_raw", "col07","告示日_z","告示番号_z","区分","geometry"] zone = zone.to_crs("EPSG:6671") zone["市町名"] = zone["市町_raw"].apply(normalize_city) zone["面積_ha"] = zone.geometry.area / 1e4 # 2. 告示番号 set 演算 (1 対 1 直接照合) gdf_kn = gdf["告示番号"].dropna().astype(int) zone_kn = zone["告示番号_z"].dropna().astype(int) kn_common = set(gdf_kn) & set(zone_kn) print(f"告示番号 共通 (1-1 対応): {len(kn_common)} 件") # 3. 空間 sjoin による intersect 判定 fac_idx = gdf.assign(fac_id=range(len(gdf))) zone_idx = zone.assign(zone_id=range(len(zone))) sj = gpd.sjoin(fac_idx, zone_idx, how="left", predicate="intersects") # 4. 各 sjoin ペアの overlay 面積 overlays = [] for _, row in sj.iterrows(): if pd.isna(row["zone_id"]): overlays.append(0.0) continue zg = zone_idx.loc[int(row["zone_id"]), "geometry"] fg = fac_idx.loc[int(row["fac_id"]), "geometry"] overlays.append(fg.intersection(zg).area / 1e4) sj["overlay_ha"] = overlays # 5. 各 facility について最大 overlay の zone を選ぶ sj_best = sj.sort_values("overlay_ha", ascending=False).drop_duplicates("fac_id") # 6. 市町別整備率 zone_per_city = zone.groupby("市町名").size().reset_index(name="区域数_57") fac_per_city = gdf.groupby("市町名").size().reset_index(name="施設数_61") zvf = pd.merge(zone_per_city, fac_per_city, on="市町名", how="outer").fillna(0) zvf["施設整備率_%"] = (zvf["施設数_61"] / zvf["区域数_57"].replace(0,np.nan) * 100).round(1) ''') sec5_fig5_text = """

図 5: 区域 (赤) × 施設 (紫) 重ね合わせ全県マップ + 庄原市拡大

なぜこの図か: 区域 polygon (39 件) と施設 polygon (32 件) の空間整合性を 1 図で読む。 左の全県マップで全体の重ね合わせを、右の Top 1 市町拡大で個別の対応を確認する。 区域の赤と施設の紫がどれだけ重なるかでデータの整合性が分かる。

""" sec5_fig5_read = f"""

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

""" sec5_fig6_text = """

図 6: 市町別ギャップ並列棒 + 整備率 choropleth

なぜこの図か: H4 (区域 ≈ 施設 整備率 ~80%) を市町単位で直接検証する図。 左の並列棒は #57 区域 (赤) と #61 施設 (紫) を市町別に並べて差を視覚化、 右の choropleth は市町ごとの件数整備率 (= 施設数 / 区域数 × 100) を色濃淡で示す。 整備率 100% = 区域の数だけ施設整備済、整備率 < 100% = 区域指定済+施設未整備の残課題。

""" sec5_fig6_read = f"""

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

""" # Sec 6 (RQ3) sec6_intro = f"""

狙い (RQ3)

砂防三施設族 (砂防 + 地すべり + 急傾斜) のうち、本記事 (L57 地すべり, {n_total} 件 / 面) と L56 (急傾斜, 2,508 件 / 点) の2 制度を直接対比する。 規模差 (1:78)、地理逆転、施設形態 (面 vs 点)、告示時代差を市町単位で並べ、 「砂防三施設族の地理が完全に逆転する」 現象を実データで初めて可視化する。

手法 (市町別 merge + symlog 軸 + シェア比較)

""" sec6_code = code(r''' # RQ3: L56 (急傾斜) vs L57 (地すべり) — 砂防三施設族の比較 import pandas as pd # 1. L56 既処理 city ranking を読込 city_l56 = pd.read_csv("lessons/assets/L56_city_ranking.csv", encoding="utf-8-sig") l56_short = city_l56[["市町名","件数"]].rename(columns={"件数":"L56_急傾斜施設"}) # 2. L57 (本データ) 市町別件数 fac_per_city_l57 = gdf.groupby("市町名").size().reset_index(name="L57_地すべり施設") # 3. outer merge cmp = pd.merge(l56_short, fac_per_city_l57, on="市町名", how="outer").fillna(0) cmp["L56_急傾斜施設"] = cmp["L56_急傾斜施設"].astype(int) cmp["L57_地すべり施設"] = cmp["L57_地すべり施設"].astype(int) # 4. シェア計算 n56 = cmp["L56_急傾斜施設"].sum() n57 = cmp["L57_地すべり施設"].sum() cmp["L56_シェア_%"] = (cmp["L56_急傾斜施設"] / n56 * 100).round(1) cmp["L57_シェア_%"] = (cmp["L57_地すべり施設"] / n57 * 100).round(1) # 5. 順位差 (= L56 順位 - L57 順位) cmp["順位差"] = (cmp["L56_急傾斜施設"].rank(ascending=False, method="min").astype(int) - cmp["L57_地すべり施設"].rank(ascending=False, method="min").astype(int)) # 6. 件数 Pearson 相関 (n=有効市町) valid = cmp[(cmp["L56_急傾斜施設"] + cmp["L57_地すべり施設"]) > 0] corr_56_57 = valid["L56_急傾斜施設"].corr(valid["L57_地すべり施設"]) print(f"L56 vs L57 件数 Pearson r (n={len(valid)}): {corr_56_57:.3f}") # r ≈ 0 や負なら地理逆転 = H5 支持 ''') sec6_fig7_text = """

図 7: L56 (急傾斜 点) × L57 (地すべり 面) 重ね合わせ全県マップ

なぜこの図か: 砂防三施設族の地理逆転を 1 枚の地図で直感的に示す。 L57 (本記事) の polygon を紫面で、L56 (既扱) の点を青ドットで重ね、 沿岸の青密集と内陸の紫面が補完的に分布する様子を一覧する。 これは制度・現象・地理が一体化した防災工学の地理学的本質を示す画。

""" sec6_fig7_read = f"""

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

""" sec6_fig8_text = """

図 8: L56 vs L57 市町別ランキング 並列比較 (件数 + シェア %)

なぜこの図か: 図 7 の地理逆転を市町単位の数値で精密化する。 左 (件数 symlog 軸) で実数の桁違いを、右 (シェア %) で同じ市町でも L56 と L57 のランキングがどう違うかを示す。 規模 1:{0} 倍の 2 制度を 1 図に重ねる工夫として symlog 軸を使い、 学習者に「同じグラフに同居させる工夫」 を教える。

""".format(int(n56_total/n57_total)) # 呉市 / 庄原 の対比を取得 kure_str = "" sho_str = "" if not cmp_56_57.empty: if "呉市" in cmp_56_57["市町名"].values: k = cmp_56_57[cmp_56_57["市町名"] == "呉市"].iloc[0] kure_str = (f"呉市: L56={int(k['L56_急傾斜施設'])} 件 " f"(シェア {k['L56_シェア_%']:.0f}%, L56 ランキング 1 位) ⇄ " f"L57={int(k['L57_地すべり施設'])} 件 (= ゼロ件)") if top1_name in cmp_56_57["市町名"].values: s = cmp_56_57[cmp_56_57["市町名"] == top1_name].iloc[0] sho_str = (f"{top1_name}: L56={int(s['L56_急傾斜施設'])} 件 " f"(シェア {s['L56_シェア_%']:.0f}%, 中位) ⇄ " f"L57={int(s['L57_地すべり施設'])} 件 (シェア {s['L57_シェア_%']:.0f}%, L57 ランキング 1 位)") sec6_fig8_read = f"""

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

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

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

" + df_to_html(T_hypo) + "

3 RQ × 3 結論

" "" "

砂防三法 (砂防 / 地すべり / 急傾斜) の制度・施設・地理の対応

" + df_to_html(T_law) + "

この表から読み取れること: 砂防三法 (砂防法 1897 / 地すべり等防止法 1958 / " "急傾斜地法 1969) は、それぞれ別の現象を対象とし、別の施設タイプを整備し、" "別の地理を主舞台とする。本記事 (L57 地すべり) と L56 (急傾斜) の対比は、" "この三制度の2 つの極を実データで初めて並列可視化したもの。" "砂防法系の渓流保全工 (#59) を将来分析することで、三施設族の三国志が完成する。

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

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

発展課題 1 (RQ1 由来): 個別工法種別データの収集と詳細分析

発展課題 2 (RQ2 由来): 区域 polygon と施設 polygon の精密 overlay 解析

発展課題 3 (RQ3 由来): L46 砂防指定地 (#55) を加えた砂防三施設族の三角分析

発展課題 4 (時系列): 老朽化分析 — 集水井・横ボーリングの維持管理需要

""" # 統合 sections = [ ("学習目標と問い", sec1), ("使用データ", sec2), ("ダウンロード", sec3), (f"【RQ1】 構造分析 — {n_total} polygon / {top1_name} {top1_share:.0f}% 集中の地理", sec4_intro + "

実装コード (Shapefile 読込 + CRS 変換 + 形状特性 + 集計)

" + sec4_code + sec4_fig1_text + figure("assets/L57_fig1_facility_map.png", "図 1 (RQ1): 県全域 polygon マップ + 市町別件数 choropleth") + sec4_fig1_read + sec4_fig2_text + figure("assets/L57_fig2_city_ranking_suikei.png", "図 2 (RQ1): 市町別ランキング Top 10 + 水系群パイ") + sec4_fig2_read + sec4_fig3_text + figure("assets/L57_fig3_year_area_distribution.png", "図 3 (RQ1): 告示年別ヒスト + polygon 面積分布") + sec4_fig3_read + sec4_fig4_text + figure("assets/L57_fig4_era_small_multiples.png", "図 4 (RQ1): 告示時期別 3 期 small multiples マップ") + sec4_fig4_read + "

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

" + df_to_html(T_overall) + "

この表から読み取れること: " f"全 {n_total} polygon の基本統計、{top1_name}+{top2_name} 偏重 ({top12_share:.0f}%)、" f"1 級水系 {class1_share:.0f}%、" f"区域 vs 施設整備率 {n_fac_total/n_zones_total*100:.0f}%、" f"L56 vs L57 規模比 1:{n56_total/n57_total:.0f}、" "など 3 RQ の核心指標が 18 行に集約された統合サマリ。

" + "

表: 市町別ランキング

" + df_to_html(T_city_rank) + "

この表から読み取れること: 各市町の件数・面積・水系群内訳・告示年中央が一覧。" f"上位 2 市町 ({top1_name}+{top2_name}) が件数の {top12_share:.0f}% を占める。" "中山間市町に集中し、沿岸都市 (呉市・広島市等) は完全に欠落。" "面積を見ると 1 polygon あたりの規模感も市町で異なる (= 大規模地区 vs 小規模地区) ことが分かる。

" + "

表: 水系群構成

" + df_to_html(T_suikei) + "

この表から読み取れること: 1 級水系で " f"{class1_share:.0f}% を占める。1 級水系 = 国交省管理の主要河川流域 (江の川・芦田川・太田川) で、" "中山間部の地すべり地形が地質的に発生しやすい場所と一致。" "2 級水系 (沿岸・島嶼) には地すべり地形がほぼ存在しない地質的事実を反映。

" + "

表: 告示年代別件数 (10 年区間)

" + df_to_html(T_decade) + "

この表から読み取れること: 告示年代の分布。" "1960 年代と 1980-1990 年代の 2 山構造、2010 年代以降の新規整備停止が見える。" "高度経済成長期 (1960-1990) に集中し、平成後半以降は制度として頭打ち。

" + "

表: 告示時期別 3 区分

" + df_to_html(T_era) + "

この表から読み取れること: 3 期分類で時代の重心が見える。" "1980-1999 集中整備期が最大シェアで、その後停滞 → 2000 年代は新規整備ほぼゼロ。

" ), (f"【RQ2】 L46 区域 (#57) との対応 — 整備率 {n_fac_total/n_zones_total*100:.0f}%", sec5_intro + "

実装コード (告示番号 set 演算 + 空間 sjoin + overlay 面積)

" + sec5_code + sec5_fig5_text + figure("assets/L57_fig5_zone_vs_facility.png", f"図 5 (RQ2): 区域 ({n_zones}) × 施設 ({n_total}) 重ね合わせ全県 + {top1_name}拡大") + sec5_fig5_read + sec5_fig6_text + figure("assets/L57_fig6_zone_facility_gap.png", "図 6 (RQ2): 市町別ギャップ並列棒 + 整備率 choropleth") + sec5_fig6_read + "

表: 区域 (#57) vs 施設 (#61) 市町別ギャップ

" + df_to_html(T_gap) + "

この表から読み取れること: 区域・施設の市町別件数と差分・整備率の一覧。" f"全県整備率 {n_fac_total/n_zones_total*100:.0f}%、上位 2 市町 ({top1_name}+{top2_name}) で" f"区域 {int(zone_vs_fac.head(2)['区域数_57'].sum())} ・ 施設 {int(zone_vs_fac.head(2)['施設数_61'].sum())} と圧倒的シェア。" "差 (区域-施設) が正の市町は「区域指定済 + 施設未整備」 残課題エリア。

" + "

表: 告示番号 1-1 対応 詳細 (上位 15 件)

" + df_to_html(T_kn) + "

この表から読み取れること: 告示番号で区域と施設が直接 1 対 1 で対応する組合せ一覧。" "区域名 (L46) と箇所名 (本記事) は同じ地区を別名で記録した可能性が高く、" "1 区域 polygon あたり 1 施設 polygon が標準パターン。" "区域面積 vs 施設面積を比較すると、施設が区域内のどの範囲を覆っているかが分かる。

" ), (f"【RQ3】 L56 急傾斜地崩壊防止施設との比較 — 砂防三施設族の地理逆転", sec6_intro + "

実装コード (市町別 outer merge + symlog 軸 + シェア比較)

" + sec6_code + sec6_fig7_text + figure("assets/L57_fig7_l56_vs_l57_overlay.png", "図 7 (RQ3): L56 (青点) × L57 (紫面) 重ね合わせ全県マップ") + sec6_fig7_read + sec6_fig8_text + figure("assets/L57_fig8_l56_vs_l57_ranking.png", "図 8 (RQ3): L56 vs L57 市町別 並列比較 (件数 symlog + シェア %)") + sec6_fig8_read + "

表: L56 (急傾斜) vs L57 (地すべり) 市町別比較 Top 15

" + df_to_html(T_compare) + "

この表から読み取れること: 各市町の L56 件数・シェア・L57 件数・シェア・順位差。" "順位差 = L56 順位 - L57 順位。正の値が大きい市町ほどL56 では下位だが L57 では上位 (= 中山間特化)、" "負の値が大きい市町ほどL56 では上位だが L57 では下位 (= 沿岸都市特化)。" f"呉市は順位差 {-(int(cmp_56_57[cmp_56_57['市町名']=='呉市']['順位差'].iloc[0])) if (cmp_56_57['市町名']=='呉市').any() else '?'} 程度の差で典型的な L56 偏重 (沿岸都市)、" f"{top1_name}は順位差 +{abs(int(cmp_56_57[cmp_56_57['市町名']==top1_name]['順位差'].iloc[0])) if (cmp_56_57['市町名']==top1_name).any() else '?'} 程度の差で L57 偏重 (中山間)。" "制度ごとに別の市町ランキングが存在する事実を 1 表で実証。

" ), ("仮説検証総合", sec7), ("発展課題", sec8), ] html = render_lesson( num=57, title="地すべり防止施設 単独 3 研究例分析 — 地すべり等防止法 1958 の polygon 台帳から砂防三施設族の地理を読む", tags=["L57", "地すべり等防止法", "地すべり防止施設", "RQ×3", "Format B", "防災施設", "中山間偏重", "庄原市", "L46連携", "L56比較", "砂防三施設族", "polygon 解析", "地理逆転"], time="35 分", level="中級", data_label=f"DoBoX dataset 61 (Shapefile 1 セット, {n_total} polygon) + L46/L56 連携", sections=sections, script_filename="L57_landslide_facility.py", ) OUT_HTML = LESSONS / "L57_landslide_facility.html" OUT_HTML.write_text(html, encoding="utf-8") print(f" HTML written: {OUT_HTML}") print(f" ({time.time()-t10:.1f}s)", flush=True) print(f"\n=== L57 DONE in {time.time()-t_all:.1f}s ===")