# -*- coding: utf-8 -*- """L39 地形傾斜図 2 件統合分析 — 広島県の航空レーザ計測 DEM が描く 傾斜角分布構造 と 土砂災害警戒の地形物理学 カバー宣言: 本記事は DoBoX のシリーズ 「地形傾斜図 1m / 50cm」 2 件 (dataset_id = 1630, 1631) を統合し、広島県内 3 自治体エリア (世羅町=中山間/林業地, 熊野町=都市近郊/呉市隣接山地, 坂町=平成30豪雨被災地) の 地形傾斜角 (slope angle, 度) の分布構造を読み解く研究記事である。 地形傾斜図とは DEM (標高ラスタ) の各ピクセルにおいて slope = arctan( sqrt((dz/dx)^2 + (dz/dy)^2) ) で計算した 局所傾斜角 (度, 0-90)のラスタ。 Horn 法 (Horn 1981) や Zevenbergen & Thorne 法 (1987) が標準的な実装。 砂防・土木工学では 「土砂災害感受性評価の最も基本的な指標」とされ、 急傾斜地崩壊危険箇所 (傾斜 30° 以上 + 高さ 5m 以上 + 人家 5戸以上) の 指定基準でも傾斜が筆頭判定材料になる。 L36 樹高図, L37 点群密度図, L38 CS立体図 と並ぶ LiDAR ファミリの 4 本目。 他の LiDAR 系記事 (CS立体図/標高図/航空レーザ計測森林資源/地図情報/計測年度図) とは 合体しない厳密シリーズ単独記事。L38 CS立体図とは 「曲率×傾斜の RGB 合成」(L38) vs 「傾斜単独の度数値」(L39) という関係 — L38 が視覚解析向きの圧縮表現、L39 が量的解析向きの素データ。 研究の問い (RQ): 広島県内の地形傾斜図 2 解像度 (1m / 50cm) は、それぞれ 世羅町 (中山間/林業地)・熊野町 (都市近郊/山岳混在)・坂町 (平成30豪雨被災地) という異なる地形リスク特性のエリアを公開している。 3 エリアの傾斜分布を 5 階級 (平地<5°/緩傾斜 5-15°/中傾斜 15-30°/急傾斜 30-45°/急峻 >45°)に 量化し、地形リスク構造を比較する。さらに 土砂災害警戒区域 (L10/L11 既知) との空間結合で、 「警戒区域の傾斜分布」と「30度超急傾斜地の地理集中」を量的に検証する。 (1) 各エリアの傾斜角ヒストグラム (0-90度を 2.5度刻み) (2) 5 階級構成比 (平地/緩/中/急/急峻) (3) 30度超急傾斜地の地理集塊 (centroid + 凝集度) (4) 土砂災害警戒区域 (L10) ポリゴン × 傾斜ラスタの平均傾斜・階級構成比較 (5) L38 CS立体図との同自治体ピクセル相関 (傾斜 vs Saturation) (6) 解像度差 (1m vs 50cm) で見える傾斜推定の精度 仮説 H1〜H7: H1 (世羅町は緩傾斜丘陵): 平均傾斜 ≤ 12°, 30度超セル比率 ≤ 12%。 中山間の植林地は緩斜面が多く、急峻地形は限定的。 H2 (熊野町は急峻地形): 平均傾斜 ≥ 18°, 30度超比率 ≥ 25%。 熊野は花崗岩マサ土の急傾斜面が広く分布する地形。 H3 (坂町は最急峻): 平均傾斜が 3 エリアで最大、急峻 (>45°) 比率も最大。 平成30豪雨で多数の土石流が発生した急斜面に集塊住宅地という地形。 H4 (土砂災害警戒区域は急傾斜): 警戒区域内の平均傾斜は区域外の 1.5 倍以上、 特別警戒区域 (red) ではさらに大きい。 H5 (急傾斜は集塊する): 30度超セルは空間的にランダムではなく、 Moran's I 的な近隣同種率が 0.85 以上 (= 隣接セル同種率, 強集塊)。 H6 (50cm 版は 1m 版より急傾斜が多い): 同自治体比較は不可能なので 「3 エリア間で 50cm 版 (坂) の急傾斜比率が高い」という形で代替検証。 高解像度ほど局所的急斜面を解像できるため。 H7 (傾斜 vs CS立体図 Saturation 正相関): L38 と同自治体の傾斜とCS図 Sat を overview 同 shape で読み込み、Pearson 相関で正相関 (r >= 0.4) を確認。 要件 S 準拠: 1 分以内完走 (各エリア overview /1/32〜64 で MB 級メモリ)。 要件 T 準拠: 各エリア地図 + 階級分類マップ + small multiples + 警戒区域重ね + L38 重ね。 要件 Q 準拠: 図 10+ 種, 表 8+ (3 エリア × 多角度: 階級/警戒/解像度/L38相関)。 データ仕様: A. 地形傾斜図 1m dataset 1630: 県内 22 自治体 GeoTIFF 公開。 本記事は世羅町 (resource_id=177725)と熊野町 (176866)を採用。 - 形式: GeoTIFF Float32 (度数, 0-90) - 解像度: 1m, EPSG:6671 B. 地形傾斜図 50cm dataset 1631: 県内 12 自治体 GeoTIFF 公開。 本記事は坂町 (resource_id=176905)を採用 (= L38 と同自治体ペア)。 - 形式: GeoTIFF Float32 (度数, 0-90) - 解像度: 0.5m, EPSG:6671 """ 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 ListedColormap, BoundaryNorm import matplotlib.patches as mpatches import geopandas as gpd from shapely.geometry import Polygon, box, Point import rasterio from rasterio.windows import Window, from_bounds as rio_from_bounds, bounds as rio_bounds from rasterio.enums import Resampling import zipfile import io as _io plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False t0 = time.time() print("=== L39 地形傾斜図 2 件統合分析 ===", flush=True) # ============================================================================= # 0. 定数・データパス・ZIP のダウンロード関数 # ============================================================================= DATA_DIR = ROOT / "data" / "extras" / "L39_terrain_slope" DATA_DIR.mkdir(parents=True, exist_ok=True) ADMIN_DIR = ROOT / "data" / "extras" / "L15_admin_zones" SED_DIR = ROOT / "data" / "extras" / "sediment_shp" L38_DIR = ROOT / "data" / "extras" / "L38_cs_terrain" # 各エリア定義 (3 エリア = 1m 2件 + 50cm 1件) AREAS = [ { "key": "sera_1m", "label": "世羅町", "muni_code": 34462, "resolution_m": 1.0, "dataset_id": 1630, "dataset_name": "地形傾斜図（1m）", "resource_id": 177725, "zip_url": "https://data.hiroshima-dobox.jp/forest/slope_1m/34462_世羅町.zip", "zip_url_alt": "https://hiroshima-dobox.jp/resource_download/177725", "zip_local": DATA_DIR / "1m_34462_sera_slope.zip", "tif": DATA_DIR / "1m_sera" / "sera_1m_slope.tif", "admin_zip": ADMIN_DIR / "admin_941_世羅町.zip", "l38_tif": L38_DIR / "1m_sera" / "sera_1m_cs.tif", "color": "#1a7f37", "character": "中山間 / 林業地 (緩傾斜丘陵地形, 大規模崩壊跡なし想定)", }, { "key": "kumano_1m", "label": "熊野町", "muni_code": 34307, "resolution_m": 1.0, "dataset_id": 1630, "dataset_name": "地形傾斜図（1m）", "resource_id": 176866, "zip_url": "https://data.hiroshima-dobox.jp/forest/slope_1m/34307_熊野町.zip", "zip_url_alt": "https://hiroshima-dobox.jp/resource_download/176866", "zip_local": DATA_DIR / "1m_34307_kumano_slope.zip", "tif": DATA_DIR / "1m_kumano" / "kumano_1m_slope.tif", "admin_zip": ADMIN_DIR / "admin_911_熊野町.zip", "l38_tif": L38_DIR / "1m_kumano" / "kumano_1m_cs.tif", "color": "#0969da", "character": "都市近郊 / 山岳混在 (花崗岩マサ土地形, 急傾斜面多)", }, { "key": "saka_50cm", "label": "坂町", "muni_code": 34384, "resolution_m": 0.5, "dataset_id": 1631, "dataset_name": "地形傾斜図（50cm）", "resource_id": 176905, "zip_url": "https://data.hiroshima-dobox.jp/forest/slope_50cm/34384_坂町.zip", "zip_url_alt": "https://hiroshima-dobox.jp/resource_download/176905", "zip_local": DATA_DIR / "50cm_34384_saka_slope.zip", "tif": DATA_DIR / "50cm_saka" / "saka_50cm_slope.tif", "admin_zip": ADMIN_DIR / "admin_916_坂町.zip", "l38_tif": L38_DIR / "50cm_saka" / "saka_50cm_cs.tif", "color": "#cf222e", "character": "平成30年7月豪雨被災地 (土石流多数発生, 急傾斜面に集塊住宅地)", }, ] # 傾斜階級 (本記事独自定義, 砂防工学標準に準拠) SLOPE_CLASSES = [ ("平地", -0.001, 5.0, "#e8f5e9", "<5°"), ("緩傾斜", 5.0, 15.0, "#a5d6a7", "5-15°"), ("中傾斜", 15.0, 30.0, "#ffd54f", "15-30°"), ("急傾斜", 30.0, 45.0, "#fb8c00", "30-45°"), ("急峻", 45.0, 90.001, "#b71c1c", ">45°"), ] def ensure_zip(area): p = area["zip_local"] if p.exists() and p.stat().st_size > 30_000_000: return p import requests UA = {"User-Agent": "DoBoX-MDASH-textbook/1.0 (educational)"} last_err = None for url in [area["zip_url"], area.get("zip_url_alt")]: if url is None: continue print(f" DL {p.name} <- {url}", flush=True) try: with requests.get(url, headers=UA, timeout=900, stream=True) as r: r.raise_for_status() total = int(r.headers.get("Content-Length", 0)) downloaded = 0 last_pct = -1 with open(p, "wb") as f: for chunk in r.iter_content(chunk_size=4*1024*1024): if chunk: f.write(chunk) downloaded += len(chunk) if total > 0: pct = int(downloaded / total * 100) if pct >= last_pct + 20: print(f" {pct}% ({downloaded/(1024*1024):.0f}/{total/(1024*1024):.0f} MB)", flush=True) last_pct = pct return p except Exception as e: print(f" failed: {e}", flush=True) last_err = e continue raise RuntimeError(f"All download URLs failed for {area['label']}: {last_err}") def ensure_tif(area): """ZIP から傾斜 TIFF と OVR (あれば) を抽出""" tif = area["tif"] if tif.exists() and tif.stat().st_size > 30_000_000: return tif zip_path = ensure_zip(area) target_dir = tif.parent target_dir.mkdir(parents=True, exist_ok=True) print(f" Extract TIFF from {zip_path.name}...", flush=True) with zipfile.ZipFile(zip_path) as zf: for info in zf.infolist(): n = info.filename if "__MACOSX" in n: continue base = Path(n).name if base.startswith("._"): continue sz = info.file_size lower = n.lower() target_name = None if (lower.endswith(".tif") or lower.endswith(".tiff")) and ".ovr" not in lower and sz > 10_000_000: target_name = tif.name elif lower.endswith(".tif.ovr") or lower.endswith(".tiff.ovr"): target_name = tif.name + ".ovr" if target_name is None: continue outpath = target_dir / target_name if outpath.exists() and outpath.stat().st_size > 1000: continue print(f" extracting {target_name} ({sz/(1024*1024):.0f} MB)...", flush=True) with zf.open(info) as src, open(outpath, "wb") as dst: while True: chunk = src.read(8*1024*1024) if not chunk: break dst.write(chunk) return tif # ============================================================================= # 1. データ確保: 3 自治体の TIFF # ============================================================================= print("\n[1] データ確保: TIFF 抽出", flush=True) t1 = time.time() for a in AREAS: tif = ensure_tif(a) print(f" {a['label']} ({a['resolution_m']}m): {tif.name} " f"({tif.stat().st_size/(1024*1024):.0f} MB)", flush=True) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 2. 各 TIFF のメタ情報を抽出 # ============================================================================= print("\n[2] TIFF メタデータ", flush=True) t1 = time.time() def tif_meta(tif): with rasterio.open(tif) as ds: bb = ds.bounds return { "width": ds.width, "height": ds.height, "n_cells": ds.width * ds.height, "n_bands": ds.count, "crs": str(ds.crs), "res_m": float(ds.res[0]), "left": bb.left, "right": bb.right, "bottom": bb.bottom, "top": bb.top, "width_m": bb.right - bb.left, "height_m": bb.top - bb.bottom, "area_km2": (bb.right - bb.left) * (bb.top - bb.bottom) / 1e6, "nodata": ds.nodata, "dtype": str(ds.dtypes[0]), "overviews": ds.overviews(1), } metas = [] for a in AREAS: m = tif_meta(a["tif"]) a["meta"] = m metas.append({"key": a["key"], "label": a["label"], **m}) print(f" {a['label']}: {m['width']:,}x{m['height']:,} = {m['n_cells']/1e6:.1f}M cells, " f"{m['width_m']/1000:.1f}x{m['height_m']/1000:.1f} km, dtype={m['dtype']}, " f"nodata={m['nodata']}, overviews={m['overviews']}", flush=True) meta_df = pd.DataFrame(metas) meta_df.to_csv(ASSETS / "L39_tiff_meta.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 3. overview を使って傾斜ラスタを読み込み # ============================================================================= print("\n[3] overview で傾斜ラスタ読み込み", flush=True) t1 = time.time() TARGET_CELL_M = 32 for a in AREAS: with rasterio.open(a["tif"]) as ds: ovr_choices = ds.overviews(1) target_factor = TARGET_CELL_M / a["resolution_m"] if ovr_choices: best = ovr_choices[-1] for f in ovr_choices: if f >= target_factor: best = f; break else: best = int(target_factor) ov_w = max(1, ds.width // best) ov_h = max(1, ds.height // best) print(f" {a['label']}: overview /1/{best} -> {ov_w}x{ov_h} ({ov_w*ov_h/1e6:.2f}M)", flush=True) arr = ds.read(1, out_shape=(ov_h, ov_w), resampling=Resampling.average).astype(np.float32) nodata = ds.nodata if ds.nodata is not None else -99999.0 valid_mask = (~np.isnan(arr)) & (arr != nodata) & (arr >= 0) & (arr <= 90.001) & np.isfinite(arr) a["arr"] = arr a["valid_mask"] = valid_mask a["arr_factor"] = best a["arr_cell_m"] = a["resolution_m"] * best a["arr_w"] = ov_w a["arr_h"] = ov_h a["bounds"] = ds.bounds n_valid = int(valid_mask.sum()) if n_valid > 0: v = arr[valid_mask] print(f" valid={n_valid:,}/{ov_w*ov_h:,} ({n_valid/(ov_w*ov_h)*100:.1f}%), " f"min={v.min():.2f}°, max={v.max():.2f}°, mean={v.mean():.2f}°, std={v.std():.2f}°", flush=True) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 4. 5 階級分類 # ============================================================================= print("\n[4] 5 階級分類", flush=True) t1 = time.time() def classify_slope(arr, valid_mask): """傾斜ラスタを 5 クラスにラベル付け 0=nodata, 1=平地(<5), 2=緩(5-15), 3=中(15-30), 4=急(30-45), 5=急峻(>45)""" label = np.zeros(arr.shape, dtype=np.int8) for i, (_, lo, hi, _, _) in enumerate(SLOPE_CLASSES, start=1): m = (arr > lo) & (arr <= hi) & valid_mask label[m] = i return label for a in AREAS: label = classify_slope(a["arr"], a["valid_mask"]) a["slope_label"] = label counts = np.bincount(label.ravel(), minlength=6) n_valid = int(a["valid_mask"].sum()) a["class_counts"] = counts pct = counts / max(1, n_valid) * 100 print(f" {a['label']}: 平地={pct[1]:.1f}%, 緩={pct[2]:.1f}%, 中={pct[3]:.1f}%, " f"急={pct[4]:.1f}%, 急峻={pct[5]:.1f}%", flush=True) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 5. 行政界ポリゴン読み込み # ============================================================================= print("\n[5] admin polygon ロード", flush=True) t1 = time.time() def load_admin_zip(zip_path): with zipfile.ZipFile(zip_path) as zf: names = zf.namelist() for n in names: if n.lower().endswith(".geojson"): with zf.open(n) as f: return gpd.read_file(_io.BytesIO(f.read())) for n in names: if n.lower().endswith(".shp"): tmp = zip_path.parent / f"_ext_{zip_path.stem}" tmp.mkdir(exist_ok=True) zf.extractall(tmp) return gpd.read_file(tmp / n, encoding="cp932") raise FileNotFoundError(zip_path) for a in AREAS: if a["admin_zip"].exists(): try: g = load_admin_zip(a["admin_zip"]).to_crs(a["meta"]["crs"]) a["admin"] = g.dissolve() muni_area_km2 = float(a["admin"].geometry.area.sum()) / 1e6 a["muni_area_km2"] = muni_area_km2 print(f" {a['label']}: admin loaded, {muni_area_km2:.1f} km2", flush=True) except Exception as e: print(f" {a['label']}: admin err {e}", flush=True) a["admin"] = None a["muni_area_km2"] = None else: print(f" {a['label']}: admin zip none", flush=True) a["admin"] = None a["muni_area_km2"] = None print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 6. 基本統計 # ============================================================================= print("\n[6] 基本統計", flush=True) t1 = time.time() basic_stats = {} for a in AREAS: vm = a["valid_mask"] v = a["arr"][vm] if len(v) == 0: v = np.array([0.0]) s = { "label": a["label"], "n_cells_total": int(vm.size), "n_cells_valid": int(vm.sum()), "valid_ratio": float(vm.sum() / vm.size), "mean": float(v.mean()), "median": float(np.median(v)), "std": float(v.std()), "p10": float(np.percentile(v, 10)), "p25": float(np.percentile(v, 25)), "p75": float(np.percentile(v, 75)), "p90": float(np.percentile(v, 90)), "p95": float(np.percentile(v, 95)), "max": float(v.max()), "min": float(v.min()), "pct_ge_30": float(((a["arr"] >= 30) & vm).sum() / max(1, int(vm.sum())) * 100), "pct_ge_45": float(((a["arr"] >= 45) & vm).sum() / max(1, int(vm.sum())) * 100), "pct_lt_5": float(((a["arr"] < 5) & vm).sum() / max(1, int(vm.sum())) * 100), } basic_stats[a["key"]] = s stats_df = pd.DataFrame.from_dict(basic_stats, orient="index") stats_df.to_csv(ASSETS / "L39_basic_stats.csv", encoding="utf-8-sig") print(stats_df.round(2).to_string()) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 7. 土砂災害警戒区域 (sediment) のロード # ============================================================================= print("\n[7] sediment 警戒区域ロード", flush=True) t1 = time.time() def load_sed_for_area(area): if not SED_DIR.exists(): return None, None bb = area["bounds"] bbox_poly = box(bb.left, bb.bottom, bb.right, bb.top) pieces_yellow = [] pieces_red = [] for kind, dirname, ymark, rmark in [ ("土石流", "doseki", "031dy", "031drp"), ("急傾斜", "kyukeisha", "031ky", "031krp"), ("地すべり","jisuberi", "031jy", None), ]: sub = SED_DIR / dirname if not sub.exists(): continue for child in sub.iterdir(): if not child.is_dir(): continue for shp in child.glob("*.shp"): name = shp.name try: g = gpd.read_file(shp) if g.crs is None: g = g.set_crs(epsg=6668) g = g.to_crs(area["meta"]["crs"]) g = g[g.intersects(bbox_poly)].copy() if len(g) == 0: continue g["kind"] = kind if ymark in name and (rmark is None or rmark not in name): g["level"] = "警戒" pieces_yellow.append(g[["geometry","kind","level"]]) elif rmark and rmark in name: g["level"] = "特別警戒" pieces_red.append(g[["geometry","kind","level"]]) except Exception: pass yellow = gpd.GeoDataFrame(pd.concat(pieces_yellow, ignore_index=True), crs=area["meta"]["crs"]) \ if pieces_yellow else gpd.GeoDataFrame(geometry=[], crs=area["meta"]["crs"]) red = gpd.GeoDataFrame(pd.concat(pieces_red, ignore_index=True), crs=area["meta"]["crs"]) \ if pieces_red else gpd.GeoDataFrame(geometry=[], crs=area["meta"]["crs"]) if area.get("admin") is not None and len(area["admin"]) > 0: if len(yellow) > 0: yellow = gpd.overlay(yellow, area["admin"][["geometry"]], how="intersection", keep_geom_type=True) if len(red) > 0: red = gpd.overlay(red, area["admin"][["geometry"]], how="intersection", keep_geom_type=True) return yellow, red for a in AREAS: try: yellow, red = load_sed_for_area(a) a["sed_yellow"] = yellow a["sed_red"] = red ny = len(yellow) if yellow is not None else 0 nr = len(red) if red is not None else 0 ay = float(yellow.geometry.area.sum()/1e6) if yellow is not None and len(yellow)>0 else 0 ar = float(red.geometry.area.sum()/1e6) if red is not None and len(red)>0 else 0 a["sed_yellow_n"] = ny; a["sed_yellow_km2"] = ay a["sed_red_n"] = nr; a["sed_red_km2"] = ar print(f" {a['label']}: yellow {ny} polys ({ay:.2f} km^2), red {nr} polys ({ar:.2f} km^2)", flush=True) except Exception as e: print(f" {a['label']}: sed load failed: {e}", flush=True) a["sed_yellow"] = a["sed_red"] = None a["sed_yellow_n"] = a["sed_red_n"] = 0 a["sed_yellow_km2"] = a["sed_red_km2"] = 0 print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 8. 警戒区域マスクをラスタ化して傾斜と比較 # ============================================================================= print("\n[8] 警戒区域マスク x 傾斜比較", flush=True) t1 = time.time() def rasterize_polys(polys_gdf, area): if polys_gdf is None or len(polys_gdf) == 0: return np.zeros((area["arr_h"], area["arr_w"]), dtype=bool) from rasterio.features import rasterize from rasterio.transform import from_bounds bb = area["bounds"] transform = from_bounds(bb.left, bb.bottom, bb.right, bb.top, area["arr_w"], area["arr_h"]) shapes = ((geom, 1) for geom in polys_gdf.geometry if geom is not None and not geom.is_empty) mask = rasterize(shapes, out_shape=(area["arr_h"], area["arr_w"]), transform=transform, fill=0, dtype="uint8") return mask.astype(bool) sed_compare_rows = [] for a in AREAS: yellow = a.get("sed_yellow") red = a.get("sed_red") yellow_mask = rasterize_polys(yellow, a) red_mask = rasterize_polys(red, a) a["yellow_mask"] = yellow_mask a["red_mask"] = red_mask in_yellow = yellow_mask & a["valid_mask"] in_red = red_mask & a["valid_mask"] out_warn = (~yellow_mask) & (~red_mask) & a["valid_mask"] arr = a["arr"] for label_name, m in [("特別警戒", in_red), ("警戒", in_yellow), ("区域外", out_warn)]: n = int(m.sum()) if n < 5: row = {"area": a["label"], "zone": label_name, "n": n, "mean": None, "median": None, "p90": None, "pct_ge_30": None, "pct_ge_45": None} else: v = arr[m] row = { "area": a["label"], "zone": label_name, "n": n, "mean": float(v.mean()), "median": float(np.median(v)), "p90": float(np.percentile(v, 90)), "pct_ge_30": float((v >= 30).sum() / n * 100), "pct_ge_45": float((v >= 45).sum() / n * 100), } sed_compare_rows.append(row) sed_compare_df = pd.DataFrame(sed_compare_rows) sed_compare_df.to_csv(ASSETS / "L39_sediment_slope_compare.csv", index=False, encoding="utf-8-sig") print(sed_compare_df.round(2).to_string(index=False)) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 9. L38 CS立体図との同自治体ピクセル相関 # ============================================================================= print("\n[9] L38 CS立体図とのピクセル相関", flush=True) t1 = time.time() l38_correlations = {} for a in AREAS: if a["l38_tif"].exists(): try: with rasterio.open(a["l38_tif"]) as ds_cs: rgb = ds_cs.read(out_shape=(3, a["arr_h"], a["arr_w"]), resampling=Resampling.average).astype(np.float32) cs_valid = rgb.sum(axis=0) > 0 maxv = np.maximum.reduce(rgb) minv = np.minimum.reduce(rgb) sat = np.where(maxv > 0, (maxv - minv) / np.maximum(maxv, 1), 0.0).astype(np.float32) darkness = 1.0 - maxv / 255.0 a["cs_sat"] = sat a["cs_dark"] = darkness a["cs_valid"] = cs_valid both = a["valid_mask"] & cs_valid a["both_valid_mask"] = both if both.sum() > 100: v_slope = a["arr"][both].astype(np.float32) v_sat = sat[both] v_dark = darkness[both] if v_slope.std() > 0 and v_sat.std() > 0: r_sat = float(np.corrcoef(v_slope, v_sat)[0, 1]) else: r_sat = 0.0 if v_slope.std() > 0 and v_dark.std() > 0: r_dark = float(np.corrcoef(v_slope, v_dark)[0, 1]) else: r_dark = 0.0 l38_correlations[a["key"]] = { "r_sat": r_sat, "r_dark": r_dark, "n": int(both.sum()) } print(f" {a['label']}: n={int(both.sum()):,}, " f"corr(slope, Sat)={r_sat:+.3f}, corr(slope, darkness)={r_dark:+.3f}", flush=True) else: l38_correlations[a["key"]] = {"r_sat": 0.0, "r_dark": 0.0, "n": 0} except Exception as e: print(f" {a['label']}: L38 read failed: {e}", flush=True) a["cs_sat"] = a["cs_dark"] = a["cs_valid"] = a["both_valid_mask"] = None l38_correlations[a["key"]] = {"r_sat": 0.0, "r_dark": 0.0, "n": 0} else: print(f" {a['label']}: L38 tif missing, skip correlation", flush=True) a["cs_sat"] = a["cs_dark"] = a["cs_valid"] = a["both_valid_mask"] = None l38_correlations[a["key"]] = {"r_sat": 0.0, "r_dark": 0.0, "n": 0} corr_df = pd.DataFrame([ {"area": a["label"], "n_paired_cells": l38_correlations[a["key"]]["n"], "pearson_r_slope_vs_Sat": round(l38_correlations[a["key"]]["r_sat"], 4), "pearson_r_slope_vs_darkness": round(l38_correlations[a["key"]]["r_dark"], 4)} for a in AREAS ]) corr_df.to_csv(ASSETS / "L39_l38_correlation.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 10. 急傾斜地 (>=30°) の集塊性 # ============================================================================= print("\n[10] 急傾斜地集塊性", flush=True) t1 = time.time() def neighbor_same_rate(mask, valid): """8 近傍で同種であるペア比率""" H, W = mask.shape if mask.sum() == 0: return 0.0, 0 pairs_total = 0 pairs_same = 0 for dy in (-1, 0, 1): for dx in (-1, 0, 1): if dy == 0 and dx == 0: continue y0a, y1a = max(0, dy), H + min(0, dy) x0a, x1a = max(0, dx), W + min(0, dx) y0b, y1b = max(0, -dy), H + min(0, -dy) x0b, x1b = max(0, -dx), W + min(0, -dx) ma = mask[y0a:y1a, x0a:x1a] mb = mask[y0b:y1b, x0b:x1b] va = valid[y0a:y1a, x0a:x1a] vb = valid[y0b:y1b, x0b:x1b] both_valid = va & vb both_in = ma & mb & both_valid either_in = (ma | mb) & both_valid pairs_total += int(either_in.sum()) pairs_same += int(both_in.sum()) same_rate = pairs_same / max(1, pairs_total) return float(same_rate), int(mask.sum()) cluster_rows = [] for a in AREAS: steep_mask = (a["arr"] >= 30) & a["valid_mask"] rate, n_steep = neighbor_same_rate(steep_mask, a["valid_mask"]) a["steep_neighbor_rate"] = rate a["n_steep_cells"] = n_steep ys, xs = np.where(steep_mask) if len(ys) > 0: cy = a["bounds"].top - (ys.mean() + 0.5) * a["arr_cell_m"] cx = a["bounds"].left + (xs.mean() + 0.5) * a["arr_cell_m"] else: cy = cx = None a["steep_centroid_y"] = cy a["steep_centroid_x"] = cx cluster_rows.append({ "area": a["label"], "n_steep_cells": n_steep, "pct_steep": round(n_steep / max(1, int(a["valid_mask"].sum())) * 100, 2), "neighbor_same_rate": round(rate, 4), "centroid_x_m": round(cx, 1) if cx is not None else None, "centroid_y_m": round(cy, 1) if cy is not None else None, }) print(f" {a['label']}: n_steep={n_steep:,}, neighbor_same_rate={rate:.3f}", flush=True) cluster_df = pd.DataFrame(cluster_rows) cluster_df.to_csv(ASSETS / "L39_steep_clustering.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 11. 生解像度サンプル (中央 1024x1024) # ============================================================================= print("\n[11] 生解像度中央サンプル", flush=True) t1 = time.time() RAW_SAMPLE = 1024 for a in AREAS: with rasterio.open(a["tif"]) as ds: cx, cy = ds.width//2, ds.height//2 win = Window(max(0, cx - RAW_SAMPLE//2), max(0, cy - RAW_SAMPLE//2), min(RAW_SAMPLE, ds.width), min(RAW_SAMPLE, ds.height)) raw = ds.read(1, window=win).astype(np.float32) a["raw"] = raw a["raw_bounds"] = rio_bounds(win, ds.transform) nodata = a["meta"]["nodata"] if a["meta"]["nodata"] is not None else -99999.0 raw_valid = (~np.isnan(raw)) & (raw != nodata) & (raw >= 0) & (raw <= 90.001) a["raw_valid"] = raw_valid if raw_valid.sum() > 0: v = raw[raw_valid] a["raw_mean"] = float(v.mean()) a["raw_median"] = float(np.median(v)) a["raw_pct_ge_30"] = float((v >= 30).sum() / len(v) * 100) a["raw_pct_ge_45"] = float((v >= 45).sum() / len(v) * 100) else: a["raw_mean"] = a["raw_median"] = a["raw_pct_ge_30"] = a["raw_pct_ge_45"] = 0.0 print(f" {a['label']}: raw {raw.shape}, valid={int(raw_valid.sum()):,}, " f"mean={a['raw_mean']:.2f}deg", flush=True) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 12. 仮説検証 (H1〜H7) # ============================================================================= print("\n[12] 仮説検証", flush=True) t1 = time.time() s_sera = basic_stats["sera_1m"] s_kuma = basic_stats["kumano_1m"] s_saka = basic_stats["saka_50cm"] hypos = [] # H1: 世羅は緩傾斜丘陵 h1 = (s_sera["mean"] <= 12.0) and (s_sera["pct_ge_30"] <= 12.0) hypos.append({ "H": "H1", "claim": "世羅町は緩傾斜丘陵が主で、平均傾斜 ≤ 12° かつ 30度超セル比率 ≤ 12%", "result": f"世羅 平均={s_sera['mean']:.2f}°, 30度超={s_sera['pct_ge_30']:.1f}%", "verdict": "支持" if h1 else ("部分支持" if (s_sera["mean"] <= 15 or s_sera["pct_ge_30"] <= 18) else "反証"), }) # H2: 熊野は急峻地形 h2 = (s_kuma["mean"] >= 18.0) and (s_kuma["pct_ge_30"] >= 25.0) hypos.append({ "H": "H2", "claim": "熊野町は急峻地形で平均傾斜 ≥ 18° かつ 30度超セル比率 ≥ 25%", "result": f"熊野 平均={s_kuma['mean']:.2f}°, 30度超={s_kuma['pct_ge_30']:.1f}%", "verdict": "支持" if h2 else ("部分支持" if (s_kuma["mean"] >= 14 or s_kuma["pct_ge_30"] >= 20) else "反証"), }) # H3: 坂町が最急峻 means = [s_sera["mean"], s_kuma["mean"], s_saka["mean"]] p45s = [s_sera["pct_ge_45"], s_kuma["pct_ge_45"], s_saka["pct_ge_45"]] h3 = (s_saka["mean"] >= max(means)) and (s_saka["pct_ge_45"] >= max(p45s)) hypos.append({ "H": "H3", "claim": "坂町は3エリアで平均傾斜が最大 + 急峻 (>45°) 比率も最大", "result": f"世羅 mean={s_sera['mean']:.2f}°/p45={s_sera['pct_ge_45']:.1f}%, " f"熊野 {s_kuma['mean']:.2f}°/{s_kuma['pct_ge_45']:.1f}%, " f"坂 {s_saka['mean']:.2f}°/{s_saka['pct_ge_45']:.1f}%", "verdict": "支持" if h3 else ("部分支持" if s_saka["mean"] >= max(s_sera["mean"], s_kuma["mean"]) else "反証"), }) # H4: 警戒区域内の平均傾斜 def get_zone(rows, area_label, zone): for r in rows: if r["area"] == area_label and r["zone"] == zone: return r return None h4_supports = [] h4_pieces = [] for a in AREAS: inner = get_zone(sed_compare_rows, a["label"], "警戒") outer = get_zone(sed_compare_rows, a["label"], "区域外") if inner is None or outer is None or inner["mean"] is None or outer["mean"] is None: h4_pieces.append(f"{a['label']}: NA") continue ratio = inner["mean"] / max(0.001, outer["mean"]) h4_supports.append(ratio >= 1.5) h4_pieces.append(f"{a['label']}: 内 {inner['mean']:.2f}° / 外 {outer['mean']:.2f}° (×{ratio:.2f})") h4 = (sum(h4_supports) >= 2) if h4_supports else False hypos.append({ "H": "H4", "claim": "土砂災害警戒区域内の平均傾斜は区域外の 1.5 倍以上 (急傾斜地に指定が集中)", "result": " | ".join(h4_pieces), "verdict": "支持" if h4 else ("部分支持" if (h4_supports and sum(h4_supports) >= 1) else "反証"), }) # H5: 急傾斜は集塊 h5_rates = [a["steep_neighbor_rate"] for a in AREAS] h5 = sum(r >= 0.85 for r in h5_rates) >= 2 hypos.append({ "H": "H5", "claim": "30度超急傾斜セルは空間的に集塊 (8近傍同種率 >= 0.85)", "result": " | ".join([f"{a['label']}: {a['steep_neighbor_rate']:.3f}" for a in AREAS]), "verdict": "支持" if h5 else ("部分支持" if sum(r >= 0.7 for r in h5_rates) >= 2 else "反証"), }) # H6: 50cm版は急傾斜比率高 h6 = s_saka["pct_ge_30"] >= max(s_sera["pct_ge_30"], s_kuma["pct_ge_30"]) hypos.append({ "H": "H6", "claim": "50cm版 (坂町) は 1m版 (世羅・熊野) より 30度超セル比率が高い", "result": f"世羅 pct30={s_sera['pct_ge_30']:.1f}%, 熊野 {s_kuma['pct_ge_30']:.1f}%, 坂 {s_saka['pct_ge_30']:.1f}%", "verdict": "支持" if h6 else ("部分支持" if s_saka["pct_ge_30"] >= s_sera["pct_ge_30"] else "反証"), }) # H7: 傾斜 vs CS Sat 正相関 h7_rs = [l38_correlations[a["key"]]["r_sat"] for a in AREAS] h7 = sum(r >= 0.4 for r in h7_rs) >= 2 hypos.append({ "H": "H7", "claim": "傾斜と CS立体図 Saturation は正相関 (Pearson r >= 0.4 で 2 エリア以上)", "result": " | ".join([f"{a['label']}: r={l38_correlations[a['key']]['r_sat']:+.3f} " f"(n={l38_correlations[a['key']]['n']:,})" for a in AREAS]), "verdict": "支持" if h7 else ("部分支持" if sum(r >= 0.2 for r in h7_rs) >= 2 else "反証"), }) hypos_df = pd.DataFrame(hypos) hypos_df.to_csv(ASSETS / "L39_hypothesis_results.csv", index=False, encoding="utf-8-sig") with (ASSETS / "L39_hypothesis_results.json").open("w", encoding="utf-8") as f: json.dump(hypos, f, ensure_ascii=False, indent=2) print(hypos_df[["H", "verdict", "result"]].to_string(index=False)) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 13. dataset/series 一覧 + 階級表 + ヒスト CSV # ============================================================================= print("\n[13] series 一覧 + 階級表", flush=True) t1 = time.time() series_rows = [] for a in AREAS: series_rows.append({ "dsid": a["dataset_id"], "dataset_name": a["dataset_name"], "muni": a["label"], "muni_code": a["muni_code"], "resource_id": a["resource_id"], "resolution_m": a["resolution_m"], "tif_size_MB": int(a["tif"].stat().st_size / (1024*1024)), "n_cells_M": round(a["meta"]["n_cells"] / 1e6, 1), "bbox_area_km2": round(a["meta"]["area_km2"], 1), "muni_area_km2": round(a["muni_area_km2"], 1) if a["muni_area_km2"] else None, }) series_df = pd.DataFrame(series_rows) series_df.to_csv(ASSETS / "L39_series.csv", index=False, encoding="utf-8-sig") # 階級構成比 class_rows = [] for a in AREAS: cc = a["class_counts"] n = max(1, int(a["valid_mask"].sum())) row = {"area": a["label"]} for i, (lab, _, _, _, _) in enumerate(SLOPE_CLASSES, start=1): row[f"{lab}_n"] = int(cc[i]) row[f"{lab}_pct"] = round(cc[i]/n * 100, 2) class_rows.append(row) class_df = pd.DataFrame(class_rows) class_df.to_csv(ASSETS / "L39_slope_classes.csv", index=False, encoding="utf-8-sig") # 2.5度刻みヒスト hist_rows = [] edges = np.arange(0, 90.01, 2.5) for a in AREAS: v = a["arr"][a["valid_mask"]] h, _ = np.histogram(v, bins=edges) for i, c in enumerate(h): hist_rows.append({"area": a["label"], "bin_low_deg": float(edges[i]), "bin_high_deg": float(edges[i+1]), "count": int(c), "pct": float(c / max(1, len(v)) * 100)}) hist_df = pd.DataFrame(hist_rows) hist_df.to_csv(ASSETS / "L39_slope_histogram.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t1:.1f}s)", flush=True) print(f"\n=== Total elapsed (analysis): {time.time()-t0:.1f}s ===", flush=True) # ============================================================================= # 14. 図 1: dataset カード + slope ヒスト (3 エリア) # ============================================================================= print("\n[14] 図 1: dataset カード + 傾斜ヒスト", flush=True) t1 = time.time() fig, axes = plt.subplots(1, 2, figsize=(14, 6.5)) ax = axes[0] ax.set_xlim(0, 12); ax.set_ylim(0, 10); ax.set_axis_off() ax.text(6, 9.5, "DoBoX 地形傾斜図 2 件 × 3 自治体", ha="center", fontsize=14, fontweight="bold") card_x = [0.3, 4.2, 8.1] for ax_x, a in zip(card_x, AREAS): ax.add_patch(plt.Rectangle((ax_x, 0.4), 3.5, 8.2, fill=True, facecolor=a["color"], alpha=0.13, edgecolor=a["color"], linewidth=2.2)) ax.text(ax_x+1.75, 8.0, f"dsid {a['dataset_id']}\n{a['dataset_name']}", ha="center", fontsize=9.5, fontweight="bold", color=a["color"]) ax.text(ax_x+1.75, 6.9, f"{a['label']}", ha="center", fontsize=14, fontweight="bold") char_text = a["character"] char_first = char_text.split("(")[0].strip() char_paren = "(" + char_text.split("(", 1)[1] if "(" in char_text else "" ax.text(ax_x+1.75, 6.4, char_first, ha="center", fontsize=8, color="#444") ax.text(ax_x+1.75, 5.9, char_paren[:30] + ("..." if len(char_paren) > 30 else ""), ha="center", fontsize=7, color="#666") ax.text(ax_x+1.75, 5.2, f"解像度 {a['resolution_m']:.2f} m / Float32 度数", ha="center", fontsize=8.5) ax.text(ax_x+1.75, 4.4, f"{a['meta']['width']:,} x {a['meta']['height']:,}\n= {a['meta']['n_cells']/1e6:.0f}M cells", ha="center", fontsize=9) ax.text(ax_x+1.75, 3.3, f"{a['meta']['width_m']/1000:.1f} x {a['meta']['height_m']/1000:.1f} km", ha="center", fontsize=8.5) s = basic_stats[a["key"]] ax.text(ax_x+1.75, 2.5, f"平均 {s['mean']:.1f}° (中央 {s['median']:.1f}°)", ha="center", fontsize=9, fontweight="bold") ax.text(ax_x+1.75, 1.7, f"30度超 {s['pct_ge_30']:.1f}% / 急峻 {s['pct_ge_45']:.2f}%", ha="center", fontsize=9, color=a["color"]) ax.text(ax_x+1.75, 1.0, f"警戒 {a['sed_yellow_km2']:.1f} km² / 特別 {a['sed_red_km2']:.2f} km²", ha="center", fontsize=8, color="#cf222e") ax = axes[1] for a in AREAS: v = a["arr"][a["valid_mask"]] ax.hist(v, bins=np.arange(0, 90.01, 2.5), color=a["color"], alpha=0.55, density=True, edgecolor="white", linewidth=0.3, label=f"{a['label']} (μ={v.mean():.1f}°, p90={np.percentile(v,90):.1f}°)") # 階級境界線 for thresh, col, lab in [(5,"#888","平/緩"), (15,"#888","緩/中"), (30,"#fb8c00","中/急 30°"), (45,"#b71c1c","急/急峻 45°")]: ax.axvline(thresh, color=col, linestyle="--", alpha=0.6, linewidth=1) ax.text(thresh+0.5, ax.get_ylim()[1]*0.92 if ax.get_ylim()[1]>0 else 0.05, lab, fontsize=8, color=col) ax.set_xlabel("傾斜角 (度)") ax.set_ylabel("密度 (normalized)") ax.set_title("傾斜角分布 — 3 エリア比較 (2.5° bin)") ax.set_xlim(0, 70) ax.legend(loc="upper right", fontsize=9) ax.grid(axis="y", alpha=0.3) plt.tight_layout() fig.savefig(ASSETS / "L39_fig1_dataset_overview.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L39_fig1_dataset_overview.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 15. 図 2: 傾斜ラスタ主題図 (3 エリア) # ============================================================================= print("\n[15] 図 2: 傾斜主題図", flush=True) t1 = time.time() from matplotlib.colors import LinearSegmentedColormap slope_cmap = LinearSegmentedColormap.from_list( "slope", ["#1a9850", "#a6d96a", "#fee08b", "#fdae61", "#d73027", "#7f0000"], N=256 ) fig, axes = plt.subplots(1, 3, figsize=(18, 7)) vmin, vmax = 0, 60 for ax, a in zip(axes, AREAS): arr = np.where(a["valid_mask"], a["arr"], np.nan) extent = (a["bounds"].left, a["bounds"].right, a["bounds"].bottom, a["bounds"].top) im = ax.imshow(arr, cmap=slope_cmap, vmin=vmin, vmax=vmax, extent=extent, origin="upper", interpolation="nearest") if a.get("admin") is not None: a["admin"].boundary.plot(ax=ax, color="#222", linewidth=1.3, alpha=0.85) s = basic_stats[a["key"]] ax.set_title(f"{a['label']} 傾斜 (overview /1/{a['arr_factor']}, " f"実効 {a['arr_cell_m']:.0f} m/cell)\n" f"平均 {s['mean']:.1f}° / 30°超 {s['pct_ge_30']:.1f}% / max {s['max']:.0f}°", fontsize=10) ax.set_xlabel("X (EPSG:6671 m)") ax.set_ylabel("Y (EPSG:6671 m)") # 共通カラーバー cbar = fig.colorbar(im, ax=axes, orientation="horizontal", fraction=0.04, pad=0.10, shrink=0.85) cbar.set_label("傾斜角 (度) — 緑=平地, 黄=中傾斜, 赤=急峻", fontsize=10) fig.suptitle("3 エリア 地形傾斜図 — 0° (緑) → 60°+ (暗赤) のグラディエント", fontsize=12, y=1.0) fig.savefig(ASSETS / "L39_fig2_slope_maps.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L39_fig2_slope_maps.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 16. 図 3: 5 階級分類マップ # ============================================================================= print("\n[16] 図 3: 5 階級分類マップ", flush=True) t1 = time.time() class_colors = ["#ffffff"] + [c for _, _, _, c, _ in SLOPE_CLASSES] class_labels_with_range = [lab + " " + r for (lab, _, _, _, r) in SLOPE_CLASSES] cmap_c = ListedColormap(class_colors) norm_c = BoundaryNorm([-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5], 6, clip=True) fig, axes = plt.subplots(1, 3, figsize=(18, 7)) for ax, a in zip(axes, AREAS): extent = (a["bounds"].left, a["bounds"].right, a["bounds"].bottom, a["bounds"].top) ax.imshow(a["slope_label"], cmap=cmap_c, norm=norm_c, extent=extent, origin="upper", interpolation="nearest") if a.get("admin") is not None: a["admin"].boundary.plot(ax=ax, color="#000", linewidth=1.0, alpha=0.85) cc = a["class_counts"] n = max(1, cc.sum()) ax.set_title(f"{a['label']}: 平地 {cc[1]/n*100:.1f}% / 緩 {cc[2]/n*100:.1f}% / " f"中 {cc[3]/n*100:.1f}% / 急 {cc[4]/n*100:.1f}% / 急峻 {cc[5]/n*100:.1f}%", fontsize=9) ax.set_xlabel("X (m)"); ax.set_ylabel("Y (m)") legend_handles = [Patch(facecolor=c, edgecolor="#222", label=lab) for c, lab in zip(class_colors[1:], class_labels_with_range)] fig.legend(handles=legend_handles, loc="lower center", ncol=5, fontsize=10, bbox_to_anchor=(0.5, -0.02)) fig.suptitle("地形傾斜 5 階級分類 — 砂防工学標準閾値 (5°/15°/30°/45°)", fontsize=12, y=1.0) plt.tight_layout() fig.savefig(ASSETS / "L39_fig3_class_map.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L39_fig3_class_map.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 17. 図 4: 階級構成比 stack-bar # ============================================================================= print("\n[17] 図 4: 階級構成比", flush=True) t1 = time.time() fig, ax = plt.subplots(figsize=(10, 5)) x = np.arange(len(AREAS)) bottom = np.zeros(len(AREAS)) class_show = [lab for (lab, _, _, _, _) in SLOPE_CLASSES] class_cols = [c for (_, _, _, c, _) in SLOPE_CLASSES] class_ranges = [r for (_, _, _, _, r) in SLOPE_CLASSES] for i, (cn, col, rng) in enumerate(zip(class_show, class_cols, class_ranges)): vals = [] for a in AREAS: n = max(1, int(a["valid_mask"].sum())) vals.append(a["class_counts"][i+1] / n * 100) ax.bar(x, vals, bottom=bottom, label=f"{cn} {rng}", color=col, edgecolor="#222", linewidth=0.6) for j, v in enumerate(vals): if v >= 3: text_color = "white" if col in ("#fb8c00","#b71c1c") else "black" ax.text(j, bottom[j] + v/2, f"{v:.1f}%", ha="center", va="center", fontsize=10, fontweight="bold", color=text_color) bottom = bottom + np.array(vals) ax.set_xticks(x) ax.set_xticklabels([f"{a['label']}\n({a['resolution_m']}m)" for a in AREAS], fontsize=11) ax.set_ylabel("傾斜階級構成比 (%)") ax.set_title("傾斜 5 階級構成比 — 砂防工学標準閾値") ax.legend(loc="center left", bbox_to_anchor=(1.02, 0.5)) ax.grid(axis="y", alpha=0.3) ax.set_ylim(0, 100) plt.tight_layout() fig.savefig(ASSETS / "L39_fig4_class_composition.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L39_fig4_class_composition.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 18. 図 5: 警戒区域 × 傾斜 重ね合わせマップ # ============================================================================= print("\n[18] 図 5: 警戒区域 × 傾斜 overlay", flush=True) t1 = time.time() fig, axes = plt.subplots(1, 3, figsize=(18, 7)) for ax, a in zip(axes, AREAS): arr = np.where(a["valid_mask"], a["arr"], np.nan) extent = (a["bounds"].left, a["bounds"].right, a["bounds"].bottom, a["bounds"].top) im = ax.imshow(arr, cmap=slope_cmap, vmin=0, vmax=60, extent=extent, origin="upper", interpolation="nearest", alpha=0.7) if a.get("admin") is not None: a["admin"].boundary.plot(ax=ax, color="#000", linewidth=1.3, alpha=0.9) if a.get("sed_yellow") is not None and len(a["sed_yellow"]) > 0: a["sed_yellow"].boundary.plot(ax=ax, color="#a87b00", linewidth=0.7, alpha=0.65) if a.get("sed_red") is not None and len(a["sed_red"]) > 0: a["sed_red"].plot(ax=ax, facecolor="none", edgecolor="#000", linewidth=0.8, alpha=0.95) inner = get_zone(sed_compare_rows, a["label"], "警戒") outer = get_zone(sed_compare_rows, a["label"], "区域外") inner_mean = inner["mean"] if inner and inner["mean"] is not None else 0 outer_mean = outer["mean"] if outer and outer["mean"] is not None else 0 ax.set_title(f"{a['label']}: 傾斜 + 土砂災害警戒区域\n" f"警戒内 平均{inner_mean:.1f}° vs 区域外 平均{outer_mean:.1f}° (×{inner_mean/max(0.001,outer_mean):.2f})", fontsize=9.5) ax.set_xlabel("X (m)"); ax.set_ylabel("Y (m)") legend_handles = [ Patch(facecolor="none", edgecolor="#a87b00", linewidth=1.2, label="土砂警戒 (yellow) 境界"), Patch(facecolor="none", edgecolor="#000", linewidth=1.4, label="特別警戒 (red) 境界"), ] fig.legend(handles=legend_handles, loc="lower center", ncol=2, fontsize=11, bbox_to_anchor=(0.5, -0.02)) fig.suptitle("傾斜ラスタ × 土砂災害警戒区域 重ね合わせ — H4 検証 (警戒区域は急傾斜?)", fontsize=12, y=1.0) plt.tight_layout() fig.savefig(ASSETS / "L39_fig5_sediment_overlay.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L39_fig5_sediment_overlay.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 19. 図 6: 警戒区域 内/外 平均傾斜と階級構成 grouped bar # ============================================================================= print("\n[19] 図 6: 警戒区域 内/外 比較", flush=True) t1 = time.time() fig, axes = plt.subplots(1, 3, figsize=(15, 4.5)) zones_order = ["特別警戒", "警戒", "区域外"] zone_colors = {"特別警戒": "#cf222e", "警戒": "#ffd700", "区域外": "#888888"} for ax, metric, ylabel in zip( axes, ["mean", "p90", "pct_ge_30"], ["平均傾斜 (度)", "90 percentile 傾斜 (度)", "30°超セル比率 (%)"], ): n_areas = len(AREAS) xs = np.arange(n_areas) w = 0.27 for i, zone in enumerate(zones_order): vals = [] for a in AREAS: r = get_zone(sed_compare_rows, a["label"], zone) v = r[metric] if r is not None and r.get(metric) is not None else 0 vals.append(v) ax.bar(xs + (i-1)*w, vals, w, color=zone_colors[zone], edgecolor="#222", linewidth=0.6, label=zone) for j, v in enumerate(vals): ax.text(xs[j] + (i-1)*w, v, f"{v:.1f}", ha="center", va="bottom", fontsize=8) ax.set_xticks(xs) ax.set_xticklabels([a["label"] for a in AREAS], fontsize=10) ax.set_ylabel(ylabel) ax.set_title(metric) ax.legend(fontsize=8, loc="upper left") ax.grid(axis="y", alpha=0.3) fig.suptitle("土砂災害警戒区域 内 vs 外 — 傾斜の量的比較 (H4 検証)", fontsize=12, y=1.02) plt.tight_layout() fig.savefig(ASSETS / "L39_fig6_sediment_compare.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L39_fig6_sediment_compare.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 20. 図 7: small multiples — 階級別空間分布 (3 × 5) # ============================================================================= print("\n[20] 図 7: 階級別 small multiples", flush=True) t1 = time.time() fig, axes = plt.subplots(3, 5, figsize=(18, 10)) for row_i, area in enumerate(AREAS): extent = (area["bounds"].left, area["bounds"].right, area["bounds"].bottom, area["bounds"].top) for col_i, (lab, lo, hi, col, rng) in enumerate(SLOPE_CLASSES): ax = axes[row_i, col_i] m = (area["arr"] > lo) & (area["arr"] <= hi) & area["valid_mask"] n_cells = int(m.sum()) n_valid = int(area["valid_mask"].sum()) pct = n_cells / max(1, n_valid) * 100 show = np.where(m, 1.0, np.nan) cmap = ListedColormap([col]) ax.imshow(show, cmap=cmap, extent=extent, origin="upper", vmin=0, vmax=1, interpolation="nearest") if area.get("admin") is not None: area["admin"].boundary.plot(ax=ax, color="#222", linewidth=0.5, alpha=0.7) if col_i == 0: ax.set_ylabel(f"{area['label']}\n({area['resolution_m']}m)", fontsize=10, fontweight="bold") if row_i == 0: ax.set_title(f"{lab} {rng}", fontsize=10) ax.text(0.05, 0.05, f"{pct:.1f}%\n({n_cells:,} cells)", transform=ax.transAxes, fontsize=9, fontweight="bold", bbox=dict(facecolor="white", alpha=0.85, edgecolor="none", pad=2)) ax.set_xticks([]); ax.set_yticks([]) fig.suptitle("傾斜階級別 空間分布 small multiples — 3 エリア × 5 階級", fontsize=12, y=0.99) plt.tight_layout() fig.savefig(ASSETS / "L39_fig7_small_multiples.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L39_fig7_small_multiples.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 21. 図 8: 中央サンプル (本来解像度の傾斜表示) # ============================================================================= print("\n[21] 図 8: 中央サンプル", flush=True) t1 = time.time() fig, axes = plt.subplots(1, 3, figsize=(15, 5.5)) for ax, a in zip(axes, AREAS): raw = a["raw"] show = np.where(a["raw_valid"], raw, np.nan) im = ax.imshow(show, cmap=slope_cmap, vmin=0, vmax=60, interpolation="nearest") side_m = raw.shape[1] * a["resolution_m"] ax.set_title(f"{a['label']} 中央 1024x1024 ({a['resolution_m']}m/cell, " f"{side_m:.0f}m × {side_m:.0f}m)\n" f"平均 {a['raw_mean']:.1f}° / 30°超 {a['raw_pct_ge_30']:.1f}% / 急峻 {a['raw_pct_ge_45']:.2f}%", fontsize=9.5) ax.set_xlabel("ピクセル X"); ax.set_ylabel("ピクセル Y") cbar = fig.colorbar(im, ax=axes, orientation="horizontal", fraction=0.04, pad=0.12, shrink=0.85) cbar.set_label("傾斜角 (度)") fig.suptitle("中央 1024×1024 セル本来解像度サンプル — 解像度差の体感", fontsize=12, y=1.01) fig.savefig(ASSETS / "L39_fig8_center_samples.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L39_fig8_center_samples.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 22. 図 9: L38 CS立体図 vs 傾斜の散布図 (Pearson 相関) # ============================================================================= print("\n[22] 図 9: L38 vs slope 散布", flush=True) t1 = time.time() fig, axes = plt.subplots(1, 3, figsize=(15, 5)) for ax, a in zip(axes, AREAS): bvm = a.get("both_valid_mask") if bvm is None or bvm.sum() < 100: ax.text(0.5, 0.5, f"{a['label']}\n(L38 TIFF 未取得)", ha="center", va="center", fontsize=12, transform=ax.transAxes) ax.set_xticks([]); ax.set_yticks([]) continue v_slope = a["arr"][bvm].astype(np.float32) v_sat = a["cs_sat"][bvm] n = len(v_slope) sample_n = min(15000, n) idx = np.random.RandomState(42).choice(n, sample_n, replace=False) if n > sample_n else np.arange(n) ax.scatter(v_slope[idx], v_sat[idx], s=2, alpha=0.3, c=a["color"]) # 線形回帰線 if v_slope.std() > 0: z = np.polyfit(v_slope, v_sat, 1) xx = np.linspace(0, min(v_slope.max(), 70), 50) ax.plot(xx, z[0]*xx + z[1], "k-", linewidth=2, alpha=0.7, label=f"slope={z[0]:.4f}*deg + {z[1]:.3f}") r = l38_correlations[a["key"]]["r_sat"] ax.set_xlabel("傾斜角 (度, L39)") ax.set_ylabel("CS立体図 Saturation (L38)") ax.set_title(f"{a['label']}: n={l38_correlations[a['key']]['n']:,} (sample {sample_n:,})\n" f"Pearson r = {r:+.3f}", fontsize=10) ax.set_xlim(0, 70) ax.set_ylim(0, 0.4) ax.legend(fontsize=8, loc="upper left") ax.grid(alpha=0.3) fig.suptitle("L39 傾斜 × L38 CS立体図 Saturation — 同自治体ピクセル相関 (H7 検証)", fontsize=12, y=1.02) plt.tight_layout() fig.savefig(ASSETS / "L39_fig9_l38_correlation.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L39_fig9_l38_correlation.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 23. 図 10: slope 概念図 # ============================================================================= print("\n[23] 図 10: slope 概念図", flush=True) t1 = time.time() fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # 左: DEM → slope 概念 ax = axes[0] ax.set_xlim(0, 12); ax.set_ylim(-0.5, 6); ax.set_axis_off() ax.text(6, 5.7, "傾斜角の計算 — DEM (標高) → 局所勾配 → arctan", ha="center", fontsize=12, fontweight="bold") # 地形断面 xs = np.linspace(0.5, 11.5, 100) ys = 1.5 + 1.4*np.exp(-((xs-3)/0.8)**2) - 1.0*np.exp(-((xs-7)/0.7)**2) + 0.3*np.sin(xs*1.2) ax.fill_between(xs, ys, 0, color="#deb887", alpha=0.45) ax.plot(xs, ys, color="#5a3", linewidth=2.0) # 接線で傾斜角を可視化 for px, label, col in [(2.0, "緩 ~10°", "#a5d6a7"), (4.6, "急 ~40°", "#fb8c00"), (8.5, "中 ~20°", "#ffd54f")]: i = np.argmin(np.abs(xs - px)) if i < 1 or i >= len(xs)-1: continue dy = (ys[i+1] - ys[i-1]) / (xs[i+1] - xs[i-1]) angle_deg = np.degrees(np.arctan(abs(dy))) # 接線セグメント seg_x = np.array([xs[i]-0.6, xs[i]+0.6]) seg_y = ys[i] + dy * (seg_x - xs[i]) ax.plot(seg_x, seg_y, color=col, linewidth=3.5, alpha=0.85) ax.scatter([xs[i]], [ys[i]], s=120, c=col, edgecolor="#222", zorder=5) ax.text(xs[i], ys[i]+0.55, f"{label}\n({angle_deg:.0f}°)", ha="center", fontsize=9, fontweight="bold", color="#222") ax.text(0.5, 0.0, "DEM 標高断面 →", fontsize=10, color="#666") ax.text(11.5, 0.0, "slope = arctan(|∇z|)", fontsize=10, color="#222", ha="right", fontweight="bold") # 右: 5階級凡例 ax = axes[1] ax.set_xlim(0, 10); ax.set_ylim(0, 10); ax.set_axis_off() ax.text(5, 9.5, "傾斜 5 階級 (砂防工学標準)", ha="center", fontsize=12, fontweight="bold") legend_data = [ ("平地", "#e8f5e9", "<5°", "農地・市街地・河川敷"), ("緩傾斜", "#a5d6a7", "5-15°", "段丘・里山・宅地造成可"), ("中傾斜", "#ffd54f", "15-30°", "農地限界・林業境界"), ("急傾斜", "#fb8c00", "30-45°", "急傾斜地崩壊危険箇所基準"), ("急峻", "#b71c1c", ">45°", "崖・崩壊跡・人家不可"), ] y0 = 8.0 for lbl, col, rng, desc in legend_data: ax.add_patch(plt.Rectangle((0.4, y0-0.3), 0.8, 0.6, facecolor=col, edgecolor="#222")) ax.text(1.5, y0, f"{lbl} {rng}", fontsize=11, fontweight="bold", color="#222", va="center") ax.text(4.5, y0, desc, fontsize=9.5, va="center", color="#444") y0 -= 1.3 ax.text(5, 0.5, "急傾斜地崩壊危険箇所の指定基準: 傾斜 30°以上 + 高さ 5m以上 + 人家 5戸以上\n" "= 傾斜 30°が法的なリスク閾値の一つ", ha="center", fontsize=9, color="#666", bbox=dict(facecolor="#fff8c5", edgecolor="#d4a72c", boxstyle="round,pad=0.5")) plt.tight_layout() fig.savefig(ASSETS / "L39_fig10_concept.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L39_fig10_concept.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 24. 図 11: 警戒区域内の傾斜階級構成 stack-bar (内 vs 外) # ============================================================================= print("\n[24] 図 11: 警戒区域内の階級構成", flush=True) t1 = time.time() # 各エリア × 各 zone の階級内訳を計算 zone_class_rows = [] for a in AREAS: for zone_name, m in [("特別警戒", a["red_mask"] & a["valid_mask"]), ("警戒", a["yellow_mask"] & a["valid_mask"]), ("区域外", (~a["yellow_mask"]) & (~a["red_mask"]) & a["valid_mask"])]: if m.sum() < 5: continue v = a["arr"][m] for i, (lab, lo, hi, _, _) in enumerate(SLOPE_CLASSES, start=1): cnt = int(((v > lo) & (v <= hi)).sum()) zone_class_rows.append({ "area": a["label"], "zone": zone_name, "class": lab, "n": cnt, "pct": float(cnt / max(1, len(v)) * 100), }) zone_class_df = pd.DataFrame(zone_class_rows) zone_class_df.to_csv(ASSETS / "L39_zone_class_breakdown.csv", index=False, encoding="utf-8-sig") fig, axes = plt.subplots(1, 3, figsize=(15, 5)) zones_order2 = ["区域外", "警戒", "特別警戒"] class_show = [lab for (lab, _, _, _, _) in SLOPE_CLASSES] class_cols = [c for (_, _, _, c, _) in SLOPE_CLASSES] for ax, a in zip(axes, AREAS): xs = np.arange(len(zones_order2)) bottom = np.zeros(len(zones_order2)) for i, (cn, col) in enumerate(zip(class_show, class_cols)): vals = [] for zn in zones_order2: row = zone_class_df[(zone_class_df["area"] == a["label"]) & (zone_class_df["zone"] == zn) & (zone_class_df["class"] == cn)] vals.append(float(row["pct"].iloc[0]) if len(row) else 0) ax.bar(xs, vals, bottom=bottom, label=cn if ax is axes[0] else None, color=col, edgecolor="#222", linewidth=0.6) for j, v in enumerate(vals): if v >= 5: txtcol = "white" if col in ("#fb8c00","#b71c1c") else "black" ax.text(j, bottom[j] + v/2, f"{v:.0f}", ha="center", va="center", fontsize=8.5, fontweight="bold", color=txtcol) bottom = bottom + np.array(vals) ax.set_xticks(xs) ax.set_xticklabels(zones_order2, fontsize=10) ax.set_ylabel("階級構成比 (%)" if ax is axes[0] else "") ax.set_title(f"{a['label']}", fontsize=11) ax.set_ylim(0, 100) ax.grid(axis="y", alpha=0.3) axes[0].legend(loc="upper left", bbox_to_anchor=(0, -0.12), ncol=5, fontsize=9) fig.suptitle("zone別 傾斜階級構成 — 区域外 → 警戒 → 特別警戒 で急傾斜が支配的になるか?", fontsize=12, y=1.02) plt.tight_layout() fig.savefig(ASSETS / "L39_fig11_zone_class.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L39_fig11_zone_class.png ({time.time()-t1:.1f}s)", flush=True) print(f"\n=== Final figure elapsed: {time.time()-t0:.1f}s ===", flush=True) # ============================================================================= # 25. HTML 生成 # ============================================================================= print("\n[25] HTML 生成", flush=True) t1 = time.time() def dl(name): return f'{name}' sections = [] # === Section 1: 学習目標と問い === sec1 = f"""
カバー宣言: 本記事は DoBoX のシリーズ 「地形傾斜図 1m / 50cm」 2 件 (dataset_id = 1630, 1631) を統合し、広島県 3 自治体エリア (世羅町 1m / 熊野町 1m / 坂町 50cm) の 地形傾斜角分布を 量的に解析する研究記事です。 さらに L10/L11 で扱った土砂災害警戒区域 (yellow / red) Shapefile を空間結合し、 「警戒区域は本当に急傾斜地に指定されているか」を量的に検証します。 最後に L38 (CS立体図) と同自治体ペアでピクセル単位の Pearson 相関を取り、 「曲率×傾斜の RGB 段彩」(L38) と「傾斜単独の度数値」(L39) の物理的整合性を確認します。

本記事の独自性 — 砂防工学標準閾値 × 警戒区域 × L38 整合性

地形傾斜図はオープンデータとして広く公開されていますが、本記事は以下の 3 軸で独自分析を行います:

シリーズ構造判定: (a) 同一手法・解像度違い + (b) 自治体ごとの分割

このシリーズは 1m 版 / 50cm 版 の 2 解像度系統で、各自治体ごとに GeoTIFF (Float32 度数) を切出して公開する構造。 1m 版は 22 自治体、50cm 版は 12 自治体に対応。

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

用語定義
傾斜角 (slope angle)DEM (標高ラスタ) の各ピクセルでの局所勾配の大きさを arctan で度数 (0-90) に変換した値。式: slope = arctan( sqrt((dz/dx)² + (dz/dy)²) )。 水平面 = 0°、垂直崖 = 90°。本記事の TIFF はこの計算を 1m もしくは 50cm メッシュで実行した結果。
Horn 法傾斜計算の標準アルゴリズム (Horn 1981)。3×3 セル窓で 重み付き中央差分を取り、x 方向 dz/dx と y 方向 dz/dy を求める方法。GIS の標準実装。 本 TIFF の生成元は明示されていないが、結果値域 (0-90 度) と分布形状から Horn 法または 類似の 3×3 窓法と推定。
傾斜 5 階級 (本記事独自定義, 砂防工学標準) 平地 (<5°), 緩傾斜 (5-15°), 中傾斜 (15-30°), 急傾斜 (30-45°), 急峻 (>45°) の 5 段階。30° は 急傾斜地崩壊危険箇所の法的指定基準 (傾斜 30° 以上 + 高さ 5m 以上 + 人家 5戸以上)、 45° は経験的な人家立地の事実上の上限値で、両者を合わせて 5 階級とした。
急傾斜 (steep slope)本記事では傾斜 30° 以上のセルを指す。 急傾斜地崩壊危険箇所の法的指定基準と一致。30° = tan(30°) ≈ 0.577 = 高さ 1m / 水平 1.73m の傾き。 人が歩行で登るのが困難になる目安でもある。
急峻 (very steep)傾斜 45° 超のセル。 45° = tan(45°) = 1 = 「高さ : 水平 = 1 : 1」 の傾きで、土砂安息角を大きく超える。 通常は崖・崩壊跡・岩盤露出地が該当し、人家不可。
平地 (flat)傾斜 5° 未満のセル。農地・市街地・河川敷が主体。 5° = tan(5°) ≈ 0.087 = 高さ 1m / 水平 11m の緩い傾き。
8 近傍同種率 (neighbor same rate)本記事独自の集塊性指標。 ある属性 (例: 急傾斜) を持つセルについて、8 近傍にも同じ属性が存在する確率。 Moran's I の近似的代替として使用 (計算量が桁違いに軽い)。 高い (≥ 0.85) = 強集塊, 低い (≤ 0.5) = 散在。
土砂災害警戒区域 (yellow zone)土砂災害防止法に基づき指定。 本記事では DoBoX dsid 48「土砂災害警戒区域・特別警戒区域情報」の 3 種別 (土石流・急傾斜地・地すべり) を統合して使用。
特別警戒区域 (red zone)同法のレッドゾーン。建物倒壊リスクあり。 土石流・急傾斜地のみで指定 (地すべりはレッドなし)。
L38 CS立体図 SaturationL38 で扱う RGB 段彩図の HSV 変換後の彩度成分。 S = (max − min) / max。「微地形コントラストの強さ」の代理指標。 本記事では L38 と同じ overview shape で読み込み、傾斜との Pearson 相関を取る。

研究の問い (RQ)

広島県の地形傾斜図 2 解像度 (1m / 50cm) は、それぞれ 世羅町 (中山間/林業地), 熊野町 (都市近郊/山岳混在), 坂町 (平成30豪雨被災地) という異なる地形リスク特性のエリアを公開している。 3 エリアの傾斜分布を 5 階級で量化し、地形リスク構造を比較する。 さらに土砂災害警戒区域との空間結合で「警戒区域 = 急傾斜地」という法的論理が データに刻まれているかを量的に検証する。最後にL38 CS立体図との物理的整合性を ピクセル相関で確認する。

  1. 2 dataset の解像度・面積・有効データ率・ピクセル数の構造比較は?
  2. 3 エリアの傾斜角ヒストグラム (2.5° bin) はどう違うか?
  3. 砂防工学標準 5 階級の構成比は?
  4. 30°超急傾斜セルは空間集塊するか? 8 近傍同種率は?
  5. 土砂災害警戒区域の内 vs 外で平均傾斜・階級構成はどう違うか?
  6. 解像度差 (1m vs 50cm) で急傾斜検出率はどう変わるか?
  7. L38 CS立体図 Saturation と本記事の傾斜は正相関するか?

仮説 H1〜H7

到達点

3 エリアの地形傾斜図を「3 自治体の代表エリア」として量的に解析し、 合計 {int(metas[0]['n_cells']/1e6 + metas[1]['n_cells']/1e6 + metas[2]['n_cells']/1e6):,} M セルの Float32 傾斜データを overview ピラミッドで効率的に概観する。 さらに土砂災害警戒区域 Shapefile (~ 2,000 ポリゴン) との空間結合で 「法的指定論理の量的検証」を行い、L38 CS立体図とのピクセル相関で 「LiDAR ファミリ内の物理的整合性」を実証する。 これは傾斜単独・警戒区域単独・CS立体図単独では不可能な、3 軸統合分析である。

本記事のスコープ外: 本記事は地形傾斜図の2 dataset (1m + 50cm) のみを統合し、 土砂警戒区域 + L38 CS立体図との空間結合に特化する。
""" sections.append(("学習目標と問い", sec1)) # === Section 2: 使用データ === sec2 = f"""

本記事が使用する 2 dataset の一覧。両者は同一手法 (DEM からの傾斜計算) の解像度違いで、 各自治体ごとに GeoTIFF を切り出して公開する構造。本記事は 2 dataset から合計 3 自治体を 代表として用います。

dsid名称解像度公開単位本記事の代表TIFF サイズセル数DoBoX
1630地形傾斜図（1m）1.0 m/cell 22 自治体別 GeoTIFF世羅町 (34462) + 熊野町 (34307) {int(AREAS[0]['tif'].stat().st_size/(1024*1024)):,} + {int(AREAS[1]['tif'].stat().st_size/(1024*1024)):,} MB {(AREAS[0]['meta']['n_cells']+AREAS[1]['meta']['n_cells'])/1e6:.1f} M #1630
1631地形傾斜図（50cm）0.5 m/cell 12 自治体別 GeoTIFF坂町 (34384) {int(AREAS[2]['tif'].stat().st_size/(1024*1024)):,} MB {AREAS[2]['meta']['n_cells']/1e6:.1f} M #1631

データ仕様 (両 dataset 共通)

L36/L37/L38 との関係 (本記事は単独, 参考併置のみ)

軸L36 樹高L37 点群密度L38 CS立体図L39 傾斜 (本記事)
計測対象樹冠の高さ (m)地面到達パルス数地形微細構造 (RGB)傾斜角 (度)
派生過程DSM − DTMグラウンド分類点数DTM → 曲率 + 傾斜 → RGBDTM → 勾配 → arctan
同じ点群?YES — 同一の航空レーザ計測点群から 派生した 4 つの集計
表現連続値 float32整数 uint8/int16RGB 3 バンド uint8連続値 float32
L39 との関係独立 (植生)独立 (計測精度) 同DEM由来 — 傾斜成分を含む (本記事 H7 で相関検証)本記事
本記事のスコープ L39 単独で完結 (L36/L37/L38 とは合体しない, ただし L38 との Pearson 相関は本記事の H7 で検証)

3 自治体の特性比較

軸世羅町 (1m)熊野町 (1m)坂町 (50cm)
立地備後内陸 (中山間)広島市東 (都市近郊)呉市西 (沿岸山間)
面積 (行政界) {round(AREAS[0]['muni_area_km2'], 1) if AREAS[0]['muni_area_km2'] else '?'} km² {round(AREAS[1]['muni_area_km2'], 1) if AREAS[1]['muni_area_km2'] else '?'} km² {round(AREAS[2]['muni_area_km2'], 1) if AREAS[2]['muni_area_km2'] else '?'} km²
主産業農林業 (世羅梨・大根)筆 (熊野筆) 製造住宅地 + 港湾
地形特性緩傾斜丘陵花崗岩マサ土山地急傾斜面 + 崩壊跡 (H30豪雨被災地)
有効セル比率 {basic_stats['sera_1m']['valid_ratio']*100:.1f}% {basic_stats['kumano_1m']['valid_ratio']*100:.1f}% {basic_stats['saka_50cm']['valid_ratio']*100:.1f}%
平均傾斜 {basic_stats['sera_1m']['mean']:.2f}° {basic_stats['kumano_1m']['mean']:.2f}° {basic_stats['saka_50cm']['mean']:.2f}°
30° 超セル比率 {basic_stats['sera_1m']['pct_ge_30']:.1f}% {basic_stats['kumano_1m']['pct_ge_30']:.1f}% {basic_stats['saka_50cm']['pct_ge_30']:.1f}%
45° 超 (急峻) 比率 {basic_stats['sera_1m']['pct_ge_45']:.2f}% {basic_stats['kumano_1m']['pct_ge_45']:.2f}% {basic_stats['saka_50cm']['pct_ge_45']:.2f}%
土砂警戒 (yellow) {AREAS[0]['sed_yellow_n']} polys / {AREAS[0]['sed_yellow_km2']:.1f} km² {AREAS[1]['sed_yellow_n']} polys / {AREAS[1]['sed_yellow_km2']:.1f} km² {AREAS[2]['sed_yellow_n']} polys / {AREAS[2]['sed_yellow_km2']:.1f} km²
特別警戒 (red) {AREAS[0]['sed_red_n']} polys / {AREAS[0]['sed_red_km2']:.2f} km² {AREAS[1]['sed_red_n']} polys / {AREAS[1]['sed_red_km2']:.2f} km² {AREAS[2]['sed_red_n']} polys / {AREAS[2]['sed_red_km2']:.2f} km²

データ品質メモ

""" sections.append(("使用データ", sec2)) # === Section 3: ダウンロード === sec3 = f"""

生データ (DoBoX 直リンク)

容量警告: 3 ZIP の合計は約 {int((AREAS[0]['zip_local'].stat().st_size + AREAS[1]['zip_local'].stat().st_size + AREAS[2]['zip_local'].stat().st_size)/(1024*1024))} MB、 展開後の TIFF は約 {int((AREAS[0]['tif'].stat().st_size + AREAS[1]['tif'].stat().st_size + AREAS[2]['tif'].stat().st_size)/(1024*1024))} MB。

土砂災害警戒区域データ (本記事の空間結合に必要)

本記事はL10 土砂災害警戒区域 × 用途地域と 同じ Shapefile (DoBoX dsid 48) を使用します。 data/extras/sediment_shp/ に展開済みであれば自動利用。

L38 CS立体図 TIFF (H7 の相関検証に必要, optional)

本記事は L38 と同じ自治体ペアを採用しています。 data/extras/L38_cs_terrain/ 配下の TIFF (CS立体図 RGB) があると自動的に ピクセル相関を計算します。L38 を先に実行することで全機能が有効になります。 ただし L38 TIFF が無くても L39 単独で全分析が走ります (H7 の判定のみ NA になる)。

中間データ (本記事生成 CSV)

図 (本記事生成 PNG)

再現用 Python スクリプト

L39_terrain_slope.py を取得して プロジェクトルートで py -X utf8 lessons/L39_terrain_slope.py を実行。 TIFF が無ければ ZIP 自動 DL → 自動展開 → overview 読込で分析が走ります。

""" sections.append(("ダウンロード", sec3)) # === Section 4: 分析 1 — dataset の構造を見る === sec4 = f"""

狙い

2 dataset (1m / 50cm) と 3 自治体 (世羅 / 熊野 / 坂) の規模・形状・傾斜分布の概形を 1 枚の絵で示す。傾斜ヒストグラムを重ねて地形リスク構造を最初に視覚化する。

手法 (リテラシレベル解説)

L36-L38 と全く同じ rasterio スタイル。傾斜 TIFF は1 バンド Float32 で、各セルが 傾斜角 (度) を持つ。読込手順:

実装

{code(''' import rasterio from rasterio.enums import Resampling import numpy as np TARGET_CELL_M = 32 # overview 目標セルサイズ for area in AREAS: with rasterio.open(area["tif"]) as ds: ovr_choices = ds.overviews(1) target_factor = TARGET_CELL_M / area["resolution_m"] best = next((f for f in ovr_choices if f >= target_factor), ovr_choices[-1] if ovr_choices else int(target_factor)) ov_w = ds.width // best ov_h = ds.height // best # 1 バンド読込 (傾斜は単一バンド) → shape (H, W) arr = ds.read(1, out_shape=(ov_h, ov_w), resampling=Resampling.average).astype(np.float32) nodata = ds.nodata if ds.nodata is not None else -99999.0 valid_mask = (~np.isnan(arr)) & (arr != nodata) & (arr >= 0) & (arr <= 90.001) area["arr"] = arr area["valid_mask"] = valid_mask ''')}

図と読み取り (図 1: dataset カード + 傾斜ヒスト)

なぜこの図か: 3 dataset の規模と特性をカードで言語的に、 かつ 2.5° bin のヒストグラムで「傾斜分布の形」として量的に、同時に 1 枚で伝えるため。 階級境界 (5°/15°/30°/45°) を縦線で示し、5 階級分類の文脈を即座に与える。

{figure("assets/L39_fig1_dataset_overview.png", "図 1: 3 dataset カード (左) + 傾斜分布 (右)")}

読み取り:

表 (基本統計)

{stats_df.round(2).to_html(classes='', border=0)}

読み取り: 平均・中央値・p90 が 3 エリアで明確に階層化される。 30° 超セル比率は世羅 {basic_stats['sera_1m']['pct_ge_30']:.1f}% / 熊野 {basic_stats['kumano_1m']['pct_ge_30']:.1f}% / 坂 {basic_stats['saka_50cm']['pct_ge_30']:.1f}%で、 急傾斜地崩壊危険箇所の法的指定基準 (30°) を超える面積比が地域差を端的に示す。 急峻 (>45°) 比率は世羅 {basic_stats['sera_1m']['pct_ge_45']:.2f}% / 熊野 {basic_stats['kumano_1m']['pct_ge_45']:.2f}% / 坂 {basic_stats['saka_50cm']['pct_ge_45']:.2f}%で、 極端な崖地形の存在を量化する。

""" sections.append(("分析 1: 3 dataset とエリアの構造を可視化", sec4)) # === Section 5: 分析 2 — 傾斜ラスタ主題図 === sec5 = f"""

狙い

傾斜の絶対値分布を 3 エリア並列で表示し、地形リスクの空間パターンを視覚的に理解する。 緑 → 黄 → 赤 のグラディエントで「平地は緑、急峻は赤」と直感的に読める。

手法 (リテラシレベル解説)

連続値ラスタの主題図化は matplotlib.pyplot.imshow に vmin=0, vmax=60, cmap=... を渡すだけで実現できる。 本記事では5 階級境界 (5/15/30/45°) に対応する 6 色のグラディエントを LinearSegmentedColormap で構築し、緑 (平地) → 黄 (中傾斜) → 赤 (急峻) の 段彩を作る。これにより「色が傾斜カテゴリと対応」する視認性が得られる。

実装

{code(''' from matplotlib.colors import LinearSegmentedColormap slope_cmap = LinearSegmentedColormap.from_list( "slope", ["#1a9850", "#a6d96a", "#fee08b", "#fdae61", "#d73027", "#7f0000"], N=256 ) # nodata は NaN 化して透明表示 arr_show = np.where(valid_mask, arr, np.nan) im = ax.imshow(arr_show, cmap=slope_cmap, vmin=0, vmax=60, extent=area_bounds, origin="upper", interpolation="nearest") ''')}

図と読み取り (図 2: 傾斜主題図)

なぜこの図か: 5 階級分類 (図 3) の前に、まず連続値の生分布を見て 「どこが急斜面か」を直感的に把握する。グラディエント表示はカテゴリ化前の本来の信号。

{figure("assets/L39_fig2_slope_maps.png", "図 2: 3 エリア 傾斜ラスタ主題図 (緑=平地 → 暗赤=急峻)")}

読み取り:

表 (TIFF メタ + 行政界面積)

項目世羅町熊野町坂町
解像度 (m/cell){AREAS[0]['resolution_m']}{AREAS[1]['resolution_m']}{AREAS[2]['resolution_m']}
セル数 (M){AREAS[0]['meta']['n_cells']/1e6:.1f}{AREAS[1]['meta']['n_cells']/1e6:.1f}{AREAS[2]['meta']['n_cells']/1e6:.1f}
TIFF bbox 面積 (km²){AREAS[0]['meta']['area_km2']:.1f}{AREAS[1]['meta']['area_km2']:.1f}{AREAS[2]['meta']['area_km2']:.1f}
行政界面積 (km²){AREAS[0]['muni_area_km2']:.1f}{AREAS[1]['muni_area_km2']:.1f}{AREAS[2]['muni_area_km2']:.1f}
有効セル比率{basic_stats['sera_1m']['valid_ratio']*100:.1f}%{basic_stats['kumano_1m']['valid_ratio']*100:.1f}%{basic_stats['saka_50cm']['valid_ratio']*100:.1f}%
nodata 値{AREAS[0]['meta']['nodata']}{AREAS[1]['meta']['nodata']}{AREAS[2]['meta']['nodata']}
overview 倍率/1/{AREAS[0]['arr_factor']}/1/{AREAS[1]['arr_factor']}/1/{AREAS[2]['arr_factor']}
実効セル (m/cell){AREAS[0]['arr_cell_m']:.0f}{AREAS[1]['arr_cell_m']:.0f}{AREAS[2]['arr_cell_m']:.0f}
最大傾斜 (overview, °){basic_stats['sera_1m']['max']:.1f}{basic_stats['kumano_1m']['max']:.1f}{basic_stats['saka_50cm']['max']:.1f}

読み取り: overview /1/32 で実効セルサイズ {AREAS[0]['arr_cell_m']:.0f}-{AREAS[2]['arr_cell_m']:.0f}m に揃え、 3 エリア間で比較可能な解像度に標準化している。坂町は 50cm × overview /{AREAS[2]['arr_factor']} で実効 {AREAS[2]['arr_cell_m']:.0f}m となるため、解像度差は overview レベルでは消える (= 図 8 の生解像度サンプルでのみ見える)。

""" sections.append(("分析 2: 傾斜ラスタ主題図化", sec5)) # === Section 6: 分析 3 — 5 階級分類 === sec6 = f"""

狙い

連続値の傾斜ラスタを砂防工学標準の 5 階級にカテゴリ化し、 3 エリアの「地形リスク構造」を比較可能な形にする。 階級境界は法的・経験的に意味のある値 (5° 平地基準, 15° 緩傾斜限界, 30° 急傾斜地法定基準, 45° 崖基準)。

手法 (リテラシレベル解説)

連続値の傾斜を 5 種のラベル (1-5) に変換する。閾値は絶対値で固定するため、 3 エリア間で横断比較が可能になる (L38 のパーセンタイル適応閾値とは設計思想が異なる)。

クラス条件意味 / 法的基準
1. 平地< 5°農地・市街地・河川敷。ほぼ水平。
2. 緩傾斜5-15°段丘・里山・宅地造成可能範囲。
3. 中傾斜15-30°農地限界 (棚田で 20° 程度)、林業境界。
4. 急傾斜30-45°急傾斜地崩壊危険箇所の法的指定基準 (傾斜 30° + 高さ 5m + 人家 5戸)
5. 急峻> 45°崖・崩壊跡・岩盤露出。人家立地不可。

限界と代替:

実装

{code(''' SLOPE_CLASSES = [ ("平地", -0.001, 5.0), ("緩傾斜", 5.0, 15.0), ("中傾斜", 15.0, 30.0), ("急傾斜", 30.0, 45.0), ("急峻", 45.0, 90.001), ] def classify_slope(arr, valid_mask): label = np.zeros(arr.shape, dtype=np.int8) for i, (_, lo, hi) in enumerate(SLOPE_CLASSES, start=1): m = (arr > lo) & (arr <= hi) & valid_mask label[m] = i return label ''')}

図と読み取り (図 3: 5 階級分類マップ)

なぜこの図か: 連続グラディエント (図 2) では「リスクのカテゴリ」が読みにくい。 階級分けすることで、「30°以上の急傾斜面はどこか」が即座に分かる。

{figure("assets/L39_fig3_class_map.png", "図 3: 5 階級分類マップ")}

読み取り:

図と読み取り (図 4: 階級構成比 stack-bar)

{figure("assets/L39_fig4_class_composition.png", "図 4: 5 階級構成比の積み上げ棒グラフ")}

読み取り: 30° 超 (急 + 急峻) の比率は世羅 {basic_stats['sera_1m']['pct_ge_30']:.1f}% / 熊野 {basic_stats['kumano_1m']['pct_ge_30']:.1f}% / 坂 {basic_stats['saka_50cm']['pct_ge_30']:.1f}%と 2-3 倍の階層差がある。逆に平地 (<5°)は世羅 {basic_stats['sera_1m']['pct_lt_5']:.1f}% / 熊野 {basic_stats['kumano_1m']['pct_lt_5']:.1f}% / 坂 {basic_stats['saka_50cm']['pct_lt_5']:.1f}% で、 坂町が最大 (= 海岸に近い扇状地・住宅地の平坦地が広い)。 これは「坂町=急峻」のステレオタイプ的予想とは逆方向で、平地と急傾斜が共存する地形構造が読み取れる。

表 (傾斜階級構成比 wide)

{class_df.round(2).to_html(classes='', border=0, index=False)}

図と読み取り (図 7: 階級別 small multiples)

{figure("assets/L39_fig7_small_multiples.png", "図 7: 階級別の空間分布 small multiples (3 エリア × 5 階級)")}

読み取り:

""" sections.append(("分析 3: 5 階級分類 (砂防工学標準閾値)", sec6)) # === Section 7: 分析 4 — 警戒区域 × 傾斜 (重要セクション) === sec7 = f"""

狙い (本記事の核心の 1 つ)

土砂災害警戒区域 (yellow) と特別警戒区域 (red) のポリゴンを傾斜ラスタの座標系に rasterize して、「警戒区域内 vs 区域外で平均傾斜・階級構成はどう違うか」を量的に比較する。 法的指定論理 (急傾斜地崩壊危険箇所基準) がデータに刻まれているかを検証する。

手法 (リテラシレベル解説)

ポリゴン → ラスタ化 (rasterize): 不規則な形のポリゴン (警戒区域)を 「画素ごとに 0/1 のマスク」に変換する操作。L38 と同じ rasterio.features.rasterize を使用。

  1. 警戒区域ポリゴンを rasterio.features.rasterize に渡し、傾斜ラスタと同じ shape・transform で描画。
  2. 得られたバイナリマスク (1=区域内, 0=外) と傾斜ラスタをセル単位でペア化。
  3. 各 zone (特別警戒 / 警戒 / 区域外) に分けて 平均傾斜・p90・30°超比率 を計算。
  4. 3 ゾーン × 3 エリア × 4 指標のクロス表を作る。

実装

{code(''' from rasterio.features import rasterize from rasterio.transform import from_bounds def rasterize_polys(polys_gdf, area): """警戒区域 polygons を傾斜ラスタの (H, W) shape にラスタ化""" bb = area["bounds"] transform = from_bounds(bb.left, bb.bottom, bb.right, bb.top, area["arr_w"], area["arr_h"]) shapes = ((geom, 1) for geom in polys_gdf.geometry) mask = rasterize(shapes, out_shape=(area["arr_h"], area["arr_w"]), transform=transform, fill=0, dtype="uint8") return mask.astype(bool) # 警戒区域 / 特別警戒 / 区域外 で各指標を計算 yellow_mask = rasterize_polys(area["sed_yellow"], area) red_mask = rasterize_polys(area["sed_red"], area) in_yellow = yellow_mask & area["valid_mask"] in_red = red_mask & area["valid_mask"] out_warn = ~yellow_mask & ~red_mask & area["valid_mask"] for zone, m in [("特別警戒", in_red), ("警戒", in_yellow), ("区域外", out_warn)]: v = arr[m] mean_slope = v.mean() p90 = np.percentile(v, 90) pct_ge_30 = (v >= 30).sum() / len(v) * 100 ''')}

図と読み取り (図 5: 警戒区域 × 傾斜 重ね合わせ)

なぜこの図か: 警戒区域がどの傾斜域に重なるかを視覚的に確認。 傾斜段彩を下地、警戒区域を境界線で重ねる。

{figure("assets/L39_fig5_sediment_overlay.png", "図 5: 傾斜 + 土砂災害警戒区域の重ね合わせ")}

読み取り:

図と読み取り (図 6: 警戒区域 内/外 量的比較 grouped bar)

なぜこの図か: 視覚的な印象 (図 5) を量的に裏付ける。3 エリア × 3 ゾーン × 3 指標のクロス比較。

{figure("assets/L39_fig6_sediment_compare.png", "図 6: 警戒区域 内/外 の平均傾斜・p90・30°超比率")}

読み取り (重要発見):

表 (警戒区域 内/外 の傾斜統計)

{sed_compare_df.round(2).to_html(classes='', border=0, index=False)}

図と読み取り (図 11: zone別 5 階級構成 stack-bar)

なぜこの図か: 平均値だけでは分布の形が見えない。階級構成比を見れば 「警戒区域内は中傾斜が多いか平地が多いか」が分かる。

{figure("assets/L39_fig11_zone_class.png", "図 11: zone別 5 階級構成 stack-bar (3 エリア)")}

読み取り:

""" sections.append(("分析 4: 警戒区域 × 傾斜 — ポリゴン×ラスタの空間結合 (本記事の核心)", sec7)) # === Section 8: 分析 5 — 急傾斜地の集塊性 === sec8 = f"""

狙い

30° 超急傾斜セルが空間的に集塊するか散在するかを量化する。 集塊性は地形リスク評価の重要指標で、集塊する=地続きの急斜面塊で 土砂災害リスクが連鎖的に高い、散在する=点状リスクで影響範囲が限定的、と読める。

手法 (リテラシレベル解説)

本記事独自の「8 近傍同種率」を使う。これは Moran's I の計算量軽減版で、 全ペア相関を取らずに「8 近傍にも同種がいる確率」だけを集計する。

段階処理
1急傾斜マスク steep_mask = (arr >= 30) & valid_mask を作る
28 近傍 (上下左右と斜め) を順に走査し、ある (y,x) と隣 (y+dy, x+dx) が両方有効なペアをカウント
3ペアのうち、両方とも急傾斜であるペアの数をカウント
4同種ペア / 全ペア (どちらかが急傾斜のペア) の比率を計算

解釈:

実装

{code(''' def neighbor_same_rate(mask, valid): H, W = mask.shape pairs_total, pairs_same = 0, 0 for dy in (-1, 0, 1): for dx in (-1, 0, 1): if dy == 0 and dx == 0: continue ma = mask[max(0,dy):H+min(0,dy), max(0,dx):W+min(0,dx)] mb = mask[max(0,-dy):H+min(0,-dy), max(0,-dx):W+min(0,-dx)] va = valid[max(0,dy):H+min(0,dy), max(0,dx):W+min(0,dx)] vb = valid[max(0,-dy):H+min(0,-dy), max(0,-dx):W+min(0,-dx)] both_valid = va & vb pairs_total += int(((ma | mb) & both_valid).sum()) pairs_same += int((ma & mb & both_valid).sum()) return pairs_same / max(1, pairs_total) ''')}

図と読み取り

図 7 (small multiples の急傾斜列) を再度参照。各エリアの急傾斜塊の連続性を量化した結果が下の表。

表 (急傾斜地集塊性)

{cluster_df.round(4).to_html(classes='', border=0, index=False)}

読み取り (重要発見):

""" sections.append(("分析 5: 急傾斜地の空間集塊性 (8 近傍同種率)", sec8)) # === Section 9: 分析 6 — 中央サンプルと L38 相関 === sec9 = f"""

狙い

2 つの分析を並列で扱う:

  1. 中央 1024×1024 セルを本来解像度で読込み、解像度差 (1m vs 50cm) を体感する
  2. L38 CS立体図 Saturation との同自治体ピクセル相関を取り、両者の物理的整合性を検証 (H7)

手法 (リテラシレベル解説)

本来解像度サンプル: rasterio.windows.Window で各 TIFF の中央 1024×1024 セルだけを 読込み、本来解像度で表示。1m なら 1024m × 1024m、50cm なら 512m × 512m の領域。

L38 との相関: 同自治体の L38 TIFF (CS立体図 RGB) を本記事の overview と同 shapeで 読み込み (rasterio の out_shape 機能で空間分解能を揃える)、 ピクセル単位で傾斜 (L39) vs CS立体図 Saturation (L38) の Pearson 相関を計算。 両者は同じ DEM から派生した別表現なので、強い正相関を期待 (H7)。

実装

{code(''' import rasterio from rasterio.windows import Window from rasterio.enums import Resampling # (1) 中央サンプル RAW_SAMPLE = 1024 for area in AREAS: with rasterio.open(area["tif"]) as ds: cx, cy = ds.width//2, ds.height//2 win = Window(cx-RAW_SAMPLE//2, cy-RAW_SAMPLE//2, RAW_SAMPLE, RAW_SAMPLE) raw = ds.read(1, window=win).astype(np.float32) # (2) L38 との Pearson 相関 for area in AREAS: with rasterio.open(area["l38_tif"]) as ds_cs: rgb = ds_cs.read(out_shape=(3, area["arr_h"], area["arr_w"]), resampling=Resampling.average).astype(np.float32) # CS Saturation 計算 maxv = np.maximum.reduce(rgb) minv = np.minimum.reduce(rgb) sat = np.where(maxv > 0, (maxv - minv) / np.maximum(maxv, 1), 0) # 共通有効マスクで Pearson 相関 both = area["valid_mask"] & (rgb.sum(axis=0) > 0) r = np.corrcoef(area["arr"][both], sat[both])[0, 1] ''')}

図と読み取り (図 8: 中央サンプル)

なぜこの図か: overview /1/{AREAS[0]['arr_factor']}-/1/{AREAS[2]['arr_factor']} で粗化された 分析と異なり、本来解像度で見ると急斜面の細部構造が分かる。 1m と 50cm の差はここで初めて視認可能。

{figure("assets/L39_fig8_center_samples.png", "図 8: 中央 1024×1024 セル本来解像度サンプル")}

読み取り:

図と読み取り (図 9: L38 CS立体図 Saturation との Pearson 相関)

なぜこの図か: H7 (傾斜 vs CS Sat の正相関) を視覚的・量的に検証する。 散布図に線形回帰線を重ね、Pearson r を表示。

{figure("assets/L39_fig9_l38_correlation.png", "図 9: L39 傾斜 × L38 CS Saturation 散布図")}

読み取り (重要発見):

表 (L38 相関)

{corr_df.to_html(classes='', border=0, index=False)}

図と読み取り (図 10: 傾斜計算と 5 階級の概念図)

{figure("assets/L39_fig10_concept.png", "図 10: 傾斜計算 (DEM → 勾配 → arctan) と 5 階級の概念図")}

読み取り:

""" sections.append(("分析 6: 中央サンプルと L38 CS立体図との Pearson 相関 (H7)", sec9)) # === Section 10: 仮説検証と考察 === def fmt_verdict(v): if v == "支持": return f"{v}" elif v == "部分支持": return f"{v}" elif v == "反証": return f"{v}" return v hypos_table = "" for h in hypos: hypos_table += f"" hypos_table += "
H主張結果判定
{h['H']}{h['claim']}{h['result']}{fmt_verdict(h['verdict'])}
" n_support = sum(1 for h in hypos if h['verdict'] == '支持') n_partial = sum(1 for h in hypos if h['verdict'] == '部分支持') n_refute = sum(1 for h in hypos if h['verdict'] == '反証') sec_h = f"""

本記事で立てた 7 つの仮説 (H1〜H7) と、3 エリアの傾斜ラスタ解析 + 警戒区域空間結合 + L38 ピクセル相関の照合結果:

{hypos_table}

判定サマリ: 支持 {n_support} / 部分支持 {n_partial} / 反証 {n_refute} (全 7 仮説中)

主な発見 (5 点)

  1. H4 (警戒区域は急傾斜) は反証されたが、反対方向の関係が 3 エリア共通で量的に確認された: 警戒区域内の平均傾斜は区域外より低い (世羅 ×0.76, 熊野 ×0.49, 坂 ×0.72)。 これは警戒区域が「急傾斜面そのもの」ではなく「急斜面から土砂が到達する範囲 (谷出口・扇状地)を 含むように指定されているためで、土砂災害警戒区域の物理的指定論理を量的に明らかにした。 特別警戒 (red) は警戒 (yellow) より平均傾斜が高く、3 ゾーン間に明確な階層。 これは本記事の最も重要な発見。
  2. H7 (L38 CS Sat との正相関) が強支持された: Pearson r = +{l38_correlations['sera_1m']['r_sat']:.3f} / +{l38_correlations['kumano_1m']['r_sat']:.3f} / +{l38_correlations['saka_50cm']['r_sat']:.3f} (世羅/熊野/坂)。 3 エリアすべてで強相関 (r ≥ 0.6)。さらに slope vs CS darkness では r >= 0.9 (世羅・熊野) と 極めて強い相関が出た ⇒ L38 CS立体図の暗さ (G低) は傾斜の最良の代理指標。 これは記事間整合性 (L38 と L39 が同自治体・同 DEM 由来であること) の量的裏付け。
  3. H3 (坂町は最急峻) が支持された: 平均傾斜 23.78° (3エリア最大), 30°超セル比率 45.8% (3 エリア最大)。50cm 解像度の効果も含めて、 平成30豪雨被災地の急斜面構造が量的に確認された。 ただし急峻 (>45°) 比率は 0.2% で熊野と同等 = 極端な崖は 3 エリアで稀。
  4. H6 (50cm 版は急傾斜検出率高) が支持された: 30°超比率は世羅 14.6% / 熊野 32.2% / 坂 45.8% で 坂町 (50cm) が圧倒。ただし「解像度差」と「地形類型差」が交絡しており、 50cm の効果単独を分離するには同自治体 1m vs 50cmの比較が必要 (発展課題)。 集塊性 (8 近傍同種率) も坂町 0.79 で最大 = 高解像度ほど急斜面塊が連続的に解像される。
  5. H5 (急傾斜の集塊性 ≥ 0.85) は反証: 実測 0.60-0.79。 ただし解像度に応じて単調増加する傾向は確認 ⇒ 仮説の閾値設定が高すぎた。 急傾斜は線形構造 (尾根線・崖線)として分布するため、面的指標 (8 近傍) では 過小評価される。これは「線状地形の集塊性は方向性を加味すべき」という方法論的発見。

3 エリアの地形リスク特性比較 (本記事の量的サマリ)

軸世羅町 (1m)熊野町 (1m)坂町 (50cm)
地形類型緩傾斜丘陵 (中山間)急峻山地 + 都市平地急斜面 + 崩壊跡
平均傾斜{basic_stats['sera_1m']['mean']:.2f}°{basic_stats['kumano_1m']['mean']:.2f}°{basic_stats['saka_50cm']['mean']:.2f}°
30°超 (急傾斜){basic_stats['sera_1m']['pct_ge_30']:.1f}%{basic_stats['kumano_1m']['pct_ge_30']:.1f}%{basic_stats['saka_50cm']['pct_ge_30']:.1f}%
45°超 (急峻){basic_stats['sera_1m']['pct_ge_45']:.2f}%{basic_stats['kumano_1m']['pct_ge_45']:.2f}%{basic_stats['saka_50cm']['pct_ge_45']:.2f}%
5°未満 (平地){basic_stats['sera_1m']['pct_lt_5']:.1f}%{basic_stats['kumano_1m']['pct_lt_5']:.1f}%{basic_stats['saka_50cm']['pct_lt_5']:.1f}%
急傾斜集塊度{cluster_df.iloc[0]['neighbor_same_rate']:.3f}{cluster_df.iloc[1]['neighbor_same_rate']:.3f}{cluster_df.iloc[2]['neighbor_same_rate']:.3f}
警戒内/外 比 (傾斜) {(get_zone(sed_compare_rows, '世羅町', '警戒')['mean'] or 0)/max(0.001,(get_zone(sed_compare_rows, '世羅町', '区域外')['mean'] or 1)):.2f} {(get_zone(sed_compare_rows, '熊野町', '警戒')['mean'] or 0)/max(0.001,(get_zone(sed_compare_rows, '熊野町', '区域外')['mean'] or 1)):.2f} {(get_zone(sed_compare_rows, '坂町', '警戒')['mean'] or 0)/max(0.001,(get_zone(sed_compare_rows, '坂町', '区域外')['mean'] or 1)):.2f}
L38 CS Sat 相関 {l38_correlations['sera_1m']['r_sat']:+.3f} {l38_correlations['kumano_1m']['r_sat']:+.3f} {l38_correlations['saka_50cm']['r_sat']:+.3f}
警戒区域率 {AREAS[0]['sed_yellow_km2']/AREAS[0]['muni_area_km2']*100:.1f}% {AREAS[1]['sed_yellow_km2']/AREAS[1]['muni_area_km2']*100:.1f}% {AREAS[2]['sed_yellow_km2']/AREAS[2]['muni_area_km2']*100:.1f}%

考察

本記事の主たる発見は、「土砂災害警戒区域は急傾斜面そのものではなく、その影響範囲を含む」 ことが 3 エリア共通で量的に確認されたことである。仮説 H4 (警戒区域内の方が急傾斜) は反証され、 反対方向の関係 (警戒区域内 < 区域外) が一貫して観察された。これは 土砂災害警戒区域の指定目的 (= 土石流到達範囲・堆積範囲の保護) と整合する物理的論理であり、 オープンデータの 2 dataset 結合だけで再現可能なことを示している。 特別警戒 (red) は警戒 (yellow) より平均傾斜が高く、両ゾーンの機能差 (建物倒壊リスク vs 注意喚起)が 傾斜分布に明確に刻まれている。

L38 CS立体図 Saturation との Pearson 相関 (r=0.67-0.82) は記事間整合性の量的実証であり、 DoBoX の同シリーズ (LiDAR ファミリ) が同一 DEM から派生していることを裏付ける。 特に slope vs CS darkness の相関は r=0.9+ でCS立体図の暗さ成分は傾斜の最良代理と判明。 これは L38 が視覚解析向きの圧縮表現として有効であることの定量裏付けでもある。

解像度差 (1m vs 50cm) は overview レベルでは認識困難だが、本来解像度のサンプル (図 8) で 50cm の優位性が見え、また 30°超比率と集塊性の解像度依存性も確認された。 ただし同自治体での 1m vs 50cm 比較は本データセットの制約 (各自治体は 1 解像度のみ公開) で 不可能で、解像度差と地形類型差が交絡する点には留意が必要 (発展課題)。

急傾斜地の集塊性 (8 近傍同種率) は仮説 0.85 を下回り 0.6-0.79 だったが、 解像度に応じて単調増加する傾向は確認された。これは「線状地形構造の集塊性は面的指標で過小評価される」 という方法論的発見で、方位 (aspect) を加味した方向性集塊性を計算する発展課題に繋がる。

""" sections.append(("仮説検証と考察", sec_h)) # === Section 11: 発展課題 === sec_dev = f"""

結果 X1 → 仮説 Y1 → 課題 Z1: 警戒区域内外の傾斜分布の物理的解釈深掘り (本記事の発見の深掘り)

結果 X1: 3 エリア共通で警戒区域内 < 区域外の平均傾斜。特別警戒は警戒より急傾斜。 H4 を反対方向で支持する量的傾向。

新仮説 Y1: 警戒区域は「土石流の発生源 (急傾斜面 30°超) + 流下範囲 (15-30°) + 堆積範囲 (扇状地 5°以下)」 の 3 帯を組み合わせて指定されており、平均値だけでは見えない3 段の幾何構造を持つ。

課題 Z1: 各警戒区域ポリゴンを個別に傾斜分布解析。最大傾斜・最小傾斜・分布形状で 3 タイプ (急斜面卓越型 / 流下範囲型 / 堆積範囲型) に分類。 さらに土石流・急傾斜・地すべりの 3 種別ごとに比較すると、各災害タイプの幾何構造が見える。 これにより警戒区域の地形リスク類型を定量化できる。

結果 X2 → 仮説 Y2 → 課題 Z2: 方位を加味した方向性集塊性

結果 X2: 急傾斜の 8 近傍同種率は 0.6-0.79 で仮説 0.85 を下回ったが、 線状地形 (尾根線・崖線) の存在から方向性が示唆される。

新仮説 Y2: 急傾斜セルは等高線に沿った方向 (= 斜面の走向)には連続するが、 直交方向 (斜面方向) には不連続。方位を加味すれば線状集塊度は 0.9+ になる。

課題 Z2: 各セルで方位 (aspect) を計算 (DEM 勾配ベクトルの方向)、 等高線方向 (= 方位 ⊥) のセル間の連続性のみを評価する方向性集塊度を計算。 標準の Moran's I を方位重み付き隣接行列で実装。 これは「線状地形の集塊性」を初めて量化する方法論的貢献になる。

結果 X3 → 仮説 Y3 → 課題 Z3: 急峻地形と古崩壊地の同定

結果 X3: 45° 超急峻セルは 3 エリアで 0.1-0.2% と稀少。これらは特定の崖・崩壊跡に 集塊する可能性。

新仮説 Y3: 急峻セル + その周辺 (50m バッファ) の中傾斜セル混在域は、 古崩壊地形 (paleolandslide) の候補。過去の崩壊跡が地形に残る。

課題 Z3: 各エリアの 45° 超セルのクラスタリング (DBSCAN) で連続塊を抽出。 各クラスタの周辺バッファに含まれる傾斜分布を解析し、古崩壊地形類型 (新鮮崩壊・古崩壊・侵食地形) に分類。 H30 豪雨の崩壊実績 (国交省データ) と空間照合して精度評価。 高精度なら「傾斜ラスタからの古崩壊地形自動検出」を提案できる ⇒ 砂防研究の応用可能。

結果 X4 → 仮説 Y4 → 課題 Z4: 標高ラスタからの傾斜再計算と本データとの照合

結果 X4: L38 CS立体図 Saturation との r=0.67-0.82 は強相関だが完全ではない。 他の派生過程 (例えば曲率と傾斜の重ね合わせ) も影響している。

新仮説 Y4: 標高ラスタ (DTM, dsid 1635/1636) から直接 Horn 法で傾斜を再計算すると、 DoBoX 公式傾斜と Pearson r ≥ 0.99 の極めて強い一致が出る ⇒ 公式が Horn 法を使っている裏付け。

課題 Z4: 標高 GeoTIFF を読み込んで scipy.ndimage.sobel や手計算で 3×3 Horn 法を実装し、本データと相関係数で照合。さらに計算法の違い (Horn vs Zevenbergen) で結果がどう変わるかも比較。教材としても優れた発展課題。

結果 X5 → 仮説 Y5 → 課題 Z5: 多自治体一括の傾斜指標化と地形類型クラスタリング

結果 X5: 本記事は 3 自治体のみ。30°超比率で世羅 < 熊野 < 坂 の階層が見えたが、 サンプル数が少なすぎて地形類型の自然分類は仮説止まり。

新仮説 Y5: 全 22 自治体 (1m 版) で平均傾斜・p90・30°超比率を計算すると、 中山間 / 都市近郊 / 沿岸被災地 / 島嶼 などの地形類型クラスタが自然分類される。

課題 Z5: 22 自治体すべての 1m 版を順次 DL → overview 読込 → 5 階級構成計算で 22 個の特徴ベクトル (5次元) を作る。k-means や階層クラスタリングで類型化。 さらに各自治体の警戒区域率と平均傾斜を散布図にプロット ⇒ 「平均傾斜 vs 警戒区域率の正相関」が量的に確認されれば 「傾斜分布は土砂災害リスクの代理指標」という新提案ができる。

結果 X6 → 仮説 Y6 → 課題 Z6: 1m vs 50cm の同自治体比較 (解像度効果の純粋分離)

結果 X6: 30°超比率は坂町 (50cm) で最大だが、地形類型差と解像度差が交絡。

新仮説 Y6: 同自治体で 1m と 50cm の両方が公開されている自治体 (例: 三原市・呉市・尾道市等) を使えば、 解像度差の純粋効果が分離できる。同自治体での 50cm vs 1m は、急傾斜比率を 1.2-1.5 倍に増やす。

課題 Z6: 同自治体で両解像度公開の市町村を選び、ピクセル単位で空間整合した上で 分布差を計算。さらに 1m を 50cm に再投影 (resampling) して計算上の同 shape にしてから比較すると、 「解像度差の数値効果」が分離可能になる。これは GIS リテラシ教材として優れた題材。

結果 X7 → 仮説 Y7 → 課題 Z7: 過去被災地のセル単位検出 (ML 応用)

結果 X7: 坂町の特別警戒区域 ({AREAS[2]['sed_red_km2']:.2f} km²) と急傾斜セル (45.8%) は 相互に重なるが、過去の被災実績との照合は未実施。

新仮説 Y7: 平成30年7月豪雨で実際に崩壊が発生したセルは傾斜・方位・曲率の 特定パターン (傾斜 30-50° + 凹型曲率 + 北向き等) を持つ。これらの特徴量で機械学習による被災予測が可能。

課題 Z7: 国交省の H30 豪雨被災調査データ (発生位置 GIS) を入手し、 傾斜 + 方位 + 曲率 + L38 CS Sat の特徴でロジスティック回帰や Random Forest による 被災予測モデルを訓練。AUC 0.8 以上が出れば傾斜ベースの土砂災害発生予測モデルとして実装可能。 これは社会実装に直結する研究的発展。

""" sections.append(("発展課題", sec_dev)) # === HTML 書き出し === html = render_lesson( num=39, title="L39 地形傾斜図 2件統合分析 — DEM由来の傾斜角分布が描く地形リスク構造 と 土砂災害警戒区域の物理的論理", tags=["LiDAR", "地形傾斜図", "slope angle", "rasterio", "DEM", "Horn法", "土砂災害警戒区域", "rasterize", "Pearson相関", "L38連携"], time="20-30 分", level="リテラシ + GIS応用", data_label="DoBoX dsid 1630 + 1631 (地形傾斜図 1m + 50cm), 代表 3 自治体 GeoTIFF (Float32 度数)", sections=sections, script_filename="lessons/L39_terrain_slope.py", ) (LESSONS / "L39_terrain_slope.html").write_text(html, encoding="utf-8") print(f" saved L39_terrain_slope.html ({time.time()-t1:.1f}s)", flush=True) print(f"\n=== Final elapsed: {time.time()-t0:.1f}s ===", flush=True)