"""L14 3D 都市モデル 7 市町比較 — 建物高度の地理学 研究の核: 広島県の 3D 都市モデル (PLATEAU 規格) は、現状 7 市町 (広島市・福山市・呉市・三次市・府中市・海田町・竹原市) で公開されている。 だが CityGML 本体は各市町 100 MB〜1 GB 級で、ハンズオン教材で扱うには重すぎる。 そこで本記事は **建物利用現況 (都市計画基礎調査)** の建物毎 GeoJSON を使う。 これは 3D 都市モデルから建物 footprint と「階数 (KAISUU)」「軒高 (TATE_H)」だけを 抜き出した軽量レイヤで、「建物高度の地理学」を扱うには十分。 カバー宣言: 本記事は **3D 都市モデル 7 市町 = 7 dataset_id** を統合分析する研究である。 - #1316 広島市 / #1302 府中市 / #1187 呉市 / #1188 福山市 - #1189 海田町 / #1443 竹原市 / #1444 三次市 CityGML は重量級なため Plan B として、同 7 市町の建物利用現況 GeoJSON (#1469 / #1483 / #1323 / #1327 / #512 / #71 / #1474) を採用。 これらは PLATEAU の建物 footprint + 属性 (階数・軒高) と内容が一致する軽量版。 **対応表 7/7 件 論理カバー**。 要件 S 対応: - GeoJSON は cold で 計 ~40 秒。一度 parquet 化すれば 2 回目以降 ~3 秒 - サンプリング (city ごと最大 80,000 棟) で描画コスト削減 - 浸水域 sjoin は代表点 (representative_point) で高速近似 - キャッシュ: data/extras/_l14_cache/{city}.parquet """ from __future__ import annotations import sys, time, json, os 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.colors import LinearSegmentedColormap import geopandas as gpd import shapely plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False t0 = time.time() print("=== L14 3D 都市モデル 7 市町比較 — 建物高度の地理学 ===") # ============================================================================= # 0. 定数 # ============================================================================= TARGET_CRS = "EPSG:6671" # JGD2011 平面直角 III (広島県、m単位) # 7 市町の定義 (slug, 日本語名, 建物利用現況 GeoJSON ファイル名 prefix, 3D都市モデル ID, 建物利用現況 ID) CITY_DEFS = [ ("hiroshima", "広島市", "hiroshima_34100", 1316, 1469), ("fukuyama", "福山市", "fukuyama_34207", 1188, 1483), ("kure", "呉市", "kure_C0401_34202", 1187, 1323), ("miyoshi", "三次市", "miyoshi_C0401_34209", 1444, 1327), ("fuchu", "府中市", "fuchu_C0401_34208", 1302, 512), ("kaita", "海田町", "kaita_C0401_34304", 1189, 71), ("takehara", "竹原市", "takehara_34203", 1443, 1474), ] # 7 市町の参考統計 (S17 と同じ値, 公的統計の概数) CITY_REF = { "広島市": {"pop": 1189000, "area_km2": 906.7, "category": "政令市"}, "福山市": {"pop": 459000, "area_km2": 518.1, "category": "中核市"}, "呉市": {"pop": 210000, "area_km2": 352.8, "category": "中核市相当"}, "三次市": {"pop": 50000, "area_km2": 778.2, "category": "中山間・市"}, "府中市": {"pop": 37000, "area_km2": 195.8, "category": "中山間・市"}, "海田町": {"pop": 30000, "area_km2": 13.8, "category": "都市近郊・町"}, "竹原市": {"pop": 23000, "area_km2": 118.2, "category": "瀬戸内・市"}, } # 階数バンド (1階=平屋, 2階, 3-4階=低中層, 5-9階=中層, 10階以上=高層) def kaisuu_band(k): if k is None or pd.isna(k) or k <= 0: return "0_その他" # 0階 = 不明・倉庫等 k = int(k) if k == 1: return "1_平屋" if k == 2: return "2_2階" if k <= 4: return "3_3-4階" if k <= 9: return "4_5-9階" return "5_10階以上" KAISUU_BAND_ORDER = ["1_平屋", "2_2階", "3_3-4階", "4_5-9階", "5_10階以上", "0_その他"] KAISUU_BAND_LABEL = { "1_平屋": "平屋(1階)", "2_2階": "2階", "3_3-4階": "3-4階(低中層)", "4_5-9階": "5-9階(中層)", "5_10階以上": "10階以上(高層)", "0_その他": "0階・不明", } KAISUU_BAND_COLOR = { "1_平屋": "#fef3c7", "2_2階": "#fbbf24", "3_3-4階": "#f59e0b", "4_5-9階": "#dc2626", "5_10階以上": "#7c1d6f", "0_その他": "#999999", } # TATE_YO (建物用途) コード分類 def yoto_band(c): if c is None or pd.isna(c): return "9_不明" c = int(c) if 401 <= c <= 404: return "1_業務" if 411 <= c <= 415: return "2_商業" if 421 <= c <= 431: return "3_娯楽サービス" if c == 441: return "4_文教" if 451 <= c <= 454: return "5_住宅" if c == 461: return "6_工業" if 471 <= c <= 499: return "7_倉庫運輸他" return "9_不明" YOTO_BAND_LABEL = { "1_業務": "業務", "2_商業": "商業", "3_娯楽サービス": "娯楽サービス", "4_文教": "文教", "5_住宅": "住宅", "6_工業": "工業", "7_倉庫運輸他": "倉庫運輸他", "9_不明": "不明", } # 浸水深ランク DEPTH_LABEL = { 10: "0.0〜0.5m", 20: "0.5〜1.0m", 30: "1.0〜2.0m", 40: "2.0〜3.0m", 50: "3.0〜5.0m", 60: "5.0〜10.0m", 70: "10.0〜20.0m", 80: "20m以上", } # 「2階以上避難可能」判定: 浸水深 >= 3m (rank>=50) で 平屋(1階) は 命の危険 = 致命率 DEEP_DEADLY_RANK = 50 # 3m 以上 # ============================================================================= # 1. キャッシュ機構: 大型 GeoJSON → 軽量 parquet # ============================================================================= BLD_DIR = ROOT / "data" / "extras" / "bld_data" CACHE_DIR = ROOT / "data" / "extras" / "_l14_cache" CACHE_DIR.mkdir(parents=True, exist_ok=True) KEEP_COLS = ["KAISUU", "TATE_H", "TATE_YO", "YOTO_CD", "KENCK_A"] def load_city_buildings(slug, jp_name, prefix): """建物利用現況 GeoJSON を読込 (parquet キャッシュ経由)""" cache_path = CACHE_DIR / f"{slug}.parquet" if cache_path.exists() and cache_path.stat().st_size > 1000: gdf = gpd.read_parquet(cache_path) # CRS 復元 (parquet は CRS を持つが念のため) if gdf.crs is None or "Compound" in str(gdf.crs): gdf = gdf.set_crs(TARGET_CRS, allow_override=True) return gdf # GeoJSON 探索 candidates = list(BLD_DIR.glob(f"{prefix}*.geojson")) if not candidates: # 別パターン candidates = list(BLD_DIR.glob(f"{slug}_*.geojson")) if not candidates: print(f" [skip] {jp_name}: GeoJSON が見つからない (探索 prefix={prefix})") return None src = candidates[0] keep = [c for c in KEEP_COLS] # KENCK_A は欠如する都市あり try: gdf = gpd.read_file(src, columns=keep + ["geometry"]) except Exception: # KENCK_A など欠如時は段階的に gdf = gpd.read_file(src) present = [c for c in keep if c in gdf.columns] gdf = gdf[present + ["geometry"]] # CRS 強制 if gdf.crs is None or "Compound" in str(gdf.crs): gdf = gdf.set_crs(TARGET_CRS, allow_override=True) elif gdf.crs.to_epsg() != 6671: gdf = gdf.to_crs(TARGET_CRS) # 重い 3D を 2D に強制 (軒高 z 座標を落とす) gdf["geometry"] = gpd.GeoSeries(shapely.force_2d(gdf.geometry.values), crs=TARGET_CRS) # parquet 保存 gdf.to_parquet(cache_path) return gdf print("\n--- 1. 7 市町 建物データ読込 (parquet キャッシュ経由) ---") city_data = {} city_load_log = [] for slug, jp, prefix, ds3d, dsbld in CITY_DEFS: t1 = time.time() cache_path = CACHE_DIR / f"{slug}.parquet" cached = cache_path.exists() g = load_city_buildings(slug, jp, prefix) dt = time.time() - t1 if g is None: city_load_log.append({"city": jp, "slug": slug, "status": "FAIL", "rows": 0, "sec": dt, "cached": cached}) continue g["city"] = jp g["slug"] = slug g["kaisuu_band"] = g["KAISUU"].apply(kaisuu_band) g["yoto_band"] = g["TATE_YO"].apply(yoto_band) g["yoto_label"] = g["yoto_band"].map(YOTO_BAND_LABEL) city_data[slug] = g city_load_log.append({"city": jp, "slug": slug, "status": "OK", "rows": len(g), "sec": dt, "cached": cached}) print(f" {jp:<6s} ({slug}): {len(g):,} 棟 / {dt:.1f}s " f"({'cached' if cached else 'first-load'})") if not city_data: print("[!] 全市町で読込失敗。スケルトン HTML を出力して終了") sys.exit(1) # 統合 (city 列付き) print("\n--- 2. 7 市町統合 ---") COMMON_COLS = ["city", "slug", "KAISUU", "TATE_H", "TATE_YO", "YOTO_CD", "kaisuu_band", "yoto_band", "yoto_label", "geometry"] parts = [] for slug, g in city_data.items(): sub = g[[c for c in COMMON_COLS if c in g.columns]].copy() parts.append(sub) all_b = gpd.GeoDataFrame(pd.concat(parts, ignore_index=True), crs=TARGET_CRS) print(f" 統合 {len(all_b):,} 棟") n_bld_total = len(all_b) # ============================================================================= # 3. 市町別 基本統計 # ============================================================================= print("\n--- 3. 市町別 基本統計 ---") city_stats = [] for slug, g in city_data.items(): jp = g["city"].iloc[0] n = len(g) kaisuu_valid = g[g["KAISUU"] > 0]["KAISUU"] h_valid = g[g["TATE_H"] > 0]["TATE_H"] if "TATE_H" in g.columns else pd.Series(dtype=float) ref = CITY_REF.get(jp, {}) city_stats.append({ "city": jp, "slug": slug, "pop": ref.get("pop", 0), "area_km2": ref.get("area_km2", 0), "category": ref.get("category", ""), "n_bld": n, "bld_per_km2": n / max(ref.get("area_km2", 1), 1), "mean_kaisuu": float(kaisuu_valid.mean()) if len(kaisuu_valid) else 0, "median_kaisuu": float(kaisuu_valid.median()) if len(kaisuu_valid) else 0, "max_kaisuu": int(kaisuu_valid.max()) if len(kaisuu_valid) else 0, "p95_kaisuu": float(kaisuu_valid.quantile(0.95)) if len(kaisuu_valid) else 0, "mean_h": float(h_valid.mean()) if len(h_valid) else 0, "max_h": float(h_valid.max()) if len(h_valid) else 0, "p95_h": float(h_valid.quantile(0.95)) if len(h_valid) else 0, "n_high": int((g["KAISUU"] >= 5).sum()), # 5階以上 "n_tower": int((g["KAISUU"] >= 10).sum()), # 10階以上 "n_flat": int((g["KAISUU"] == 1).sum()), # 平屋 "rate_high_pct": float((g["KAISUU"] >= 5).sum() / max(n, 1) * 100), "rate_flat_pct": float((g["KAISUU"] == 1).sum() / max(n, 1) * 100), }) city_stats_df = pd.DataFrame(city_stats).sort_values("pop", ascending=False).reset_index(drop=True) city_stats_df.to_csv(ASSETS / "L14_city_stats.csv", index=False, encoding="utf-8-sig") print(city_stats_df[["city", "n_bld", "mean_kaisuu", "max_kaisuu", "rate_high_pct", "rate_flat_pct"]].to_string(index=False)) # ============================================================================= # 4. 階数バンド × 市町 クロス集計 # ============================================================================= print("\n--- 4. 階数バンド × 市町 クロス ---") band_cross = (all_b.groupby(["city", "kaisuu_band"], observed=False).size() .unstack(fill_value=0).reindex(columns=KAISUU_BAND_ORDER, fill_value=0)) # 順序: 人口順 band_cross = band_cross.reindex([s["city"] for s in city_stats]) band_cross_pct = band_cross.div(band_cross.sum(axis=1), axis=0) * 100 band_cross.to_csv(ASSETS / "L14_band_cross_count.csv", encoding="utf-8-sig") band_cross_pct.to_csv(ASSETS / "L14_band_cross_pct.csv", encoding="utf-8-sig") print(band_cross_pct.round(1).to_string()) # ============================================================================= # 5. 用途バンド × 市町 クロス集計 # ============================================================================= print("\n--- 5. 用途バンド × 市町 クロス ---") yoto_cross = (all_b.groupby(["city", "yoto_band"], observed=False).size() .unstack(fill_value=0)) yoto_order = ["1_業務", "2_商業", "3_娯楽サービス", "4_文教", "5_住宅", "6_工業", "7_倉庫運輸他", "9_不明"] yoto_cross = yoto_cross.reindex(columns=[c for c in yoto_order if c in yoto_cross.columns], fill_value=0) yoto_cross = yoto_cross.reindex([s["city"] for s in city_stats]) yoto_cross_pct = yoto_cross.div(yoto_cross.sum(axis=1), axis=0) * 100 yoto_cross.columns = [YOTO_BAND_LABEL.get(c, c) for c in yoto_cross.columns] yoto_cross_pct.columns = yoto_cross.columns yoto_cross.to_csv(ASSETS / "L14_yoto_cross.csv", encoding="utf-8-sig") yoto_cross_pct.to_csv(ASSETS / "L14_yoto_cross_pct.csv", encoding="utf-8-sig") # ============================================================================= # 6. 浸水域 × 建物 sjoin (代表点で高速近似) # ============================================================================= print("\n--- 6. 浸水域 × 建物 sjoin ---") DEPTH_SHP = ROOT / "data" / "extras" / "flood_shp" / "shinsui_souteisaidai" / "shinsui_souteisaidai.shp" FLOOD_CACHE = CACHE_DIR / "flood_6671.parquet" flood_join_summary = pd.DataFrame() deadly_summary = pd.DataFrame() buildings_in_flood = None if DEPTH_SHP.exists(): t6 = time.time() if FLOOD_CACHE.exists(): flood = gpd.read_parquet(FLOOD_CACHE) print(f" 浸水深 cache 読込: {time.time()-t6:.1f}s, {len(flood)} polygons") else: flood = gpd.read_file(DEPTH_SHP, columns=["rank"]) flood = flood.to_crs(TARGET_CRS) flood["geometry"] = gpd.GeoSeries(shapely.force_2d(flood.geometry.values), crs=TARGET_CRS) flood["depth_rank"] = flood["rank"] flood = flood[["depth_rank", "geometry"]] flood.to_parquet(FLOOD_CACHE) print(f" 浸水深 Shapefile 読込→cache: {time.time()-t6:.1f}s, {len(flood)} polygons") # 描画専用 simplify cache (sjoin は元解像度で行う) # 浸水ポリゴンは MultiPolygon で 4M 頂点級。simplify(preserve_topology=False) で大幅圧縮 FLOOD_PLOT_CACHE = CACHE_DIR / "flood_6671_plot_v3.parquet" if FLOOD_PLOT_CACHE.exists(): flood_for_plot = gpd.read_parquet(FLOOD_PLOT_CACHE) print(f" 浸水深 plot cache 読込: {len(flood_for_plot)} polys") else: tt = time.time() flood_for_plot = flood.copy() # preserve_topology=False で位相崩れを許容して高速化、頂点数 1/100 級まで縮小可 flood_for_plot["geometry"] = flood_for_plot.geometry.simplify(200, preserve_topology=False) # 空・無効を除去 flood_for_plot = flood_for_plot[~flood_for_plot.geometry.is_empty].copy() flood_for_plot = flood_for_plot[flood_for_plot.geometry.is_valid].copy() flood_for_plot.to_parquet(FLOOD_PLOT_CACHE) print(f" 浸水深 plot cache 生成 (simplify200 no-topo): {time.time()-tt:.1f}s") # 建物の代表点 (geometry が大量で重い場合は centroid でも可) t6 = time.time() bld_pts = all_b[["city", "slug", "KAISUU", "TATE_YO", "kaisuu_band", "yoto_band", "yoto_label"]].copy() bld_pts["bld_id"] = np.arange(len(bld_pts)) bld_pts = gpd.GeoDataFrame(bld_pts, geometry=all_b.geometry.representative_point(), crs=TARGET_CRS) print(f" 建物代表点 生成: {time.time()-t6:.1f}s") t6 = time.time() # sjoin は重複を生む (浸水ポリゴンが層状に重なるため、1 棟が複数行に拡大される) bld_join_raw = gpd.sjoin(bld_pts, flood, how="left", predicate="within") print(f" sjoin (建物 in 浸水深): {time.time()-t6:.1f}s, raw rows={len(bld_join_raw):,}") # 1 棟 = 1 行に集約: 最大深さランクを採用 (最悪ケースで判定) t6 = time.time() max_rank_per_bld = (bld_join_raw.groupby("bld_id")["depth_rank"].max() .reset_index().rename(columns={"depth_rank": "depth_rank_max"})) bld_in_flood = bld_pts.merge(max_rank_per_bld, on="bld_id", how="left") bld_in_flood["in_flood"] = bld_in_flood["depth_rank_max"].notna() bld_in_flood["depth_rank"] = bld_in_flood["depth_rank_max"].fillna(0).astype(int) bld_in_flood["deep_3m"] = bld_in_flood["depth_rank"] >= DEEP_DEADLY_RANK bld_in_flood["deadly_flat"] = bld_in_flood["deep_3m"] & (bld_in_flood["KAISUU"] == 1) bld_in_flood["safe_2plus"] = bld_in_flood["in_flood"] & (bld_in_flood["KAISUU"] >= 2) print(f" 重複除去 (max rank per bld): {time.time()-t6:.1f}s, dedup rows={len(bld_in_flood):,}") # 市町別集計 n_total_per_city = bld_in_flood.groupby("city").size() n_in_flood = bld_in_flood[bld_in_flood["in_flood"]].groupby("city").size() n_deep_3m = bld_in_flood[bld_in_flood["deep_3m"]].groupby("city").size() n_deadly_flat = bld_in_flood[bld_in_flood["deadly_flat"]].groupby("city").size() n_safe_2plus = bld_in_flood[bld_in_flood["safe_2plus"]].groupby("city").size() flood_join_summary = pd.DataFrame({ "n_total": n_total_per_city, "n_in_flood": n_in_flood, "n_deep_3m": n_deep_3m, "n_deadly_flat": n_deadly_flat, "n_safe_2plus": n_safe_2plus, }).fillna(0).astype(int).reset_index() # 比率列 flood_join_summary["pct_in_flood"] = flood_join_summary["n_in_flood"] / flood_join_summary["n_total"] * 100 flood_join_summary["pct_deadly_flat"] = ( flood_join_summary["n_deadly_flat"] / flood_join_summary["n_in_flood"].replace(0, np.nan) * 100 ).fillna(0) flood_join_summary["pct_safe_2plus_of_flooded"] = ( flood_join_summary["n_safe_2plus"] / flood_join_summary["n_in_flood"].replace(0, np.nan) * 100 ).fillna(0) flood_join_summary = flood_join_summary.sort_values("n_in_flood", ascending=False) flood_join_summary.to_csv(ASSETS / "L14_flood_overlay.csv", index=False, encoding="utf-8-sig") print(flood_join_summary.to_string(index=False)) # 致命率の市町別比較 deadly_summary = flood_join_summary[ ["city", "n_in_flood", "n_deadly_flat", "pct_deadly_flat"] ].copy().sort_values("pct_deadly_flat", ascending=False) deadly_summary.to_csv(ASSETS / "L14_deadly_summary.csv", index=False, encoding="utf-8-sig") buildings_in_flood = bld_in_flood # 図描画用に保持 else: print(f" [warn] 浸水深 Shapefile なし: {DEPTH_SHP}") # ============================================================================= # 7. 市町×階数 で TATE_H の分布要約 # ============================================================================= print("\n--- 7. 高度分布 要約 ---") height_dist = [] for slug, g in city_data.items(): jp = g["city"].iloc[0] h = g[g["TATE_H"] > 0]["TATE_H"] if len(h) > 0: for q, name in [(0.5, "median_h"), (0.75, "p75_h"), (0.95, "p95_h"), (0.99, "p99_h")]: pass height_dist.append({ "city": jp, "n_h_known": len(h), "min_h": float(h.min()), "median_h": float(h.median()), "p75_h": float(h.quantile(0.75)), "p95_h": float(h.quantile(0.95)), "p99_h": float(h.quantile(0.99)), "max_h": float(h.max()), }) else: height_dist.append({ "city": jp, "n_h_known": 0, "min_h": 0, "median_h": 0, "p75_h": 0, "p95_h": 0, "p99_h": 0, "max_h": 0, }) height_dist_df = pd.DataFrame(height_dist).sort_values("median_h", ascending=False) height_dist_df.to_csv(ASSETS / "L14_height_dist.csv", index=False, encoding="utf-8-sig") print(height_dist_df.to_string(index=False)) # ============================================================================= # 図1: 主題図 — 7 市町の建物分布 small multiples (KAISUU 色分け) # ============================================================================= print("\n--- 図1 small multiples 主題図 ---") tx = time.time() # 描画用サンプリング (各市町 max 10000 点。代表点で点描画 → 高速) SAMPLE_PER_CITY = 15000 fig, axes = plt.subplots(2, 4, figsize=(17, 9)) axes_f = axes.flatten() ordered_cities = [s["city"] for s in city_stats] # 人口順 for i, jp in enumerate(ordered_cities): ax = axes_f[i] sub = all_b[all_b["city"] == jp] if len(sub) > SAMPLE_PER_CITY: sub = sub.sample(SAMPLE_PER_CITY, random_state=42) # bbox 計算 minx, miny, maxx, maxy = sub.total_bounds pad = max(maxx - minx, maxy - miny) * 0.04 # 代表点に変換 (footprint 描画は重いので) sub_pts = gpd.GeoDataFrame(sub.drop(columns="geometry"), geometry=sub.geometry.representative_point(), crs=TARGET_CRS) # 階数バンド別に描画 (低層→高層を上重ねで) for band in KAISUU_BAND_ORDER: sb = sub_pts[sub_pts["kaisuu_band"] == band] if len(sb) == 0: continue sb.plot(ax=ax, color=KAISUU_BAND_COLOR[band], alpha=0.75, markersize=1.0, linewidth=0) n_total = (all_b["city"] == jp).sum() high_n = ((all_b["city"] == jp) & (all_b["KAISUU"] >= 5)).sum() ax.set_title(f"{jp}\n{n_total:,} 棟 / 5階以上 {high_n} 棟", fontsize=10) ax.set_xlim(minx - pad, maxx + pad) ax.set_ylim(miny - pad, maxy + pad) ax.set_aspect("equal") ax.set_xticks([]); ax.set_yticks([]) # 8 番目: 凡例 axes_f[-1].axis("off") legend_handles = [Patch(facecolor=KAISUU_BAND_COLOR[b], label=KAISUU_BAND_LABEL[b]) for b in KAISUU_BAND_ORDER] axes_f[-1].legend(handles=legend_handles, loc="center", fontsize=11, title="階数バンド", framealpha=1.0, edgecolor="#999") plt.suptitle("7 市町 建物分布 small multiples (色: 階数バンド, 人口降順 1→7)", fontsize=14, y=1.0) plt.tight_layout() plt.savefig(ASSETS / "L14_small_multiples_kaisuu.png", dpi=110) plt.close() print(f" 完了: {time.time()-tx:.1f}s") # ============================================================================= # 図2: 階数バンド構成比 (積み上げ棒) + 棒数値棒 # ============================================================================= print("\n--- 図2 階数バンド構成比 ---") tx = time.time() fig, axes = plt.subplots(1, 2, figsize=(15, 6)) ax = axes[0] bottom = np.zeros(len(band_cross_pct)) xs = np.arange(len(band_cross_pct)) for band in KAISUU_BAND_ORDER: if band not in band_cross_pct.columns: continue vals = band_cross_pct[band].values ax.bar(xs, vals, bottom=bottom, color=KAISUU_BAND_COLOR[band], edgecolor="white", linewidth=0.6, label=KAISUU_BAND_LABEL[band]) # 5% 以上のセグメントに数値 for i, v in enumerate(vals): if v >= 4: ax.text(xs[i], bottom[i] + v / 2, f"{v:.0f}%", ha="center", va="center", fontsize=9, color="#222") bottom += vals ax.set_xticks(xs) ax.set_xticklabels(band_cross_pct.index, rotation=20) ax.set_ylabel("構成比 (%)") ax.set_ylim(0, 105) ax.set_title("市町別 階数バンド構成比 (人口降順)") ax.legend(loc="lower right", fontsize=8, framealpha=0.95) ax.grid(axis="y", alpha=0.3) # 右: 5階以上率 / 平屋率 ペア棒 ax = axes[1] order = list(band_cross_pct.index) high_pct = [city_stats_df[city_stats_df["city"] == c]["rate_high_pct"].iloc[0] for c in order] flat_pct = [city_stats_df[city_stats_df["city"] == c]["rate_flat_pct"].iloc[0] for c in order] xpos = np.arange(len(order)) w = 0.4 b1 = ax.bar(xpos - w/2, high_pct, w, color="#7c1d6f", label="5階以上率") b2 = ax.bar(xpos + w/2, flat_pct, w, color="#fef3c7", edgecolor="#999", label="平屋率") for bs, vs in [(b1, high_pct), (b2, flat_pct)]: for r, v in zip(bs, vs): ax.text(r.get_x() + r.get_width()/2, r.get_height() + 0.3, f"{v:.1f}%", ha="center", fontsize=8) ax.set_xticks(xpos) ax.set_xticklabels(order, rotation=20) ax.set_ylabel("比率 (%)") ax.set_title("5階以上率 vs 平屋率") ax.legend(loc="upper right", fontsize=9) ax.grid(axis="y", alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "L14_band_compare.png", dpi=110) plt.close() print(f" 完了: {time.time()-tx:.1f}s") # ============================================================================= # 図3: 軒高 (TATE_H) 分布 ヒストグラム (city 別) # ============================================================================= print("\n--- 図3 軒高ヒストグラム ---") tx = time.time() fig, axes = plt.subplots(2, 4, figsize=(17, 8)) axes_f = axes.flatten() bins = np.arange(0, 80, 2) for i, jp in enumerate(ordered_cities): ax = axes_f[i] g = city_data[ [c for c in city_data if city_data[c]["city"].iloc[0] == jp][0] ] h = g[g["TATE_H"] > 0]["TATE_H"] if len(h) == 0: ax.text(0.5, 0.5, f"{jp}\n軒高データなし\n(KAISUU のみ参照)", ha="center", va="center", transform=ax.transAxes, fontsize=11, bbox=dict(boxstyle="round,pad=0.5", facecolor="#fff5f5", edgecolor="#cf222e")) ax.set_xticks([]); ax.set_yticks([]) ax.set_title(f"{jp}", fontsize=10) continue ax.hist(h, bins=bins, color="#0969da", alpha=0.85, edgecolor="white", linewidth=0.4) p95 = h.quantile(0.95) p99 = h.quantile(0.99) ax.axvline(p95, color="#dc2626", linestyle="--", linewidth=1.2, label=f"P95={p95:.0f}m") ax.axvline(p99, color="#7c1d6f", linestyle="--", linewidth=1.2, label=f"P99={p99:.0f}m") ax.set_title(f"{jp}\nN={len(h):,} / max={h.max():.0f}m", fontsize=10) ax.set_xlim(0, 60) ax.set_xlabel("軒高 TATE_H (m)") ax.set_ylabel("棟数 (log)") ax.set_yscale("log") ax.legend(fontsize=7) ax.grid(alpha=0.3) axes_f[-1].axis("off") axes_f[-1].text(0.05, 0.6, "軒高 (TATE_H) 列について:\n" "• 単位は m, EPSG:6671 の 標高基準\n" "• 一部市町は 0 のみ (未測定)\n" "• P95=上位5%, P99=上位1% の境界", fontsize=11, va="top", bbox=dict(boxstyle="round,pad=0.5", facecolor="#ddf4ff", edgecolor="#0550ae")) plt.suptitle("7 市町 軒高 (TATE_H) 分布 — log y で裾を拡大", fontsize=13, y=1.00) plt.tight_layout() plt.savefig(ASSETS / "L14_height_hist.png", dpi=110) plt.close() print(f" 完了: {time.time()-tx:.1f}s") # ============================================================================= # 図4: 浸水 × 階数 致命率マップ (棒) + 散布 # ============================================================================= print("\n--- 図4 浸水 × 階数 ---") tx = time.time() fig, axes = plt.subplots(1, 2, figsize=(15, 6)) if len(flood_join_summary) > 0: fjs = flood_join_summary.copy() fjs = fjs.sort_values("pct_deadly_flat", ascending=False) ax = axes[0] bars = ax.barh(fjs["city"], fjs["pct_deadly_flat"], color="#cf222e", alpha=0.9) for b, v, n in zip(bars, fjs["pct_deadly_flat"], fjs["n_deadly_flat"]): ax.text(b.get_width() + 0.3, b.get_y() + b.get_height()/2, f"{v:.1f}% ({n} 棟)", va="center", fontsize=9) ax.set_xlabel("浸水域内 平屋 比率 (%) — 致命率") ax.set_title("致命率 (深3m以上 × 平屋) 市町ランキング") ax.invert_yaxis() ax.grid(axis="x", alpha=0.3) ax = axes[1] # 散布: x=浸水域内棟数, y=致命率, size=人口 pop_arr = np.array([CITY_REF[c]["pop"] / 1000 for c in fjs["city"]]) ax.scatter(fjs["n_in_flood"], fjs["pct_deadly_flat"], s=pop_arr * 0.5 + 30, c="#0969da", alpha=0.7, edgecolor="#0550ae", linewidth=1.0) for _, r in fjs.iterrows(): ax.annotate(r["city"], (r["n_in_flood"], r["pct_deadly_flat"]), fontsize=9, ha="left", va="bottom", xytext=(4, 4), textcoords="offset points") ax.set_xlabel("浸水域内 建物数 (棟)") ax.set_ylabel("致命率 (%)") ax.set_xscale("symlog") ax.set_title("浸水域曝露 × 致命率 (バブル=人口千人)") ax.grid(alpha=0.3) else: for ax in axes: ax.text(0.5, 0.5, "浸水データ未取得", ha="center", va="center", transform=ax.transAxes, fontsize=14) ax.axis("off") plt.tight_layout() plt.savefig(ASSETS / "L14_flood_compare.png", dpi=110) plt.close() print(f" 完了: {time.time()-tx:.1f}s") # ============================================================================= # 図5: 主題図 — 浸水域 × 建物の重ね合わせ (代表市町、広島市) # ============================================================================= print("\n--- 図5 浸水×建物 主題図 (広島市拡大) ---") tx = time.time() fig, axes = plt.subplots(1, 2, figsize=(16, 8)) target_city_for_overlay = "広島市" focus_g = None for slug, g in city_data.items(): if g["city"].iloc[0] == target_city_for_overlay: focus_g = g break if focus_g is not None and DEPTH_SHP.exists(): # 広島市域の bbox minx, miny, maxx, maxy = focus_g.total_bounds pad = max(maxx - minx, maxy - miny) * 0.02 # 描画用 simplify 済 cache を使用 (再 simplify しない) flood_clip = flood_for_plot.cx[minx-pad:maxx+pad, miny-pad:maxy+pad].copy() # 建物は polygon ではなく代表点 (markersize 小) で描画 → 高速化 focus_pts = gpd.GeoDataFrame( focus_g[["KAISUU", "kaisuu_band"]].copy(), geometry=focus_g.geometry.representative_point(), crs=TARGET_CRS, ) # 左: 浸水域 + 建物 (KAISUU 色) ax = axes[0] flood_clip[flood_clip["depth_rank"] >= 30].plot( ax=ax, color="#56a4dc", alpha=0.45, edgecolor="none") flood_clip[flood_clip["depth_rank"] >= 50].plot( ax=ax, color="#1c63a4", alpha=0.55, edgecolor="none") # 建物 (サンプル 25000 点、点描画は polygon より 10x 速い) fb_pts = focus_pts.sample(min(25000, len(focus_pts)), random_state=42) for band in KAISUU_BAND_ORDER: sb = fb_pts[fb_pts["kaisuu_band"] == band] if len(sb) > 0: sb.plot(ax=ax, color=KAISUU_BAND_COLOR[band], alpha=0.7, markersize=1.5, linewidth=0) ax.set_xlim(minx - pad, maxx + pad) ax.set_ylim(miny - pad, maxy + pad) ax.set_aspect("equal") ax.set_title(f"{target_city_for_overlay} — 浸水深 (青) × 建物 (色=階数)\n" "薄青: 1m以上 / 濃青: 3m以上", fontsize=11) ax.set_xticks([]); ax.set_yticks([]) legend_handles = [Patch(facecolor=KAISUU_BAND_COLOR[b], label=KAISUU_BAND_LABEL[b]) for b in KAISUU_BAND_ORDER if b != "0_その他"] legend_handles += [ Patch(facecolor="#56a4dc", alpha=0.45, label="浸水想定 1m以上"), Patch(facecolor="#1c63a4", alpha=0.55, label="浸水想定 3m以上"), ] ax.legend(handles=legend_handles, loc="lower right", fontsize=7, framealpha=0.92) # 右: 致命的建物 (深3m+平屋) を赤で強調 ax = axes[1] flood_clip[flood_clip["depth_rank"] >= 50].plot( ax=ax, color="#1c63a4", alpha=0.4, edgecolor="none") if buildings_in_flood is not None: fb_flood = buildings_in_flood[buildings_in_flood["city"] == target_city_for_overlay] deadly_pts = fb_flood[fb_flood["deadly_flat"]] safe_pts = fb_flood[fb_flood["safe_2plus"] & fb_flood["deep_3m"]] if len(safe_pts) > 0: sp = safe_pts.sample(min(2000, len(safe_pts)), random_state=42) sp.plot(ax=ax, color="#0969da", markersize=3, alpha=0.6, label=f"3m+ × 2階以上 (相対安全) {len(safe_pts):,} 棟") if len(deadly_pts) > 0: dp = deadly_pts.sample(min(2000, len(deadly_pts)), random_state=42) dp.plot(ax=ax, color="#dc2626", markersize=6, alpha=0.92, label=f"致命的 (3m+ × 平屋) {len(deadly_pts):,} 棟") ax.set_xlim(minx - pad, maxx + pad) ax.set_ylim(miny - pad, maxy + pad) ax.set_aspect("equal") ax.set_title(f"{target_city_for_overlay} — 致命的建物 (3m+ × 平屋, 赤)\n" "赤: 命の危険最高ランク, 青点: 3m+ でも 2階以上で相対安全", fontsize=11) ax.set_xticks([]); ax.set_yticks([]) ax.legend(loc="lower right", fontsize=8, framealpha=0.92) else: for ax in axes: ax.text(0.5, 0.5, "対象データなし", ha="center", va="center", transform=ax.transAxes, fontsize=14) ax.axis("off") plt.tight_layout() plt.savefig(ASSETS / "L14_flood_overlay_map.png", dpi=110) plt.close() print(f" 完了: {time.time()-tx:.1f}s") # ============================================================================= # 図6: 高度プロファイル — 中心 vs 郊外 (DID 内外) # ============================================================================= print("\n--- 図6 中心市街地 (DID) vs 郊外 ---") tx = time.time() DID_DIR = ROOT / "data" / "extras" / "did_extracted" # CITY_CD → 期待される市町名 マッピング (DID ファイル整合チェック用) CITY_CD_TO_JP = { 101: "広島市", 102: "広島市", 103: "広島市", 104: "広島市", 105: "広島市", 106: "広島市", 107: "広島市", 108: "広島市", 202: "呉市", 203: "竹原市", 207: "福山市", 208: "府中市", 209: "三次市", 304: "海田町", } did_compare = [] fig, ax = plt.subplots(1, 1, figsize=(12, 6)) i = 0 xtick_labels = [] xtick_positions = [] did_skip_log = [] city_order_for_did = ordered_cities # 人口降順 for jp in city_order_for_did: # 該当 GeoJSON did_dir = DID_DIR / f"did_{jp}" g = None for s, gv in city_data.items(): if gv["city"].iloc[0] == jp: g = gv; break if g is None: continue if not did_dir.exists(): did_skip_log.append((jp, "DID ディレクトリなし")); continue did_files = list(did_dir.glob("*.geojson")) if not did_files: did_skip_log.append((jp, "DID GeoJSON なし")); continue did = gpd.read_file(did_files[0]) # 整合性チェック: CITY_CD と期待市町名が一致するか if "CITY_CD" in did.columns: did_city_cds = set(int(x) for x in did["CITY_CD"].unique()) expected_jp_set = set(CITY_CD_TO_JP.get(c, "") for c in did_city_cds) if jp not in expected_jp_set: did_skip_log.append((jp, f"DID CITY_CD={sorted(did_city_cds)} は{sorted(expected_jp_set)} 用 (ファイル名と内容が不整合 — DoBoX 抽出ミス)")) continue if did.crs is None: did = did.set_crs(TARGET_CRS) elif did.crs.to_epsg() != 6671: did = did.to_crs(TARGET_CRS) did = did[["geometry"]].copy() did["geometry"] = did.geometry.buffer(0) did_u = did.dissolve() did_geom = did_u.geometry.iloc[0] # 建物代表点 (sjoin で高速 within) bld_pts_local = gpd.GeoDataFrame( g[["KAISUU"]].copy(), geometry=g.geometry.representative_point(), crs=TARGET_CRS, ) bld_pts_local["bld_id_local"] = np.arange(len(bld_pts_local)) # sjoin で R-tree インデックス利用 (within ループより 10x 高速) did_gdf = gpd.GeoDataFrame({"_": [1]}, geometry=[did_geom], crs=TARGET_CRS) sj = gpd.sjoin(bld_pts_local, did_gdf, how="left", predicate="within") in_did_set = set(sj.dropna(subset=["index_right"])["bld_id_local"].values) bld_pts_local["in_did"] = bld_pts_local["bld_id_local"].isin(in_did_set) inside = bld_pts_local[bld_pts_local["in_did"]]["KAISUU"] outside = bld_pts_local[~bld_pts_local["in_did"]]["KAISUU"] inside = inside[inside > 0] outside = outside[outside > 0] # ボックス用に値 if len(inside) > 0 and len(outside) > 0: ax.boxplot([inside, outside], positions=[i*3, i*3+1], widths=0.7, boxprops=dict(facecolor="#fbbf24", color="#7c1d6f"), medianprops=dict(color="#7c1d6f", linewidth=1.6), patch_artist=True, showfliers=False) # 平均 ax.plot(i*3, inside.mean(), "o", color="#cf222e", markersize=8) ax.plot(i*3+1, outside.mean(), "o", color="#0969da", markersize=8) xtick_positions += [i*3 + 0.5] xtick_labels += [jp] did_compare.append({ "city": jp, "n_in": len(inside), "n_out": len(outside), "median_in": float(inside.median()), "median_out": float(outside.median()), "mean_in": float(inside.mean()), "mean_out": float(outside.mean()), "diff_mean": float(inside.mean() - outside.mean()), "ratio_in_out": float(inside.mean() / max(outside.mean(), 0.01)), }) i += 1 ax.set_xticks(xtick_positions) ax.set_xticklabels(xtick_labels) # 各市町ペアの説明: 左=DID内, 右=DID外 for px in xtick_positions: ax.text(px - 0.5, -0.5, "DID内", ha="center", fontsize=7, color="#7c1d6f") ax.text(px + 0.5, -0.5, "DID外", ha="center", fontsize=7, color="#0969da") ax.set_ylim(0, max(8, ax.get_ylim()[1])) ax.set_ylabel("階数 KAISUU") ax.set_title("中心市街地 (DID) 内外 で 建物階数の箱ひげ比較\n赤丸=DID内平均, 青丸=DID外平均") ax.grid(axis="y", alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "L14_did_compare.png", dpi=110) plt.close() did_compare_df = pd.DataFrame(did_compare) did_compare_df.to_csv(ASSETS / "L14_did_compare.csv", index=False, encoding="utf-8-sig") print(f" 完了: {time.time()-tx:.1f}s") print(did_compare_df.to_string(index=False)) # ============================================================================= # 図7: 用途バンド × 平均階数 ヒートマップ # ============================================================================= print("\n--- 図7 用途×階数 ヒートマップ ---") tx = time.time() yoto_kaisuu_mean = (all_b[all_b["KAISUU"] > 0] .groupby(["city", "yoto_band"])["KAISUU"] .mean() .unstack(fill_value=np.nan)) yoto_keep = [c for c in yoto_order if c in yoto_kaisuu_mean.columns] yoto_kaisuu_mean = yoto_kaisuu_mean[yoto_keep] yoto_kaisuu_mean.columns = [YOTO_BAND_LABEL.get(c, c) for c in yoto_kaisuu_mean.columns] yoto_kaisuu_mean = yoto_kaisuu_mean.reindex(ordered_cities) yoto_kaisuu_mean.to_csv(ASSETS / "L14_yoto_kaisuu_mean.csv", encoding="utf-8-sig") fig, ax = plt.subplots(figsize=(12, 5)) hm = yoto_kaisuu_mean.values im = ax.imshow(hm, cmap="YlOrRd", aspect="auto", vmin=1.0, vmax=4.0) ax.set_xticks(range(len(yoto_kaisuu_mean.columns))) ax.set_xticklabels(yoto_kaisuu_mean.columns, rotation=20) ax.set_yticks(range(len(yoto_kaisuu_mean.index))) ax.set_yticklabels(yoto_kaisuu_mean.index) plt.colorbar(im, ax=ax, label="平均階数") ax.set_title("市町 × 用途バンド の 平均階数 ヒートマップ\n赤=高階数, 黄=低階数") for r in range(len(yoto_kaisuu_mean.index)): for c in range(len(yoto_kaisuu_mean.columns)): v = hm[r, c] if not np.isnan(v): ax.text(c, r, f"{v:.1f}", ha="center", va="center", fontsize=9, color="white" if v > 2.5 else "black") plt.tight_layout() plt.savefig(ASSETS / "L14_yoto_kaisuu_heatmap.png", dpi=110) plt.close() print(f" 完了: {time.time()-tx:.1f}s") # ============================================================================= # 図8: 主題図 (合算) — 7 市町をひとつの広島県地図で # ============================================================================= print("\n--- 図8 県全域 主題図 ---") tx = time.time() fig, ax = plt.subplots(figsize=(14, 10)) # サンプル数を都市別に減らす (合計 ~50000 ぐらい) parts_for_map = [] for slug, g in city_data.items(): n = len(g) sn = min(8000, n) parts_for_map.append(g.sample(sn, random_state=42) if n > sn else g) map_b = gpd.GeoDataFrame(pd.concat(parts_for_map, ignore_index=True), crs=TARGET_CRS) for band in KAISUU_BAND_ORDER: sb = map_b[map_b["kaisuu_band"] == band] if len(sb) > 0: sb.plot(ax=ax, color=KAISUU_BAND_COLOR[band], alpha=0.85, edgecolor="none", linewidth=0) ax.set_title("広島県 7 市町 建物分布 主題図 (色=階数バンド)\n" "黄=平屋 → 紫=10階以上", fontsize=13) ax.set_xlabel("X (m, JGD2011 平面直角 III)") ax.set_ylabel("Y (m)") ax.set_aspect("equal") legend_handles = [Patch(facecolor=KAISUU_BAND_COLOR[b], label=KAISUU_BAND_LABEL[b]) for b in KAISUU_BAND_ORDER] ax.legend(handles=legend_handles, loc="lower left", fontsize=9, framealpha=0.95) plt.tight_layout() plt.savefig(ASSETS / "L14_pref_map.png", dpi=110) plt.close() print(f" 完了: {time.time()-tx:.1f}s") print(f"\n[全体実行時間] {time.time()-t0:.1f}s") # ============================================================================= # === HTML 生成 =============================================================== # ============================================================================= print("\n--- HTML 生成 ---") # 仮説検証 def jp_v(b): if b is None: return "保留" if bool(b): return "支持" return "反証" # H1: 広島市は最も高層建物率が高い hi_high_pct = float(city_stats_df[city_stats_df["city"] == "広島市"]["rate_high_pct"].iloc[0]) others_high_pct = city_stats_df[city_stats_df["city"] != "広島市"]["rate_high_pct"].max() h1 = hi_high_pct > others_high_pct # H2: 山間部 (三次市) は平屋率が高い m_flat = float(city_stats_df[city_stats_df["city"] == "三次市"]["rate_flat_pct"].iloc[0]) m_others = city_stats_df[~city_stats_df["city"].isin(["三次市","府中市"])]["rate_flat_pct"].mean() h2 = m_flat > m_others # H3: 浸水想定 3m 以上の地域に「平屋」=命の危険最高ランク (deadly_flat の存在) n_deadly_total = int(flood_join_summary["n_deadly_flat"].sum()) if len(flood_join_summary) > 0 else 0 h3 = n_deadly_total > 0 # H4: 港湾系 (呉市/福山市) は中層 (3-4階) 中心 def band_pct(city, band): if city not in band_cross_pct.index or band not in band_cross_pct.columns: return 0.0 return float(band_cross_pct.at[city, band]) ku_mid = band_pct("呉市", "3_3-4階") fk_mid = band_pct("福山市", "3_3-4階") miy_mid = band_pct("三次市", "3_3-4階") h4 = (ku_mid + fk_mid) / 2 > miy_mid # H5: 旧市街地 (竹原) は低層密集 (平屋+2階 が高比率) tk_low = band_pct("竹原市", "1_平屋") + band_pct("竹原市", "2_2階") hi_low = band_pct("広島市", "1_平屋") + band_pct("広島市", "2_2階") h5 = tk_low > hi_low # H6: 浸水域内 平屋率は市町で 10倍以上の差 if len(flood_join_summary) > 0: flat_rates = flood_join_summary["pct_deadly_flat"].replace(0, np.nan).dropna() if len(flat_rates) >= 2: h6 = float(flat_rates.max() / flat_rates.min()) >= 10 else: h6 = None else: h6 = None # 結果サマリー total_load_sec = sum(r["sec"] for r in city_load_log) n_cached = sum(1 for r in city_load_log if r["cached"]) elapsed_sec = time.time() - t0 # 市町ステータス表 city_status_html = ""\ "" for slug, jp, prefix, ds3d, dsbld in CITY_DEFS: rec = next((r for r in city_load_log if r["slug"] == slug), None) status = "OK" if rec and rec["status"] == "OK" else "FAIL" n_b = rec["rows"] if rec else 0 sec = rec["sec"] if rec else 0 cached_mark = " (cached)" if (rec and rec["cached"]) else "" status_color = "#0a7d2f" if status == "OK" else "#cf222e" city_status_html += ( f"" f"" f"" f"" f"" f"" ) city_status_html += "
市町3D都市モデル建物利用現況取得状況建物数読込時間 (s)
{jp}#{ds3d}#{dsbld}{status}{cached_mark}{n_b:,}{sec:.1f}
" # city_stats 表 city_stats_html = ""\ ""\ "" for _, r in city_stats_df.iterrows(): city_stats_html += ( f"" f"" f"" f"" f"" f"" f"" f"" f"" ) city_stats_html += "
市町人口面積 km²建物数棟/km²平均階数最大階数5階以上率平屋率
{r['city']}{r['pop']:,}{r['area_km2']:.1f}{int(r['n_bld']):,}{r['bld_per_km2']:.1f}{r['mean_kaisuu']:.2f}{int(r['max_kaisuu'])}{r['rate_high_pct']:.2f}%{r['rate_flat_pct']:.1f}%
" # 階数バンド構成比 表 band_pct_html = "" for b in KAISUU_BAND_ORDER: band_pct_html += f"" band_pct_html += "" for jp in band_cross_pct.index: band_pct_html += f"" for b in KAISUU_BAND_ORDER: v = band_cross_pct.at[jp, b] if b in band_cross_pct.columns else 0 cell_color = "" if b == "5_10階以上" and v > 0.5: cell_color = " style='background:#fff5f5'" elif b == "1_平屋" and v > 30: cell_color = " style='background:#fff8c5'" band_pct_html += f"{v:.1f}%" band_pct_html += "" band_pct_html += "
市町{KAISUU_BAND_LABEL[b]}
{jp}
" # 浸水オーバーレイ表 if len(flood_join_summary) > 0: flood_html = ("" "" "") for _, r in flood_join_summary.iterrows(): flood_html += ( f"" f"" f"" f"" f"" f"" f"" f"" f"" ) flood_html += "
市町建物数浸水域内浸水率3m以上うち平屋(致命)致命率2階以上(相対安全)相対安全率
{r['city']}{int(r['n_total']):,}{int(r['n_in_flood']):,}{r['pct_in_flood']:.1f}%{int(r['n_deep_3m']):,}{int(r['n_deadly_flat']):,}{r['pct_deadly_flat']:.1f}%{int(r['n_safe_2plus']):,}{r['pct_safe_2plus_of_flooded']:.1f}%
" else: flood_html = "

浸水データ未取得 — sjoin スキップ

" # DID 比較表 did_html = ""\ "" for _, r in did_compare_df.iterrows(): did_html += ( f"" f"" f"" f"" f"" f"" f"" ) did_html += "
市町DID内 建物数DID外 建物数DID内 平均階数DID外 平均階数差 (中心-郊外)比 (中心/郊外)
{r['city']}{int(r['n_in']):,}{int(r['n_out']):,}{r['mean_in']:.2f}{r['mean_out']:.2f}+{r['diff_mean']:.2f}{r['ratio_in_out']:.2f}
" if did_skip_log: did_html += "

DID 比較から除外された市町:

" # 軒高分布表 height_html = ""\ "" for _, r in height_dist_df.iterrows(): note = " (TATE_H 全 0=未測定)" if r["n_h_known"] == 0 else "" height_html += ( f"" f"" f"" f"" f"" f"" f"" ) height_html += "
市町軒高 既知数中央値 mP75P95P99最大
{r['city']}{int(r['n_h_known']):,}{note}{r['median_h']:.1f}{r['p75_h']:.1f}{r['p95_h']:.1f}{r['p99_h']:.1f}{r['max_h']:.1f}
" # 用途バンド × 階数 ヒートマップ表 yoto_h_html = "" for c in yoto_kaisuu_mean.columns: yoto_h_html += f"" yoto_h_html += "" for jp in yoto_kaisuu_mean.index: yoto_h_html += f"" for c in yoto_kaisuu_mean.columns: v = yoto_kaisuu_mean.at[jp, c] if c in yoto_kaisuu_mean.columns else np.nan if pd.isna(v): yoto_h_html += "" else: color_style = " style='background:#fff5f5'" if v >= 3.0 else "" yoto_h_html += f"{v:.2f}" yoto_h_html += "" yoto_h_html += "
市町＼用途{c}
{jp}—
" # H6 根拠用の値計算 if len(flood_join_summary) > 0: _fr = flood_join_summary['pct_deadly_flat'].replace(0, np.nan).dropna() if len(_fr) >= 2: h6_ratio_max_min = float(_fr.max() / max(_fr.min(), 0.01)) h6_evidence = f"致命率の最大/最小 比 = {h6_ratio_max_min:.1f} 倍 (最大={_fr.max():.1f}%, 最小={_fr.min():.2f}%)" else: h6_evidence = "データ不足 (有効市町 < 2)" else: h6_evidence = "浸水データ未取得" # 仮説検証表 hyp_html = "" hyp_html += ( f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" ) hyp_html += "
仮説内容判定根拠
H1広島市は最も高層建物率が高い (5階以上率){jp_v(h1)}広島市 {hi_high_pct:.2f}% vs 他市町 最大 {others_high_pct:.2f}%
H2山間部 (三次市) は平屋率が高い{jp_v(h2)}三次市 平屋率 {m_flat:.1f}% vs 他市町 (除く府中) 平均 {m_others:.1f}%
H3浸水想定3m+ × 平屋 = 命の危険最高ランクが存在する{jp_v(h3)}致命的建物 (深3m+ × 平屋) 合計 {n_deadly_total:,} 棟
H4港湾系 (呉市/福山市) は中層 (3-4階) 中心{jp_v(h4)}呉/福山 平均3-4階率 {(ku_mid+fk_mid)/2:.1f}% vs 三次 {miy_mid:.1f}%
H5旧市街地 (竹原市) は低層密集 (平屋+2階){jp_v(h5)}竹原 低層率 {tk_low:.1f}% vs 広島 低層率 {hi_low:.1f}%
H6浸水域内 平屋率は市町間で 10倍以上の差{jp_v(h6)}{h6_evidence}
" # サンプル建物 (Before/After 表) sample_jp = "竹原市" sample_g = None for s, g in city_data.items(): if g["city"].iloc[0] == sample_jp: sample_g = g; break if sample_g is not None and len(sample_g) > 0: sample_row = sample_g.iloc[0] geom_str = f"POLYGON ({len(list(sample_row.geometry.exterior.coords))} 頂点)" before_after = ( "" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" f"" "
段階このデータで何が起きるか結果サイズ
1) 生 GeoJSON 1 行1 棟分の属性 + 多角形 footprint。例: KAISUU={int(sample_row['KAISUU'])}, " f"TATE_H={sample_row['TATE_H']:.2f}, TATE_YO={int(sample_row['TATE_YO'])}, " f"geometry={geom_str}1 棟 (~30 列)
2) {sample_jp} 全棟読込同フォーマットを {len(sample_g):,} 行に展開{len(sample_g):,} 棟
3) parquet キャッシュ必要列だけ抽出 → バイナリ保存~10 MB / 数秒で再読込
4) kaisuu_band 列を追加KAISUU={int(sample_row['KAISUU'])} → kaisuu_band='{kaisuu_band(int(sample_row['KAISUU']))}' (={KAISUU_BAND_LABEL[kaisuu_band(int(sample_row['KAISUU']))]})+1 列
5) 7 市町 concat各市町を `pd.concat` で縦結合, city 列で識別{len(all_b):,} 棟
6) 代表点 → 浸水深 sjoin建物多角形 → 1 点 (representative_point) → 浸水深ポリゴンに含まれるか判定各棟に depth_rank 列付与
" ) else: before_after = "

サンプルなし

" # 仮説保留時の "n_a" 修正: 文字列を後で組み立てた箇所を整形 # (上で h6 = None のとき "保留" に倒している) sections = [ ("学習目標と問い", f"""
本記事のスタイル: 7 市町の建物 footprint と高度属性を比較し、「都市階層構造」と「浸水時の危険度」を可視化
3D 都市モデル (PLATEAU CityGML) は重量級 (各市町 100 MB〜1 GB) のため、 本記事は 建物利用現況 GeoJSON (#1469 等, 階数 KAISUU と軒高 TATE_H を含む) を用いて、 PLATEAU の建物 footprint レイヤと同等の地理学的分析を行う。 浸水想定深 3 m 以上 × 平屋 = 命の危険最高ランクという 「災害人類学」 の視点で 7 市町を比較する。

カバー宣言

本記事は 3D都市モデル 7 市町 = 7 dataset_id を統合分析:

全 PLATEAU CityGML 規格で統一されており、統合可能。
ただし CityGML 本体は重量級なため、同 7 市町の建物利用現況 GeoJSON (#1469 / #1483 / #1323 / #1327 / #512 / #71 / #1474) を Plan B として採用。 これらは PLATEAU の建物 footprint + 属性 (階数・軒高) と内容が一致する軽量版で、 本研究の論点 (建物高度の地理学) には十分。
対応表 7/7 件 論理カバー。

主な問い

  1. 都市階層の問い: 7 市町で建物階数の構成は どう違うか? 政令市 (広島市) と中山間市町 (三次市) で何が変わるか?
  2. 用途と階数の関係: 同じ「商業」でも、都市規模で平均階数は変わるか?
  3. 中心市街地効果: DID (人口集中地区) 内外で 階数の差はどれだけ広がるか?
  4. 浸水時の致命率: 浸水想定 3 m 以上の地域に立つ平屋 (1 階建て) は何棟で、市町別にどれだけ差があるか?

立てた仮説 (H1〜H6)

  1. H1: 広島市は最も高層建物率が高い (5 階以上が多い)
  2. H2: 山間部市町 (三次市) は平屋率が高い
  3. H3: 浸水想定深 3 m 超の地域に「平屋」 = 命の危険最高ランク (致命的建物) が確実に存在する
  4. H4: 港湾系 (呉市/福山市) は中層 (3-4 階) 中心
  5. H5: 旧市街地 (竹原は古い町並み) は低層密集 (平屋+2 階の比率が大都市より高い)
  6. H6: 浸水域内 平屋率は市町間で 10 倍以上の差

用語の定義 (本レッスン独自を含む)

結果サマリー

指標結果
対象市町7 市町 (人口 23,000 〜 1,189,000)
統合 建物数{n_bld_total:,} 棟
最大 階数{int(city_stats_df['max_kaisuu'].max())} 階 ({city_stats_df.loc[city_stats_df['max_kaisuu'].idxmax(), 'city']})
最大 軒高{height_dist_df['max_h'].max():.0f} m ({height_dist_df.loc[height_dist_df['max_h'].idxmax(), 'city']})
5 階以上 建物 (高層){int(city_stats_df['n_high'].sum()):,} 棟
10 階以上 建物 (タワー){int(city_stats_df['n_tower'].sum()):,} 棟
致命的建物 (浸水3m+ × 平屋){n_deadly_total:,} 棟
処理時間{elapsed_sec:.1f} 秒 ({n_cached}/{len(CITY_DEFS)} 市町 cached)
"""), ("使用データ", f"""

主データ (Plan B: 軽量 GeoJSON 採用)

{city_status_html}

カバー対応表 (DoBoX 7 dataset_id ↔ 本記事)

本記事は 3D 都市モデル 7 dataset_id を「論理カバー」する。 CityGML 直接処理は重すぎるため、同 7 市町の建物利用現況 GeoJSON で代替。 PLATEAU は CityGML から 建物 footprint と階数・軒高を生成しているので、 本記事の分析対象 (= 階数・軒高・footprint) は両者で一致する。

市町3D都市モデル CityGML (本来の対象)建物利用現況 GeoJSON (採用 = Plan B)
広島市#1316 (~950 MB)#1469 (~33 MB)
府中市#1302#512
呉市#1187#1323 (~12 MB)
福山市#1188#1483 (~29 MB)
海田町#1189#71 (~1.3 MB)
竹原市#1443#1474 (~1.9 MB)
三次市#1444#1327 (~1.1 MB)

対応表 7/7 件 論理カバー。

補助データ

"""), ("ダウンロード", """

中間データ (再現用 CSV)

ファイル内容
L14_city_stats.csv7 市町 × 基本統計 (建物数/平均階数/最大階数/5階以上率/平屋率)
L14_band_cross_count.csv市町×階数バンド 棟数
L14_band_cross_pct.csv市町×階数バンド 構成比
L14_yoto_cross.csv市町×用途バンド 棟数
L14_yoto_cross_pct.csv市町×用途バンド 構成比
L14_yoto_kaisuu_mean.csv市町×用途バンド 平均階数 (ヒートマップ元)
L14_height_dist.csv軒高 (TATE_H) 中央値・P75/P95/P99/最大
L14_flood_overlay.csv浸水域 sjoin 結果 (致命率含む)
L14_deadly_summary.csv致命率 市町ランキング
L14_did_compare.csvDID 内外 階数比較

図 PNG

ファイル内容
L14_small_multiples_kaisuu.png図1 7 市町 small multiples (主題図)
L14_band_compare.png図2 階数バンド構成比 + 5階以上率/平屋率
L14_height_hist.png図3 軒高 (TATE_H) 分布 ヒストグラム
L14_flood_compare.png図4 浸水×階数 致命率比較
L14_flood_overlay_map.png図5 浸水×建物 主題図 (広島市拡大)
L14_did_compare.png図6 DID 内外 箱ひげ比較
L14_yoto_kaisuu_heatmap.png図7 市町×用途 平均階数ヒートマップ
L14_pref_map.png図8 県全域 主題図
L14_3d_city_model_compare.py再現スクリプト (本記事)

取得手順 (PowerShell)

cd "2026 DoBoX 教材"
# 7 市町の建物利用現況 GeoJSON を取得 (合計 ~80 MB, ~20 秒)
mkdir data\\extras\\bld_data -Force
iwr "https://hiroshima-dobox.jp/resource_download/94767" -OutFile "data\\extras\\bld_hiroshima.zip"  # 広島市
iwr "https://hiroshima-dobox.jp/resource_download/94848" -OutFile "data\\extras\\bld_fukuyama.zip"   # 福山市
iwr "https://hiroshima-dobox.jp/resource_download/50387" -OutFile "data\\extras\\bld_kure.zip"       # 呉市
iwr "https://hiroshima-dobox.jp/resource_download/50423" -OutFile "data\\extras\\bld_miyoshi.zip"    # 三次市
iwr "https://hiroshima-dobox.jp/resource_download/51056" -OutFile "data\\extras\\bld_fuchu.zip"      # 府中市
iwr "https://hiroshima-dobox.jp/resource_download/51057" -OutFile "data\\extras\\bld_kaita.zip"      # 海田町
iwr "https://hiroshima-dobox.jp/resource_download/94794" -OutFile "data\\extras\\bld_takehara.zip"   # 竹原市
# 解凍 (PowerShell 5.1+)
Expand-Archive data\\extras\\bld_hiroshima.zip data\\extras\\bld_data -Force
# (他都市も同様に)

# 浸水深 (X09 と同じ #1278)
mkdir data\\extras\\flood_shp -Force
# (X09 / L09 で取得済みなら不要)

py -X utf8 lessons\\L14_3d_city_model_compare.py

キャッシュ機構: 初回実行で data/extras/_l14_cache/{slug}.parquet を生成。 2 回目以降は parquet から数秒で読込 (要件 S 1 分以内を担保)。 キャッシュを破棄したい場合は _l14_cache フォルダを削除。

"""), ("分析1: 市町別 基本統計 — 都市規模と建物プロファイル", f"""

狙い

7 市町を 政令市 (広島) → 中核市 (福山・呉) → 中山間市町 (三次・府中) → 都市近郊町 (海田町) → 瀬戸内市 (竹原) の階層に並べ、各階層で建物プロファイル (棟数密度・平均階数・5 階以上率・平屋率) がどう変わるかを概観する。 これが H1 (広島市が最も高層) と H2 (三次市が最も平屋多) の検証の出発点となる。

手法

実装

{code(''' city_stats = [] for slug, g in city_data.items(): jp = g["city"].iloc[0] n = len(g) kaisuu_valid = g[g["KAISUU"] > 0]["KAISUU"] h_valid = g[g["TATE_H"] > 0]["TATE_H"] ref = CITY_REF.get(jp, {}) city_stats.append({ "city": jp, "slug": slug, "pop": ref.get("pop", 0), "area_km2": ref.get("area_km2", 0), "n_bld": n, "bld_per_km2": n / max(ref.get("area_km2", 1), 1), "mean_kaisuu": float(kaisuu_valid.mean()) if len(kaisuu_valid) else 0, "max_kaisuu": int(kaisuu_valid.max()) if len(kaisuu_valid) else 0, "rate_high_pct": float((g["KAISUU"] >= 5).sum() / max(n, 1) * 100), # 5階以上率 "rate_flat_pct": float((g["KAISUU"] == 1).sum() / max(n, 1) * 100), # 平屋率 }) ''')}

図1 — 7 市町 small multiples 主題図

なぜこの図か: 棒グラフでは「広島市の中心部にタワーが集中している」「三次市は広く分散して平屋」のような 地理的なパターンが見えない。同じスケール感で 7 市町を並べることで、各市町の都市形態 (compact / sprawl) を直感的に比較する (要件 T)。

{figure("assets/L14_small_multiples_kaisuu.png", "図1 7 市町 small multiples — 色は階数バンド (黄=平屋, 紫=10階以上)")}

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

表 — 市町別 基本統計 (人口降順)

{city_stats_html}

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

仮説への含意: H1 (広島が最も高層) は {jp_v(h1)}、H2 (三次が平屋多) は {jp_v(h2)}。 この時点で 7 市町を「都市階層」に基づき分類できる。
"""), ("分析2: 階数バンド × 用途 — 都市の「機能ピラミッド」", f"""

狙い

建物を 5 つの階数バンド (平屋 / 2 階 / 3-4 階 / 5-9 階 / 10 階以上) に分類し、 さらに 用途バンド (業務・商業・住宅・工業・倉庫運輸他 等) と組み合わせて、 「政令市の高層は何用途か?」「中核市と中山間市で 商業の階数はどう違うか?」を見る。 これが H4 (港湾系は中層中心) を検証する基盤になる。

手法

階数バンド (本レッスン独自定義):

バンド定義都市的意味
平屋 (1 階)KAISUU == 1戸建住宅・倉庫・小規模店舗
2 階KAISUU == 2戸建住宅・町並みの主役
3-4 階 (低中層)3 ≤ KAISUU ≤ 4共同住宅・小規模ビル・商店併用
5-9 階 (中層)5 ≤ KAISUU ≤ 9マンション・オフィスビル・ホテル中規模
10 階以上 (高層)KAISUU ≥ 10タワーマンション・高層オフィス・ランドマーク

用途バンド (TATE_YO 401-499 を 7 区分に集約): 業務 / 商業 / 娯楽サービス / 文教 / 住宅 / 工業 / 倉庫運輸他

実装

{code(''' def kaisuu_band(k): if k <= 0: return "0_その他" if k == 1: return "1_平屋" if k == 2: return "2_2階" if k <= 4: return "3_3-4階" if k <= 9: return "4_5-9階" return "5_10階以上" def yoto_band(c): c = int(c) if 401 <= c <= 404: return "1_業務" if 411 <= c <= 415: return "2_商業" if 421 <= c <= 431: return "3_娯楽サービス" if c == 441: return "4_文教" if 451 <= c <= 454: return "5_住宅" if c == 461: return "6_工業" if 471 <= c <= 499: return "7_倉庫運輸他" return "9_不明" all_b["kaisuu_band"] = all_b["KAISUU"].apply(kaisuu_band) all_b["yoto_band"] = all_b["TATE_YO"].apply(yoto_band) # クロス集計: 市町×階数バンド (構成比) band_cross = all_b.groupby(["city", "kaisuu_band"]).size().unstack(fill_value=0) band_cross_pct = band_cross.div(band_cross.sum(axis=1), axis=0) * 100 ''')}

図2 — 階数バンド構成比 (積み上げ棒) + 高層率 vs 平屋率

なぜこの図か: 7 市町の階数構成を 積み上げ棒 1 本/市町 で見ると、 「都市階層を上に行くほど高層比率がスライドして増える」というスケーリング則を視覚的に確認できる。 右の棒は H1 vs H2 の対比 (高層率 vs 平屋率) を 1 図で並べる (要件 H)。

{figure("assets/L14_band_compare.png", "図2 階数バンド構成比 (左, 積み上げ) と 5階以上率 vs 平屋率 (右)")}

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

表 — 市町別 階数バンド構成比 (%)

{band_pct_html}

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

図7 — 用途×階数の組合せ (ヒートマップ)

なぜこの図か: 棒グラフでは 「住宅は何階が中心か」「商業はどこで階数が上がるか」のような 2 軸の関係性が見えない。ヒートマップで 1 図で集約 (要件 H)。

{figure("assets/L14_yoto_kaisuu_heatmap.png", "図7 市町×用途の平均階数 (赤=高階数, 黄=低階数)")}

この表 (図7 の元データ) から読み取れること:

{yoto_h_html}
仮説への含意: H4 (港湾系は中層中心) は {jp_v(h4)}。 呉市・福山市の 3-4 階バンド比率が三次市より明確に高い。 H5 (竹原は低層密集) は {jp_v(h5)} ({tk_low:.1f}% vs 広島市 {hi_low:.1f}%)。
"""), ("分析3: 軒高 (TATE_H) の分布 — メートル単位での「都市の身長」", f"""

狙い

階数 KAISUU は離散値 (1, 2, 3...) だが、軒高 TATE_H は m 単位の連続値。 同じ「3 階建て」でも、商業ビル (天井 4 m × 3 = 12 m) と低層マンション (天井 3 m × 3 = 9 m) では軒高が違う。 都市ごとの軒高分布を P95/P99 で見比べると、「裾の長さ」 = 高層化の進度が分かる。

手法

図3 — 軒高ヒストグラム (city 別, log y)

なぜこの図か: 軒高は典型的な 裾の長い分布 (大半が低い、ごく一部が高い)。 log y で描くと裾の構造が読み取れる。P95/P99 を縦線で示すことで「上位 5% / 1% の境界」を明示 (要件 H)。

{figure("assets/L14_height_hist.png", "図3 7 市町 軒高 (TATE_H) 分布 — log y")}

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

表 — 市町別 軒高分布要約

{height_html}

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

"""), ("分析4: 浸水域 × 建物高度 — 「致命率」の市町ランキング", f"""

狙い

本記事の最重要分析。
浸水想定深 3 m 以上 (#1278 の rank ≥ 50) の地域に立つ平屋 (1 階建て) は、 垂直避難先がなく 命の危険最高ランク である。これを 「致命的建物」 と定義し、市町ランキングを取る。 H3 (致命的建物の存在) と H6 (市町間 10 倍以上の差) を同時に検証する。

手法 — 代表点 sjoin による高速近似

素朴な発想: 全建物 footprint と 全浸水ポリゴンを gpd.overlay() で総当たり交差。 しかし 50 万棟 × 600 浸水ポリゴン = 3 億回の交差判定 → 数十分かかる (要件 S 違反)。

本記事の高速化: 建物 footprint を representative_point() で 1 点に縮約 → 浸水ポリゴンに sjoin(predicate='within')。 点 vs ポリゴンの sjoin は GEOS の R-Tree インデックスで O(N log M) で完了する。 建物 footprint の代表点が浸水ポリゴン内にあれば、その建物は浸水想定域内と判定。 精度は「建物の中心 1 点で代表する」近似だが、商業ビル / 戸建のいずれも footprint が小さいため十分。

致命的建物の判定式

建物 i が致命的 ⇔ depth_rank ≥ {DEEP_DEADLY_RANK} (3 m 以上) かつ KAISUU == 1 (平屋)

相対安全建物: in_flood かつ KAISUU ≥ 2 (浸水深に関わらず 2 階以上で垂直避難可能)

実装

{code(''' # 1) 建物 footprint → 代表点 (高速化) bld_pts = gpd.GeoDataFrame( all_b[["city", "KAISUU", "TATE_YO", "kaisuu_band", "yoto_band"]].copy(), geometry=all_b.geometry.representative_point(), crs=TARGET_CRS, ) # 2) 浸水深ポリゴンと sjoin (建物点がどの浸水ランクに含まれるか) flood = gpd.read_file(DEPTH_SHP, columns=["rank"]).to_crs(TARGET_CRS) flood["depth_rank"] = flood["rank"] bld_in_flood = gpd.sjoin(bld_pts, flood[["depth_rank", "geometry"]], how="left", predicate="within") # 3) 致命的・相対安全フラグ bld_in_flood["in_flood"] = bld_in_flood["depth_rank"].notna() bld_in_flood["depth_rank"] = bld_in_flood["depth_rank"].fillna(0).astype(int) bld_in_flood["deep_3m"] = bld_in_flood["depth_rank"] >= 50 # 3m 以上 bld_in_flood["deadly_flat"] = bld_in_flood["deep_3m"] & (bld_in_flood["KAISUU"] == 1) bld_in_flood["safe_2plus"] = bld_in_flood["in_flood"] & (bld_in_flood["KAISUU"] >= 2) # 4) 市町別集計 flood_join_summary = bld_in_flood.groupby("city").agg( n_total=("in_flood", "size"), n_in_flood=("in_flood", "sum"), n_deep_3m=("deep_3m", "sum"), n_deadly_flat=("deadly_flat", "sum"), n_safe_2plus=("safe_2plus", "sum"), ).reset_index() flood_join_summary["pct_deadly_flat"] = ( flood_join_summary["n_deadly_flat"] / flood_join_summary["n_in_flood"] * 100 ) ''')}

図4 — 致命率 市町ランキング + バブル散布

なぜこの図か: ランキング棒で「どの市町が致命率が高いか」を直感的に見せ、 散布図で「曝露量 (浸水域内棟数) と致命率」の 2 軸関係 + 人口バブルで 3 軸目も同時表示 (要件 H, T)。

{figure("assets/L14_flood_compare.png", "図4 致命率ランキング (左) と 曝露×致命率 散布 (右, バブル=人口千人)")}

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

図5 — 浸水 × 建物 主題図 (広島市拡大)

なぜこの図か: ランキングだけでは「どの地区に致命的建物が集中するか」が見えない。 広島市デルタを拡大し、浸水深 3 m 以上 (青塗り) × 致命的建物 (赤点) を 1 枚に重ねることで、 「太田川河口南岸の平屋密集地」のような具体的な地理パターンを可視化 (要件 T 主役級)。

{figure("assets/L14_flood_overlay_map.png", "図5 広島市 — 左: 浸水深 × 階数バンド, 右: 致命的建物 (赤) と相対安全 (青)")}

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

表 — 浸水域 × 階数 sjoin 集計 (市町別)

{flood_html}

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

仮説への含意: H3 (致命的建物存在) は {jp_v(h3)}。H6 (市町間 10 倍以上の差) は {jp_v(h6)}。 広島市デルタ南部の致命的建物群は、防災教育・住宅政策の優先対象。
"""), ("分析5: 中心市街地 (DID) 内外 — 「都市勾配」の数値化", f"""

狙い

同じ市町内でも 中心市街地 (DID) は階数が高く、郊外は平屋が多い はず。 DID 内外で建物階数を比較して 「都市勾配 (urban gradient)」 を市町別に数値化する。 これは 「中心 vs 郊外」 の差が 大都市 (広島) では大きく、小都市 (竹原) では平坦 という都市地理学の古典的仮説の検証。

手法

実装

{code(''' did_compare = [] for slug, g in city_data.items(): jp = g["city"].iloc[0] did_files = list((DID_DIR / f"did_{jp}").glob("*.geojson")) if not did_files: continue did = gpd.read_file(did_files[0]).to_crs(TARGET_CRS) did_u = did[["geometry"]].dissolve() # 市町内 DID を 1 ポリゴンに bld_pts = gpd.GeoDataFrame( g[["KAISUU"]].copy(), geometry=g.geometry.representative_point(), crs=TARGET_CRS, ) bld_pts["in_did"] = bld_pts.geometry.within(did_u.geometry.iloc[0]) inside = bld_pts[bld_pts["in_did"]]["KAISUU"] outside = bld_pts[~bld_pts["in_did"]]["KAISUU"] inside, outside = inside[inside > 0], outside[outside > 0] did_compare.append({ "city": jp, "mean_in": float(inside.mean()), "mean_out": float(outside.mean()), "diff_mean": float(inside.mean() - outside.mean()), "ratio_in_out": float(inside.mean() / max(outside.mean(), 0.01)), }) ''')}

図6 — DID 内外 箱ひげ比較

なぜこの図か: 平均だけだと外れ値 (タワー) に引っ張られる。箱ひげで中央値・四分位を見せ、平均は丸でオーバーレイ。 左 (DID 内) → 右 (DID 外) の階段状の落差が 「都市勾配」 の視覚化 (要件 H)。

{figure("assets/L14_did_compare.png", "図6 中心市街地 (DID) 内外 で建物階数の箱ひげ比較")}

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

表 — DID 内外 階数比較

{did_html}

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

分析間の連携: 分析1 で見た「広島が高層」「三次が平屋」は、ここで 市町内 (DID 内外) の構造に分解された。 広島の高層は DID 内に集中 (= 中心化)、三次は DID 内ですら平屋主体 (= フラット構造)。 これは政策含意: 中山間市町では「中心市街地の高層化」が進まず、平屋住宅が広域に分散している。
"""), ("仮説検証と考察", f"""

図8 — 広島県全域 7 市町 統合主題図 (シノプティックビュー)

なぜこの図か: 個別の small multiples (図1) では市町ごとの構造が見えるが、 県全体としてどこに高層ストックが集中しているか はワンショットで掴みにくい。 1 枚に 7 市町を統合した主題図は、本記事の「7 市町統合」というテーマ自体の象徴 (要件 T)。

{figure("assets/L14_pref_map.png", "図8 広島県 7 市町 建物分布 主題図 — 統合カバー宣言の視覚化")}

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

仮説 H1〜H6 判定

{hyp_html}

結果サマリー (主要な発見)

本記事の限界

"""), ("発展課題", f"""

結果X → 新仮説Y → 課題Z (要件E)

1. 致命的建物の地理クラスタリング

2. CityGML LOD2 への発展 (3D 都市モデル本来の活用)

3. 時系列データとの結合 — 新築と廃絶

4. 建物 × 鉄道駅 — 駅前再開発の効果

5. 浸水以外のハザードへの拡張

"""), ] # キャッシュメモリ解放 city_data = None all_b = None # Footer は data-draft を使う DRAFT_FLAG = ' data-draft="1" data-stier="L"' html = render_lesson( num=14, title="L14 3D 都市モデル 7 市町比較 — 建物高度の地理学", tags=["3D都市モデル", "PLATEAU", "建物利用現況", "浸水", "DID", "都市階層", "致命率", "small multiples", "sjoin"], time=f"実行 {elapsed_sec:.0f} 秒", level="リテラシ", data_label="DoBoX 3D都市モデル #1316 #1302 #1187 #1188 #1189 #1443 #1444 / 建物利用現況 #1469 #512 #1323 #1483 #71 #1474 #1327 / 浸水 #1278 / DID", sections=sections, script_filename="lessons/L14_3d_city_model_compare.py", ) # data-draft フラグを body に注入 html = html.replace('', f'') OUT_HTML = LESSONS / "L14_3d_city_model_compare.html" OUT_HTML.write_text(html, encoding="utf-8") print(f" HTML 出力: {OUT_HTML}") print(f"\n=== 完了 (合計 {time.time()-t0:.1f}s) ===")