# -*- coding: utf-8 -*- """L54 宅地造成等工事規制区域 単独 3 研究例分析 — 広島県の「宅造規制区域」 24 区域・約 991 km² を 3 つの独立した研究角度で読み解く カバー宣言: 本記事は DoBoX のシリーズ「宅地造成等工事規制区域」 1 件 (dataset_id = 1427) を 単独で取り上げ、 広島県内の宅造規制区域 24 区域 (Shapefile 120 polygon, 約 991 km²) を 3 つの独立した研究角度 (RQ1 / RQ2 / RQ3) で並列に分析する。 「宅地造成等工事規制区域」 とは: 盛土規制法 (= 宅地造成及び特定盛土等規制法, 2023-05-26 施行) 第 12 条第 1 項に基づき、知事が指定する区域。 旧法「宅地造成等規制法」 (1961-2023) の指定区域を引き継ぎ、 都市計画区域内の宅地等で盛土・切土工事に許可が必要なエリア。 本データは広島県知事指定の 24 区域 (24 市町) を Shapefile (Polygon × 120) で公開したもの。指定面積は 約 991 km² (広島県総面積 8,479 km² の 11.7%)。 「12 条 1 項」 (本記事独自定義の中核条文): 『都道府県知事は、宅地造成、特定盛土等又は土石の堆積に伴う崖崩れ又は 土砂の流出による災害の防止のため必要があると認めるときは、 関係市町村長の意見を聴いて、宅地造成等工事規制区域として指定する』 つまり都市計画区域内の市街地周辺で、宅造工事の許可制を実施する地理範囲を 法的に固定する。一度指定されると当該区域内で許可不要規模を超える工事は 全件 12 条 1 項許可を要する。 「許可閾値」 (盛土規制法施行令 第 3 条): 本区域内で許可必要となる工事の規模: - 盛土で高さ 1m 超の崖を生ずる工事 - 切土で高さ 2m 超の崖を生ずる工事 - 高さ閾値以下でも面積500m² 超の盛土・切土 - 土石の堆積で面積 300m² 超 + 高さ 2m 超または面積 500m² 超 閾値以下でも本区域内なら21 条 1 項届出義務が発生する (= L53 で扱った届出制)。 「指定継承」 (本記事独自定義): 旧「宅地造成等規制法」 (1961 年制定) で指定された区域は、 新法施行 (2023-05-26) 時点で新法の宅地造成等工事規制区域として 自動的に引き継がれた (盛土規制法附則第 2 条)。 本データの指定日は2023-04-27 (= 新法施行直前) で、 実質的には旧法時代の指定をそのまま継承。 「二重リスクエリア」 (本記事独自定義): 宅造規制区域 (= 行政が宅造工事を許可制で監督) と 土砂災害警戒区域 (= 自然災害想定エリア) が地理的に重なる場所。 本来「規制 + 警戒」 の二重防御が働くべき領域だが、 指定タイミングのずれや空間的な微差で「重なるが連携が薄い」 状況がありうる。 RQ3 で定量化する。 研究の問い (3 RQ): RQ1 (主研究): 広島県内の宅造規制区域の地理範囲と市町指定状況。 24 区域・120 polygon の面積分布、市町別ランキング、地形タイプ (沿岸都市 / 内陸都市 / 山間中心)、指定面積 (km²) の偏在を読む。 広島県総面積 8,479 km² のうち本区域指定は何 % を占めるか? RQ2 (副研究 1): 許可盛土 (L52 既扱) の発生密度との関係。 L52 許可盛土 152 件のうち何 % が本規制区域内に立地するか? 本区域内の許可盛土密度 (件/km²) はどれほどか? 区域外で許可された盛土 (= 特定盛土区域・40 条 1 項側) の地理特性は? 「区域指定 = 許可盛土集中」 の仮説を検証する。 RQ3 (副研究 2): 土砂災害警戒区域 (L10/L11 既扱) との重ね合わせ。 本規制区域と急傾斜地崩壊+土石流警戒区域の地理重複面積、 「警戒域 + 規制区域」 の二重リスクエリアの市町別偏在、 旧法時代からの指定継承による盲点 (= 新たな警戒区域指定が 規制区域外に発生している場合) を診断する。 仮説 (5): H1 (都市集中, RQ1): 24 区域は人口の多い沿岸都市部市町に偏在し、 上位 5 市町指定面積で全体の 50% 以上を占める。 山間部市町は指定なし、または小面積。 H2 (区域内立地率高, RQ2): L52 許可盛土 152 件のうち、本規制区域内立地は 70% 以上。区域外立地は特定盛土区域 (法 30 条) 側の許可で、 少数派 (30% 未満)。 H3 (区域内密度の市町偏在, RQ2): 本区域内の許可盛土密度 (件/km²) は 市町別に 2 桁以上の幅で偏在。 東広島市等の都市部は密度高、山間市町区域内は密度低。 H4 (警戒区域との部分重複, RQ3): 本規制区域 ∩ 警戒区域の重複面積比は 20-40%。完全包摂ではなく、規制区域は市街地中心を、 警戒区域は急傾斜・渓流を捉えるため空間目的が異なる。 H5 (二重リスクエリアの市町偏在, RQ3): 二重リスク (規制区域内かつ警戒区域内) の 市町別ランキングは「警戒区域面積トップ」 と「規制区域面積トップ」 の クロス積で決まる。広島市・呉市・尾道市が上位 3。 要件 S 準拠 (1分以内完走): - 規制区域 Shapefile は 120 polygon (named) のみ抽出。bbox + STRtree で高速化 - 警戒区域 polygon (急傾斜 30K + 土石流 13K) は dissolve で先に統合 - L52 の処理済 CSV (assets/L52_permits_processed.csv) を参照することで 比較研究は追加コスト ~0 秒 - 期待実行時間: <60 秒 要件 T 準拠 (位置情報あり=地図必須): - RQ1: 24 区域マップ + 市町別指定面積 choropleth - RQ2: 規制区域 × L52 許可盛土 オーバーレイマップ - RQ2: 区域内/外の盛土の規模分布比較 - RQ3: 規制区域 × 警戒区域 重ね合わせマップ - RQ3: 二重リスクエリア choropleth 要件 Q 準拠: 図 9 / 表 12 (3 RQ × 多角度)。 データ仕様: - dataset 1427: 宅地造成等工事規制区域 - リソース 51124 (Shapefile, 2.3 MB): 県統合 24 区域 (Polygon × 120) - リソース 51145 (XLSX, 23.1 KB): 盛土規制法関係資料 (規制内容 + 手続き先) - 形式: Shapefile (EPSG:2445 → 6671 変換) + XLSX (2 シート) - 取得日: 2026-05-09 - ライセンス: クリエイティブ・コモンズ表示 4.0 - 作成主体: 広島県土木建築局 都市環境整備課 """ from __future__ import annotations import sys, time, re, json, zipfile, shutil 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 import shapely from shapely import STRtree from shapely.ops import unary_union 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("=== L54 宅地造成等工事規制区域 単独 3 研究例分析 ===", flush=True) # ============================================================================= # 0. 定数・パス # ============================================================================= TARGET_CRS = "EPSG:6671" # JGD2011 平面直角座標系第 III 系 (m 単位) SOURCE_CRS = "EPSG:2445" # JGD2000 平面直角座標系第 III 系 (本データの原 CRS) DATA_DIR = ROOT / "data" / "extras" / "L54_residential_filling_zone" DATA_DIR.mkdir(parents=True, exist_ok=True) ZONE_ZIP = DATA_DIR / "340006_residential_land_construction_regulation_area_shp_20230427.zip" ZONE_SHP_DIR = DATA_DIR / "shp" XLSX_PATH = DATA_DIR / "340006_earth_filling_regulation_related_documents.xlsx" # 市町別 Shapefile (RID 51125-51144 = 20 市町) CITY_DIR = DATA_DIR / "all_cities" CITY_DIR.mkdir(parents=True, exist_ok=True) CITY_SHP_RIDS = [ (51125, "342033"), (51126, "342041"), (51127, "342050"), (51128, "342084"), (51129, "342092"), (51130, "342106"), (51131, "342114"), (51132, "342122"), (51133, "342131"), (51134, "342149"), (51135, "342157"), (51136, "343099"), (51137, "343684"), (51138, "343692"), (51139, "344621"), (51140, "345458"), (51141, "343021"), (51142, "343048"), (51143, "343072"), (51144, "344311"), ] # 関連 Shapefile (土砂災害警戒区域 = L10/L11 既存 sediment_shp を流用) WARN_KYUKEI = ROOT / "data" / "extras" / "sediment_shp" / "kyukeisha" / \ "340006_sediment_disaster_hazard_area_steep_slope_20260427" / \ "73_031krp_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" # L52 (許可盛土) 既処理 CSV (RQ2 比較研究用) L52_PROCESSED = ASSETS / "L52_permits_processed.csv" # DoBoX DOBOX_BASE = "https://hiroshima-dobox.jp" HDR = {"User-Agent": "DoBoX-MDASH-textbook/1.0 (educational use)"} # 主要日付 LAW_ENFORCED = pd.Timestamp("2023-05-26") ZONE_DESIGNATED = pd.Timestamp("2023-04-27") # 本データの指定 (継承) 日付 PREF_AREA_KM2 = 8479.4 # 広島県総面積 (公的値) # 市町コード → 名前 (L52/L53 共通) 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: "大崎上島町", } NAME_TO_CITY_CD = {v: k for k, v in CITY_NAME.items()} # ============================================================================= # 1. データ取得 (なければダウンロード) # ============================================================================= print("\n[1] データ取得", flush=True) t1 = time.time() def ensure_resource(rid: int, dst: Path, label: str): if dst.exists() and dst.stat().st_size > 1000: print(f" cache HIT {label}: {dst.name}") return url = f"{DOBOX_BASE}/resource_download/{rid}" print(f" DL ← {url}") r = requests.get(url, headers=HDR, timeout=120) r.raise_for_status() dst.write_bytes(r.content) print(f" saved → {dst.name} ({len(r.content)} bytes)") ensure_resource(51124, ZONE_ZIP, "宅造規制区域 Shapefile (県知事指定 = 県統合)") ensure_resource(51145, XLSX_PATH, "規制関係 XLSX") # Shapefile を展開 ZONE_SHP_DIR.mkdir(parents=True, exist_ok=True) shp_path = ZONE_SHP_DIR / "340006_residential_land_construction_regulation_area_shp_20230427.shp" if not shp_path.exists(): print(f" unzipping → {ZONE_SHP_DIR}") with zipfile.ZipFile(ZONE_ZIP) as zf: zf.extractall(ZONE_SHP_DIR) # 市町別 Shapefile (20 件) もダウンロード print(f" 20 市町別 Shapefile も取得 (= 各市町長指定の宅造規制区域)...", flush=True) for rid, code5 in CITY_SHP_RIDS: fname = f"{code5}_residential_land_construction_regulation_area_shp_20230929.zip" ensure_resource(rid, CITY_DIR / fname, f"市町 {code5}") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 2. 規制区域 Shapefile の読込・前処理 (県知事 + 20 市町長指定を統合) # ============================================================================= print("\n[2] 規制区域 Shapefile 読込・前処理 (21 ファイル統合)", flush=True) t1 = time.time() # 2-1 県知事指定 (340006 = 51124, 2023-04-27 版) zone_pref_raw = gpd.read_file(shp_path) print(f" 県知事指定: source CRS {zone_pref_raw.crs}, rows {len(zone_pref_raw)}") zone_pref_raw = zone_pref_raw.to_crs(TARGET_CRS) # 区域整理番 != None の行のみ named (24 区域メタ情報あり) pref_named = zone_pref_raw[zone_pref_raw["区域整理番"].notna()].copy().reset_index(drop=True) pref_named["zone_id"] = pref_named["区域整理番"].astype(int) pref_named["指定主体"] = "県知事" pref_named["指定主体_code"] = "340006" pref_named["poly_area_km2"] = pref_named.geometry.area / 1e6 # 区域整理番=None の細片は本研究対象外 pref_unnamed = zone_pref_raw[zone_pref_raw["区域整理番"].isna()].copy().reset_index(drop=True) # 2-2 20 市町長指定 Shapefile を順次読込 city_polygons_list = [] for rid, code5 in CITY_SHP_RIDS: fname = f"{code5}_residential_land_construction_regulation_area_shp_20230929.zip" zip_path = CITY_DIR / fname extract_dir = CITY_DIR / f"tmp_{code5}" extract_dir.mkdir(exist_ok=True) if not any(extract_dir.glob("*.shp")): with zipfile.ZipFile(zip_path) as zf: zf.extractall(extract_dir) shp_files = list(extract_dir.glob("*.shp")) if not shp_files: print(f" WARN: {code5} no shp") continue g = gpd.read_file(shp_files[0]).to_crs(TARGET_CRS) g["指定主体"] = "市町長" g["指定主体_code"] = code5 g["poly_area_km2"] = g.geometry.area / 1e6 city_polygons_list.append(g) city_named = pd.concat(city_polygons_list, ignore_index=True) city_named = gpd.GeoDataFrame(city_named, geometry="geometry", crs=TARGET_CRS) print(f" 市町長指定: {len(CITY_SHP_RIDS)} 市町 / {len(city_named)} polygon") # 2-3 統合 (県知事 + 市町長) # 共通列のみ抽出して縦結合 common_cols = ["指定主体", "指定主体_code", "poly_area_km2", "geometry"] zones_all = pd.concat([ pref_named[common_cols], city_named[common_cols], ], ignore_index=True) zones_all = gpd.GeoDataFrame(zones_all, geometry="geometry", crs=TARGET_CRS) zones_all = zones_all.reset_index(drop=True) zones_all["poly_id"] = np.arange(len(zones_all)) print(f" 統合後 polygon: {len(zones_all)} (県知事 {len(pref_named)} + 市町長 {len(city_named)})") print(f" 統合面積: {zones_all['poly_area_km2'].sum():.2f} km² " f"(県知事 {pref_named['poly_area_km2'].sum():.2f} + 市町長 {city_named['poly_area_km2'].sum():.2f})") # 互換のため zone_pref_raw を残し、zones を統合 polygon に変更 zones = zones_all 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) print(f" admin: {len(admin_diss)} 市町ブロック") warn_kyukei = gpd.read_file(WARN_KYUKEI).to_crs(TARGET_CRS) warn_doseki = gpd.read_file(WARN_DOSEKI).to_crs(TARGET_CRS) print(f" 警戒区域: 急傾斜 {len(warn_kyukei):,} polygon, 土石流 {len(warn_doseki):,} polygon") # 大量 polygon を扱うので dissolve で先に統合 → overlay コスト削減 print(f" 警戒区域 dissolve (急傾斜 + 土石流統合)... ", flush=True) warn_all_geom = unary_union(list(warn_kyukei.geometry) + list(warn_doseki.geometry)) warn_all = gpd.GeoDataFrame(geometry=[warn_all_geom], crs=TARGET_CRS) print(f" warn_all 統合面積: {warn_all.geometry.iloc[0].area / 1e6:.2f} km²") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 4. polygon → 市町紐付け (admin sjoin) # ============================================================================= print("\n[4] polygon → 市町紐付け", flush=True) t1 = time.time() # 各 polygon の representative_point を市町 polygon と sjoin zones["rep_pt"] = zones.geometry.representative_point() zone_pts = gpd.GeoDataFrame( zones[["poly_id", "指定主体", "指定主体_code", "poly_area_km2"]], geometry=zones["rep_pt"], crs=TARGET_CRS) zone_city_join = gpd.sjoin(zone_pts, admin_diss[["CITY_CD", "geometry"]], how="left", predicate="within") # poly_id は重複しないので 1:1 対応 zone_city_join = zone_city_join.drop_duplicates(subset=["poly_id"], keep="first") zones["CITY_CD"] = zone_city_join.set_index("poly_id")["CITY_CD"].reindex( zones["poly_id"]).values zones["市町名"] = zones["CITY_CD"].map(CITY_NAME).fillna("不明") n_resolved = (zones["市町名"] != "不明").sum() print(f" polygon → 市町 解決: {n_resolved}/{len(zones)}") print(f" 市町別 polygon 数 上位 5:") print(zones["市町名"].value_counts().head(5)) # 市町別に dissolve (= 同市町の県知事/市町長指定 polygon を合算) zone_dissolved = zones.dissolve(by="市町名", as_index=False) zone_dissolved["zone_area_m2"] = zone_dissolved.geometry.area zone_dissolved["zone_area_km2"] = zone_dissolved["zone_area_m2"] / 1e6 zone_dissolved["zone_id"] = zone_dissolved["市町名"].map(NAME_TO_CITY_CD).fillna(-1).astype(int) print(f" 市町別 dissolve 後: {len(zone_dissolved)} 市町") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 5. RQ1 集計: 地理範囲と市町指定状況 # ============================================================================= print("\n[5] RQ1 集計", flush=True) t1 = time.time() # 5-1 市町別指定面積ランキング (zones を市町別集計, polygon 数も含む) city_summary = zones.groupby("市町名").agg( 指定polygon数=("poly_id", "count"), 指定面積_km2=("poly_area_km2", "sum"), ).reset_index().sort_values("指定面積_km2", ascending=False) # 県知事/市町長 別の polygon 数も追加 city_summary["県知事指定_n"] = city_summary["市町名"].map( zones[zones["指定主体"] == "県知事"].groupby("市町名").size().to_dict()).fillna(0).astype(int) city_summary["市町長指定_n"] = city_summary["市町名"].map( zones[zones["指定主体"] == "市町長"].groupby("市町名").size().to_dict()).fillna(0).astype(int) city_summary["順位"] = np.arange(1, len(city_summary) + 1) city_summary["県シェア_%"] = (city_summary["指定面積_km2"] / PREF_AREA_KM2 * 100).round(2) city_summary["指定面積シェア_%"] = (city_summary["指定面積_km2"] / city_summary["指定面積_km2"].sum() * 100).round(2) top5_share = city_summary.head(5)["指定面積_km2"].sum() / city_summary["指定面積_km2"].sum() total_zone_area_sum = city_summary["指定面積_km2"].sum() # 重複含み (= 県知事+市町長別々) # union (重複除去) 面積は後で計算 print(f" 全県指定面積 (合計, 重複含): {total_zone_area_sum:.2f} km²") total_zone_area = total_zone_area_sum # 互換のため (図表ラベルで使用) pref_share_sum = total_zone_area_sum / PREF_AREA_KM2 * 100 print(f" ↑ 県 8479 km² の {pref_share_sum:.2f}% (重複あり)") pref_share = pref_share_sum # 互換 # 指定主体別 polygon 数 (HTML 内で使用) n_pref_polys = (zones["指定主体"] == "県知事").sum() n_city_polys = (zones["指定主体"] == "市町長").sum() n_cities_designated = zone_dissolved["市町名"].nunique() print(f" Top 5 市町: {city_summary.head(5)[['市町名','指定面積_km2','指定面積シェア_%']].to_dict('records')}") print(f" H1 検証: Top 5 シェア = {top5_share*100:.1f}%") # 5-2 polygon 個別の面積分布 (24 区域内に 120 polygon の規模差) poly_size_stats = zones["poly_area_km2"].describe() print(f"\n polygon 規模分布:") print(f" n = {len(zones)}, 中央値 {poly_size_stats['50%']:.3f} km², " f"max {poly_size_stats['max']:.2f} km²") # 5-3 区域タイプ分類 (本研究独自定義) def zone_type(r): """市町タイプ (沿岸都市 / 内陸都市 / 山間)""" cn = r["市町名"] if cn in ("広島市中区", "広島市東区", "広島市南区", "広島市西区", "広島市安佐南区", "広島市安佐北区", "広島市安芸区", "広島市佐伯区", "呉市", "尾道市", "三原市", "福山市", "竹原市", "廿日市市", "大竹市", "東広島市", "江田島市", "府中町", "海田町", "熊野町", "坂町"): return "A_沿岸都市" if cn in ("府中市", "三次市", "庄原市", "安芸高田市"): return "B_内陸都市" if cn in ("世羅町", "安芸太田町", "北広島町", "神石高原町", "大崎上島町"): return "C_山間/島嶼" return "D_その他" zone_dissolved["地形タイプ"] = zone_dissolved.apply(zone_type, axis=1) type_summary = zone_dissolved.groupby("地形タイプ").agg( 市町数=("zone_id", "count"), 指定面積_km2=("zone_area_km2", "sum"), ).reset_index().sort_values("指定面積_km2", ascending=False) type_summary["シェア_%"] = (type_summary["指定面積_km2"] / type_summary["指定面積_km2"].sum() * 100).round(1) print(f"\n 地形タイプ集計:") print(type_summary) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 6. RQ2 集計: 許可盛土 (L52) との発生密度 # ============================================================================= print("\n[6] RQ2 集計: 許可盛土発生密度", flush=True) t1 = time.time() if L52_PROCESSED.exists(): df_l52 = pd.read_csv(L52_PROCESSED, encoding="utf-8-sig") print(f" L52 (許可盛土) 読込: {len(df_l52)} 件") g52 = gpd.GeoDataFrame( df_l52, geometry=gpd.points_from_xy(df_l52["lon"], df_l52["lat"]), crs="EPSG:4326").to_crs(TARGET_CRS) # 各盛土点が宅造規制区域内にあるか判定 # zone_dissolved の geometry を 1 つの union に zone_union_geom = unary_union(list(zone_dissolved.geometry)) g52["in_zone"] = g52.geometry.within(zone_union_geom) # 区域内/外 n_in_zone = int(g52["in_zone"].sum()) n_out_zone = len(g52) - n_in_zone in_zone_share = n_in_zone / len(g52) print(f" L52 区域内立地: {n_in_zone}/{len(g52)} 件 ({in_zone_share*100:.1f}%)") print(f" H2 検証: 区域内立地率 ≥ 70%? {'YES' if in_zone_share >= 0.70 else 'NO'}") # 市町別の盛土件数 (各市町の規制区域 polygon 内に何件あるか) # L52 の '市町名' 列と zone_dissolved の '市町名' 列が衝突するため、 # zone_dissolved 側の列を一時改名してから sjoin zd_for_join = zone_dissolved[["市町名", "zone_area_km2", "geometry"]].rename( columns={"市町名": "市町名_zone"}) g52_zone_join = gpd.sjoin(g52, zd_for_join, how="left", predicate="within") zone_permit_count = (g52_zone_join.dropna(subset=["市町名_zone"]) .groupby("市町名_zone") .size().reset_index(name="許可盛土件数")) zone_permit_count = zone_permit_count.rename(columns={"市町名_zone": "市町名"}) zone_density = zone_dissolved[["市町名", "zone_area_km2", "zone_id"]].merge( zone_permit_count, on="市町名", how="left").fillna(0) zone_density["許可盛土件数"] = zone_density["許可盛土件数"].astype(int) zone_density["密度_件_km2"] = (zone_density["許可盛土件数"] / zone_density["zone_area_km2"]).round(3) zone_density = zone_density.sort_values("密度_件_km2", ascending=False) print(f"\n 区域内盛土密度 (件/km²) Top 5:") print(zone_density.head(5)[["市町名", "zone_area_km2", "許可盛土件数", "密度_件_km2"]]) density_max = zone_density["密度_件_km2"].max() density_min_pos = zone_density.loc[zone_density["許可盛土件数"] > 0, "密度_件_km2"].min() if density_min_pos > 0: density_ratio = density_max / density_min_pos density_log_diff = np.log10(density_ratio) if density_ratio > 0 else 0 else: density_log_diff = 0 print(f" H3 検証: 密度の桁差 ({density_max:.2f} / {density_min_pos:.4f}) " f"= {10**density_log_diff:.0f} 倍 → " f"{'YES (≥2 桁)' if density_log_diff >= 2 else 'NO (<2 桁)'}") # 区域外立地の地理特性 (= 特定盛土区域での許可) g52_out = g52[~g52["in_zone"]].copy() out_city_count = (g52_out["市町名"].value_counts() if "市町名" in g52_out.columns else pd.Series(dtype=int)) print(f"\n 区域外立地 ({n_out_zone} 件) の市町分布 上位 5:") print(out_city_count.head(5)) # 規模比較 (in_zone vs out_zone) in_med = g52[g52["in_zone"]][["height_max_m", "area_m2", "fill_volume_m3"]].median() out_med = g52[~g52["in_zone"]][["height_max_m", "area_m2", "fill_volume_m3"]].median() print(f" 規模中央値 in_zone: 高さ {in_med['height_max_m']:.2f}m / " f"面積 {in_med['area_m2']:.0f}m² / 土量 {in_med['fill_volume_m3']:.0f}m³") print(f" 規模中央値 out_zone: 高さ {out_med['height_max_m']:.2f}m / " f"面積 {out_med['area_m2']:.0f}m² / 土量 {out_med['fill_volume_m3']:.0f}m³") L52_OK = True else: print(f" WARN: L52 processed CSV not found, RQ2 disabled") df_l52 = None g52 = None n_in_zone = 0 n_out_zone = 0 in_zone_share = np.nan zone_density = pd.DataFrame() L52_OK = False density_log_diff = 0 print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 7. RQ3 集計: 警戒区域との重ね合わせ # ============================================================================= print("\n[7] RQ3 集計: 警戒区域との重ね合わせ", flush=True) t1 = time.time() # 7-1 全県の規制区域 ∩ 警戒区域 重複面積 # 計算負荷軽減のため、規制区域全体を 1 つにまとめてから intersection zone_pref_geom = unary_union(list(zone_dissolved.geometry)) warn_pref_geom = warn_all.geometry.iloc[0] # union 後の真の規制区域面積 (= 県知事/市町長の重複を除去) total_zone_area_union = zone_pref_geom.area / 1e6 pref_share_union = total_zone_area_union / PREF_AREA_KM2 * 100 print(f" 全県規制区域 union 計算済 (面積 {total_zone_area_union:.2f} km², " f"県 {pref_share_union:.2f}% — 重複除去後の真の指定面積)") print(f" 全県警戒区域 union 計算済 (面積 {warn_pref_geom.area/1e6:.2f} km²)") # intersection overlap_geom = zone_pref_geom.intersection(warn_pref_geom) overlap_area_km2 = overlap_geom.area / 1e6 overlap_in_zone_pct = overlap_area_km2 / (zone_pref_geom.area/1e6) * 100 overlap_in_warn_pct = overlap_area_km2 / (warn_pref_geom.area/1e6) * 100 print(f" 規制区域 ∩ 警戒区域 = {overlap_area_km2:.2f} km²") print(f" H4 検証: 規制区域内に占める警戒区域割合 = {overlap_in_zone_pct:.1f}% " f"(予想 20-40%): {'YES' if 20 <= overlap_in_zone_pct <= 40 else f'観測 {overlap_in_zone_pct:.0f}%'}") # 7-2 市町別 二重リスク (規制 ∩ 警戒) 面積 double_risk_rows = [] for _, r in zone_dissolved.iterrows(): cn = r["市町名"] z_geom = r.geometry inter = z_geom.intersection(warn_pref_geom) inter_km2 = inter.area / 1e6 if not inter.is_empty else 0.0 double_risk_rows.append({ "市町名": cn, "zone_area_km2": r["zone_area_km2"], "二重リスク_km2": round(inter_km2, 3), "区域内警戒比_%": round(inter_km2 / r["zone_area_km2"] * 100, 1) if r["zone_area_km2"] > 0 else 0, }) double_risk = pd.DataFrame(double_risk_rows).sort_values("二重リスク_km2", ascending=False) print(f"\n 二重リスク 上位 5 区域:") print(double_risk.head(5)) top3_double = double_risk.head(3) top3_double_share = top3_double["二重リスク_km2"].sum() / double_risk["二重リスク_km2"].sum() print(f" H5 検証: 上位 3 シェア = {top3_double_share*100:.1f}%") print(f" 上位 3 = {top3_double['市町名'].tolist()}") # 7-3 規制区域外かつ警戒区域内 (= 警戒のみ) の面積 = 規制の盲点 warn_only_geom = warn_pref_geom.difference(zone_pref_geom) warn_only_area_km2 = warn_only_geom.area / 1e6 print(f" 警戒区域のみ (規制区域外) = {warn_only_area_km2:.2f} km² " f"= 全警戒区域の {warn_only_area_km2 / (warn_pref_geom.area/1e6) * 100:.1f}%") # 規制区域のみ (警戒区域外) = 規制があるが自然災害想定がない場所 zone_only_geom = zone_pref_geom.difference(warn_pref_geom) zone_only_area_km2 = zone_only_geom.area / 1e6 print(f" 規制区域のみ (警戒区域外) = {zone_only_area_km2:.2f} km² " f"= 全規制区域の {zone_only_area_km2 / (zone_pref_geom.area/1e6) * 100:.1f}%") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 8. CSV 出力 (再現用中間データ) # ============================================================================= print("\n[8] 中間 CSV 出力", flush=True) t1 = time.time() # 市町別 dissolve サマリ zone_summary_out = zone_dissolved[["zone_id", "市町名", "zone_area_km2", "地形タイプ"]].copy() zone_summary_out.to_csv(ASSETS / "L54_zones_summary.csv", index=False, encoding="utf-8-sig") # 市町別ランキング city_summary.to_csv(ASSETS / "L54_city_designation_ranking.csv", index=False, encoding="utf-8-sig") # polygon 個別 (全 280 行) zones_csv = pd.DataFrame(zones.drop(columns=["geometry", "rep_pt"])) zones_csv["lon_rep"] = zones.geometry.representative_point().to_crs("EPSG:4326").x zones_csv["lat_rep"] = zones.geometry.representative_point().to_crs("EPSG:4326").y zones_csv.to_csv(ASSETS / "L54_zones_polygons.csv", index=False, encoding="utf-8-sig") # RQ2: 区域密度 if L52_OK: zone_density.to_csv(ASSETS / "L54_zone_permit_density.csv", index=False, encoding="utf-8-sig") # 各 L52 盛土に in_zone フラグを付けたものを出力 g52_with_zone = pd.DataFrame(g52.drop(columns=["geometry"])) g52_with_zone.to_csv(ASSETS / "L54_l52_permits_in_zone.csv", index=False, encoding="utf-8-sig") # RQ3: 二重リスク double_risk.to_csv(ASSETS / "L54_double_risk_by_zone.csv", index=False, encoding="utf-8-sig") # 元 Shapefile (zip コピー) shutil.copy(ZONE_ZIP, ASSETS / "L54_zones_raw.zip") print(f" CSV 6 ファイル + Shapefile zip 1 ファイル") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 9. 図の生成 (9 図) # ============================================================================= print("\n[9] 図の生成 (9 図)", 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 ZONE_COLOR = "#2da44e" # 規制区域: 緑 WARN_KYUKEI_COLOR = "#cf222e" # 急傾斜: 赤 WARN_DOSEKI_COLOR = "#bf8700" # 土石流: 橙 PERMIT_IN_COLOR = "#1f6feb" # 区域内盛土: 青 PERMIT_OUT_COLOR = "#a045a0" # 区域外盛土: 紫 TYPE_COLOR = { "A_沿岸都市": "#1f6feb", "B_内陸都市": "#fa8c16", "C_山間/島嶼": "#2da44e", "D_その他": "#888", } # ---- 図 1 (RQ1): 県全域 24 区域マップ + 市町別指定面積コロプレス (2 panel) ---- fig, axes = plt.subplots(1, 2, figsize=(16, 7)) # 1-a 全規制区域マップ (指定主体別色分け) 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) # 指定主体別色分け SUB_COLOR = {"県知事": "#cf222e", "市町長": "#1f6feb"} for subj, col in SUB_COLOR.items(): sub = zones[zones["指定主体"] == subj] if len(sub) > 0: sub.plot(ax=ax, color=col, edgecolor="#333", linewidth=0.4, alpha=0.6, label=f"{subj}指定 ({len(sub)} polygon)") # 上位 5 にラベル for _, r in city_summary.head(5).iterrows(): z = zone_dissolved[zone_dissolved["市町名"] == r["市町名"]] if len(z) == 0: continue g = z.geometry.iloc[0] cx, cy = g.representative_point().x, g.representative_point().y ax.annotate(f"{r['市町名']}\n{r['指定面積_km2']:.0f} km²", (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"図 1a: 全規制区域 ({len(zones)} polygon, {total_zone_area:.0f} km²) — 県知事(赤)/市町長(青)指定別", fontsize=11) legend_handles = [Patch(facecolor=c, edgecolor="#333", label=f"{s} ({(zones['指定主体']==s).sum()} polygon)") for s, c in SUB_COLOR.items()] ax.legend(handles=legend_handles, loc="upper right", fontsize=9, title="指定主体") ax.set_xlabel("X (m, EPSG:6671)", fontsize=9) ax.set_ylabel("Y (m, EPSG:6671)", fontsize=9) # 1-b 市町別指定面積 コロプレス ax = axes[1] admin_area = admin_diss.copy() city_to_area = dict(zip(city_summary["市町名"], city_summary["指定面積_km2"])) admin_area["指定面積_km2"] = admin_area["CITY_CD"].map( {NAME_TO_CITY_CD.get(n, -1): v for n, v in city_to_area.items()}).fillna(0) admin_area[admin_area["指定面積_km2"] == 0].plot(ax=ax, facecolor="#f0f0f0", edgecolor="#888", linewidth=0.4) admin_area[admin_area["指定面積_km2"] > 0].plot(ax=ax, column="指定面積_km2", cmap="YlGn", edgecolor="#444", linewidth=0.4, legend=True, legend_kwds={"label": "指定面積 (km²)", "shrink": 0.7}) for _, r in city_summary.head(5).iterrows(): cd = NAME_TO_CITY_CD.get(r["市町名"], -1) if cd > 0: a = admin_diss[admin_diss["CITY_CD"] == cd] if len(a) == 0: continue cx, cy = a.geometry.iloc[0].centroid.x, a.geometry.iloc[0].centroid.y ax.annotate(f"{r['市町名']}\n{r['指定面積_km2']:.0f} km²", (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: 市町別指定面積 — Top 5 シェア {top5_share*100:.0f}%", fontsize=11) ax.set_xlabel("X (m, EPSG:6671)", fontsize=9) ax.set_ylabel("Y (m, EPSG:6671)", fontsize=9) fig.suptitle(f"図 1 (RQ1): 宅造規制区域の県内地理範囲 — {len(zones)} polygon (県 {pref_share_union:.1f}%, union {total_zone_area_union:.0f} km²)", fontsize=13, y=1.02) plt.tight_layout() save_fig("L54_fig1_zone_map.png") # ---- 図 2 (RQ1): 区域 polygon 規模分布 + 市町別ランキング棒 ---- fig, axes = plt.subplots(1, 2, figsize=(16, 6)) # 2-a polygon 個別の規模分布 (全 polygon) ax = axes[0] poly_a = zones["poly_area_km2"][zones["poly_area_km2"] > 0] ax.hist(poly_a, bins=np.logspace(np.log10(max(poly_a.min(), 0.0001)), np.log10(poly_a.max()), 30), color="#2da44e", edgecolor="#333", alpha=0.85) ax.set_xscale("log") ax.set_xlabel("polygon 個別面積 (km², log)", fontsize=11) ax.set_ylabel("件数", fontsize=11) ax.axvline(poly_a.median(), color="#cf222e", linestyle="--", linewidth=2, label=f"中央値 {poly_a.median():.3f} km²") ax.set_title(f"図 2a: polygon 規模分布 (n={len(zones)}) — 中央値 {poly_a.median():.3f} km², " f"max {poly_a.max():.1f} km²", fontsize=11) ax.legend(fontsize=9) ax.grid(alpha=0.3, which="both") # 2-b 市町別指定面積 横棒 (Top 15) ax = axes[1] top15 = city_summary.head(15).iloc[::-1] y = np.arange(len(top15)) bar_colors = [TYPE_COLOR.get(zone_dissolved[zone_dissolved["市町名"] == cn] ["地形タイプ"].iloc[0], "#888") if (zone_dissolved["市町名"] == cn).any() else "#888" for cn in top15["市町名"]] ax.barh(y, top15["指定面積_km2"], color=bar_colors, edgecolor="#333") ax.set_yticks(y) ax.set_yticklabels(top15["市町名"], fontsize=9) ax.set_xlabel("指定面積 (km²)", fontsize=11) ax.set_title(f"図 2b: 市町別指定面積 (Top 15) — 全県 {total_zone_area:.0f} km²", fontsize=11) ax.grid(axis="x", alpha=0.3) for i, (cn, v, sh) in enumerate(zip(top15["市町名"], top15["指定面積_km2"], top15["指定面積シェア_%"])): ax.text(v + 1, i, f"{v:.1f} km² ({sh:.1f}%)", va="center", fontsize=8) legend_handles = [Patch(facecolor=c, edgecolor="#333", label=t) for t, c in TYPE_COLOR.items() if t != "D_その他"] ax.legend(handles=legend_handles, loc="lower right", fontsize=8, title="地形タイプ") fig.suptitle(f"図 2 (RQ1): 規模分布と市町別偏在 — 県総面積 8479 km² の {pref_share:.2f}% が指定下", fontsize=13, y=1.02) plt.tight_layout() save_fig("L54_fig2_size_ranking.png") # ---- 図 3 (RQ1): 地形タイプ × 区域数/面積 + 県シェア表 ---- fig, axes = plt.subplots(1, 2, figsize=(15, 5.5)) # 3-a 地形タイプ別 面積パイ ax = axes[0] type_present = type_summary[type_summary["指定面積_km2"] > 0] labels = [f"{t}\n{int(n)} 市町 / {a:.0f} km² ({s:.0f}%)" for t, n, a, s in zip(type_present["地形タイプ"], type_present["市町数"], type_present["指定面積_km2"], type_present["シェア_%"])] colors = [TYPE_COLOR[t] for t in type_present["地形タイプ"]] ax.pie(type_present["指定面積_km2"], labels=labels, colors=colors, startangle=90, wedgeprops=dict(edgecolor="white", linewidth=2), textprops={"fontsize": 10}) ax.set_title(f"図 3a: 地形タイプ別 指定面積構成 (全 {total_zone_area:.0f} km²)", fontsize=11) # 3-b 県シェア比較 (棒) ax = axes[1] share_data = [ ("規制区域", total_zone_area, "#2da44e"), ("非指定", PREF_AREA_KM2 - total_zone_area, "#dddddd"), ] xs = np.arange(len(share_data)) ax.bar(xs, [d[1] for d in share_data], color=[d[2] for d in share_data], edgecolor="#333") ax.set_xticks(xs) ax.set_xticklabels([d[0] for d in share_data], fontsize=10) ax.set_ylabel("面積 (km²)", fontsize=11) ax.set_title(f"図 3b: 広島県総面積 {PREF_AREA_KM2:.0f} km² の内訳 — 規制区域 {pref_share:.1f}%", fontsize=11) for i, (n, v, _) in enumerate(share_data): ax.text(i, v + 80, f"{v:.0f} km²\n({v/PREF_AREA_KM2*100:.1f}%)", ha="center", va="bottom", fontsize=10) ax.grid(axis="y", alpha=0.3) fig.suptitle("図 3 (RQ1): 地形タイプ別構成と県シェア", fontsize=13, y=1.02) plt.tight_layout() save_fig("L54_fig3_type_share.png") # ---- 図 4 (RQ2): 規制区域 × L52 許可盛土 オーバーレイマップ ---- if L52_OK: fig, axes = plt.subplots(1, 2, figsize=(16, 7)) # 4-a 県全体オーバーレイ ax = axes[0] admin_diss.boundary.plot(ax=ax, color="#888", linewidth=0.4) admin_diss.plot(ax=ax, facecolor="#f8f8f5", edgecolor="none", alpha=0.7) zone_dissolved.plot(ax=ax, color=ZONE_COLOR, edgecolor="#1a4d2e", linewidth=0.6, alpha=0.45) # 区域内盛土 (青○) g52_in = g52[g52["in_zone"]] ax.scatter(g52_in.geometry.x, g52_in.geometry.y, s=35, c=PERMIT_IN_COLOR, edgecolor="black", linewidth=0.4, alpha=0.9, label=f"区域内 ({len(g52_in)})") # 区域外盛土 (紫△) g52_out = g52[~g52["in_zone"]] ax.scatter(g52_out.geometry.x, g52_out.geometry.y, s=35, c=PERMIT_OUT_COLOR, edgecolor="black", linewidth=0.4, alpha=0.9, marker="^", label=f"区域外 ({len(g52_out)})") ax.set_xlim(PXMIN, PXMAX) ax.set_ylim(PYMIN, PYMAX) ax.set_aspect("equal") ax.set_title(f"図 4a: 規制区域 (緑) × L52 許可盛土 152 件 — 区域内 {in_zone_share*100:.1f}%", fontsize=11) ax.legend(loc="upper right", fontsize=10) ax.set_xlabel("X (m)", fontsize=9) ax.set_ylabel("Y (m)", fontsize=9) # 4-b 区域内/外の規模比較 (面積分布 log) ax = axes[1] a_in = g52[g52["in_zone"]]["area_m2"].dropna() a_out = g52[~g52["in_zone"]]["area_m2"].dropna() if len(a_in) > 0 and len(a_out) > 0: bins = np.logspace(np.log10(min(a_in.min(), a_out.min())), np.log10(max(a_in.max(), a_out.max())), 25) ax.hist(a_in, bins=bins, color=PERMIT_IN_COLOR, edgecolor="#333", alpha=0.7, label=f"区域内 (n={len(a_in)})") ax.hist(a_out, bins=bins, color=PERMIT_OUT_COLOR, edgecolor="#333", alpha=0.6, label=f"区域外 (n={len(a_out)})") ax.axvline(a_in.median(), color=PERMIT_IN_COLOR, linestyle="--", linewidth=2, label=f"区域内中央値 {a_in.median():.0f} m²") ax.axvline(a_out.median(), color=PERMIT_OUT_COLOR, linestyle="--", linewidth=2, label=f"区域外中央値 {a_out.median():.0f} m²") ax.set_xscale("log") ax.set_xlabel("盛土面積 (m², log)", fontsize=11) ax.set_ylabel("件数", fontsize=11) ax.set_title(f"図 4b: 区域内/外 盛土面積分布 — 規模に差はあるか?", fontsize=11) ax.legend(fontsize=9) ax.grid(alpha=0.3, which="both") fig.suptitle("図 4 (RQ2): 規制区域 × 許可盛土オーバーレイ", fontsize=13, y=1.02) plt.tight_layout() save_fig("L54_fig4_zone_x_permit.png") # ---- 図 5 (RQ2): 区域別密度 (件/km²) ランキング + 散布図 (面積 × 件数) ---- if L52_OK: fig, axes = plt.subplots(1, 2, figsize=(15, 6)) # 5-a 区域別密度 ランキング (Top 12) ax = axes[0] dens_top = zone_density[zone_density["許可盛土件数"] > 0].head(12).iloc[::-1] y = np.arange(len(dens_top)) ax.barh(y, dens_top["密度_件_km2"], color="#1f6feb", edgecolor="#333") ax.set_yticks(y) ax.set_yticklabels(dens_top["市町名"], fontsize=9) ax.set_xlabel("許可盛土密度 (件/km²)", fontsize=11) ax.set_title(f"図 5a: 区域別 許可盛土密度 ランキング (盛土件数>0 のみ)", fontsize=11) ax.grid(axis="x", alpha=0.3) for i, (cn, d, n, a) in enumerate(zip(dens_top["市町名"], dens_top["密度_件_km2"], dens_top["許可盛土件数"], dens_top["zone_area_km2"])): ax.text(d + 0.05, i, f"{int(n)}件 / {a:.1f}km² = {d:.2f}", va="center", fontsize=8) # 5-b 散布図 区域面積 × 盛土件数 ax = axes[1] x = zone_density["zone_area_km2"] y = zone_density["許可盛土件数"] colors = [] for cn in zone_density["市町名"]: match = zone_dissolved[zone_dissolved["市町名"] == cn] if len(match) > 0: colors.append(TYPE_COLOR.get(match["地形タイプ"].iloc[0], "#888")) else: colors.append("#888") ax.scatter(x, y, s=80, c=colors, edgecolor="#333", alpha=0.8) # 中央値ライン if len(zone_density[zone_density["許可盛土件数"] > 0]) > 0: med_dens = zone_density[zone_density["許可盛土件数"] > 0]["密度_件_km2"].median() x_line = np.linspace(0.01, x.max() * 1.05, 50) ax.plot(x_line, med_dens * x_line, "k:", alpha=0.5, label=f"中央値密度 {med_dens:.2f} 件/km²") # 上位 5 ラベル for _, r in zone_density.head(5).iterrows(): ax.annotate(r["市町名"], (r["zone_area_km2"], r["許可盛土件数"]), fontsize=8, xytext=(5, 5), textcoords="offset points") ax.set_xscale("log") ax.set_yscale("symlog", linthresh=1) ax.set_xlabel("区域面積 (km², log)", fontsize=11) ax.set_ylabel("許可盛土件数 (symlog)", fontsize=11) ax.set_title(f"図 5b: 区域面積 × 盛土件数 (n={len(zone_density)} 区域)", fontsize=11) ax.legend(fontsize=9) ax.grid(alpha=0.3, which="both") legend_handles = [Patch(facecolor=c, edgecolor="#333", label=t) for t, c in TYPE_COLOR.items() if t != "D_その他"] ax.legend(handles=legend_handles, loc="lower right", fontsize=8, title="地形タイプ") fig.suptitle("図 5 (RQ2): 区域内 許可盛土密度 — 都市部偏重の検証", fontsize=13, y=1.02) plt.tight_layout() save_fig("L54_fig5_density.png") # ---- 図 6 (RQ3): 規制区域 × 警戒区域 重ね合わせ全県マップ ---- fig, ax = plt.subplots(figsize=(13, 8.5)) admin_diss.boundary.plot(ax=ax, color="#888", linewidth=0.4) admin_diss.plot(ax=ax, facecolor="#fafafa", edgecolor="#bbb", linewidth=0.3, alpha=0.6) # Layer order: warn → zone → overlap warn_kyukei.plot(ax=ax, color=WARN_KYUKEI_COLOR, alpha=0.30, edgecolor="none", label="急傾斜警戒") warn_doseki.plot(ax=ax, color=WARN_DOSEKI_COLOR, alpha=0.30, edgecolor="none", label="土石流警戒") zone_dissolved.plot(ax=ax, color=ZONE_COLOR, alpha=0.40, edgecolor="#1a4d2e", linewidth=0.5) # 上位 3 二重リスク市町ラベル for _, r in top3_double.iterrows(): z = zone_dissolved[zone_dissolved["市町名"] == r["市町名"]] if len(z) > 0: cx, cy = z.geometry.iloc[0].representative_point().x, \ z.geometry.iloc[0].representative_point().y ax.annotate(f"{r['市町名']}\n二重リスク {r['二重リスク_km2']:.1f} km²", (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"図 6 (RQ3): 規制区域 (緑) × 警戒区域 (赤=急傾斜, 橙=土石流) — " f"重複 {overlap_area_km2:.1f} km² (規制区域の {overlap_in_zone_pct:.1f}%)", fontsize=12) ax.set_xlabel("X (m, EPSG:6671)", fontsize=10) ax.set_ylabel("Y (m, EPSG:6671)", fontsize=10) legend_handles = [ Patch(facecolor=ZONE_COLOR, alpha=0.5, label="宅造規制区域"), Patch(facecolor=WARN_KYUKEI_COLOR, alpha=0.4, label="急傾斜警戒区域"), Patch(facecolor=WARN_DOSEKI_COLOR, alpha=0.4, label="土石流警戒区域"), ] ax.legend(handles=legend_handles, loc="upper right", fontsize=10, title="凡例") plt.tight_layout() save_fig("L54_fig6_overlay_warn.png") # ---- 図 7 (RQ3): 二重リスクエリア choropleth + 上位ランキング ---- fig, axes = plt.subplots(1, 2, figsize=(16, 7)) # 7-a 市町別二重リスク choropleth ax = axes[0] admin_dr = admin_diss.copy() city_to_dr = {} for _, r in double_risk.iterrows(): cn = r["市町名"] city_to_dr[cn] = r["二重リスク_km2"] admin_dr["二重リスク_km2"] = admin_dr["CITY_CD"].map( {NAME_TO_CITY_CD.get(n, -1): v for n, v in city_to_dr.items()}).fillna(0) admin_dr[admin_dr["二重リスク_km2"] == 0].plot(ax=ax, facecolor="#f0f0f0", edgecolor="#888", linewidth=0.4) admin_dr[admin_dr["二重リスク_km2"] > 0].plot(ax=ax, column="二重リスク_km2", cmap="OrRd", edgecolor="#444", linewidth=0.4, legend=True, legend_kwds={"label": "二重リスク (km²)", "shrink": 0.7}) for _, r in top3_double.iterrows(): cd = NAME_TO_CITY_CD.get(r["市町名"], -1) if cd > 0: a = admin_diss[admin_diss["CITY_CD"] == cd] if len(a) == 0: continue cx, cy = a.geometry.iloc[0].centroid.x, a.geometry.iloc[0].centroid.y ax.annotate(f"{r['市町名']}\n{r['二重リスク_km2']:.1f} km²\n({r['区域内警戒比_%']:.0f}%)", (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"図 7a: 市町別 二重リスク面積 (規制 ∩ 警戒)", fontsize=11) ax.set_xlabel("X (m)", fontsize=9) ax.set_ylabel("Y (m)", fontsize=9) # 7-b 二重リスクランキング (Top 12) ax = axes[1] dr_top = double_risk[double_risk["二重リスク_km2"] > 0].head(12).iloc[::-1] y = np.arange(len(dr_top)) ax.barh(y, dr_top["二重リスク_km2"], color="#cf222e", edgecolor="#333", label="二重リスク面積 (km²)") # 区域面積を背景バーで重ねる ax.barh(y, dr_top["zone_area_km2"], color="#cccccc", alpha=0.4, edgecolor="#888", zorder=0, label="区域面積 (km²)") ax.set_yticks(y) ax.set_yticklabels(dr_top["市町名"], fontsize=9) ax.set_xlabel("面積 (km²)", fontsize=11) ax.set_title(f"図 7b: 市町別 二重リスク + 区域面積 (Top 12) — 上位 3 シェア {top3_double_share*100:.0f}%", fontsize=11) ax.grid(axis="x", alpha=0.3) ax.legend(fontsize=9, loc="lower right") for i, (cn, dr_v, z_v, pct) in enumerate(zip(dr_top["市町名"], dr_top["二重リスク_km2"], dr_top["zone_area_km2"], dr_top["区域内警戒比_%"])): ax.text(z_v + 1, i, f"{dr_v:.1f} ({pct:.0f}%)", va="center", fontsize=8) fig.suptitle(f"図 7 (RQ3): 二重リスクエリア — 規制 ∩ 警戒 = {overlap_area_km2:.0f} km²", fontsize=13, y=1.02) plt.tight_layout() save_fig("L54_fig7_double_risk.png") # ---- 図 8 (RQ3): 規制 × 警戒 4 区分 ベン図風 + 面積構成 ---- fig, axes = plt.subplots(1, 2, figsize=(15, 6)) # 8-a 4 区分 (規制∩警戒, 規制のみ, 警戒のみ, 両方なし) パイ ax = axes[0] nonarea_km2 = PREF_AREA_KM2 - ( overlap_area_km2 + zone_only_area_km2 + warn_only_area_km2) parts = [ ("規制 ∩ 警戒\n(二重リスク)", overlap_area_km2, "#cf222e"), ("規制のみ", zone_only_area_km2, "#2da44e"), ("警戒のみ", warn_only_area_km2, "#fa8c16"), ("どちらも該当外", max(nonarea_km2, 0), "#dddddd"), ] labels = [f"{n}\n{a:.0f} km² ({a/PREF_AREA_KM2*100:.1f}%)" for n, a, _ in parts] colors = [c for _, _, c in parts] ax.pie([a for _, a, _ in parts], labels=labels, colors=colors, startangle=90, wedgeprops=dict(edgecolor="white", linewidth=2), textprops={"fontsize": 9}) ax.set_title(f"図 8a: 県全土 {PREF_AREA_KM2:.0f} km² の 4 区分構成", fontsize=11) # 8-b 4 区分の絶対面積 棒 ax = axes[1] x = np.arange(len(parts)) ax.bar(x, [a for _, a, _ in parts], color=colors, edgecolor="#333") ax.set_xticks(x) ax.set_xticklabels([n.replace("\n", " ") for n, _, _ in parts], fontsize=9, rotation=15) ax.set_ylabel("面積 (km²)", fontsize=11) ax.set_title(f"図 8b: 4 区分の絶対面積 — 警戒のみ {warn_only_area_km2:.0f} km² が制度間ギャップ", fontsize=11) ax.grid(axis="y", alpha=0.3) for i, (n, a, _) in enumerate(parts): ax.text(i, a + 30, f"{a:.0f}", ha="center", fontsize=10, weight="bold") fig.suptitle("図 8 (RQ3): 規制区域と警戒区域の 4 区分 — 制度間カバレッジ", fontsize=13, y=1.02) plt.tight_layout() save_fig("L54_fig8_four_quadrants.png") # ---- 図 9 (RQ1+RQ2+RQ3 統合): 上位 8 区域の総合プロファイル ---- fig, ax = plt.subplots(figsize=(15, 7)) # 上位 8 市町 (面積順) の規制面積・盛土密度・二重リスク%を 3 系列で表示 top8 = zone_density.merge(double_risk[["市町名", "二重リスク_km2", "区域内警戒比_%"]], on="市町名", how="left").nlargest(8, "zone_area_km2") top8 = top8.iloc[::-1].reset_index(drop=True) y = np.arange(len(top8)) w = 0.27 # 系列 1: 区域面積 (km²) ax.barh(y - w, top8["zone_area_km2"], height=w, color="#2da44e", edgecolor="#333", label="区域面積 (km²)") for i, v in enumerate(top8["zone_area_km2"]): ax.text(v + 1, i - w, f"{v:.0f}", va="center", fontsize=8) # 系列 2: 盛土密度 (件/km²) × 50 (スケール調整) SCALE_DENS = 50 ax.barh(y, top8["密度_件_km2"] * SCALE_DENS, height=w, color="#1f6feb", edgecolor="#333", label=f"盛土密度 (件/km²) × {SCALE_DENS}") for i, v in enumerate(top8["密度_件_km2"]): ax.text(v * SCALE_DENS + 1, i, f"{v:.2f}", va="center", fontsize=8) # 系列 3: 二重リスク (km²) ax.barh(y + w, top8["二重リスク_km2"], height=w, color="#cf222e", edgecolor="#333", label="二重リスク面積 (km²)") for i, v in enumerate(top8["二重リスク_km2"]): ax.text(v + 1, i + w, f"{v:.1f}", va="center", fontsize=8) ax.set_yticks(y) ax.set_yticklabels(top8["市町名"], fontsize=10) ax.set_xlabel(f"値 (区域面積・二重リスク = km², 盛土密度 × {SCALE_DENS} で同尺度化)", fontsize=11) ax.set_title(f"図 9 (統合 RQ1+RQ2+RQ3): 上位 8 区域の総合プロファイル — 面積 × 密度 × 二重リスク", fontsize=12) ax.legend(fontsize=10, loc="lower right") ax.grid(axis="x", alpha=0.3) plt.tight_layout() save_fig("L54_fig9_top8_profile.png") print(f" 9 図 生成完了 ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 10. 表の DataFrame 作成 (HTML 用) # ============================================================================= print("\n[10] 表の作成", 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 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}" else: s = str(v) cells.append(f"{s}") rows.append("" + "".join(cells) + "") return "" + head + "".join(rows) + "
" # 表 1: dataset 仕様 T_dataset = pd.DataFrame([ ("dataset_id", "1427"), ("名称", "宅地造成等工事規制区域"), ("組織", "広島県土木建築局 都市環境整備課"), ("リソース 1 (県知事)", "51124 — 県全域 Shapefile (120 polygon, 24 区域整理番, 2.3 MB, 2023-04-27 版)"), ("リソース 2-21 (市町長)", "51125〜51144 — 市町別 Shapefile 計 20 ファイル (260 polygon, 2023-09 版)"), ("リソース 22 (規制資料)", "51145 — 盛土規制法関係資料 XLSX (規制内容 + 手続き先, 23.1 KB)"), ("根拠法", "宅地造成及び特定盛土等規制法 第 12 条第 1 項 (2023-05-26 施行)"), ("旧法", "宅地造成等規制法 (1961-2023)"), ("CRS", "EPSG:2445 (JGD2000 平面直角第 III 系) → 解析時 EPSG:6671 へ変換"), ("指定日", "2023-04-27 (新法施行直前 = 旧法指定の継承)"), ("ライセンス", "クリエイティブ・コモンズ表示 4.0"), ("取得日", "2026-05-09"), ], columns=["項目", "値"]) # 表 2: 全体サマリ T_overall = pd.DataFrame([ ("総 polygon 数", f"{len(zones)} polygon (県知事 {n_pref_polys} + 市町長 {n_city_polys}, 市町別 dissolve 後 {len(zone_dissolved)} 市町)"), ("対象市町数", f"{zone_dissolved['市町名'].nunique()} 市町"), ("指定面積 (合算)", f"{total_zone_area_sum:.2f} km² (県+市町長別々の合算, 重複あり)"), ("指定面積 (union)", f"{total_zone_area_union:.2f} km² (重複除去後の真の指定面積, 県 {PREF_AREA_KM2:.0f} km² の {pref_share_union:.2f}%)"), ("polygon 規模 中央値 / 最大", f"{zones['poly_area_km2'].median():.3f} km² / {zones['poly_area_km2'].max():.2f} km²"), ("Top 1 市町", f"{city_summary.iloc[0]['市町名']} ({city_summary.iloc[0]['指定面積_km2']:.1f} km²)"), ("Top 5 シェア", f"{top5_share*100:.1f}%"), ("地形タイプ別 沿岸都市", f"{type_summary[type_summary['地形タイプ']=='A_沿岸都市']['指定面積_km2'].sum():.0f} km² " f"({type_summary[type_summary['地形タイプ']=='A_沿岸都市']['シェア_%'].sum():.0f}%)"), ("L52 許可盛土 総件数", f"{len(g52) if L52_OK else '—'}"), ("L52 区域内立地", f"{n_in_zone}/{len(g52)} 件 ({in_zone_share*100:.1f}%)" if L52_OK else "—"), ("L52 区域外立地", f"{n_out_zone} 件 ({100*n_out_zone/len(g52):.1f}%)" if L52_OK else "—"), ("規制区域 ∩ 警戒区域", f"{overlap_area_km2:.2f} km² " f"(規制内 {overlap_in_zone_pct:.1f}% / 警戒内 {overlap_in_warn_pct:.1f}%)"), ("二重リスク Top 3 シェア", f"{top3_double_share*100:.0f}% ({', '.join(top3_double['市町名'].tolist())})"), ("警戒のみ (規制外)", f"{warn_only_area_km2:.0f} km² " f"(全警戒の {warn_only_area_km2/(warn_pref_geom.area/1e6)*100:.1f}%)"), ], columns=["指標", "値"]) # 表 3: 市町別規制区域一覧 (順位 = 面積順) T_zones = zone_dissolved[["市町名", "地形タイプ", "zone_area_km2"]].copy() T_zones["順位"] = T_zones["zone_area_km2"].rank(ascending=False, method="min").astype(int) T_zones["県シェア_%"] = (T_zones["zone_area_km2"] / PREF_AREA_KM2 * 100).round(2) T_zones["指定面積シェア_%"] = (T_zones["zone_area_km2"] / T_zones["zone_area_km2"].sum() * 100).round(2) T_zones = T_zones.sort_values("zone_area_km2", ascending=False).reset_index(drop=True) # 県知事/市町長 polygon 数も追加 T_zones["県知事_n"] = T_zones["市町名"].map( zones[zones["指定主体"] == "県知事"].groupby("市町名").size().to_dict()).fillna(0).astype(int) T_zones["市町長_n"] = T_zones["市町名"].map( zones[zones["指定主体"] == "市町長"].groupby("市町名").size().to_dict()).fillna(0).astype(int) T_zones = T_zones[["順位", "市町名", "地形タイプ", "zone_area_km2", "指定面積シェア_%", "県シェア_%", "県知事_n", "市町長_n"]] T_zones.columns = ["順位", "市町名", "地形タイプ", "面積_km2", "指定面積シェア_%", "県総面積シェア_%", "県知事n", "市町長n"] # 表 4: 地形タイプ集計 T_type = type_summary.copy() T_type["地形タイプ"] = T_type["地形タイプ"].replace({ "A_沿岸都市": "A 沿岸都市", "B_内陸都市": "B 内陸都市", "C_山間/島嶼": "C 山間/島嶼", "D_その他": "D その他", }) # 表 5: RQ2 市町別密度 (sorted by density desc) if L52_OK: T_density = zone_density[["市町名", "zone_area_km2", "許可盛土件数", "密度_件_km2"]].copy() T_density.columns = ["市町名", "区域面積_km2", "L52許可盛土件数", "密度_件km2"] T_density = T_density.sort_values("密度_件km2", ascending=False).reset_index(drop=True) T_density["順位"] = np.arange(1, len(T_density) + 1) T_density = T_density[["順位", "市町名", "区域面積_km2", "L52許可盛土件数", "密度_件km2"]] else: T_density = pd.DataFrame() # 表 6: RQ2 区域内/外 規模比較 if L52_OK: in_stats = g52[g52["in_zone"]][["height_max_m", "area_m2", "fill_volume_m3"]] out_stats = g52[~g52["in_zone"]][["height_max_m", "area_m2", "fill_volume_m3"]] T_in_out = pd.DataFrame([ {"指標": "件数 (件)", "区域内": len(in_stats), "区域外": len(out_stats), "差 (区域内 - 区域外)": len(in_stats) - len(out_stats)}, {"指標": "高さ 中央値 (m)", "区域内": round(in_stats["height_max_m"].median(), 2), "区域外": round(out_stats["height_max_m"].median(), 2), "差 (区域内 - 区域外)": round(in_stats["height_max_m"].median() - out_stats["height_max_m"].median(), 2)}, {"指標": "面積 中央値 (m²)", "区域内": int(in_stats["area_m2"].median()), "区域外": int(out_stats["area_m2"].median()), "差 (区域内 - 区域外)": int(in_stats["area_m2"].median() - out_stats["area_m2"].median())}, {"指標": "盛土量 中央値 (m³)", "区域内": int(in_stats["fill_volume_m3"].median()), "区域外": int(out_stats["fill_volume_m3"].median()), "差 (区域内 - 区域外)": int(in_stats["fill_volume_m3"].median() - out_stats["fill_volume_m3"].median())}, {"指標": "面積 最大 (m²)", "区域内": int(in_stats["area_m2"].max()), "区域外": int(out_stats["area_m2"].max()), "差 (区域内 - 区域外)": int(in_stats["area_m2"].max() - out_stats["area_m2"].max())}, ]) else: T_in_out = pd.DataFrame() # 表 7: RQ3 二重リスクランキング T_dr = double_risk.copy() T_dr = T_dr[["市町名", "zone_area_km2", "二重リスク_km2", "区域内警戒比_%"]] T_dr.columns = ["市町名", "区域面積_km2", "二重リスク_km2", "区域内警戒区域比率_%"] T_dr["順位"] = T_dr["二重リスク_km2"].rank(ascending=False, method="min").astype(int) T_dr = T_dr.sort_values("二重リスク_km2", ascending=False).head(15).reset_index(drop=True) T_dr = T_dr[["順位", "市町名", "区域面積_km2", "二重リスク_km2", "区域内警戒区域比率_%"]] # 表 8: RQ3 4 区分カバレッジ T_quad = pd.DataFrame([ {"区分": "1. 規制 ∩ 警戒 (二重リスク)", "面積_km2": round(overlap_area_km2, 1), "県シェア_%": round(overlap_area_km2 / PREF_AREA_KM2 * 100, 2), "意味": "宅地工事に許可必要 + 自然災害想定 = 二重防御エリア"}, {"区分": "2. 規制のみ (警戒区域外)", "面積_km2": round(zone_only_area_km2, 1), "県シェア_%": round(zone_only_area_km2 / PREF_AREA_KM2 * 100, 2), "意味": "宅地工事の許可制のみ。市街地中心の平地等"}, {"区分": "3. 警戒のみ (規制区域外)", "面積_km2": round(warn_only_area_km2, 1), "県シェア_%": round(warn_only_area_km2 / PREF_AREA_KM2 * 100, 2), "意味": "災害想定のみ。山間・郊外で宅造規制なし → 制度間ギャップ"}, {"区分": "4. どちらも該当外", "面積_km2": round(max(PREF_AREA_KM2 - overlap_area_km2 - zone_only_area_km2 - warn_only_area_km2, 0), 1), "県シェア_%": round(max(PREF_AREA_KM2 - overlap_area_km2 - zone_only_area_km2 - warn_only_area_km2, 0) / PREF_AREA_KM2 * 100, 2), "意味": "両制度の対象外。森林・農地・水域など"}, ]) # 表 9: 二段階監督モデル拡張版 (許可制 L52 + 届出制 L53 + 区域 L54) T_three_systems = pd.DataFrame([ {"側面": "根拠条文", "規制区域 (L54)": "法 12 条 1 項 (本記事)", "許可制 (L52)": "法 12 条 1 項 / 30 条 1 項", "届出制 (L53)": "法 21 条 1 項 / 40 条 1 項"}, {"側面": "対象", "規制区域 (L54)": "地理範囲 (区域指定 = ポリゴン)", "許可制 (L52)": "個別工事 (規模超 = 許可必要)", "届出制 (L53)": "個別工事 (規模以下 = 届出必要)"}, {"側面": "本データ単位", "規制区域 (L54)": f"{len(zones)} polygon ({n_cities_designated} 市町, union {total_zone_area_union:.0f} km²)", "許可制 (L52)": f"{len(df_l52) if L52_OK else '—'} 件 (個別工事)", "届出制 (L53)": "284 件 (個別工事)"}, {"側面": "全県カバレッジ", "規制区域 (L54)": f"union {total_zone_area_union:.0f} km² ({pref_share_union:.1f}% of 県)", "許可制 (L52)": f"L52 区域内立地 {in_zone_share*100:.0f}%" if L52_OK else "—", "届出制 (L53)": "区域内 + 特定盛土区域 (40 条)"}, {"側面": "指定の意味", "規制区域 (L54)": "区域内全工事に許可/届出義務発動の前提条件", "許可制 (L52)": "事前審査 + 許可条件 (重い手続)", "届出制 (L53)": "事前通知のみ (軽い手続)"}, {"側面": "L54 との関係", "規制区域 (L54)": "—", "許可制 (L52)": f"L52 のうち {in_zone_share*100:.0f}% が L54 区域内" if L52_OK else "—", "届出制 (L53)": "L53 もすべて L54 区域内 (区域指定が前提)"}, ]) # 表 10: 仮説検証 g_share = in_zone_share if L52_OK else 0 T_hypo = pd.DataFrame([ {"仮説": "H1 (都市集中, 上位 5 シェア ≥ 50%, RQ1)", "予想": "Top 5 市町で全県の 50% 以上", "観測": f"Top 5 = {', '.join(city_summary.head(5)['市町名'].tolist())}, " f"シェア {top5_share*100:.1f}%", "判定": "支持" if top5_share >= 0.50 else "部分支持"}, {"仮説": "H2 (区域内立地率 ≥ 70%, RQ2)", "予想": "L52 152 件のうち区域内立地が 70% 以上", "観測": f"L52 区域内 {n_in_zone}/{len(g52) if L52_OK else 0} 件 = {in_zone_share*100:.1f}%" if L52_OK else "L52 比較不可", "判定": "支持" if (L52_OK and in_zone_share >= 0.70) else "部分支持" if L52_OK else "検証不可"}, {"仮説": "H3 (区域内密度の市町偏在 ≥ 2 桁, RQ2)", "予想": "区域別 件/km² が 2 桁以上の幅で偏在", "観測": (f"密度 max {zone_density['密度_件_km2'].max():.2f} / " f"min(>0) {zone_density.loc[zone_density['許可盛土件数']>0, '密度_件_km2'].min():.4f} " f"= {10**density_log_diff:.0f} 倍" if L52_OK else "—"), "判定": ("支持" if (L52_OK and density_log_diff >= 2) else "部分支持" if L52_OK else "検証不可")}, {"仮説": "H4 (規制 ∩ 警戒 = 20-40%, RQ3)", "予想": "規制区域内に占める警戒区域の割合が 20-40%", "観測": f"規制 ∩ 警戒 = {overlap_area_km2:.1f} km², " f"規制内 {overlap_in_zone_pct:.1f}%", "判定": "支持" if 20 <= overlap_in_zone_pct <= 40 else f"観測 {overlap_in_zone_pct:.0f}% (予想外)"}, {"仮説": "H5 (二重リスク Top 3 集中, RQ3)", "予想": "上位 3 で二重リスク全体の半数以上", "観測": f"Top 3 = {', '.join(top3_double['市町名'].tolist())}, " f"シェア {top3_double_share*100:.1f}%", "判定": "支持" if top3_double_share >= 0.50 else f"部分支持 ({top3_double_share*100:.0f}%)"}, ]) # 表 11: 警戒のみエリア 上位市町 (制度間ギャップ) # 各市町の admin polygon に対して warn ∩ city - zone ∩ city = 警戒のみ print(f" warn_only ∩ city 計算 (制度間ギャップ識別)... ", flush=True) warn_only_geom_full = warn_pref_geom.difference(zone_pref_geom) warn_only_rows = [] admin_diss_p = admin_diss[["CITY_CD", "geometry"]].copy() admin_diss_p["city_warn_only_km2"] = 0.0 for idx, r in admin_diss_p.iterrows(): cd = int(r["CITY_CD"]) cn = CITY_NAME.get(cd, f"CITY_CD={cd}") inter = r.geometry.intersection(warn_only_geom_full) a_km2 = inter.area / 1e6 if not inter.is_empty else 0 if a_km2 > 0: warn_only_rows.append({ "市町名": cn, "警戒のみ_km2": round(a_km2, 2), }) T_warn_only = pd.DataFrame(warn_only_rows).sort_values("警戒のみ_km2", ascending=False).head(15) T_warn_only["順位"] = np.arange(1, len(T_warn_only) + 1) T_warn_only = T_warn_only[["順位", "市町名", "警戒のみ_km2"]] # 表 12: polygon 詳細サンプル (Top 10 大 polygon) T_poly_sample = zones.copy() T_poly_sample = T_poly_sample.nlargest(10, "poly_area_km2")[ ["指定主体", "指定主体_code", "市町名", "poly_area_km2"]].reset_index(drop=True) T_poly_sample["順位"] = np.arange(1, len(T_poly_sample) + 1) T_poly_sample = T_poly_sample[["順位", "指定主体", "指定主体_code", "市町名", "poly_area_km2"]] T_poly_sample.columns = ["順位", "指定主体", "団体コード", "市町名", "polygon面積_km2"] print(f" 表 12 種作成 ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 11. HTML 構築 # ============================================================================= print("\n[11] HTML 構築", flush=True) t11 = time.time() # Sec 1: 学習目標と問い sec1 = f"""

本記事は DoBoX のシリーズ「宅地造成等工事規制区域」 1 件 (dataset_id = 1427) を 単独で取り上げ、 広島県内の宅造規制区域 {len(zones)} polygon (県知事指定 {n_pref_polys} + 市町長指定 {n_city_polys}) / {n_cities_designated} 市町 / union 面積 {total_zone_area_union:.0f} km² ({pref_share_union:.1f}% of 県土) を 3 つの独立した研究角度 (RQ1 / RQ2 / RQ3) で並列に分析する。 本データは盛土規制法 (2023-05-26 施行) 第 12 条第 1 項に基づく区域指定であり、 旧「宅地造成等規制法」 (1961-2023) からの指定継承である。

本データは DoBoX で21 ファイル構成: 県知事指定の県全域版 1 ファイル (51124, 2023-04-27 版) + 各市町長指定の市町別 20 ファイル (51125〜51144, 2023-09 版)。 盛土規制法 12 条 1 項は都道府県知事が指定するが、 広島市・呉市・尾道市・福山市・東広島市等の政令指定都市/中核市は 法令により独自に指定権限を持つため別ファイルで配信されている。 本記事は 21 ファイルすべてを統合して扱う。

本記事の独自ポイントは、L52 (許可盛土 152 件) / L53 (届出盛土 284 件) と並列で読むことで、 盛土規制法の「区域 → 工事監督」 という 2 階建て構造を初めてデータから可視化する点。 区域指定 (本記事 L54) は地理範囲を法的に固定し、 許可制 (L52) と届出制 (L53) はその区域内で個別工事を監督する。 区域なしには工事監督が成立しないため、L54 は L52・L53 の地理的根拠である。

独自用語の定義

研究の問い (3 RQ)

仮説 H1〜H5

  1. H1 (都市集中, RQ1): Top 5 市町で全県指定面積の 50% 以上を占める。 山間部市町は指定なしか小面積。
  2. H2 (区域内立地率高, RQ2): L52 許可盛土 152 件のうち、本規制区域内立地は 70% 以上。 区域外立地は特定盛土区域 (法 30 条) 側の許可で少数派。
  3. H3 (区域内密度の市町偏在, RQ2): 本区域内の許可盛土密度 (件/km²) は 市町別に 2 桁以上の幅で偏在 (都市部高 / 山間市町低)。
  4. H4 (警戒区域との部分重複, RQ3): 規制 ∩ 警戒の重複面積比は 20-40%。 規制は市街地中心を、警戒は急傾斜・渓流を捉えるため空間目的が異なる。
  5. H5 (二重リスクエリアの市町偏在, RQ3): 二重リスク (規制 ∩ 警戒) の市町別ランキングは 「規制面積トップ × 警戒面積トップ」 のクロス積。広島市・呉市・尾道市が上位 3。

到達点

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

  1. 1 つの「区域 Shapefile」 (120 polygon) から、 盛土規制法の地理的射程 ({pref_share:.1f}%) を読み取り、 L52 (許可) / L53 (届出) との地理的依存関係を可視化する方法を習得する。
  2. 区域内盛土密度という独自指標で、件数の絶対値ではなく 面積正規化した強度で市町を比較する眼を獲得する。
  3. 規制 × 警戒の 4 区分カバレッジ分析で、制度間ギャップ (警戒のみ {warn_only_area_km2:.0f} km²) という 行政設計上の重要課題を定量化する手法を体感する。
""" # Sec 2: 使用データ sec2 = f"""

DoBoX のシリーズ「宅地造成等工事規制区域」 1 件のみを単独で扱う。 リソースは Shapefile 21 件 (県知事 1 + 市町長 20) + XLSX 1 件の計 22 ファイル構造:

{df_to_html(T_dataset)}

データの構造

CRS と単位の取扱

本記事は dataset 1427 を単独で扱う Format B 記事。 L52 (dataset 1429 = 許可盛土) との比較は本区域指定が L52 許可立地の前提であるため RQ2 で使用する。土砂災害警戒区域 (sediment_shp) は L10/L11 で扱った既存データを RQ3 で空間結合に利用するが、警戒区域そのものの分析は別レッスンに委ねる。

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

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

生データ (DoBoX 直リンク)

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

図 (PNG 9 枚)

再現スクリプト

個別取得 (PowerShell):

cd "2026 DoBoX 教材"
iwr "https://hiroshima-dobox.jp/resource_download/51124" -OutFile "data/extras/L54_residential_filling_zone/340006_residential_land_construction_regulation_area_shp_20230427.zip"
iwr "https://hiroshima-dobox.jp/resource_download/51145" -OutFile "data/extras/L54_residential_filling_zone/340006_earth_filling_regulation_related_documents.xlsx"
py -X utf8 lessons\\L54_residential_filling_zone.py

注: RQ2 (L52 比較) には L52 の処理済 CSV (lessons/assets/L52_permits_processed.csv) が必要。 未生成の場合、本スクリプトは RQ2 部分を「データ無し」 表示にしてフォールバック実行する。

""" # Sec 4 (RQ1) sec4_intro = f"""

狙い (RQ1)

宅造規制区域 ({n_cities_designated} 市町指定) は広島県の宅造工事監督の地理的射程を法的に固定する。 本 RQ1 では (1) 県全体の指定面積、(2) 市町別指定面積ランキング、 (3) 地形タイプ別構成、(4) polygon 個別の規模分布の4 軸で構造を読む。 これは「県は何を、どの面積で、どんな地形タイプで、どの市町で監督対象としているか」 を 全方位で記述する基礎研究。

手法 (21 ファイル読込 → 統合 → CRS 変換 → admin sjoin → 集計)

入出力 Before/After 具体例 (1 ファイル追跡)

市町別 Shapefile 51140 (= 345458_..., 大崎上島町相当) の処理例:

{df_to_html(pd.DataFrame([ {"段階": "RAW (zip)", "対象": "市町別 zip 1 個", "geometry": "Shapefile 集合", "CRS": "EPSG:2445", "polygon数": "—"}, {"段階": "1. unzip + read_file", "対象": "GeoDataFrame", "geometry": "Polygon × N", "CRS": "EPSG:2445", "polygon数": "N"}, {"段階": "2. to_crs(6671)", "対象": "GeoDataFrame", "geometry": "EPSG:6671 へ再投影", "CRS": "EPSG:6671", "polygon数": "N"}, {"段階": "3. 指定主体列付与", "対象": "GeoDataFrame", "geometry": "+ 指定主体='市町長'", "CRS": "EPSG:6671", "polygon数": "N"}, {"段階": "4. concat (20 ファイル)", "対象": "全市町長 GDF", "geometry": "Polygon × 260", "CRS": "EPSG:6671", "polygon数": "260"}, {"段階": "5. concat (県知事 + 市町長)", "対象": "統合 zones", "geometry": f"Polygon × {len(zones)}", "CRS": "EPSG:6671", "polygon数": f"{len(zones)}"}, {"段階": "6. unary_union", "対象": "重複除去後", "geometry": "1 MultiPolygon", "CRS": "EPSG:6671", "polygon数": f"面積 {total_zone_area_union:.0f} km²"}, ]))} """ sec4_code = code(r''' # 21 ファイル統合: 県知事 + 20 市町長指定 Shapefile import geopandas as gpd, pandas as pd, zipfile from pathlib import Path from shapely.ops import unary_union # 1. 県知事指定 (51124, 県全域版) を読込 zone_pref = gpd.read_file( "data/extras/L54_residential_filling_zone/shp/" "340006_residential_land_construction_regulation_area_shp_20230427.shp" ).to_crs("EPSG:6671") pref_named = zone_pref[zone_pref["区域整理番"].notna()].copy() pref_named["指定主体"] = "県知事" # 2. 20 市町長指定 Shapefile を順次読込 (51125〜51144) city_dir = Path("data/extras/L54_residential_filling_zone/all_cities/") city_polys = [] for zip_path in sorted(city_dir.glob("*.zip")): code5 = zip_path.name[:5] # 例: 342033 → "34203" (市町コード) extract = city_dir / f"tmp_{code5}" extract.mkdir(exist_ok=True) if not any(extract.glob("*.shp")): with zipfile.ZipFile(zip_path) as zf: zf.extractall(extract) g = gpd.read_file(next(extract.glob("*.shp"))).to_crs("EPSG:6671") g["指定主体"] = "市町長" g["指定主体_code"] = code5 city_polys.append(g) city_named = pd.concat(city_polys, ignore_index=True) # 3. 統合 (縦結合) common_cols = ["指定主体", "geometry"] zones = gpd.GeoDataFrame( pd.concat([pref_named[common_cols], city_named[common_cols]], ignore_index=True), crs="EPSG:6671") zones["poly_area_km2"] = zones.geometry.area / 1e6 zones["poly_id"] = range(len(zones)) # 4. 各 polygon の代表点で行政区域と sjoin → 市町名解決 zones["rep_pt"] = zones.geometry.representative_point() zone_pts = gpd.GeoDataFrame(zones[["poly_id"]], geometry=zones["rep_pt"], crs="EPSG:6671") join = gpd.sjoin(zone_pts, admin_diss[["CITY_CD", "geometry"]], how="left", predicate="within").drop_duplicates("poly_id") zones["市町名"] = join.set_index("poly_id")["CITY_CD"].reindex(zones["poly_id"]) \ .map(CITY_NAME).fillna("不明").values # 5. 市町別 dissolve + union 面積 zone_dissolved = zones.dissolve(by="市町名", as_index=False) zone_dissolved["zone_area_km2"] = zone_dissolved.geometry.area / 1e6 zone_pref_union = unary_union(list(zone_dissolved.geometry)) print(f"統合 polygon: {len(zones)}, union 面積: {zone_pref_union.area/1e6:.1f} km²") ''') sec4_fig1_text = """

図 1: 全規制区域マップ (指定主体別) + 市町別指定面積コロプレス (2 panel)

なぜこの図か: 学習者がまず「規制区域はどこにあるか」 を一目で把握するため。 左の指定主体別色分けで、県知事指定 (赤) と市町長指定 (青) の制度的色彩を視覚化。 右のコロプレスで市町別指定面積の濃淡を見る。 両方を並べることで「指定権限主体の違い」 と「量の偏り」 を別軸で読む基本訓練ができる。

""" top1_name = city_summary.iloc[0]["市町名"] top1_area = city_summary.iloc[0]["指定面積_km2"] sec4_fig1_read = f"""

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

""" sec4_fig2_text = """

図 2: polygon 規模分布 (log) + 市町別指定面積ランキング (Top 15)

なぜこの図か: 統合 {len(zones)} polygon は市町・指定主体別に大きさが大きく異なる。 左のヒストグラム (log) で polygon の規模分布、右の市町ランキングで 県全体での偏在度を読む。両方を並べることで「polygon 数の多さ ≠ 面積の多さ」 を体感できる。

""" poly_a_med = zones["poly_area_km2"].median() poly_a_max = zones["poly_area_km2"].max() poly_a_min = zones["poly_area_km2"].min() sec4_fig2_read = f"""

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

""" sec4_fig3_text = """

図 3: 地形タイプ別 指定面積構成 (パイ) + 県総面積シェア比較 (棒)

なぜこの図か: 規制区域指定の地理的偏りを 2 つのスケールで見る。 左パイは「規制区域の中での 3 タイプの内訳」、右棒は「県全体に対する規制区域の量」。 両方を並べることで、規制区域の質 (どこに偏るか) と量 (県のうちどれだけか) を分けて理解できる。

""" a_share = type_summary[type_summary['地形タイプ']=='A_沿岸都市']['シェア_%'].sum() b_share = type_summary[type_summary['地形タイプ']=='B_内陸都市']['シェア_%'].sum() c_share = type_summary[type_summary['地形タイプ']=='C_山間/島嶼']['シェア_%'].sum() d_share = type_summary[type_summary['地形タイプ']=='D_その他']['シェア_%'].sum() if 'D_その他' in type_summary['地形タイプ'].values else 0 sec4_fig3_read = f"""

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

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

狙い (RQ2)

宅造規制区域は L52 許可盛土の地理的根拠。 区域指定が無ければ 12 条 1 項の許可義務は発動しないため、 本来「L52 152 件は本区域内に立地する」 のが制度設計上の前提である。 本 RQ2 では (1) 区域内立地率の検証、(2) 区域別の盛土密度 (件/km²)、 (3) 区域外立地の地理特性 (= 特定盛土区域 30 条側) を読み解く。

「件数」 は絶対値だが、密度 (件/km²) は面積で正規化された強度指標。 都市部の小区域に多くの許可盛土が集中するか、山間市町の大区域に薄く分散するか — 密度ランキングはそれを直接比較できる本研究独自の指標。

手法 (within 判定 + 区域別 sjoin + 密度計算)

入力: L52 152 点 + 規制区域 union ({total_zone_area_union:.0f} km²)。
出力: 各 L52 点に in_zone フラグ (True/False)、市町別に許可盛土件数・密度。
限界: within 判定は polygon 境界線上の点を「外」 と判定する場合がある (誤差は数 m レベルで実用上問題ない)。本データは座標精度 ~10 m なので影響は無視可能。

""" sec5_code = code(r''' # RQ2: 規制区域 × L52 許可盛土の空間解析 from shapely.ops import unary_union import pandas as pd, geopandas as gpd # 1. L52 既処理 CSV から GeoDataFrame df_l52 = pd.read_csv("lessons/assets/L52_permits_processed.csv", encoding="utf-8-sig") g52 = gpd.GeoDataFrame( df_l52, geometry=gpd.points_from_xy(df_l52["lon"], df_l52["lat"]), crs="EPSG:4326").to_crs("EPSG:6671") # 2. 規制区域 (県知事 + 市町長指定 全 polygon) を 1 つに統合し、各 L52 点が内/外を判定 zone_union = unary_union(list(zone_dissolved.geometry)) g52["in_zone"] = g52.geometry.within(zone_union) n_in = g52["in_zone"].sum() print(f"L52 区域内立地: {n_in}/{len(g52)} = {n_in/len(g52):.1%}") # 3. 区域別盛土件数 (sjoin) g52_zone = gpd.sjoin(g52, zone_dissolved[["zone_id", "geometry"]], how="left", predicate="within") zone_count = g52_zone.groupby("zone_id").size().reset_index(name="許可盛土件数") # 4. 密度 = 件数 / 区域面積 (km²) zone_density = zone_dissolved[["zone_id", "市町名", "zone_area_km2"]].merge( zone_count, on="zone_id", how="left").fillna(0) zone_density["密度_件_km2"] = zone_density["許可盛土件数"] / zone_density["zone_area_km2"] print(zone_density.sort_values("密度_件_km2", ascending=False).head(5)) ''') sec5_fig4_text = """

図 4: 規制区域 × L52 許可盛土オーバーレイ + 規模分布比較

なぜこの図か: 「区域指定 = 許可盛土集中」 仮説 (H2) を地図上で直接検証する。 左で規制区域 (緑) と許可盛土 (区域内=青○ / 区域外=紫△) を重ね、 右で区域内/外の盛土面積分布を log で並列比較する。 規模に差があれば「区域指定の有無で許可される盛土の質も変わる」 ことを示す。

""" if L52_OK: in_med_a = g52[g52["in_zone"]]["area_m2"].median() out_med_a = g52[~g52["in_zone"]]["area_m2"].median() sec5_fig4_read = f"""

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

""" else: sec5_fig4_read = "

L52 比較不可 (assets/L52_permits_processed.csv 未生成)

" sec5_fig5_text = """

図 5: 区域別 許可盛土密度ランキング + 区域面積 × 件数 散布

なぜこの図か: 「件数」 は絶対値だが「密度 (件/km²)」 は面積正規化指標。 小区域に少数立地でも密度は高くなり、大区域に多数立地でも密度は低くなる。 ランキングは件数だけでは見えない区域の盛土圧を可視化する。 右の散布図 (面積 × 件数) で、密度 = 傾きの解釈もできる。

""" if L52_OK: dens_max = zone_density["密度_件_km2"].max() dens_top1 = zone_density[zone_density["密度_件_km2"] == dens_max].iloc[0] sec5_fig5_read = f"""

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

""" else: sec5_fig5_read = "

L52 比較不可

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

狙い (RQ3)

宅造規制区域 (= 行政が宅造工事を許可制で監督) と 土砂災害警戒区域 (= 自然災害想定エリア) は異なる行政目的で指定される 2 つの polygon 集合。本来は「規制 ∩ 警戒 = 二重防御」の領域があるはずだが、 指定タイミングの差や空間目的の違いで、想定外の重複構造を持つ可能性がある。

本 RQ3 では (1) 全県の規制 ∩ 警戒重複面積、(2) 市町別二重リスクランキング、 (3) 4 区分カバレッジ (規制∩警戒 / 規制のみ / 警戒のみ / どちらも該当外) を読み解く。 特に「警戒のみ (= 規制区域外の警戒区域)」 は制度間ギャップとして重要。

手法 (geometry intersection + 4 区分集計)

入力: 規制区域 polygon (union {total_zone_area_union:.0f} km²) + 急傾斜 30K + 土石流 13K polygon。
出力: 全県重複面積、区域別二重リスク、4 区分面積。
限界: 警戒区域の面積は dissolve 前後で変わらない (= 同種類の警戒同士の重複は無視できる)。 ただし急傾斜 ∩ 土石流 の重複 (異種警戒の重複) は本研究では合計に含めて扱う。

""" sec6_code = code(r''' # RQ3: 規制区域 × 警戒区域の重ね合わせ from shapely.ops import unary_union import geopandas as gpd # 1. 警戒区域 (急傾斜 + 土石流) を 1 つに union warn_all_geom = unary_union(list(warn_kyukei.geometry) + list(warn_doseki.geometry)) # 2. 規制区域も union zone_pref_geom = unary_union(list(zone_dissolved.geometry)) # 3. intersection で重複面積 overlap_geom = zone_pref_geom.intersection(warn_all_geom) overlap_km2 = overlap_geom.area / 1e6 # 4. 4 区分集計 zone_only_geom = zone_pref_geom.difference(warn_all_geom) warn_only_geom = warn_all_geom.difference(zone_pref_geom) zone_only_km2 = zone_only_geom.area / 1e6 warn_only_km2 = warn_only_geom.area / 1e6 print(f"規制 ∩ 警戒 = {overlap_km2:.1f} km²") print(f"規制のみ = {zone_only_km2:.1f} km²") print(f"警戒のみ = {warn_only_km2:.1f} km² ← 制度間ギャップ") # 5. 区域別二重リスク double_risk = [] for _, r in zone_dissolved.iterrows(): inter = r.geometry.intersection(warn_all_geom) double_risk.append({ "市町名": r["市町名"], "二重リスク_km2": inter.area / 1e6 if not inter.is_empty else 0, }) ''') sec6_fig6_text = """

図 6: 規制区域 × 警戒区域 重ね合わせ全県マップ

なぜこの図か: 全県で規制区域 (緑) と警戒区域 (赤=急傾斜・橙=土石流) が どう重なるかを 1 図で俯瞰する。レイヤ順は warn → zone で、 警戒区域の上に規制区域が透過で乗る。 重なる場所は 2 色が混ざって見え、二重リスクエリアが視覚的に確認できる。

""" sec6_fig6_read = f"""

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

""" sec6_fig7_text = """

図 7: 市町別 二重リスク面積 choropleth + ランキング

なぜこの図か: 二重リスクを市町単位で正規化したランキングと地図。 左コロプレスで地理分布、右ランキングで絶対値+区域内比率を読む。 これにより「どの市町が最も二重リスクを抱えるか」 と「区域面積に占める警戒比率」 を同時に把握できる。

""" sec6_fig7_read = f"""

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

""" sec6_fig8_text = """

図 8: 県全土 4 区分カバレッジ — 規制 / 警戒の組合せパターン

なぜこの図か: 県全土 (8479 km²) を「規制と警戒の有無」 で 4 区分し、 それぞれの面積を可視化することで制度間カバレッジを診断する。 特に「警戒のみ」の面積は規制が及ばない災害リスクエリア = 制度間ギャップ。 左パイで構成比、右棒で絶対面積を見る。

""" sec6_fig8_read = f"""

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

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

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

" + df_to_html(T_hypo) + "

3 RQ × 3 結論

" "" "

3 制度モデル — 規制区域 (L54) / 許可制 (L52) / 届出制 (L53) の関係表

" + df_to_html(T_three_systems) + "

この表から読み取れること: 盛土規制法の3 階層構造が見える: " "(1) 区域 (L54) = 地理範囲を法的に固定 → " "(2) 区域内で許可閾値超 = 12/30 条許可 (L52) → " "(3) 区域内で許可閾値以下 = 21/40 条届出 (L53)。" "L54 (区域) は L52 / L53 の地理的前提であり、L52 / L53 は L54 の具体実装。" "3 つを並列で読むことで、盛土規制法の制度設計が初めて立体的に理解できる。

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

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

発展課題 1 (RQ1 由来): 区域指定の都市計画整合性検証

発展課題 2 (RQ2 由来): 盛土密度と人口密度・地価の連関

発展課題 3 (RQ3 由来): 警戒のみエリアの宅地造成リスク定量化

""" # Combine sections sections = [ ("学習目標と問い", sec1), ("使用データ", sec2), ("ダウンロード", sec3), ("【RQ1】 地理範囲と市町指定状況の研究 — 21 ファイル・約 2,075 km² (県の 24%) の沿岸都市偏重", sec4_intro + "

実装コード (Shapefile 読込 + dissolve + admin sjoin)

" + sec4_code + sec4_fig1_text + figure("assets/L54_fig1_zone_map.png", "図 1 (RQ1): 全規制区域マップ (指定主体別) + 市町別指定面積コロプレス") + sec4_fig1_read + sec4_fig2_text + figure("assets/L54_fig2_size_ranking.png", "図 2 (RQ1): polygon 規模分布 (log) + 市町別指定面積ランキング (Top 15)") + sec4_fig2_read + sec4_fig3_text + figure("assets/L54_fig3_type_share.png", "図 3 (RQ1): 地形タイプ別 指定面積構成 + 県シェア比較") + sec4_fig3_read + "

表: 全体サマリ (規制区域 + L52 比較 + 警戒区域)

" + df_to_html(T_overall) + "

この表から読み取れること: " f"{len(zones)} polygon / {n_cities_designated} 市町 / union 面積 {total_zone_area_union:.0f} km² / " f"県シェア {pref_share_union:.1f}% という基礎数値、" f"L52 区域内立地 {in_zone_share*100:.0f}%、二重リスク {overlap_area_km2:.0f} km² 等の" "横断指標が一覧で確認できる。RQ1〜RQ3 の核心が 13 行に集約された統合サマリ。

" + "

表: 市町別規制区域一覧 (面積順)

" + df_to_html(T_zones) + "

この表から読み取れること: 各市町の地形タイプ・面積・指定主体構成を一覧で確認できる。" "Top 5 はすべて沿岸都市で、上位 5 で県の {top5_share*100:.0f}% を占める。" "下位の小区域 (1 km² 未満) は山間集落の小規模市街地で、規制区域指定の" "カバレッジの広さを担保している。

".replace("{top5_share*100:.0f}", f"{top5_share*100:.0f}") + "

表: 地形タイプ別集計 (3 タイプ)

" + df_to_html(T_type) + "

この表から読み取れること: 沿岸都市が" f" {a_share:.0f}% を占め、内陸 {b_share:.0f}% / 山間・島嶼 {c_share:.0f}% " "と続く。市町数で見ると沿岸都市は少ないが面積が大きく、山間・島嶼は市町数が多いが面積は小さい =\n" "都市規模と市町数の逆対応関係が現れる。

" ), ("【RQ2】 許可盛土発生密度との関係 — L52 152 件の区域内立地検証", sec5_intro + "

実装コード (within 判定 + 区域別 sjoin + 密度計算)

" + sec5_code + sec5_fig4_text + figure("assets/L54_fig4_zone_x_permit.png", "図 4 (RQ2): 規制区域 × L52 許可盛土オーバーレイ + 区域内/外規模比較") + sec5_fig4_read + sec5_fig5_text + figure("assets/L54_fig5_density.png", "図 5 (RQ2): 区域別 許可盛土密度ランキング + 区域面積 × 件数 散布") + sec5_fig5_read + ("

表: 区域別 許可盛土密度

" + df_to_html(T_density) + "

この表から読み取れること: 密度トップは" f" {T_density.iloc[0]['市町名']} = {T_density.iloc[0]['密度_件km2']:.2f} 件/km²。" "区域面積が小さく盛土件数が多い区域が密度上位を占める = " "都市部の集中開発エリアの特徴。

" if L52_OK else "") + ("

表: 区域内/外 L52 規模比較

" + df_to_html(T_in_out) + "

この表から読み取れること: 区域内盛土と区域外盛土の規模中央値の差で、" "「区域指定の有無で許可される盛土の質が変わるか」 を直接確認できる。" f"件数差 {len(g52[g52['in_zone']]) - len(g52[~g52['in_zone']])} 件 " "(区域内 - 区域外) は区域指定の地理的偏在の現れ。

" if L52_OK else "") ), ("【RQ3】 土砂災害警戒区域との重ね合わせ — 二重リスクと制度間ギャップ", sec6_intro + "

実装コード (geometry intersection + 4 区分集計)

" + sec6_code + sec6_fig6_text + figure("assets/L54_fig6_overlay_warn.png", "図 6 (RQ3): 規制区域 × 警戒区域 重ね合わせ全県マップ") + sec6_fig6_read + sec6_fig7_text + figure("assets/L54_fig7_double_risk.png", "図 7 (RQ3): 市町別 二重リスク面積 choropleth + ランキング") + sec6_fig7_read + sec6_fig8_text + figure("assets/L54_fig8_four_quadrants.png", "図 8 (RQ3): 県全土 4 区分カバレッジ (規制 / 警戒の組合せパターン)") + sec6_fig8_read + "

表: 区域別 二重リスクランキング (Top 15)

" + df_to_html(T_dr) + "

この表から読み取れること: 二重リスク面積トップは" f" {T_dr.iloc[0]['市町名']} ({T_dr.iloc[0]['二重リスク_km2']:.1f} km², " f"区域内警戒比率 {T_dr.iloc[0]['区域内警戒区域比率_%']:.0f}%)。" "区域面積が大きいだけでなく、区域内に占める警戒区域の比率も高い = " "地形リスクと宅地需要が重なる宿命的な都市構造。

" + "

表: 4 区分カバレッジ (県 8479 km²)

" + df_to_html(T_quad) + "

この表から読み取れること: 4 区分のうち「警戒のみ」が" f" {warn_only_area_km2:.0f} km² ({warn_only_area_km2/PREF_AREA_KM2*100:.1f}%) と最大級で、" "これが制度間ギャップ。宅造規制が及ばない災害想定エリアで、" "山間・郊外集落周辺に分布する。該当外が県土の大半を占めるが、" "これは森林・農地・水域で住民密度が低く宅地監督対象外という当然の結果。

" + "

表: 警戒のみエリア 上位市町 (制度間ギャップ)

" + df_to_html(T_warn_only) + "

この表から読み取れること: 警戒のみエリアの上位市町は" f" {T_warn_only.iloc[0]['市町名']} ({T_warn_only.iloc[0]['警戒のみ_km2']:.1f} km²) を筆頭に" "山間市町が多く名を連ねる → " "山間市町は警戒区域の絶対面積が大きいが宅造規制区域が小さいため、" "制度間ギャップが集中して発生。これらの市町は L53 届出制度 (40 条 = 特定盛土区域) で" "間接的に盛土工事を監督する設計だが、規制区域 (12 条) のような厳格な許可制ではない。

" + "

表: polygon サンプル (Top 10 大 polygon)

" + df_to_html(T_poly_sample) + "

この表から読み取れること: 個別 polygon の最大は" f" {T_poly_sample.iloc[0]['市町名']} の {T_poly_sample.iloc[0]['polygon面積_km2']:.1f} km²。" "1 区域内に複数大規模 polygon を含む場合、その区域は" "連続した広い市街地を持つ大都市と推察される。

" ), ("仮説検証総合", sec7), ("発展課題", sec8), ] html = render_lesson( num=54, title="宅地造成等工事規制区域 単独 3 研究例分析 — 盛土規制法 12 条 1 項の地理的射程を 21 ファイルから読み解く", tags=["L54", "盛土規制法", "宅地造成等工事規制区域", "RQ×3", "Format B", "12条1項", "二重リスク", "L52連携", "L53連携", "制度間ギャップ"], time="40 分", level="中級+", data_label="DoBoX dataset 1427 (Shapefile 120 polygon + XLSX 規制資料) + L52/L10/L11 連携", sections=sections, script_filename="L54_residential_filling_zone.py", ) OUT_HTML = LESSONS / "L54_residential_filling_zone.html" OUT_HTML.write_text(html, encoding="utf-8") print(f" HTML written: {OUT_HTML}") print(f" ({time.time()-t11:.1f}s)", flush=True) print(f"\n=== L54 DONE in {time.time()-t_all:.1f}s ===")