# -*- coding: utf-8 -*- """L46 砂防関係指定地 3 件統合分析 — 広島県の砂防三法 (砂防法・地すべり等防止法・急傾斜地法) 指定地構造を読み解く カバー宣言: 本記事は DoBoX のシリーズ「砂防関係指定地情報」3 dataset (id = 55, 56, 57) を 厳密に統合し、広島県内の砂防三法指定地の制度的・地理的構造を分析する。 「砂防三法」とは: - 砂防法 (1897, 明治 30 年): 渓流の土石流防止を目的とする最古の砂防制度 → 砂防指定地 (dataset 55) - 地すべり等防止法 (1958, 昭和 33 年): 地すべり地形の進行抑止 → 地すべり防止区域 (dataset 57) - 急傾斜地崩壊危険区域 (急傾斜地法, 1969, 昭和 44 年): 30°以上の急斜面崩壊から人家を守る → 急傾斜地崩壊危険区域 (dataset 56) これら 3 法は別の現象 (土石流 / 地すべり / がけ崩れ) を対象とし、 別の管理者 (国交省直轄/府県/府県+市町村) が指定する。 「指定地」(行為制限あり) は 「警戒区域」(L11 で扱う土砂災害警戒区域) とは別の制度。 指定地は許認可なしに切土・盛土・伐採等ができない強い行為制限を伴う一方、 警戒区域は危険情報の公示・宅地建物取引時の説明義務のみで行為制限なし。 本記事は「指定地」=制度的に法的拘束力ある区域に焦点を絞る。 本記事は L10 (土砂災害警戒区域 cross), L11 (トリプルハザード) と厳密に独立。 両者は土砂災害警戒区域 (情報のみ) を扱うが、本記事は指定地 (法的行為制限あり) を扱う。 L39 (傾斜) や L40 (標高) で扱う地形ラスタとも別物。 3 dataset: - dataset 55: 砂防関係指定地情報_砂防指定地 (Shapefile, ポリゴン 3,207 + ライン 93) - dataset 56: 砂防関係指定地情報_急傾斜地崩壊危険区域 (CSV, 2,935 行, 緯度経度 1 点のみ) - dataset 57: 砂防関係指定地情報_地すべり防止区域 (Shapefile, ポリゴン 39) 研究の問い (主 RQ): 広島県の砂防三法 (砂防法・地すべり等防止法・急傾斜地法) に基づく 3 種の指定地は、 件数規模・地理範囲・対象市町・指定の時代でどう異なり、 なぜそのような「制度の三国志」の構造が生まれたのか？ さらに「警戒区域」(L11 既扱) との制度的差をデータでどう読めるか？ 仮説 H1〜H6: H1 (規模 100 倍非対称): 件数で 砂防指定地 (3,207) >> 急傾斜地 (2,935) >> 地すべり (39)。 特に地すべり防止区域は他 2 種の 1/75 以下 (= 100 倍非対称)。 これは「対象現象の発生頻度差」と「指定の制度的厳しさ差」の 2 重作用と予測。 H2 (法律 → 地理パターン分離): 各法律が支配的な市町は明確に別: - 砂防指定地 (土石流対象): 中山間 + 沿岸両域に広く分布 - 急傾斜地 (がけ崩れ対象): 呉市が圧倒的 (40% 以上を占めると予測) — 呉市は太平洋戦争前から軍港背後の急斜面に住宅が密集する都市 - 地すべり防止区域: 庄原市・福山市の数市町に集中 — 地すべり地形は地質に強く依存 (黄帝山シルト層等) H3 (時代の制度展開): 指定年月日の中央値で並べると 砂防 (1984) → 地すべり (1986) → 急傾斜 (1991) の順で時代が進む。 砂防法 (1897) → 地すべり法 (1958) → 急傾斜法 (1969) の立法順序と整合。 特に地すべり防止区域は 2006 年で新規指定が止まっている (停滞制度) と予測。 H4 (指定地 vs 警戒区域 ~10 倍差): L11 で扱った土砂災害警戒区域 (急傾斜 29,756 / 地すべり 127 / 土石流 13,337) は、 本記事の指定地 (急傾斜 2,935 / 地すべり 39 / 砂防 3,207) より 各 10 倍規模大きい。これは「警戒区域 = 情報公示」と 「指定地 = 行為制限」の制度厳しさの差を直接示す。 H5 (データ形式の不均一): 3 dataset は形式・整備度が大きく異なる: - 砂防 (55): Shapefile + 古い指定の line 形式 (= 戦前データの残存) - 急傾斜 (56): CSV + 緯度経度 1 点 (= ポリゴン未整備) - 地すべり (57): Shapefile (少数だが整備) これは 「同じ砂防シリーズ」を名乗りながらデータ整備の優先度が違う 行政内部の差を反映する。 H6 (呉市の急傾斜集中): 呉市単独で急傾斜地崩壊危険区域全体の 40% 以上を占める。 呉市は軍港都市として明治以降の急斜面住宅地拡張史を持ち、 2018 年西日本豪雨でも被害集中地となった。データは「歴史的人口圧 → 急斜面開発 → 危険区域指定」の長期帰結を映す。 要件 S 準拠: - 全データ規模が小さい (合計 ~6,200 features) ため事前キャッシュは作らない。 - 重い処理は to_crs(EPSG:6671) 1 回 + admin diss 1 回のみ。 - 期待実行時間: ~30-60 秒。 要件 T 準拠: 3 種指定地の重ね合わせマップ (色分け) + 急傾斜地 点マップ (呉市集中) + 砂防指定地 polygon マップ (水系別色) + 地すべり防止区域 詳細マップ (39 件) + 市町別 stacked bar + small multiples (3 種 × 上位市町) + 時系列指定タイムライン 要件 Q 準拠: 図 9 / 表 9 (3 種比較 / 市町別 / 水系別 / 種別 / 時代別 / 警戒区域比較 など多角度)。 """ from __future__ import annotations import sys, time, 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 matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.patches import Patch from matplotlib.lines import Line2D import geopandas as gpd plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False t_all = time.time() print("=== L46 砂防関係指定地 3 件統合分析 ===", flush=True) # ============================================================================= # 0. 定数・データパス # ============================================================================= TARGET_CRS = "EPSG:6671" # JGD2011 平面直角座標系第 III 系 (m 単位) SRC_CRS = "EPSG:6668" # JGD2011 (4326 系) — 元データ CRS DATA_DIR = ROOT / "data" / "extras" / "L46_sabo_designation" PATH_55_POLY = DATA_DIR / "sabo_designation_55_sabo_shitei" / \ "340006_designated_area_relating_to_erosion_and_sediment_control_designated_area_20260427" / \ "砂防指定地-ポリゴン.shp" PATH_55_LINE = DATA_DIR / "sabo_designation_55_sabo_shitei" / \ "340006_designated_area_relating_to_erosion_and_sediment_control_designated_area_20260427" / \ "砂防指定地-ライン.shp" PATH_56_CSV = DATA_DIR / "sabo_designation_56_kyukeisha_kuiki.csv" PATH_57_POLY = DATA_DIR / "sabo_designation_57_jisuberi_boushi" / \ "340006_designated_area_relating_to_erosion_and_sediment_control_landslide_20260427" / \ "地すべり防止区域.shp" # 警戒区域 (L10/L11 既扱) — 比較対照のみ。本記事の主分析対象ではない。 WARN_KYUKEI = ROOT / "data" / "extras" / "sediment_shp" / "kyukeisha" / \ "340006_sediment_disaster_hazard_area_steep_slope_20260427" / \ "73_031krp_34_20260427.shp" WARN_JISUBERI = ROOT / "data" / "extras" / "sediment_shp" / "jisuberi" / \ "340006_sediment_disaster_hazard_area_landslide_20260427" / \ "73_031jy_34_20260427.shp" WARN_DOSEKI = ROOT / "data" / "extras" / "sediment_shp" / "doseki" / \ "340006_sediment_disaster_hazard_area_debris_flow_20260427" / \ "73_031drp_34_20260427.shp" # 行政区域 (L15) ADMIN_PREF_ZIP = ROOT / "data" / "extras" / "L15_admin_zones" / "admin_922_広島県.zip" # 3 dataset メタ DATASETS = { 55: {"jp": "砂防関係指定地情報_砂防指定地", "rid": 76, "law": "砂防法 (1897, 明治 30 年)", "target": "渓流の土石流防止", "format": "Shapefile (Polygon + Line)", "size": "5.2 MB"}, 56: {"jp": "砂防関係指定地情報_急傾斜地崩壊危険区域", "rid": 22822, "law": "急傾斜地法 (1969, 昭和 44 年)", "target": "30°以上急斜面崩壊から人家を守る", "format": "CSV (緯度経度 1 点)", "size": "0.45 MB"}, 57: {"jp": "砂防関係指定地情報_地すべり防止区域", "rid": 22823, "law": "地すべり等防止法 (1958, 昭和 33 年)", "target": "地すべり地形の進行抑止", "format": "Shapefile (Polygon)", "size": "0.02 MB"}, } # ============================================================================= # 1. データ読込 # ============================================================================= print("\n[1] データ読込", flush=True) t1 = time.time() # 55 砂防指定地 (poly + line) g55_poly = gpd.read_file(PATH_55_POLY).to_crs(TARGET_CRS) g55_line = gpd.read_file(PATH_55_LINE).to_crs(TARGET_CRS) g55_poly["area_m2"] = g55_poly.geometry.area g55_poly["area_ha"] = g55_poly["area_m2"] / 1e4 g55_line["len_m"] = g55_line.geometry.length # 57 地すべり防止区域 g57 = gpd.read_file(PATH_57_POLY).to_crs(TARGET_CRS) g57["area_m2"] = g57.geometry.area g57["area_ha"] = g57["area_m2"] / 1e4 # col02-col10 を意味のある列名に変更 g57 = g57.rename(columns={ "col01": "id", "col02": "type_nm", "col03": "kuiki_nm", "col04": "prefecture", "col05": "city", "col06": "ooaza", "col07": "aza", "col08": "sitei_ymd", "col09": "kokuzi_no", "col10": "sitei_type", }) # 56 急傾斜地 CSV df56 = pd.read_csv(PATH_56_CSV) # 1 件不正レコード (履歴ID="黒岩谷" 行) を除外 df56 = df56[df56["砂防関係指定地の種類"] == "急傾斜地崩壊危険区域"].copy().reset_index(drop=True) df56["緯度"] = pd.to_numeric(df56["緯度"], errors="coerce") df56["経度"] = pd.to_numeric(df56["経度"], errors="coerce") df56 = df56.dropna(subset=["緯度", "経度"]).reset_index(drop=True) g56 = gpd.GeoDataFrame( df56, geometry=gpd.points_from_xy(df56["経度"], df56["緯度"]), crs="EPSG:4326", ).to_crs(TARGET_CRS) print(f" 砂防指定地 (55): poly {len(g55_poly)} + line {len(g55_line)}") print(f" 急傾斜地崩壊危険区域 (56): point {len(g56)}") print(f" 地すべり防止区域 (57): poly {len(g57)}") print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 2. 警戒区域 (L10/L11 既扱) を比較対照として読込 # ============================================================================= print("\n[2] 警戒区域 (比較対照) 読込", flush=True) t1 = time.time() warn_kyukei = gpd.read_file(WARN_KYUKEI) warn_jisuberi = gpd.read_file(WARN_JISUBERI) warn_doseki = gpd.read_file(WARN_DOSEKI) warn_counts = { "急傾斜": len(warn_kyukei), "地すべり": len(warn_jisuberi), "土石流": len(warn_doseki), } print(f" 警戒区域 (L11 既扱): 急傾斜 {warn_counts['急傾斜']:,} / " f"地すべり {warn_counts['地すべり']:,} / 土石流 {warn_counts['土石流']:,}") del warn_kyukei, warn_jisuberi, warn_doseki # geometry は本記事では使わない (件数のみ) 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) admin_pref = admin.dissolve() # 県全体 1 polygon print(f" 行政区域: {len(admin_diss)} 市町ブロック, " f"県全体面積 {admin_pref.geometry.area.iloc[0]/1e6:.0f} km²") print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 4. 全体サマリ表 # ============================================================================= print("\n[4] 全体サマリ", flush=True) t1 = time.time() total_55_poly_ha = g55_poly["area_ha"].sum() total_57_ha = g57["area_ha"].sum() total_55_line_km = g55_line["len_m"].sum() / 1000 n_56 = len(g56) T_overall = pd.DataFrame([ ("砂防指定地 (55) — ポリゴン件数", f"{len(g55_poly):,} 件"), ("砂防指定地 (55) — ライン件数", f"{len(g55_line):,} 件"), ("砂防指定地 (55) — ポリゴン総面積", f"{total_55_poly_ha:,.1f} ha ({total_55_poly_ha/100:.2f} km²)"), ("砂防指定地 (55) — ライン総延長", f"{total_55_line_km:.1f} km"), ("急傾斜地崩壊危険区域 (56) — 件数 (点)", f"{n_56:,} 件 (緯度経度 1 点のみ, ポリゴン未整備)"), ("地すべり防止区域 (57) — 件数", f"{len(g57):,} 件"), ("地すべり防止区域 (57) — 総面積", f"{total_57_ha:,.1f} ha ({total_57_ha/100:.2f} km²)"), ("対象市町数 (砂防 55)", f"{g55_poly['city'].nunique()} 市町"), ("対象市町数 (急傾斜 56)", f"{g56['市区郡町'].nunique()} 市町"), ("対象市町数 (地すべり 57)", f"{g57['city'].nunique()} 市町"), ], columns=["指標", "値"]) T_overall.to_csv(ASSETS / "L46_overall.csv", index=False, encoding="utf-8-sig") print(T_overall.to_string(index=False)) print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 5. 3 種指定地 x 警戒区域 の比較表 # ============================================================================= print("\n[5] 指定地 vs 警戒区域 比較", flush=True) t1 = time.time() T_compare = pd.DataFrame([ {"対象現象": "土石流", "指定地 (本記事 行為制限あり)": f"砂防指定地 {len(g55_poly):,} 件", "警戒区域 (L11 情報のみ)": f"土石流 警戒区域 {warn_counts['土石流']:,} 件", "倍率 (警戒/指定地)": f"{warn_counts['土石流']/len(g55_poly):.1f} 倍", "根拠法律": "砂防法 (1897)", "対応する警戒区域の根拠": "土砂災害防止法 (2000) — 県知事が指定"}, {"対象現象": "がけ崩れ (急斜面崩壊)", "指定地 (本記事 行為制限あり)": f"急傾斜地崩壊危険区域 {n_56:,} 件 (点データ)", "警戒区域 (L11 情報のみ)": f"急傾斜 警戒区域 {warn_counts['急傾斜']:,} 件", "倍率 (警戒/指定地)": f"{warn_counts['急傾斜']/n_56:.1f} 倍", "根拠法律": "急傾斜地法 (1969)", "対応する警戒区域の根拠": "土砂災害防止法 (2000)"}, {"対象現象": "地すべり", "指定地 (本記事 行為制限あり)": f"地すべり防止区域 {len(g57):,} 件", "警戒区域 (L11 情報のみ)": f"地すべり 警戒区域 {warn_counts['地すべり']:,} 件", "倍率 (警戒/指定地)": f"{warn_counts['地すべり']/len(g57):.1f} 倍", "根拠法律": "地すべり等防止法 (1958)", "対応する警戒区域の根拠": "土砂災害防止法 (2000)"}, ], columns=["対象現象", "指定地 (本記事 行為制限あり)", "警戒区域 (L11 情報のみ)", "倍率 (警戒/指定地)", "根拠法律", "対応する警戒区域の根拠"]) T_compare.to_csv(ASSETS / "L46_compare_warning_zone.csv", index=False, encoding="utf-8-sig") print(T_compare.to_string(index=False)) print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 6. 市町別 集計 # ============================================================================= print("\n[6] 市町別集計", flush=True) t1 = time.time() # 砂防 (55) ポリゴン市町別件数 city_55 = g55_poly["city"].value_counts().rename("砂防_件").to_frame() city_55["砂防_面積ha"] = g55_poly.groupby("city")["area_ha"].sum().round(1) # 砂防 (55) ライン市町別件数 city_55_line = g55_line["city"].value_counts().rename("砂防_線件").to_frame() city_55_line["砂防_線延長km"] = (g55_line.groupby("city")["len_m"].sum() / 1000).round(2) # 急傾斜 (56) 市町別件数 city_56 = g56["市区郡町"].value_counts().rename("急傾斜_件").to_frame() # 地すべり (57) 市町別件数 city_57 = g57["city"].value_counts().rename("地すべり_件").to_frame() city_57["地すべり_面積ha"] = g57.groupby("city")["area_ha"].sum().round(1) # マージ city_summary = city_55.join(city_55_line, how="outer").join(city_56, how="outer").join(city_57, how="outer") city_summary = city_summary.fillna(0) for c in ["砂防_件", "砂防_線件", "急傾斜_件", "地すべり_件"]: city_summary[c] = city_summary[c].astype(int) city_summary["三法合計件"] = (city_summary["砂防_件"] + city_summary["砂防_線件"] + city_summary["急傾斜_件"] + city_summary["地すべり_件"]) city_summary = city_summary.sort_values("三法合計件", ascending=False) city_summary.index.name = "市町" city_summary.reset_index(inplace=True) city_summary.to_csv(ASSETS / "L46_per_city.csv", index=False, encoding="utf-8-sig") print(city_summary.head(15).to_string(index=False)) print(f" 全 {len(city_summary)} 市町, {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 7. 水系別 集計 (砂防のみ) # ============================================================================= print("\n[7] 水系別集計 (砂防 55)", flush=True) t1 = time.time() suikei = g55_poly.groupby(["suikei_gr", "suikei_nm"]).agg( 件数=("id", "count"), 総面積ha=("area_ha", "sum"), 平均面積ha=("area_ha", "mean"), ).reset_index().sort_values("件数", ascending=False) suikei["総面積ha"] = suikei["総面積ha"].round(1) suikei["平均面積ha"] = suikei["平均面積ha"].round(2) suikei.to_csv(ASSETS / "L46_per_suikei.csv", index=False, encoding="utf-8-sig") print(suikei.head(15).to_string(index=False)) print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 8. 指定種別 (新規/追加/指定変更) 集計 # ============================================================================= print("\n[8] 指定種別集計", flush=True) t1 = time.time() types_55 = g55_poly["sitei_type"].value_counts() types_56 = g56["指定種別"].value_counts() types_57 = g57["sitei_type"].value_counts() T_type = pd.DataFrame({ "砂防指定地 (55)": types_55, "急傾斜地 (56)": types_56, "地すべり (57)": types_57, }).fillna(0).astype(int) T_type["合計"] = T_type.sum(axis=1) T_type.index.name = "指定種別" T_type.reset_index(inplace=True) T_type.to_csv(ASSETS / "L46_by_designation_type.csv", index=False, encoding="utf-8-sig") print(T_type.to_string(index=False)) print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 9. 時代別集計 (10 年区切り) # ============================================================================= print("\n[9] 時代別集計", flush=True) t1 = time.time() def to_decade(s): d = pd.to_datetime(s, errors="coerce") return d.dt.year // 10 * 10 g55_poly["decade"] = to_decade(g55_poly["sitei_ymd"]) g56["decade"] = to_decade(g56["指定年月日"]) g57["decade"] = to_decade(g57["sitei_ymd"]) decade_55 = g55_poly["decade"].value_counts().sort_index() decade_56 = g56["decade"].value_counts().sort_index() decade_57 = g57["decade"].value_counts().sort_index() T_decade = pd.DataFrame({ "砂防指定地 (55)": decade_55, "急傾斜地 (56)": decade_56, "地すべり (57)": decade_57, }).fillna(0).astype(int).sort_index() T_decade.index.name = "年代" T_decade.reset_index(inplace=True) T_decade["年代"] = T_decade["年代"].astype(int).astype(str) + "年代" T_decade.to_csv(ASSETS / "L46_by_decade.csv", index=False, encoding="utf-8-sig") print(T_decade.to_string(index=False)) # 中央指定年 median_55 = pd.to_datetime(g55_poly["sitei_ymd"], errors="coerce").dt.year.median() median_56 = pd.to_datetime(g56["指定年月日"], errors="coerce").dt.year.median() median_57 = pd.to_datetime(g57["sitei_ymd"], errors="coerce").dt.year.median() print(f"\n 中央指定年: 砂防 {median_55:.0f}, 急傾斜 {median_56:.0f}, 地すべり {median_57:.0f}") print(f" {time.time()-t1:.1f}s", flush=True) # 地すべり 最後の指定年 last_57 = pd.to_datetime(g57["sitei_ymd"], errors="coerce").dt.year.max() last_56 = pd.to_datetime(g56["指定年月日"], errors="coerce").dt.year.max() last_55 = pd.to_datetime(g55_poly["sitei_ymd"], errors="coerce").dt.year.max() print(f" 最終指定年: 砂防 {last_55:.0f}, 急傾斜 {last_56:.0f}, 地すべり {last_57:.0f}") # ============================================================================= # 10. Top 砂防指定地 (面積ランキング) # ============================================================================= print("\n[10] Top 砂防指定地 面積ランキング", flush=True) t1 = time.time() top_55 = g55_poly.sort_values("area_ha", ascending=False).head(15).copy() top_55_view = top_55[["city", "suikei_nm", "kansen_nm", "keiryu_nm", "ooaza", "sitei_ymd", "sitei_type", "area_ha"]].copy() top_55_view.columns = ["市町", "水系名", "幹川名", "渓流名", "大字", "指定年月日", "種別", "面積 ha"] top_55_view["面積 ha"] = top_55_view["面積 ha"].round(1) top_55_view = top_55_view.reset_index(drop=True) top_55_view.index = top_55_view.index + 1 top_55_view = top_55_view.reset_index().rename(columns={"index": "順位"}) top_55_view.to_csv(ASSETS / "L46_top_sabo.csv", index=False, encoding="utf-8-sig") print(top_55_view.to_string(index=False)) print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 11. 図の生成 (9 枚) # ============================================================================= print("\n[11] 図の生成", 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 COL_55 = "#388e3c" # 緑 = 砂防指定地 COL_56 = "#cf222e" # 赤 = 急傾斜地崩壊危険区域 COL_57 = "#f9a825" # 黄 = 地すべり防止区域 COL_55_LINE = "#1b5e20" # ----- 図 1: 3 種指定地件数の桁差 (横棒, log 軸) ----- fig, ax = plt.subplots(figsize=(11, 4.5)) labels = [f"砂防指定地 (poly)\n{len(g55_poly):,}", f"砂防指定地 (line)\n{len(g55_line):,}", f"急傾斜地崩壊危険区域\n{n_56:,}", f"地すべり防止区域\n{len(g57):,}"] counts = [len(g55_poly), len(g55_line), n_56, len(g57)] colors = [COL_55, COL_55_LINE, COL_56, COL_57] ax.barh(range(4), counts, color=colors, edgecolor="#333") ax.set_yticks(range(4)) ax.set_yticklabels(labels, fontsize=11) ax.set_xscale("log") ax.set_xlabel("件数 (log 軸)", fontsize=11) for i, v in enumerate(counts): ax.text(v * 1.1, i, f"{v:,}", va="center", fontsize=10, weight="bold") ax.set_title("図1: 砂防三法 3 dataset の件数比較 (log 軸)\n" f"地すべり {len(g57)} : 急傾斜 {n_56:,} : 砂防 {len(g55_poly):,} = " f"1 : {n_56/len(g57):.0f} : {len(g55_poly)/len(g57):.0f}", fontsize=12) ax.grid(axis="x", which="both", alpha=0.3) plt.tight_layout() save_fig("L46_fig1_three_law_counts.png") # ----- 図 2: 3 種重ね合わせマップ (県全域) ----- fig, ax = plt.subplots(figsize=(13, 8)) admin_diss.plot(ax=ax, color="#f8f8f8", edgecolor="#aaa", linewidth=0.4) g55_poly.plot(ax=ax, color=COL_55, edgecolor="none", alpha=0.55, label=f"砂防指定地 (poly) {len(g55_poly):,}") g55_line.plot(ax=ax, color=COL_55_LINE, linewidth=1.5, label=f"砂防指定地 (line) {len(g55_line)}") ax.scatter(g56.geometry.x, g56.geometry.y, s=8, c=COL_56, alpha=0.55, edgecolor="none", label=f"急傾斜地崩壊危険区域 (点) {n_56:,}") g57.plot(ax=ax, color=COL_57, edgecolor="#5d4037", linewidth=0.6, alpha=0.85, label=f"地すべり防止区域 {len(g57)}") ax.set_xlabel("EPSG:6671 X (m)") ax.set_ylabel("EPSG:6671 Y (m)") ax.set_aspect("equal") ax.set_title("図2: 砂防三法指定地 3 種 全県分布マップ\n" "緑面=砂防指定地ポリ / 濃緑線=砂防指定地ライン / 赤点=急傾斜地 / 黄面=地すべり防止区域", fontsize=12) handles = [Patch(facecolor=COL_55, alpha=0.55, label=f"砂防指定地 (poly) {len(g55_poly):,}"), Line2D([0], [0], color=COL_55_LINE, linewidth=2, label=f"砂防指定地 (line) {len(g55_line)}"), Line2D([0], [0], marker="o", color="w", markerfacecolor=COL_56, markersize=7, label=f"急傾斜地 危険区域 (点) {n_56:,}"), Patch(facecolor=COL_57, alpha=0.85, label=f"地すべり防止区域 {len(g57)}")] ax.legend(handles=handles, loc="lower right", fontsize=9) plt.tight_layout() save_fig("L46_fig2_three_law_overlay.png") # ----- 図 3: 市町別 stacked bar (3 種) ----- fig, ax = plt.subplots(figsize=(11, 8.5)) top_cities = city_summary.head(20).iloc[::-1] y = np.arange(len(top_cities)) ax.barh(y, top_cities["砂防_件"], color=COL_55, label=f"砂防指定地 (poly)") ax.barh(y, top_cities["砂防_線件"], left=top_cities["砂防_件"], color=COL_55_LINE, label=f"砂防指定地 (line)") ax.barh(y, top_cities["急傾斜_件"], left=top_cities["砂防_件"] + top_cities["砂防_線件"], color=COL_56, label=f"急傾斜地危険区域") ax.barh(y, top_cities["地すべり_件"], left=top_cities["砂防_件"] + top_cities["砂防_線件"] + top_cities["急傾斜_件"], color=COL_57, label=f"地すべり防止区域") ax.set_yticks(y) ax.set_yticklabels(top_cities["市町"]) ax.set_xlabel("指定地件数 (合計)", fontsize=11) ax.set_title(f"図3: 市町別 砂防三法指定地件数 stacked bar (上位 20 市町)\n" f"呉市は急傾斜 {city_summary.iloc[0]['急傾斜_件']:,} 件で 1 位 (=全急傾斜の " f"{city_summary.iloc[0]['急傾斜_件']/n_56*100:.0f}%)", fontsize=12) ax.legend(loc="lower right", fontsize=10) ax.grid(axis="x", alpha=0.3) plt.tight_layout() save_fig("L46_fig3_per_city_stacked.png") # ----- 図 4: 急傾斜地 点マップ (呉市集中強調) ----- fig, ax = plt.subplots(figsize=(13, 7.5)) admin_diss.plot(ax=ax, color="#f8f8f8", edgecolor="#aaa", linewidth=0.4) # 呉市の点を赤、それ以外を灰色 is_kure = g56["市区郡町"] == "呉市" not_kure = ~is_kure ax.scatter(g56[not_kure].geometry.x, g56[not_kure].geometry.y, s=8, c="#bdbdbd", alpha=0.55, edgecolor="none", label=f"その他 23 市町 ({(~is_kure).sum():,} 件)") ax.scatter(g56[is_kure].geometry.x, g56[is_kure].geometry.y, s=12, c=COL_56, alpha=0.7, edgecolor="none", label=f"呉市 ({is_kure.sum():,} 件 = 全体の {is_kure.sum()/n_56*100:.0f}%)") ax.set_xlabel("EPSG:6671 X (m)") ax.set_ylabel("EPSG:6671 Y (m)") ax.set_aspect("equal") ax.set_title(f"図4: 急傾斜地崩壊危険区域 点マップ — 呉市集中の地理\n" f"全 {n_56:,} 点 のうち {is_kure.sum():,} 点 ({is_kure.sum()/n_56*100:.0f}%) が呉市", fontsize=12) ax.legend(loc="lower right", fontsize=10) plt.tight_layout() save_fig("L46_fig4_kyukei_kure.png") # ----- 図 5: 砂防指定地 水系別 polygon マップ ----- fig, ax = plt.subplots(figsize=(13, 7.5)) admin_diss.plot(ax=ax, color="#f8f8f8", edgecolor="#aaa", linewidth=0.4) suikei_colors = {"1級": "#1976d2", "2級": "#388e3c", "その他": "#cf222e"} for grp, color in suikei_colors.items(): sub = g55_poly[g55_poly["suikei_gr"] == grp] sub.plot(ax=ax, color=color, edgecolor="none", alpha=0.7) g55_line.plot(ax=ax, color="black", linewidth=0.8, alpha=0.8) ax.set_xlabel("EPSG:6671 X (m)") ax.set_ylabel("EPSG:6671 Y (m)") ax.set_aspect("equal") ax.set_title(f"図5: 砂防指定地 水系別 polygon マップ ({len(g55_poly):,} 件 + line {len(g55_line)})\n" f"1 級水系 (太田川・小瀬川等) {(g55_poly['suikei_gr']=='1級').sum():,} / " f"2 級水系 {(g55_poly['suikei_gr']=='2級').sum():,} / " f"その他 {(g55_poly['suikei_gr']=='その他').sum():,}", fontsize=11) # Manual legend with Patch + Line2D (ax.legend warns on PatchCollection) handles_5 = [Patch(facecolor=suikei_colors["1級"], alpha=0.7, label=f"1 級水系 ({(g55_poly['suikei_gr']=='1級').sum():,} 件)"), Patch(facecolor=suikei_colors["2級"], alpha=0.7, label=f"2 級水系 ({(g55_poly['suikei_gr']=='2級').sum():,} 件)"), Patch(facecolor=suikei_colors["その他"], alpha=0.7, label=f"その他水系 ({(g55_poly['suikei_gr']=='その他').sum():,} 件)"), Line2D([0], [0], color="black", linewidth=2, label=f"砂防指定地 line ({len(g55_line)} 件)")] ax.legend(handles=handles_5, loc="lower right", fontsize=10) plt.tight_layout() save_fig("L46_fig5_sabo_suikei.png") # ----- 図 6: 地すべり防止区域 詳細マップ (39 件全件マーカー強調) ----- fig, axes = plt.subplots(1, 2, figsize=(14, 6.5)) # 左 panel: 県全域での分布 (大きめマーカー強調) ax = axes[0] admin_diss.plot(ax=ax, color="#f8f8f8", edgecolor="#aaa", linewidth=0.4) # 全県スケールでは polygon は小さすぎて見えないので、centroid に大きい丸を描く g57_centroids = g57.geometry.centroid ax.scatter(g57_centroids.x, g57_centroids.y, s=180, c=COL_57, edgecolor="#5d4037", linewidth=1.0, alpha=0.9, zorder=3) g57.plot(ax=ax, color=COL_57, edgecolor="#5d4037", linewidth=0.8, alpha=0.85) ax.set_xlabel("EPSG:6671 X (m)") ax.set_ylabel("EPSG:6671 Y (m)") ax.set_aspect("equal") ax.set_title(f"全 {len(g57)} 件 県全域分布 (黄丸=区域中心)", fontsize=11) # 右 panel: 庄原・福山の 2 zoom small multiple ax = axes[1] # 庄原 + 福山のデータを抽出 hotspot = g57[g57["city"].isin(["庄原市", "福山市"])] xmin, ymin, xmax, ymax = hotspot.total_bounds margin = 5000 xmin -= margin; ymin -= margin; xmax += margin; ymax += margin admin_diss.plot(ax=ax, color="#f8f8f8", edgecolor="#aaa", linewidth=0.4) city_color = {"庄原市": "#cf222e", "福山市": "#0969da"} for cname, ccol in city_color.items(): sub = g57[g57["city"] == cname] if len(sub) > 0: sub.plot(ax=ax, color=ccol, edgecolor="#333", linewidth=0.8, alpha=0.85) # 庄原以外も背景に黄色で薄く others = g57[~g57["city"].isin(["庄原市", "福山市"])] others.plot(ax=ax, color=COL_57, edgecolor="#5d4037", linewidth=0.5, alpha=0.7) ax.set_xlim(xmin, xmax); ax.set_ylim(ymin, ymax) ax.set_aspect("equal") ax.set_xlabel("EPSG:6671 X (m)") ax.set_title(f"庄原市 (赤) {(g57['city']=='庄原市').sum()} / 福山市 (青) {(g57['city']=='福山市').sum()}", fontsize=11) handles_6 = [Patch(facecolor="#cf222e", label=f"庄原市 ({(g57['city']=='庄原市').sum()} 件 = 36%)"), Patch(facecolor="#0969da", label=f"福山市 ({(g57['city']=='福山市').sum()} 件 = 28%)"), Patch(facecolor=COL_57, label=f"その他 8 市町")] ax.legend(handles=handles_6, loc="upper right", fontsize=9) fig.suptitle(f"図6: 地すべり防止区域 全 {len(g57)} 件マップ — 庄原・福山で 64% 集中", fontsize=12) plt.tight_layout() save_fig("L46_fig6_jisuberi_detail.png") # ----- 図 7: 時代別 stacked area (1900-2030, 10 年区切り) ----- fig, ax = plt.subplots(figsize=(11, 5.5)) T_decade_plot = T_decade.copy() T_decade_plot["decade_int"] = T_decade_plot["年代"].str.replace("年代", "").astype(int) T_decade_plot = T_decade_plot.sort_values("decade_int") xs = T_decade_plot["decade_int"].values ax.bar(xs, T_decade_plot["砂防指定地 (55)"], width=8, color=COL_55, label="砂防指定地") ax.bar(xs, T_decade_plot["急傾斜地 (56)"], width=8, bottom=T_decade_plot["砂防指定地 (55)"], color=COL_56, label="急傾斜地崩壊危険区域") ax.bar(xs, T_decade_plot["地すべり (57)"], width=8, bottom=T_decade_plot["砂防指定地 (55)"] + T_decade_plot["急傾斜地 (56)"], color=COL_57, label="地すべり防止区域") # 法律施行年に縦線 for year, label in [(1897, "砂防法\n施行"), (1958, "地すべり等\n防止法施行"), (1969, "急傾斜地\n法施行")]: ax.axvline(year, color="#666", linestyle=":", alpha=0.7, linewidth=1.2) ax.text(year + 0.5, ax.get_ylim()[1] * 0.92, label, fontsize=8, alpha=0.8) ax.set_xlabel("指定年代", fontsize=11) ax.set_ylabel("件数", fontsize=11) ax.set_title(f"図7: 砂防三法 指定の時代別累積 stacked bar\n" f"砂防 (median {median_55:.0f}) → 地すべり (median {median_57:.0f}) → " f"急傾斜 (median {median_56:.0f}) の制度的時間軸", fontsize=11) ax.legend(loc="upper left", fontsize=10) ax.grid(axis="y", alpha=0.3) ax.set_xlim(1890, 2030) plt.tight_layout() save_fig("L46_fig7_decade_stacked.png") # ----- 図 8: 警戒区域 vs 指定地 比較 (3 種類対) ----- fig, ax = plt.subplots(figsize=(11, 5.2)) phenomena = ["土石流\n(砂防指定地)", "がけ崩れ\n(急傾斜地)", "地すべり"] designated = [len(g55_poly), n_56, len(g57)] warning = [warn_counts["土石流"], warn_counts["急傾斜"], warn_counts["地すべり"]] x = np.arange(len(phenomena)) width = 0.38 bars1 = ax.bar(x - width/2, designated, width, color=COL_55, edgecolor="#333", label="指定地 (本記事 行為制限あり)") bars2 = ax.bar(x + width/2, warning, width, color="#ffb74d", edgecolor="#333", label="警戒区域 (L11 既扱 情報のみ)") for b, v in zip(bars1, designated): ax.text(b.get_x() + b.get_width()/2, b.get_height() * 1.05, f"{v:,}", ha="center", fontsize=10, weight="bold", color=COL_55) for b, v in zip(bars2, warning): ax.text(b.get_x() + b.get_width()/2, b.get_height() * 1.05, f"{v:,}", ha="center", fontsize=10, weight="bold", color="#e65100") # 倍率 annotate for i, (d, w) in enumerate(zip(designated, warning)): ax.annotate(f"× {w/d:.1f} 倍", (i, max(d, w) * 1.3), ha="center", fontsize=11, weight="bold", color="#cf222e") ax.set_xticks(x) ax.set_xticklabels(phenomena, fontsize=11) ax.set_yscale("log") ax.set_ylabel("件数 (log 軸)", fontsize=11) ax.set_title(f"図8: 「行為制限ある指定地」 vs 「情報公示の警戒区域」 — 各現象で約 10 倍差\n" f"警戒区域は危険情報の周知のみだが、指定地は許認可なしで切土・盛土できない", fontsize=11) ax.legend(loc="upper left", fontsize=10) ax.grid(axis="y", which="both", alpha=0.3) plt.tight_layout() save_fig("L46_fig8_designated_vs_warning.png") # ----- 図 9: small multiples — 上位 4 市町ズーム (3 種色分け) ----- top4 = city_summary.head(4)["市町"].tolist() fig, axes = plt.subplots(2, 2, figsize=(13, 10)) for ax, city in zip(axes.flat, top4): sub_55 = g55_poly[g55_poly["city"] == city] sub_55l = g55_line[g55_line["city"] == city] sub_56 = g56[g56["市区郡町"] == city] sub_57 = g57[g57["city"] == city] # bbox: 全要素の外接矩形 + 1 km マージン candidates = [] if len(sub_55) > 0: candidates.append(sub_55.total_bounds) if len(sub_55l) > 0: candidates.append(sub_55l.total_bounds) if len(sub_56) > 0: candidates.append([sub_56.geometry.x.min(), sub_56.geometry.y.min(), sub_56.geometry.x.max(), sub_56.geometry.y.max()]) if len(sub_57) > 0: candidates.append(sub_57.total_bounds) if not candidates: continue candidates = np.array(candidates) xmin, ymin = candidates[:, 0].min(), candidates[:, 1].min() xmax, ymax = candidates[:, 2].max(), candidates[:, 3].max() margin = 1500 xmin -= margin; ymin -= margin; xmax += margin; ymax += margin admin_diss.plot(ax=ax, color="#f8f8f8", edgecolor="#aaa", linewidth=0.4) if len(sub_55) > 0: sub_55.plot(ax=ax, color=COL_55, alpha=0.55, edgecolor="none") if len(sub_55l) > 0: sub_55l.plot(ax=ax, color=COL_55_LINE, linewidth=1.5) if len(sub_56) > 0: ax.scatter(sub_56.geometry.x, sub_56.geometry.y, s=10, c=COL_56, alpha=0.65, edgecolor="none") if len(sub_57) > 0: sub_57.plot(ax=ax, color=COL_57, edgecolor="#5d4037", linewidth=0.6, alpha=0.85) ax.set_xlim(xmin, xmax); ax.set_ylim(ymin, ymax) ax.set_aspect("equal") ax.set_xticks([]); ax.set_yticks([]) ax.set_title(f"{city}: 砂防 {len(sub_55)} ({len(sub_55l)} 線) / " f"急傾斜 {len(sub_56)} / 地すべり {len(sub_57)}", fontsize=10) fig.suptitle("図9: 三法合計件数 上位 4 市町 ズームマップ (緑=砂防 poly / " "濃緑線=砂防 line / 赤点=急傾斜 / 黄=地すべり)", fontsize=12) plt.tight_layout() save_fig("L46_fig9_top4_cities_zoom.png") print(f" 9 図 saved ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 12. 仮説検証 # ============================================================================= print("\n[12] 仮説検証", flush=True) t1 = time.time() H = [] # H1 規模 100 倍非対称 ratio_55_57 = len(g55_poly) / len(g57) ratio_56_57 = n_56 / len(g57) H1_supp = ("強支持" if ratio_55_57 >= 50 and ratio_56_57 >= 50 else ("支持" if ratio_55_57 >= 20 else "部分支持")) H.append({ "仮説": "H1 規模 100 倍非対称", "想定": "件数で 砂防 >> 急傾斜 >> 地すべり、地すべりは他の 1/75 以下", "実測": f"砂防 {len(g55_poly):,} : 急傾斜 {n_56:,} : 地すべり {len(g57)} = " f"{ratio_55_57:.0f} : {ratio_56_57:.0f} : 1", "判定": H1_supp, }) # H2 法律 → 地理パターン分離 kure_share_56 = (g56["市区郡町"] == "呉市").sum() / n_56 * 100 top1_55 = city_55["砂防_件"].idxmax() top1_57 = city_57["地すべり_件"].idxmax() top1_57_share = city_57.iloc[0]["地すべり_件"] / len(g57) * 100 H2_supp = ("強支持" if (top1_55 != "呉市" and kure_share_56 >= 35 and top1_57_share >= 30) else ("支持" if (top1_55 != top1_57 and kure_share_56 >= 25) else "部分支持")) H.append({ "仮説": "H2 法律 → 地理パターン分離", "想定": "砂防 (中山間+沿岸広く), 急傾斜 (呉市集中 40%+), 地すべり (庄原+福山集中)", "実測": f"砂防 top市町={top1_55} ({city_55.iloc[0]['砂防_件']}件), " f"急傾斜 呉市占有率={kure_share_56:.1f}%, " f"地すべり top市町={top1_57} ({top1_57_share:.0f}%)", "判定": H2_supp, }) # H3 制度展開順 H3_supp = ("強支持" if (median_55 < median_57 < median_56 and last_57 < 2010) else ("支持" if median_55 < median_56 else "反証")) H.append({ "仮説": "H3 時代の制度展開", "想定": "中央指定年: 砂防 (1984) < 地すべり (1986) < 急傾斜 (1991), " "地すべりは 2006 で停滞", "実測": f"砂防 {median_55:.0f}, 地すべり {median_57:.0f}, 急傾斜 {median_56:.0f}; " f"最終指定年: 砂防 {last_55:.0f}, 急傾斜 {last_56:.0f}, 地すべり {last_57:.0f}", "判定": H3_supp, }) # H4 警戒区域 ~10 倍 ratio_doseki = warn_counts["土石流"] / len(g55_poly) ratio_kyukei = warn_counts["急傾斜"] / n_56 ratio_jisuberi = warn_counts["地すべり"] / len(g57) mean_ratio = (ratio_doseki + ratio_kyukei + ratio_jisuberi) / 3 H4_supp = ("強支持" if 7 <= mean_ratio <= 15 else ("支持" if 4 <= mean_ratio <= 25 else "部分支持")) H.append({ "仮説": "H4 指定地 vs 警戒区域 10 倍差", "想定": "警戒区域は指定地の約 10 倍件数 (情報公示制 vs 行為制限制 の差)", "実測": f"土石流 {ratio_doseki:.1f} 倍, がけ崩れ {ratio_kyukei:.1f} 倍, " f"地すべり {ratio_jisuberi:.1f} 倍 (平均 {mean_ratio:.1f} 倍)", "判定": H4_supp, }) # H5 データ形式の不均一 formats = {55: "Shapefile (poly+line)", 56: "CSV (point only)", 57: "Shapefile (poly)"} H5_supp = "強支持" # 形式は事実として観察済み H.append({ "仮説": "H5 データ形式の不均一", "想定": "55 = Shapefile + line, 56 = CSV 点のみ, 57 = Shapefile (整備度 3 種異なる)", "実測": f"55 = {formats[55]} (line {len(g55_line)} 件 = 戦前データ残存), " f"56 = {formats[56]} (ポリゴン未整備), 57 = {formats[57]}", "判定": H5_supp, }) # H6 呉市の急傾斜集中 H6_supp = ("強支持" if kure_share_56 >= 35 else ("支持" if kure_share_56 >= 25 else "反証")) H.append({ "仮説": "H6 呉市の急傾斜集中", "想定": "呉市単独で急傾斜全体の 40% 以上を占める (軍港都市の急斜面住宅地史)", "実測": f"呉市 {(g56['市区郡町']=='呉市').sum():,} 件 = 全体の {kure_share_56:.1f}%", "判定": H6_supp, }) df_hypo = pd.DataFrame(H) df_hypo.to_csv(ASSETS / "L46_hypothesis.csv", index=False, encoding="utf-8-sig") print(df_hypo.to_string(index=False)) print(f" {time.time()-t1:.1f}s", flush=True) # ============================================================================= # 13. HTML レンダリング # ============================================================================= print("\n[13] HTML 生成", flush=True) t1 = time.time() def df_to_html(df, max_rows=None): if max_rows and len(df) > max_rows: df = df.head(max_rows) return df.to_html(index=False, escape=False, classes="data") # データ仕様表 T_dataspec = pd.DataFrame([ {"dataset_id": 55, "シリーズ": "砂防関係指定地情報_砂防指定地", "形式": "Shapefile (Polygon + Line)", "サイズ": "5.2 MB (ZIP)", "件数": f"polygon {len(g55_poly):,} + line {len(g55_line)}", "根拠法律": "砂防法 (1897, 明治 30 年)", "対象現象": "渓流の土石流", "管理者": "国交省 + 県 (1 級水系直轄区間は国)", "本記事での役割": "面・線で表現された行為制限区域 (戦前指定は line として残存)"}, {"dataset_id": 56, "シリーズ": "砂防関係指定地情報_急傾斜地崩壊危険区域", "形式": "CSV (緯度経度 1 点)", "サイズ": "0.45 MB", "件数": f"{n_56:,} 行", "根拠法律": "急傾斜地法 (1969, 昭和 44 年)", "対象現象": "30°以上急斜面崩壊から人家を守る", "管理者": "県 + 市町村", "本記事での役割": "ポリゴン未整備のため点データのみ。代表点で位置を特定"}, {"dataset_id": 57, "シリーズ": "砂防関係指定地情報_地すべり防止区域", "形式": "Shapefile (Polygon)", "サイズ": "0.02 MB (ZIP)", "件数": f"{len(g57)} 件", "根拠法律": "地すべり等防止法 (1958, 昭和 33 年)", "対象現象": "地すべり地形の進行抑止", "管理者": "国交省 (建設) + 林野庁 (山林) + 農水省 (農地) — 三省共管", "本記事での役割": "極少数だが高品質な polygon。三省共管制の特殊性を映す"}, ]) T_dataspec.to_csv(ASSETS / "L46_dataspec.csv", index=False, encoding="utf-8-sig") # 結合不要 (3 dataset は別の行政台帳) を明示する STEP 表 T_steps = pd.DataFrame([ {"段階": "STEP 0", "内容": "3 dataset を独立に読込", "サイズ/件数": f"55={len(g55_poly):,}+{len(g55_line)}, 56={n_56:,}, 57={len(g57)}"}, {"段階": "STEP 1", "内容": "EPSG:6668/EPSG:4326 → EPSG:6671 (m 単位) 投影変換", "サイズ/件数": "砂防/地すべり/急傾斜 すべて"}, {"段階": "STEP 2", "内容": "55: polygon area, line length 計算", "サイズ/件数": f"area sum {total_55_poly_ha:,.0f} ha, line sum {total_55_line_km:.1f} km"}, {"段階": "STEP 3", "内容": "57: polygon area 計算 + 列名 col01-col10 を意味的列名にリネーム", "サイズ/件数": f"area sum {total_57_ha:.1f} ha"}, {"段階": "STEP 4", "内容": "56: 1 件不正レコード除外 + 緯度経度→Point geometry 化", "サイズ/件数": f"2,935 → {n_56:,} 行"}, {"段階": "STEP 5", "内容": "市町別 / 水系別 / 種別 / 時代別の集計表を作成", "サイズ/件数": f"市町別 {len(city_summary)} / 水系別 {len(suikei)} / 時代別 {len(T_decade)}"}, {"段階": "STEP 6", "内容": "L11 警戒区域の件数を読み込み、3 種比較表を作成 (面積処理は行わない)", "サイズ/件数": "急傾斜 29,756, 地すべり 127, 土石流 13,337 (件数のみ)"}, ]) T_steps.to_csv(ASSETS / "L46_steps.csv", index=False, encoding="utf-8-sig") sections = [] # (1) 学習目標と問い sections.append(("学習目標と問い", f"""

本レッスンでは、広島県オープンデータ DoBoX が公開する 砂防関係指定地情報シリーズ 3 dataset (id = 55, 56, 57) を統合し、 砂防三法 (砂防法・地すべり等防止法・急傾斜地法) に基づく 3 種の指定地の 制度的・地理的構造を分析する。

独自用語の定義 (本記事内での使用)

主 RQ (研究の問い)

広島県の砂防三法 (砂防法・地すべり等防止法・急傾斜地法) に基づく 3 種の指定地は、 件数規模・地理範囲・対象市町・指定の時代でどう異なり、 なぜそのような「制度の三国志」の構造が生まれたのか？ さらに「警戒区域」(L11 既扱) との制度的差をデータでどう読めるか？

仮説 H1〜H6 (本記事で検証する)

到達点

3 dataset を市町・水系・指定種別・時代でクロスし、6 仮説を量的に検証する。 学習者は (a) 異種形式 (Shapefile + CSV) の並列読込、 (b) 「同じシリーズ」を名乗りながら制度的に独立な 3 dataset を扱う際の 分析設計、(c) 砂防三法という歴史と地理が交錯する制度群を データで読み解く視座、を体得する。

厳密性の宣言: 本記事は L10 (土砂災害警戒区域 cross), L11 (トリプルハザード) と厳密に独立。両者は警戒区域 (情報公示制) を扱うが、本記事は指定地 (行為制限制) を扱う。同じ土砂災害というドメインだが、法的位置づけが質的に違う。 L39 (傾斜) や L40 (標高) で扱う地形ラスタとも別物。 比較対照として L11 の警戒区域件数を一部参照するが、ジオメトリ処理は行わず件数比較のみ。
""")) # (2) 使用データ sections.append(("使用データ", f"""

本記事は DoBoX の「砂防関係指定地情報」関連 3 dataset を統合する:

表 (1) — 3 dataset の仕様サマリ

{df_to_html(T_dataspec)}

この表から読み取れること: (1) 形式が 3 種それぞれ違う (Shapefile poly+line / CSV / Shapefile poly のみ) — 同じ 「砂防関係指定地情報」シリーズを名乗りながらデータ整備の優先度が異なる。 (2) サイズが 250 倍差 (55: 5.2 MB ≫ 57: 0.02 MB) — これは件数差を反映するが、 57 は polygon 1 件あたりの geometry 詳細度は 55 と同等であり、件数自体が圧倒的に少ない。 (3) 根拠法律は 112 年差 (1897 vs 2009 のような新法ではない、1897 vs 1969 = 72 年差) で、 日本の砂防行政の歴史的層がデータに刻まれている。 (4) 地すべり防止区域の三省共管 (国交+林野+農水) は他 2 法にない特殊性で、 これが指定数を抑える制度的要因の 1 つと推測される。

表 (2) — 処理 STEP 段階表 (要件 K Before/After)

{df_to_html(T_steps)}

この表から読み取れること: 3 dataset は結合 (merge) 不要な独立した行政台帳である。 そこが本記事と L45 (ため池) との大きな違い: L45 は CSV と SHP が「同じ池の異なる断面」 として merge できたが、L46 の 3 dataset は「別の現象を別の法律で別に管理」 した別物。よって本記事は merge ではなく「並列分析と比較」が中心となる。

形式の詳細

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

HTML から直接以下を取得できる:

(A) DoBoX 直リンク (3 dataset)

datasetカタログresource_download (直 DL)形式サイズ
55 砂防関係指定地情報_砂防指定地 dataset 55 直 DL (rid 76) Shapefile (ZIP)5.2 MB
56 砂防関係指定地情報_急傾斜地崩壊危険区域 dataset 56 直 DL (rid 22822) CSV (UTF-8 BOM)0.45 MB
57 砂防関係指定地情報_地すべり防止区域 dataset 57 直 DL (rid 22823) Shapefile (ZIP)0.02 MB

(B) PowerShell 一括取得 (再現用)

cd "2026 DoBoX 教材"
mkdir -Force data\\extras\\L46_sabo_designation
iwr "https://hiroshima-dobox.jp/resource_download/76"    -OutFile "data\\extras\\L46_sabo_designation\\sabo_designation_55_sabo_shitei.zip"
iwr "https://hiroshima-dobox.jp/resource_download/22822" -OutFile "data\\extras\\L46_sabo_designation\\sabo_designation_56_kyukeisha_kuiki.csv"
iwr "https://hiroshima-dobox.jp/resource_download/22823" -OutFile "data\\extras\\L46_sabo_designation\\sabo_designation_57_jisuberi_boushi.zip"
# ZIP 展開 (55 と 57 のみ)
Expand-Archive -Force "data\\extras\\L46_sabo_designation\\sabo_designation_55_sabo_shitei.zip"   "data\\extras\\L46_sabo_designation\\sabo_designation_55_sabo_shitei"
Expand-Archive -Force "data\\extras\\L46_sabo_designation\\sabo_designation_57_jisuberi_boushi.zip" "data\\extras\\L46_sabo_designation\\sabo_designation_57_jisuberi_boushi"
# 本記事スクリプト (~30 秒)
py -X utf8 lessons\\L46_sabo_designation.py

(C) 中間データ・図 (本記事生成物の直リンク)

""")) # (4) 分析 1: 規模 100 倍非対称 sections.append(("分析 1: 規模 100 倍非対称 — 3 dataset の件数構造", f"""

狙い (RQ ↔ 仮説 H1)

3 dataset の件数を比較し、「同じ砂防関係指定地」シリーズを名乗りながら 規模が桁レベルで違う構造を量化する。これは砂防三法の制度的非対称の最も基礎的な指標であり、 後続分析の出発点となる。

手法 — 単純カウントと比率

各 dataset を独立に読み込み、行数 / polygon 数 / point 数を取得する。 最少件数 (地すべり 39) を分母として、他 2 種の倍率を計算する。 件数を log 軸の横棒グラフで表示することで、3 桁の差を視覚化する。

log 軸の選択理由 (要件 H, J): 線形軸でグラフを描くと、最大件数 ({len(g55_poly):,}) が最少件数 ({len(g57)}) を圧倒し、 地すべりの棒が事実上見えなくなる。log 軸では「桁」がそのまま視覚化されるので、 3 dataset の規模オーダーを一目で比較できる。

実装コード

{code(''' import geopandas as gpd import pandas as pd # 砂防指定地 (55) - Shapefile g55_poly = gpd.read_file('data/extras/L46_sabo_designation/.../砂防指定地-ポリゴン.shp') g55_line = gpd.read_file('data/extras/L46_sabo_designation/.../砂防指定地-ライン.shp') # 急傾斜地 (56) - CSV df56 = pd.read_csv('data/extras/L46_sabo_designation/sabo_designation_56_kyukeisha_kuiki.csv') # 1 件不正レコードを除外 (履歴ID 列に区域名が入っている) df56 = df56[df56["砂防関係指定地の種類"] == "急傾斜地崩壊危険区域"] # 地すべり (57) - Shapefile g57 = gpd.read_file('data/extras/L46_sabo_designation/.../地すべり防止区域.shp') # 件数比較 n_55, n_55l, n_56, n_57 = len(g55_poly), len(g55_line), len(df56), len(g57) ratio_55_57 = n_55 / n_57 ratio_56_57 = n_56 / n_57 print(f"砂防 {{n_55:,}} : 急傾斜 {{n_56:,}} : 地すべり {{n_57}} = " f"{{ratio_55_57:.0f}} : {{ratio_56_57:.0f}} : 1") ''')}

結果の図

図 1 を選んだ理由: 4 つの数値 (55-poly, 55-line, 56, 57) を比較するなら、 小数の棒グラフで十分。log 軸にすることで桁差を視覚化し、 件数を棒の右に明示することで「視覚的な差」と「実数」の両方が読める。

{figure('assets/L46_fig1_three_law_counts.png', '図1: 砂防三法 3 dataset の件数比較 (log 軸)')}

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

結果の表 (1) — 全体サマリ

{df_to_html(T_overall)}

表 (1) から読み取れること:

結果の図 (2) — 3 種指定地 全県重ね合わせマップ

図 2 を選んだ理由: 件数比較 (図 1) で3 種の規模差を見たあと、 地理空間でどう分布しているかを 1 枚で見ておくと、後続の分析 2-6 への伏線になる。 3 種を別の色 (緑面 / 濃緑線 / 赤点 / 黄面) で重ねることで、 「同じ県の同じ砂防三法でも分布パターンが異なる」ことが直感的に分かる。

{figure('assets/L46_fig2_three_law_overlay.png', '図2: 砂防三法指定地 3 種 全県分布マップ')}

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

""")) # (5) 分析 2: 法律 → 地理パターン分離 sections.append(("分析 2: 法律 → 地理パターン分離 — 3 種の市町集中構造", f"""

狙い (RQ ↔ 仮説 H2)

3 種の指定地がどの市町に集中しているかを比較する。 法律ごとに対象現象が異なる (土石流 / がけ崩れ / 地すべり) ので、 地形的・地質的に発生しやすい市町と歴史的に人口が集中した市町の 2 つの軸で分布が分かれるはずである。

手法 — 市町別 stacked bar + 占有率計算

3 dataset を市町列で value_counts() し、 外側結合 (outer join) で「ある dataset には現れるが他には現れない」市町も保持する。 スタックバーで 3 種の積み上げを表示し、上位市町の3 種混在比率を比較する。

外側結合 (outer join) の選択理由 (要件 K): 内側結合 (inner) を使うと、3 dataset すべてに該当する市町だけが残り、 たとえば地すべり防止区域がない 19 市町が消えるので公平な比較ができない。 outer join で全市町 (32 程度) を保持し、欠損は 0 で埋めるのが正しい。

実装コード

{code(''' import pandas as pd # 各 dataset から市町列の value_counts city_55 = g55_poly["city"].value_counts().rename("砂防_件").to_frame() city_55["砂防_面積ha"] = g55_poly.groupby("city")["area_ha"].sum() city_56 = g56["市区郡町"].value_counts().rename("急傾斜_件").to_frame() city_57 = g57["city"].value_counts().rename("地すべり_件").to_frame() # outer join で全市町を保持 city_summary = ( city_55.join(city_56, how="outer") .join(city_57, how="outer") .fillna(0) ) city_summary["三法合計件"] = (city_summary["砂防_件"] + city_summary["急傾斜_件"] + city_summary["地すべり_件"]) city_summary = city_summary.sort_values("三法合計件", ascending=False) print(city_summary.head(10)) ''')}

結果の図 (1) — 市町別 stacked bar

図 3 を選んだ理由: stacked bar は「合計と内訳の比率」を 1 図で示す定番手法。 市町の長い名前 (広島市安佐南区等) は横棒の方が読みやすい。

{figure('assets/L46_fig3_per_city_stacked.png', '図3: 市町別 砂防三法指定地件数 stacked bar')}

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

結果の図 (2) — 急傾斜地点マップ (呉市集中強調)

図 4 を選んだ理由: 呉市集中の地理を見るには地図がベスト。 呉市点を赤、それ以外を灰色にすることで、地理的偏りが一目で分かる。

{figure('assets/L46_fig4_kyukei_kure.png', '図4: 急傾斜地崩壊危険区域 点マップ — 呉市集中の地理')}

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

結果の表 — 上位 15 市町

{df_to_html(city_summary.head(15))}

表から読み取れること:

結果の図 (3) — 地すべり防止区域 詳細マップ

図 6 を選んだ理由: 39 件しかないので全件を地図に表示できる。 各区域に id 番号を振ると、データ集計表と地図の対応が直接取れる。

{figure('assets/L46_fig6_jisuberi_detail.png', '図6: 地すべり防止区域 全 39 件マップ')}

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

""")) # (6) 分析 3: 時代の制度展開 sections.append(("分析 3: 時代の制度展開 — 砂防三法の歴史的層", f"""

狙い (RQ ↔ 仮説 H3)

3 dataset の指定年月日から、「いつ・どの法律で・何件指定されたか」の 時代別構造を読み解く。砂防法 (1897) → 地すべり等防止法 (1958) → 急傾斜地法 (1969) の 立法順序が、データに刻まれているはずである。

手法 — 10 年区切り decade 集計 + 法律施行年マーカー

3 dataset の指定年月日を pd.to_datetime で日付化し、 decade = year // 10 * 10 で 10 年区切りに丸める。 スタックバーに 3 法の施行年 (1897/1958/1969) を縦線で重ね、 「立法 → 指定の波 → 制度的成熟・停滞」の流れを可視化する。

実装コード

{code(''' import pandas as pd def to_decade(s): d = pd.to_datetime(s, errors="coerce") return d.dt.year // 10 * 10 g55_poly["decade"] = to_decade(g55_poly["sitei_ymd"]) g56["decade"] = to_decade(g56["指定年月日"]) g57["decade"] = to_decade(g57["sitei_ymd"]) decade_55 = g55_poly["decade"].value_counts().sort_index() decade_56 = g56["decade"].value_counts().sort_index() decade_57 = g57["decade"].value_counts().sort_index() # 中央指定年 median_55 = pd.to_datetime(g55_poly["sitei_ymd"], errors="coerce").dt.year.median() median_56 = pd.to_datetime(g56["指定年月日"], errors="coerce").dt.year.median() median_57 = pd.to_datetime(g57["sitei_ymd"], errors="coerce").dt.year.median() print(f"中央年: 砂防 {{median_55:.0f}}, 急傾斜 {{median_56:.0f}}, 地すべり {{median_57:.0f}}") ''')}

結果の図

図 7 を選んだ理由: 横軸=年代 + 縦軸=件数 + 色=3 法 + 縦線=立法年、の 4 情報を 1 図で示すには stacked bar with vertical lines がベスト。

{figure('assets/L46_fig7_decade_stacked.png', '図7: 砂防三法 指定の時代別累積 stacked bar')}

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

結果の表 — 時代別件数

{df_to_html(T_decade)}

表から読み取れること:

""")) # (7) 分析 4: 指定地 vs 警戒区域 ~10 倍差 sections.append(("分析 4: 指定地 vs 警戒区域 — 制度的厳しさ 10 倍差", f"""

狙い (RQ ↔ 仮説 H4)

本記事 (砂防三法指定地) と L11 (土砂災害警戒区域) は 同じ土砂災害を扱うが法的位置づけが質的に違う。 件数比で「制度的厳しさ」を量化する。

手法 — 件数比較バー (3 種類対 × log 軸)

L11 で読み込んだ警戒区域 Shapefile の行数だけを取得し、 本記事の指定地件数と3 種類対で比較する。比率は警戒区域 ÷ 指定地。

制度的差の意味 (要件 J ツール化視点):

実装コード

{code(''' import geopandas as gpd # L11 で扱った警戒区域 Shapefile (面積処理は不要、件数のみ) warn_kyukei = gpd.read_file("data/extras/sediment_shp/kyukeisha/.../73_031krp_34_20260427.shp") warn_jisuberi = gpd.read_file("data/extras/sediment_shp/jisuberi/.../73_031jy_34_20260427.shp") warn_doseki = gpd.read_file("data/extras/sediment_shp/doseki/.../73_031drp_34_20260427.shp") # 件数比較 (gdf を保持せず長さだけ) warn_counts = {{"急傾斜": len(warn_kyukei), "地すべり": len(warn_jisuberi), "土石流": len(warn_doseki)}} # 指定地件数 (本記事) designated = {{"土石流": len(g55_poly), # 砂防指定地 "がけ崩れ": len(g56), # 急傾斜地崩壊危険区域 "地すべり": len(g57)}} # 地すべり防止区域 # 3 種類対の倍率 print(f"土石流 警戒/指定 = {{warn_counts['土石流']/designated['土石流']:.1f}} 倍") print(f"がけ崩れ 警戒/指定 = {{warn_counts['急傾斜']/designated['がけ崩れ']:.1f}} 倍") print(f"地すべり 警戒/指定 = {{warn_counts['地すべり']/designated['地すべり']:.1f}} 倍") ''')}

結果の図

図 8 を選んだ理由: 3 現象 × 2 制度 = 6 セルの比較を、 同じ x 位置にペアで並べた縦棒で表示。倍率を bar の上に annotate することで、 「警戒区域は指定地の何倍か」が即読み取れる。log 軸で件数差をさらに強調。

{figure('assets/L46_fig8_designated_vs_warning.png', '図8: 指定地 vs 警戒区域 件数比較')}

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

結果の表 — 指定地 vs 警戒区域

{df_to_html(T_compare)}

表から読み取れること:

""")) # (8) 分析 5: 砂防指定地の水系構造 sections.append(("分析 5: 砂防指定地の水系構造 — 1 級・2 級・その他", f"""

狙い (RQ ↔ 仮説 H1 補強)

砂防指定地 (55) は唯一水系メタ情報 (suikei_gr, suikei_nm, kansen_nm) を持つ。 これを使い、「どの水系に砂防指定地が集中しているか」を読み解く。 河川行政の1 級水系 = 国直轄、2 級 = 県、その他 = 市町村関与という階層と、 指定地の地理分布の関係を見る。

手法 — groupby 水系 + 水系別 polygon マップ

砂防指定地を suikei_gr (水系区分) と suikei_nm (水系名) で groupby し、件数・総面積・平均面積を集計する。 さらに 1 級/2 級/その他 で色分けした polygon マップを描き、 水系階層と地理分布の関係を可視化する。

実装コード

{code(''' # 水系別 集計 suikei = g55_poly.groupby(["suikei_gr", "suikei_nm"]).agg( 件数=("id", "count"), 総面積ha=("area_ha", "sum"), 平均面積ha=("area_ha", "mean"), ).reset_index().sort_values("件数", ascending=False) # 1 級水系トップ top_lev1 = suikei[suikei["suikei_gr"] == "1級"].head(5) print("1 級水系 トップ 5:") print(top_lev1[["suikei_nm", "件数", "総面積ha"]].to_string()) ''')}

結果の図

図 5 を選んだ理由: 水系区分 (1 級/2 級/その他) を色で分け、 全 polygon を県全域に重ねることで「県の主要水系の砂防地理」が見える。 line データも黒で重ね、戦前指定の場所も同時に表示。

{figure('assets/L46_fig5_sabo_suikei.png', '図5: 砂防指定地 水系別 polygon マップ')}

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

結果の表 — 水系別 上位 15

{df_to_html(suikei.head(15))}

表から読み取れること:

""")) # (9) 分析 6: 上位 4 市町ズーム — 三法重ね sections.append(("分析 6: 上位 4 市町ズーム — 三法重ねの空間構造", f"""

狙い (RQ 統合)

三法合計件数の上位 4 市町に絞り、3 種指定地 (+ line) を同じズーム範囲で small multiples 並置することで、「市町ごとに異なる三法の混在比率」を可視化する。 これは分析 1-3 の結果を空間的に統合する仕上げの図。

手法 — small multiples (2x2 grid)

4 市町の bbox (外接矩形) を共通計算し、各 panel で同じ凡例で 4 種要素 (砂防 poly, 砂防 line, 急傾斜点, 地すべり poly) を描く。背景は admin_diss で県全体の市町境界。

実装コード

{code(''' top4 = city_summary.head(4)["市町"].tolist() fig, axes = plt.subplots(2, 2, figsize=(13, 10)) for ax, city in zip(axes.flat, top4): sub_55 = g55_poly[g55_poly["city"] == city] sub_55l = g55_line[g55_line["city"] == city] sub_56 = g56[g56["市区郡町"] == city] sub_57 = g57[g57["city"] == city] # bbox: 全要素の外接矩形 + 1.5 km マージン candidates = [g.total_bounds for g in [sub_55, sub_55l, sub_57] if len(g)] if len(sub_56): candidates.append([sub_56.geometry.x.min(), sub_56.geometry.y.min(), sub_56.geometry.x.max(), sub_56.geometry.y.max()]) candidates = np.array(candidates) xmin, ymin = candidates[:, 0].min() - 1500, candidates[:, 1].min() - 1500 xmax, ymax = candidates[:, 2].max() + 1500, candidates[:, 3].max() + 1500 admin_diss.plot(ax=ax, color="#f8f8f8", edgecolor="#aaa", linewidth=0.4) sub_55.plot(ax=ax, color="#388e3c", alpha=0.55) sub_55l.plot(ax=ax, color="#1b5e20", linewidth=1.5) ax.scatter(sub_56.geometry.x, sub_56.geometry.y, s=10, c="#cf222e", alpha=0.65) sub_57.plot(ax=ax, color="#f9a825", alpha=0.85) ax.set_xlim(xmin, xmax); ax.set_ylim(ymin, ymax) ''')}

結果の図

図 9 を選んだ理由: small multiples は「複数の地域を同じレンズで比較する」 ための定番構図。各市町の三法混在パターンが並置で読める。

{figure('assets/L46_fig9_top4_cities_zoom.png', '図9: 三法合計件数 上位 4 市町 ズームマップ')}

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

結果の表 — Top 砂防指定地 面積ランキング

{df_to_html(top_55_view)}

表から読み取れること:

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

6 仮説の検証結果

{df_to_html(df_hypo)}

主要発見の整理

3 dataset 相互関係の構造発見

本記事の最重要発見は、3 dataset が「同じ砂防関係指定地」シリーズを名乗りながら 制度的にも形式的にも独立した別物である点である。

3 dataset は結合不能 (= merge する共通キーが無い)。 代わりに「市町別カウント + 時代別カウント + 形式比較」の並列分析で構造を読む。 これは L45 (ため池, 2 dataset を merge) との分析設計の根本的な違いである。

「指定地」と「警戒区域」の制度的補完

本記事の指定地 (砂防三法, 戦前〜1970s 法律) と L11 の警戒区域 (土砂災害防止法 2000) は、 同じ土砂災害ドメインの異なる時代の制度として相補的:

本研究の限界

""")) # (11) 発展課題 sections.append(("発展課題", f"""

本記事の結果から論理的に導かれる新仮説と、それを検証する具体的課題:

課題 1: 指定地 × 警戒区域 空間結合 — 二重保護エリアの定量化

課題 2: 急傾斜地崩壊危険区域 ポリゴン整備 — 県への提言

課題 3: 砂防指定地 line ({len(g55_line)} 件) の歴史調査

課題 4: 地すべり防止区域 × 地質図 重ね合わせ

課題 5: 呉市急傾斜地 × 西日本豪雨被害 マッチング

課題 6: 三省共管制 (地すべり防止区域) の停滞要因

""")) html = render_lesson( num=46, title="L46 砂防関係指定地 3 件統合分析 — 砂防三法 (砂防法・地すべり等防止法・急傾斜地法) 指定地構造を読み解く", tags=["砂防三法", "砂防指定地", "急傾斜地崩壊危険区域", "地すべり防止区域", "行為制限", "warning vs designated", "geopandas", "Shapefile", "CSV"], time="60〜90 分", level="リテラシ〜中級 (geopandas/Shapefile-CSV 並列分析/制度比較)", data_label="DoBoX dataset 55, 56, 57 (砂防指定地 + 急傾斜地崩壊危険区域 + 地すべり防止区域) — Shapefile 5.2 MB + CSV 0.45 MB + Shapefile 0.02 MB", sections=sections, script_filename="L46_sabo_designation.py", ) out_html = LESSONS / "L46_sabo_designation.html" out_html.write_text(html, encoding="utf-8") print(f" HTML written: {out_html} ({out_html.stat().st_size:,} B), {time.time()-t1:.1f}s", flush=True) print(f"\n=== L46 done in {time.time()-t_all:.1f}s ===", flush=True)