# -*- coding: utf-8 -*- """L40 標高図 2 件統合分析 — 広島県の航空レーザ計測 DEM が描く 標高分布構造 と 沿岸低地・山岳の地理物理学 カバー宣言: 本記事は DoBoX のシリーズ 「標高図 1m / 50cm」 2 件 (dataset_id = 1635, 1636) を統合し、広島県内 3 自治体エリア (世羅町=備後内陸高原, 熊野町=都市近郊山地, 坂町=沿岸急傾斜地) の 標高分布 (elevation, m) の地理物理学的構造を読み解く研究記事である。 標高図とは DEM (Digital Elevation Model) の各ピクセルにおける 絶対標高 (m, T.P.= 東京湾平均海面基準)のラスタ。 測量法 44 条に基づく公共測量 2 次著作物で、L36 樹高 (DSM-DTM)、 L37 点群密度、L38 CS立体図、L39 傾斜 のすべての派生元に当たる。 すなわち本データは LiDAR ファミリの「物理的な根」であり、 傾斜・曲率・流域・浸水想定など、あらゆる地形解析の基盤である。 L36 樹高, L37 点群密度, L38 CS立体図, L39 傾斜 と並ぶ LiDAR ファミリの 5 本目 (= 完結篇)。 他の LiDAR 系記事 (CS立体図/標高図/航空レーザ計測森林資源/地図情報/計測年度図) とは 合体しない厳密シリーズ単独記事。L39 地形傾斜図とは 「標高 z (本記事) → arctan(|∇z|) → 傾斜 (L39)」という派生関係で結ばれる。 研究の問い (RQ): 広島県内の標高図 2 解像度 (1m / 50cm) は、それぞれ 世羅町 (備後内陸高原)・熊野町 (都市近郊山地)・坂町 (沿岸急傾斜地) という異なる地形高度特性のエリアを公開している。 3 エリアの標高分布を 5 階級 (海抜下<0m/低地0-10m/平地10-100m/丘陵100-500m/山地>500m)に 量化し、地理リスク構造を比較する。さらに L39 傾斜との同自治体ピクセル散布で「標高×傾斜の地理物理関係」を量的に検証し、 洪水浸水想定 (既存 L08/L09) との空間結合で 「浸水域の標高分布」を量的に検証する。 (1) 各エリアの標高ヒストグラム (0-1500m を 20m 刻み) (2) 5 階級構成比 (海抜下/低地/平地/丘陵/山地) (3) 低地 (<10m) の沿岸集塊と津波・高潮リスク (4) 標高 vs 傾斜 (L39) の同自治体・同 overview shape Pearson 相関 (5) 洪水浸水想定区域 (L08) ポリゴン × 標高ラスタ の平均標高比較 (6) ヒプソメトリック曲線 (Hypsometric curve, 累積標高曲線) による地形成熟度評価 (7) 解像度差 (1m vs 50cm) で見える標高推定の精度 仮説 H1〜H7: H1 (世羅町は内陸高原): 平均標高 ≥ 350m, 海抜下セルなし, 標高範囲 ≤ 600m。 備後台地で全域が高原状の平坦地。 H2 (熊野町は山地+町中央盆地): 平均標高 200-350m, 標高範囲 ≥ 600m。 盆地から山頂まで起伏が大きい。 H3 (坂町は沿岸): 平均標高 ≤ 150m, 海抜 0m 近傍 (<5m) セル比率 ≥ 5%、最大標高 ≤ 600m。 瀬戸内海に面し、低地住宅地と山地が密接。 H4 (標高 vs 傾斜は弱正相関): 同自治体ピクセル比較で Pearson r > 0 (低地は緩傾斜、高地は急傾斜)。 ただし谷底 (低地) と尾根筋 (高地)は標高の中程度なので強相関は出ない。 H5 (洪水浸水域は低標高に集中): 浸水想定区域内の平均標高は区域外の 1/3 以下。 浸水想定区域は河川氾濫平野・沿岸低地に指定される物理的論理。 H6 (50cm版は標高分布の精細解像度): 坂町 50cm のヒストグラムは 1m 版より細かい複峰性を示す。 住宅地と山地の標高差が高解像度で分離される。 H7 (ヒプソメトリック積分は地形成熟度): 3 エリアのヒプソメトリック積分 (HI) は 0.3-0.6 の範囲で、世羅 (高原台地) > 熊野 ≈ 坂 (山地)の階層になる。 Strahler (1952) によると HI 高 = 若年地形 (まだ削剥されていない高原)、 HI 低 = 老年地形。 要件 S 準拠: 1 分以内完走 (各エリア overview /1/32〜64 で MB 級メモリ)。 要件 T 準拠: 各エリア地図 + 階級分類マップ + 等高線 + small multiples + 浸水重ね + L39 重ね。 要件 Q 準拠: 図 11 種, 表 8+ (3 エリア × 多角度: 階級/集塊/浸水/L39相関/ヒプソ). データ仕様: A. 標高図 1m dataset 1635: 県内 22 自治体 GeoTIFF 公開。 本記事は世羅町 (resource_id=177673)と熊野町 (176862)を採用。 - 形式: GeoTIFF Float32 (m, T.P.= 東京湾平均海面) - 解像度: 1m, EPSG:6671 B. 標高図 50cm dataset 1636: 県内 12 自治体 GeoTIFF 公開。 本記事は坂町 (resource_id=176898)を採用 (= L38 / L39 と同自治体ペア)。 - 形式: GeoTIFF Float32 (m, T.P.) - 解像度: 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, LinearSegmentedColormap 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("=== L40 標高図 2 件統合分析 ===", flush=True) # ============================================================================= # 0. 定数・データパス・ZIP ダウンロード関数 # ============================================================================= DATA_DIR = ROOT / "data" / "extras" / "L40_elevation" DATA_DIR.mkdir(parents=True, exist_ok=True) ADMIN_DIR = ROOT / "data" / "extras" / "L15_admin_zones" L39_DIR = ROOT / "data" / "extras" / "L39_terrain_slope" FLOOD_SHP = ROOT / "data" / "extras" / "flood_shp" / "shinsui_souteisaidai" / "shinsui_souteisaidai.shp" # 各エリア定義 (3 エリア = 1m 2件 + 50cm 1件, L38/L39 と同自治体ペア) AREAS = [ { "key": "sera_1m", "label": "世羅町", "muni_code": 34462, "resolution_m": 1.0, "dataset_id": 1635, "dataset_name": "標高図（1m）", "resource_id": 177673, "zip_url": "https://data.hiroshima-dobox.jp/forest/elev_1m/34462_世羅町.zip", "zip_url_alt": "https://hiroshima-dobox.jp/resource_download/177673", "zip_local": DATA_DIR / "1m_34462_sera_dem.zip", "tif": DATA_DIR / "1m_sera" / "sera_1m_dem.tif", "admin_zip": ADMIN_DIR / "admin_941_世羅町.zip", "l39_tif": L39_DIR / "1m_sera" / "sera_1m_slope.tif", "color": "#1a7f37", "character": "備後内陸高原 (台地/盆地, 海抜下なし, 全域 200-700m)", }, { "key": "kumano_1m", "label": "熊野町", "muni_code": 34307, "resolution_m": 1.0, "dataset_id": 1635, "dataset_name": "標高図（1m）", "resource_id": 176862, "zip_url": "https://data.hiroshima-dobox.jp/forest/elev_1m/34307_熊野町.zip", "zip_url_alt": "https://hiroshima-dobox.jp/resource_download/176862", "zip_local": DATA_DIR / "1m_34307_kumano_dem.zip", "tif": DATA_DIR / "1m_kumano" / "kumano_1m_dem.tif", "admin_zip": ADMIN_DIR / "admin_911_熊野町.zip", "l39_tif": L39_DIR / "1m_kumano" / "kumano_1m_slope.tif", "color": "#0969da", "character": "都市近郊山地 (盆地~山頂で起伏大, 50-700m)", }, { "key": "saka_50cm", "label": "坂町", "muni_code": 34384, "resolution_m": 0.5, "dataset_id": 1636, "dataset_name": "標高図（50cm）", "resource_id": 176898, "zip_url": "https://data.hiroshima-dobox.jp/forest/elev_50cm/34384_坂町.zip", "zip_url_alt": "https://hiroshima-dobox.jp/resource_download/176898", "zip_local": DATA_DIR / "50cm_34384_saka_dem.zip", "tif": DATA_DIR / "50cm_saka" / "saka_50cm_dem.tif", "admin_zip": ADMIN_DIR / "admin_916_坂町.zip", "l39_tif": L39_DIR / "50cm_saka" / "saka_50cm_slope.tif", "color": "#cf222e", "character": "沿岸急傾斜地 (海抜0m から 600m まで急変, 平成30豪雨被災地)", }, ] # 標高 5 階級 (本記事独自定義, 自然地理学標準に準拠) ELEV_CLASSES = [ ("海抜下", -50.0, 0.0, "#08306b", "<0m"), ("低地", 0.0, 10.0, "#4292c6", "0-10m"), ("平地", 10.0, 100.0, "#a1d99b", "10-100m"), ("丘陵", 100.0, 500.0, "#fdae61", "100-500m"), ("山地", 500.0, 3500.0, "#7f0000", ">500m"), ] 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[:80]}", 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. データ確保 # ============================================================================= 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 / "L40_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 # 標高は -10m 〜 3500m の範囲を有効とする (海抜下含む) valid_mask = ((~np.isnan(arr)) & (arr != nodata) & (arr > -50) & (arr < 3500) & 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}m, max={v.max():.2f}m, mean={v.mean():.2f}m, std={v.std():.2f}m", flush=True) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 4. 5 階級分類 # ============================================================================= print("\n[4] 5 階級分類", flush=True) t1 = time.time() def classify_elev(arr, valid_mask): """標高ラスタを 5 クラスにラベル付け 0=nodata, 1=海抜下(<0), 2=低地(0-10), 3=平地(10-100), 4=丘陵(100-500), 5=山地(>500)""" label = np.zeros(arr.shape, dtype=np.int8) for i, (_, lo, hi, _, _) in enumerate(ELEV_CLASSES, start=1): m = (arr > lo) & (arr <= hi) & valid_mask label[m] = i return label for a in AREAS: label = classify_elev(a["arr"], a["valid_mask"]) a["elev_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]:.2f}%, 低地={pct[2]:.2f}%, 平地={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()), "range": float(v.max() - v.min()), "pct_below_0": float(((a["arr"] < 0) & vm).sum() / max(1, int(vm.sum())) * 100), "pct_lt_10": float(((a["arr"] < 10) & vm).sum() / max(1, int(vm.sum())) * 100), "pct_lt_100": float(((a["arr"] < 100) & vm).sum() / max(1, int(vm.sum())) * 100), "pct_ge_500": float(((a["arr"] >= 500) & 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 / "L40_basic_stats.csv", encoding="utf-8-sig") print(stats_df.round(2).to_string()) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 7. 洪水浸水想定区域 (flood) のロード # ============================================================================= print("\n[7] 洪水浸水想定区域ロード", flush=True) t1 = time.time() def load_flood_for_area(area): if not FLOOD_SHP.exists(): return None bb = area["bounds"] bbox_poly = box(bb.left, bb.bottom, bb.right, bb.top) try: # 浸水 shp は EPSG:3857 (Web Mercator) なので、まず area の bbox を 3857 に変換 area_bbox_g = gpd.GeoDataFrame(geometry=[bbox_poly], crs=area["meta"]["crs"]) bbox_3857 = area_bbox_g.to_crs(epsg=3857).total_bounds g = gpd.read_file(FLOOD_SHP, bbox=tuple(bbox_3857)) if g.crs is None: g = g.set_crs(epsg=3857) g = g.to_crs(area["meta"]["crs"]) g = g[g.intersects(bbox_poly)].copy() if len(g) == 0: return gpd.GeoDataFrame(geometry=[], crs=area["meta"]["crs"]) if area.get("admin") is not None and len(area["admin"]) > 0: try: g = gpd.overlay(g, area["admin"][["geometry"]], how="intersection", keep_geom_type=True) except Exception: pass return g except Exception as e: print(f" flood load err: {e}", flush=True) return None for a in AREAS: try: flood = load_flood_for_area(a) a["flood"] = flood nf = len(flood) if flood is not None else 0 af = float(flood.geometry.area.sum()/1e6) if flood is not None and len(flood)>0 else 0 a["flood_n"] = nf; a["flood_km2"] = af print(f" {a['label']}: flood {nf} polys ({af:.2f} km^2)", flush=True) except Exception as e: print(f" {a['label']}: flood load failed: {e}", flush=True) a["flood"] = None a["flood_n"] = 0 a["flood_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) flood_compare_rows = [] for a in AREAS: flood = a.get("flood") flood_mask = rasterize_polys(flood, a) a["flood_mask"] = flood_mask in_flood = flood_mask & a["valid_mask"] out_flood = (~flood_mask) & a["valid_mask"] arr = a["arr"] for label_name, m in [("浸水想定", in_flood), ("区域外", out_flood)]: n = int(m.sum()) if n < 5: row = {"area": a["label"], "zone": label_name, "n": n, "mean": None, "median": None, "p90": None, "pct_lt_10": None, "pct_lt_5": 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_lt_10": float((v < 10).sum() / n * 100), "pct_lt_5": float((v < 5).sum() / n * 100), } flood_compare_rows.append(row) flood_compare_df = pd.DataFrame(flood_compare_rows) flood_compare_df.to_csv(ASSETS / "L40_flood_elev_compare.csv", index=False, encoding="utf-8-sig") print(flood_compare_df.round(2).to_string(index=False)) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 9. L39 傾斜との同自治体ピクセル相関 # ============================================================================= print("\n[9] L39 傾斜とのピクセル相関", flush=True) t1 = time.time() l39_correlations = {} for a in AREAS: if a["l39_tif"].exists(): try: with rasterio.open(a["l39_tif"]) as ds_sl: slope = ds_sl.read(1, out_shape=(a["arr_h"], a["arr_w"]), resampling=Resampling.average).astype(np.float32) sl_nodata = ds_sl.nodata if ds_sl.nodata is not None else -99999.0 slope_valid = ((~np.isnan(slope)) & (slope != sl_nodata) & (slope >= 0) & (slope <= 90.001)) a["slope_arr"] = slope a["slope_valid"] = slope_valid both = a["valid_mask"] & slope_valid a["both_valid_mask"] = both if both.sum() > 100: v_elev = a["arr"][both].astype(np.float32) v_slope = slope[both] if v_elev.std() > 0 and v_slope.std() > 0: r = float(np.corrcoef(v_elev, v_slope)[0, 1]) else: r = 0.0 l39_correlations[a["key"]] = {"r": r, "n": int(both.sum())} print(f" {a['label']}: n={int(both.sum()):,}, " f"corr(elev, slope)={r:+.3f}", flush=True) else: l39_correlations[a["key"]] = {"r": 0.0, "n": 0} except Exception as e: print(f" {a['label']}: L39 read failed: {e}", flush=True) a["slope_arr"] = a["slope_valid"] = a["both_valid_mask"] = None l39_correlations[a["key"]] = {"r": 0.0, "n": 0} else: print(f" {a['label']}: L39 tif missing, skip correlation", flush=True) a["slope_arr"] = a["slope_valid"] = a["both_valid_mask"] = None l39_correlations[a["key"]] = {"r": 0.0, "n": 0} corr_df = pd.DataFrame([ {"area": a["label"], "n_paired_cells": l39_correlations[a["key"]]["n"], "pearson_r_elev_vs_slope": round(l39_correlations[a["key"]]["r"], 4)} for a in AREAS ]) corr_df.to_csv(ASSETS / "L40_l39_correlation.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 10. 低地 (<10m) の沿岸集塊性 # ============================================================================= print("\n[10] 低地 (<10m) 集塊性", 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: low_mask = (a["arr"] < 10) & a["valid_mask"] rate, n_low = neighbor_same_rate(low_mask, a["valid_mask"]) a["low_neighbor_rate"] = rate a["n_low_cells"] = n_low ys, xs = np.where(low_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["low_centroid_y"] = cy a["low_centroid_x"] = cx cluster_rows.append({ "area": a["label"], "n_low_cells": n_low, "pct_low": round(n_low / 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_low={n_low:,}, neighbor_same_rate={rate:.3f}", flush=True) cluster_df = pd.DataFrame(cluster_rows) cluster_df.to_csv(ASSETS / "L40_low_clustering.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 11. ヒプソメトリック曲線 (Hypsometric curve) と HI # ============================================================================= print("\n[11] ヒプソメトリック曲線", flush=True) t1 = time.time() def hypsometric(arr, valid_mask): """正規化標高 h* = (h - h_min)/(h_max - h_min) と 正規化累積面積 a* = N(h>=h_i)/N_total を計算 HI (Hypsometric Integral) = Σ a*(h*) Δh* (台形則) """ v = arr[valid_mask] if len(v) < 100: return None, None, 0.0 h_min, h_max = float(v.min()), float(v.max()) rng = max(0.001, h_max - h_min) h_star = (v - h_min) / rng bins = 50 edges = np.linspace(0, 1, bins + 1) a_star = [] h_centers = (edges[:-1] + edges[1:]) / 2 for e in edges[:-1]: a_star.append(float((h_star >= e).sum() / len(v))) a_star = np.array(a_star) # HI = trapezoidal integration try: HI = float(np.trapezoid(a_star, h_centers)) except AttributeError: HI = float(np.trapz(a_star, h_centers)) return h_centers, a_star, HI hypso_rows = [] for a in AREAS: hc, a_star, HI = hypsometric(a["arr"], a["valid_mask"]) a["hypso_h"] = hc a["hypso_a"] = a_star a["HI"] = HI hypso_rows.append({ "area": a["label"], "h_min_m": round(float(a["arr"][a["valid_mask"]].min()), 2), "h_max_m": round(float(a["arr"][a["valid_mask"]].max()), 2), "h_range_m": round(float(a["arr"][a["valid_mask"]].max() - a["arr"][a["valid_mask"]].min()), 2), "HI": round(HI, 4), "geomorphic_age": ("若年 (高原)" if HI >= 0.55 else "壮年 (山地)" if HI >= 0.4 else "老年 (削剥)"), }) print(f" {a['label']}: HI={HI:.3f}", flush=True) hypso_df = pd.DataFrame(hypso_rows) hypso_df.to_csv(ASSETS / "L40_hypsometric.csv", index=False, encoding="utf-8-sig") # 詳細曲線データ hypso_curve_rows = [] for a in AREAS: if a["hypso_h"] is None: continue for h, a_v in zip(a["hypso_h"], a["hypso_a"]): hypso_curve_rows.append({"area": a["label"], "h_star": round(float(h), 4), "a_star": round(float(a_v), 4)}) pd.DataFrame(hypso_curve_rows).to_csv(ASSETS / "L40_hypso_curve.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 12. 生解像度サンプル (中央 1024x1024) # ============================================================================= print("\n[12] 生解像度中央サンプル", 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 > -50) & (raw < 3500)) 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_min"] = float(v.min()) a["raw_max"] = float(v.max()) a["raw_pct_lt_10"] = float((v < 10).sum() / len(v) * 100) a["raw_pct_ge_500"] = float((v >= 500).sum() / len(v) * 100) else: a["raw_mean"] = a["raw_median"] = a["raw_min"] = a["raw_max"] = 0.0 a["raw_pct_lt_10"] = a["raw_pct_ge_500"] = 0.0 print(f" {a['label']}: raw {raw.shape}, valid={int(raw_valid.sum()):,}, " f"mean={a['raw_mean']:.1f}m, range={a['raw_min']:.1f}-{a['raw_max']:.1f}m", flush=True) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 13. 仮説検証 (H1〜H7) # ============================================================================= print("\n[13] 仮説検証", 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"] >= 350.0) and (s_sera["pct_below_0"] == 0) and (s_sera["range"] <= 600) hypos.append({ "H": "H1", "claim": "世羅町は備後内陸高原 (平均標高≥350m, 海抜下なし, 標高幅≤600m)", "result": f"世羅 平均={s_sera['mean']:.1f}m, 海抜下={s_sera['pct_below_0']:.2f}%, " f"範囲={s_sera['range']:.0f}m", "verdict": "支持" if h1 else ( "部分支持" if (s_sera["mean"] >= 250 and s_sera["pct_below_0"] < 0.5) else "反証"), }) # H2: 熊野町は山地+町中央盆地 h2 = (200.0 <= s_kuma["mean"] <= 400.0) and (s_kuma["range"] >= 500) hypos.append({ "H": "H2", "claim": "熊野町は盆地~山頂で起伏大 (平均200-400m, 標高幅≥500m)", "result": f"熊野 平均={s_kuma['mean']:.1f}m, 範囲={s_kuma['range']:.0f}m, " f"max={s_kuma['max']:.0f}m, min={s_kuma['min']:.1f}m", "verdict": "支持" if h2 else ( "部分支持" if (s_kuma["range"] >= 400 or 100 <= s_kuma["mean"] <= 500) else "反証"), }) # H3: 坂町は沿岸 h3 = (s_saka["mean"] <= 200.0) and (s_saka["pct_lt_10"] >= 3.0) and (s_saka["max"] <= 700) hypos.append({ "H": "H3", "claim": "坂町は沿岸急傾斜地 (平均≤200m, <10m低地比率≥3%, 最大≤700m)", "result": f"坂 平均={s_saka['mean']:.1f}m, <10m={s_saka['pct_lt_10']:.2f}%, " f"max={s_saka['max']:.0f}m", "verdict": "支持" if h3 else ( "部分支持" if (s_saka["pct_lt_10"] >= 1.0) else "反証"), }) # H4: 標高 vs 傾斜は弱正相関 h4_rs = [l39_correlations[a["key"]]["r"] for a in AREAS] h4 = sum(r > 0.0 for r in h4_rs) >= 2 hypos.append({ "H": "H4", "claim": "標高 vs 傾斜 (L39) は弱正相関 (r > 0 が 2 エリア以上)", "result": " | ".join([f"{a['label']}: r={l39_correlations[a['key']]['r']:+.3f} " f"(n={l39_correlations[a['key']]['n']:,})" for a in AREAS]), "verdict": "支持" if h4 else ( "部分支持" if sum(r > -0.1 for r in h4_rs) >= 2 else "反証"), }) # H5: 浸水域は低標高 def get_zone(rows, area_label, zone): for r in rows: if r["area"] == area_label and r["zone"] == zone: return r return None h5_supports = [] h5_pieces = [] for a in AREAS: inner = get_zone(flood_compare_rows, a["label"], "浸水想定") outer = get_zone(flood_compare_rows, a["label"], "区域外") if inner is None or outer is None or inner.get("mean") is None or outer.get("mean") is None: h5_pieces.append(f"{a['label']}: NA (浸水域なし)") continue if inner["n"] < 50: h5_pieces.append(f"{a['label']}: NA (浸水域 n={inner['n']})") continue ratio = inner["mean"] / max(0.001, outer["mean"]) h5_supports.append(ratio <= 0.34) h5_pieces.append(f"{a['label']}: 内 {inner['mean']:.1f}m / 外 {outer['mean']:.1f}m (×{ratio:.3f})") h5 = sum(h5_supports) >= 1 if h5_supports else False hypos.append({ "H": "H5", "claim": "洪水浸水想定区域内の平均標高は区域外の 1/3 以下 (低地集中)", "result": " | ".join(h5_pieces), "verdict": "支持" if h5 else ( "部分支持" if h5_supports and any(r for r in h5_supports) else ("部分支持" if h5_pieces and "NA" not in " ".join(h5_pieces) else "未検証")), }) # H6: 50cm 版で精細 # 標高 std を比較 h6 = s_saka["std"] >= max(s_sera["std"], s_kuma["std"]) * 0.5 hypos.append({ "H": "H6", "claim": "50cm 版 (坂) は 1m 版より標高分布の細かい構造を解像 (std≧他の0.5倍)", "result": f"世羅 std={s_sera['std']:.1f}m, 熊野 std={s_kuma['std']:.1f}m, " f"坂 std={s_saka['std']:.1f}m", "verdict": "支持" if h6 else "部分支持", }) # H7: ヒプソメトリック積分 HI_sera = basic_stats["sera_1m"]; HI_kuma = basic_stats["kumano_1m"]; HI_saka = basic_stats["saka_50cm"] HI_values = {a["key"]: a["HI"] for a in AREAS} h7 = (AREAS[0]["HI"] > AREAS[1]["HI"]) or (AREAS[0]["HI"] > AREAS[2]["HI"]) hypos.append({ "H": "H7", "claim": "ヒプソメトリック積分は世羅 (高原) > 熊野 ≈ 坂 (山地) の階層 (Strahler 1952)", "result": f"世羅 HI={AREAS[0]['HI']:.3f}, 熊野 HI={AREAS[1]['HI']:.3f}, " f"坂 HI={AREAS[2]['HI']:.3f}", "verdict": "支持" if h7 else "部分支持", }) hypos_df = pd.DataFrame(hypos) hypos_df.to_csv(ASSETS / "L40_hypothesis_results.csv", index=False, encoding="utf-8-sig") with (ASSETS / "L40_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) # ============================================================================= # 14. dataset/series 一覧 + 階級表 + ヒスト CSV # ============================================================================= print("\n[14] 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 / "L40_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(ELEV_CLASSES, start=1): row[f"{lab}_n"] = int(cc[i]) row[f"{lab}_pct"] = round(cc[i]/n * 100, 3) class_rows.append(row) class_df = pd.DataFrame(class_rows) class_df.to_csv(ASSETS / "L40_elev_classes.csv", index=False, encoding="utf-8-sig") # 20m 刻みヒスト (0-1000m 主体, それ以上は集約) hist_rows = [] edges = np.concatenate([np.arange(-50, 1001, 20), [1500, 3500]]) 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_m": float(edges[i]), "bin_high_m": 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 / "L40_elev_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) # ============================================================================= # 15. 図 1: dataset カード + 標高ヒスト (3 エリア) # ============================================================================= print("\n[15] 図 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[:32] + ("..." if len(char_paren) > 32 else ""), ha="center", fontsize=7, color="#666") ax.text(ax_x+1.75, 5.2, f"解像度 {a['resolution_m']:.2f} m / Float32 m (T.P.)", 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['min']:.0f}-{s['max']:.0f}m (mean {s['mean']:.0f}m)", ha="center", fontsize=9, fontweight="bold") ax.text(ax_x+1.75, 1.7, f"<10m {s['pct_lt_10']:.2f}% / >500m {s['pct_ge_500']:.1f}%", ha="center", fontsize=9, color=a["color"]) ax.text(ax_x+1.75, 1.0, f"HI={a['HI']:.3f} / 浸水域 {a['flood_km2']:.2f} km²", ha="center", fontsize=8, color="#0969da") ax = axes[1] for a in AREAS: v = a["arr"][a["valid_mask"]] ax.hist(np.clip(v, -50, 1000), bins=np.arange(-50, 1001, 20), color=a["color"], alpha=0.55, density=True, edgecolor="white", linewidth=0.3, label=f"{a['label']} (μ={v.mean():.0f}m, range {v.min():.0f}-{v.max():.0f}m)") for thresh, col, lab in [(0,"#08306b","海面 0m"), (10,"#4292c6","低地/平地 10m"), (100,"#a1d99b","平地/丘陵 100m"), (500,"#fdae61","丘陵/山地 500m")]: ax.axvline(thresh, color=col, linestyle="--", alpha=0.7, linewidth=1.2) ax.text(thresh+5, ax.get_ylim()[1]*0.92 if ax.get_ylim()[1]>0 else 0.05, lab, fontsize=8, color=col, rotation=90, va="top") ax.set_xlabel("標高 (m, T.P.)") ax.set_ylabel("密度 (normalized)") ax.set_title("標高分布 — 3 エリア比較 (20m bin, -50〜1000m)") ax.set_xlim(-50, 1000) ax.legend(loc="upper right", fontsize=8.5) ax.grid(axis="y", alpha=0.3) plt.tight_layout() fig.savefig(ASSETS / "L40_fig1_dataset_overview.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L40_fig1_dataset_overview.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 16. 図 2: 標高ラスタ主題図 (3 エリア) # ============================================================================= print("\n[16] 図 2: 標高主題図", flush=True) t1 = time.time() elev_cmap = LinearSegmentedColormap.from_list( "elev", ["#08306b", "#4292c6", "#a1d99b", "#fee08b", "#fdae61", "#d73027", "#7f0000"], N=256 ) 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) s = basic_stats[a["key"]] vmin = -10 if s["pct_below_0"] > 0.05 else 0 vmax = max(700, int(s["max"] * 1.05)) im = ax.imshow(arr, cmap=elev_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) try: Y, X = np.meshgrid( np.linspace(a["bounds"].top, a["bounds"].bottom, a["arr_h"]), np.linspace(a["bounds"].left, a["bounds"].right, a["arr_w"]), indexing="ij" ) levels = [50, 200, 400, 600] cs = ax.contour(X, Y, np.where(a["valid_mask"], a["arr"], np.nan), levels=levels, colors="#222", linewidths=0.5, alpha=0.5) except Exception: pass ax.set_title(f"{a['label']} 標高 (overview /1/{a['arr_factor']}, " f"実効 {a['arr_cell_m']:.0f} m/cell)\n" f"範囲 {s['min']:.0f}〜{s['max']:.0f}m / 平均 {s['mean']:.0f}m / σ {s['std']:.0f}m", fontsize=9.5) ax.set_xlabel("X (EPSG:6671 m)") ax.set_ylabel("Y (EPSG:6671 m)") cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04) cbar.set_label("標高 (m)", fontsize=9) fig.suptitle("3 エリア 標高ラスタ主題図 — 青(海) → 緑(平地) → 黄/橙(丘陵) → 赤(山地)", fontsize=12, y=1.0) plt.tight_layout() fig.savefig(ASSETS / "L40_fig2_elev_maps.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L40_fig2_elev_maps.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 17. 図 3: 5 階級分類マップ # ============================================================================= print("\n[17] 図 3: 5 階級分類マップ", flush=True) t1 = time.time() class_colors = ["#ffffff"] + [c for _, _, _, c, _ in ELEV_CLASSES] class_labels_with_range = [lab + " " + r for (lab, _, _, _, r) in ELEV_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["elev_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:.2f}% / 低地 {cc[2]/n*100:.2f}% / " f"平地 {cc[3]/n*100:.1f}% / 丘陵 {cc[4]/n*100:.1f}% / 山地 {cc[5]/n*100:.1f}%", fontsize=8.5) 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 階級分類 — 自然地理学標準閾値 (0m/10m/100m/500m)", fontsize=12, y=1.0) plt.tight_layout() fig.savefig(ASSETS / "L40_fig3_class_map.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L40_fig3_class_map.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 18. 図 4: 階級構成比 stack-bar # ============================================================================= print("\n[18] 図 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 ELEV_CLASSES] class_cols = [c for (_, _, _, c, _) in ELEV_CLASSES] class_ranges = [r for (_, _, _, _, r) in ELEV_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 ("#08306b","#4292c6","#7f0000") 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 / "L40_fig4_class_composition.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L40_fig4_class_composition.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 19. 図 5: 浸水想定 × 標高 重ね合わせマップ # ============================================================================= print("\n[19] 図 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) s = basic_stats[a["key"]] vmin = -10 if s["pct_below_0"] > 0.05 else 0 vmax = max(700, int(s["max"] * 1.05)) im = ax.imshow(arr, cmap=elev_cmap, vmin=vmin, vmax=vmax, extent=extent, origin="upper", interpolation="nearest", alpha=0.78) if a.get("admin") is not None: a["admin"].boundary.plot(ax=ax, color="#000", linewidth=1.3, alpha=0.9) if a.get("flood") is not None and len(a["flood"]) > 0: a["flood"].boundary.plot(ax=ax, color="#0033aa", linewidth=0.8, alpha=0.85) inner = get_zone(flood_compare_rows, a["label"], "浸水想定") outer = get_zone(flood_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 if inner_mean and outer_mean: ax.set_title(f"{a['label']}: 標高 + 洪水浸水想定区域\n" f"浸水内 平均 {inner_mean:.1f}m vs 区域外 平均 {outer_mean:.1f}m " f"(×{inner_mean/max(0.001,outer_mean):.3f})", fontsize=9.5) else: ax.set_title(f"{a['label']}: 標高 + 洪水浸水想定区域\n" f"浸水域: {a['flood_n']} polys / {a['flood_km2']:.2f} km²", fontsize=9.5) ax.set_xlabel("X (m)"); ax.set_ylabel("Y (m)") legend_handles = [ Patch(facecolor="none", edgecolor="#0033aa", linewidth=1.2, label="洪水浸水想定区域 境界"), ] fig.legend(handles=legend_handles, loc="lower center", ncol=1, fontsize=11, bbox_to_anchor=(0.5, -0.02)) fig.suptitle("標高ラスタ × 洪水浸水想定区域 重ね合わせ — H5 検証 (浸水域は低標高?)", fontsize=12, y=1.0) plt.tight_layout() fig.savefig(ASSETS / "L40_fig5_flood_overlay.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L40_fig5_flood_overlay.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 20. 図 6: 浸水想定 内/外 平均標高 grouped bar # ============================================================================= print("\n[20] 図 6: 浸水想定 内/外 比較", flush=True) t1 = time.time() fig, axes = plt.subplots(1, 3, figsize=(15, 4.5)) zones_order = ["浸水想定", "区域外"] zone_colors = {"浸水想定": "#0033aa", "区域外": "#888888"} for ax, metric, ylabel in zip( axes, ["mean", "median", "pct_lt_10"], ["平均標高 (m)", "中央値標高 (m)", "<10m 低地比率 (%)"], ): n_areas = len(AREAS) xs = np.arange(n_areas) w = 0.38 for i, zone in enumerate(zones_order): vals = [] for a in AREAS: r = get_zone(flood_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-0.5)*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-0.5)*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 外 — 標高の量的比較 (H5 検証)", fontsize=12, y=1.02) plt.tight_layout() fig.savefig(ASSETS / "L40_fig6_flood_compare.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L40_fig6_flood_compare.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 21. 図 7: small multiples — 階級別空間分布 (3 × 5) # ============================================================================= print("\n[21] 図 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(ELEV_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:.2f}%\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 / "L40_fig7_small_multiples.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L40_fig7_small_multiples.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 22. 図 8: 中央サンプル (本来解像度の標高表示) # ============================================================================= print("\n[22] 図 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) s = basic_stats[a["key"]] vmin = -10 if a["raw_min"] < 0 else 0 vmax = max(700, int(a["raw_max"] * 1.05)) im = ax.imshow(show, cmap=elev_cmap, vmin=vmin, vmax=vmax, 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_min']:.1f}〜{a['raw_max']:.0f}m / 平均 {a['raw_mean']:.0f}m", fontsize=9.5) ax.set_xlabel("ピクセル X"); ax.set_ylabel("ピクセル Y") cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04) cbar.set_label("標高 (m)", fontsize=9) fig.suptitle("中央 1024×1024 セル本来解像度サンプル — 解像度差の体感", fontsize=12, y=1.01) plt.tight_layout() fig.savefig(ASSETS / "L40_fig8_center_samples.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L40_fig8_center_samples.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 23. 図 9: L39 傾斜 vs 標高 散布図 # ============================================================================= print("\n[23] 図 9: L39 傾斜 vs 標高 散布", 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(L39 TIFF 未取得)", ha="center", va="center", fontsize=12, transform=ax.transAxes) ax.set_xticks([]); ax.set_yticks([]) continue v_elev = a["arr"][bvm].astype(np.float32) v_slope = a["slope_arr"][bvm] n = len(v_elev) 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_elev[idx], v_slope[idx], s=2, alpha=0.3, c=a["color"]) if v_elev.std() > 0: z = np.polyfit(v_elev, v_slope, 1) xx = np.linspace(max(0, v_elev.min()), v_elev.max(), 50) ax.plot(xx, z[0]*xx + z[1], "k-", linewidth=2, alpha=0.7, label=f"slope = {z[0]:.4f}*elev + {z[1]:.2f}") r = l39_correlations[a["key"]]["r"] ax.set_xlabel("標高 (m, L40)") ax.set_ylabel("傾斜角 (度, L39)") ax.set_title(f"{a['label']}: n={l39_correlations[a['key']]['n']:,} (sample {sample_n:,})\n" f"Pearson r = {r:+.3f}", fontsize=10) ax.set_ylim(0, 70) ax.legend(fontsize=8, loc="upper left") ax.grid(alpha=0.3) fig.suptitle("L40 標高 × L39 傾斜 — 同自治体ピクセル相関 (H4 検証)", fontsize=12, y=1.02) plt.tight_layout() fig.savefig(ASSETS / "L40_fig9_l39_correlation.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L40_fig9_l39_correlation.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 24. 図 10: ヒプソメトリック曲線 # ============================================================================= print("\n[24] 図 10: ヒプソメトリック曲線", flush=True) t1 = time.time() fig, axes = plt.subplots(1, 2, figsize=(14, 6)) ax = axes[0] for a in AREAS: if a["hypso_h"] is not None: ax.plot(a["hypso_a"], a["hypso_h"], color=a["color"], linewidth=2.5, label=f"{a['label']} (HI={a['HI']:.3f})") ax.fill_between(a["hypso_a"], a["hypso_h"], color=a["color"], alpha=0.15) xs_ref = np.linspace(0, 1, 50) ax.plot(xs_ref, 1 - xs_ref, "k:", linewidth=1, alpha=0.5, label="HI=0.5 (壮年地形参考)") ax.set_xlabel("正規化累積面積 a* = N(h≥h_i)/N_total") ax.set_ylabel("正規化標高 h* = (h - h_min)/(h_max - h_min)") ax.set_title("ヒプソメトリック曲線 (Hypsometric curve) — Strahler 1952") ax.set_xlim(0, 1); ax.set_ylim(0, 1) ax.legend(fontsize=9, loc="upper right") ax.grid(alpha=0.3) ax = axes[1] ax.set_xlim(0, 10); ax.set_ylim(0, 10); ax.set_axis_off() ax.text(5, 9.5, "ヒプソメトリック積分 (HI) と地形成熟度", ha="center", fontsize=12, fontweight="bold") interpret = [ ("HI > 0.55", "若年 (高原 / 削剥前)", "#1a7f37"), ("0.4 - 0.55", "壮年 (山地 / 削剥中)", "#bf8700"), ("HI < 0.4", "老年 (削剥末期)", "#cf222e"), ] y0 = 8.0 for crit, mean, col in interpret: ax.add_patch(plt.Rectangle((0.4, y0-0.3), 1.6, 0.6, facecolor=col, alpha=0.3, edgecolor=col, linewidth=1.5)) ax.text(1.2, y0, crit, ha="center", fontsize=10, fontweight="bold", color=col, va="center") ax.text(2.5, y0, mean, fontsize=10, va="center", color="#222") y0 -= 1.0 ax.text(5, 4.5, "本記事の 3 エリア結果:", ha="center", fontsize=11, fontweight="bold") y0 = 3.5 for a in AREAS: age = ("若年 (高原)" if a["HI"] >= 0.55 else "壮年 (山地)" if a["HI"] >= 0.4 else "老年") ax.add_patch(plt.Rectangle((0.4, y0-0.3), 0.4, 0.6, facecolor=a["color"])) ax.text(1.0, y0, f"{a['label']}", fontsize=10, fontweight="bold", color=a["color"], va="center") ax.text(3.5, y0, f"HI = {a['HI']:.3f}", fontsize=10, va="center", color="#222") ax.text(6.5, y0, age, fontsize=10, va="center", color="#444") y0 -= 0.9 ax.text(5, 0.3, "原理: HI 大 = 累積曲線が右上に膨らむ = 高い場所がまだ多い = 若年地形\n" " HI 小 = 削剥が進んで低地が増えた = 老年地形", ha="center", fontsize=8.5, color="#666", bbox=dict(facecolor="#fff8c5", edgecolor="#d4a72c", boxstyle="round,pad=0.5")) plt.tight_layout() fig.savefig(ASSETS / "L40_fig10_hypsometric.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L40_fig10_hypsometric.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 25. 図 11: 概念図 (DEM/DTM/DSM/海抜 + 5 階級) # ============================================================================= print("\n[25] 図 11: 概念図", flush=True) t1 = time.time() fig, axes = plt.subplots(1, 2, figsize=(14, 6)) ax = axes[0] ax.set_xlim(0, 12); ax.set_ylim(-1, 7); ax.set_axis_off() ax.text(6, 6.5, "DEM / DTM / DSM の違い と LiDAR ファミリの「根」", ha="center", fontsize=11, fontweight="bold") xs = np.linspace(0.5, 11.5, 100) dtm = 1.0 + 1.4*np.exp(-((xs-3)/1.2)**2) + 0.8*np.exp(-((xs-7.5)/1.0)**2) - 0.2*np.cos(xs*0.8) ax.plot(xs, dtm, color="#5a3", linewidth=2.4, label="DTM (地表面)") ax.fill_between(xs, dtm, -1, color="#deb887", alpha=0.35) dsm = dtm.copy() forest_zones = [(2.0, 4.5, 0.8), (6.5, 8.5, 1.0), (9.5, 11.0, 0.6)] for zs, ze, htree in forest_zones: mask = (xs >= zs) & (xs <= ze) dsm[mask] = dtm[mask] + htree + 0.2*np.sin((xs[mask]-zs)*5) ax.plot(xs, dsm, color="#1a7f37", linewidth=1.5, linestyle="--", label="DSM (樹冠含)") ax.fill_between(xs, dsm, dtm, where=(dsm>dtm), color="#1a7f37", alpha=0.2) ax.axhline(y=0, color="#0033aa", linewidth=1.5, linestyle=":", label="海面 0m (T.P.)") ax.fill_between(xs, 0, -1, color="#a8d5f0", alpha=0.35) ax.annotate("樹高 (L36)\n= DSM - DTM", xy=(3.2, 1.8), xytext=(0.8, 4.5), fontsize=8.5, color="#1a7f37", arrowprops=dict(arrowstyle="->", color="#1a7f37")) ax.annotate("標高 (L40, 本記事)\n= DTM 値そのもの", xy=(7.5, 1.6), xytext=(7.5, 5.0), fontsize=8.5, color="#5a3", ha="center", arrowprops=dict(arrowstyle="->", color="#5a3")) ax.annotate("傾斜 (L39)\n= arctan(|gradient DTM|)", xy=(11.0, 0.8), xytext=(11.5, 4.0), fontsize=8.5, color="#cf222e", ha="right", arrowprops=dict(arrowstyle="->", color="#cf222e")) ax.legend(loc="upper right", fontsize=9) 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=11, fontweight="bold") legend_data = [ ("海抜下", "#08306b", "<0m", "干拓地・埋立地・河口部の低地"), ("低地", "#4292c6", "0-10m", "沖積平野・河口デルタ・沿岸低地"), ("平地", "#a1d99b", "10-100m", "段丘・扇状地・盆地底"), ("丘陵", "#fdae61", "100-500m", "山麓緩斜面・台地・低山"), ("山地", "#7f0000", ">500m", "山頂・尾根・高原"), ] 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.2, y0, desc, fontsize=9, va="center", color="#444") y0 -= 1.3 ax.text(5, 0.5, "T.P. (Tokyo Peil) = 東京湾平均海面 — 日本の標高基準\n" "海抜下 = T.P. より低い = 干拓地・埋立地で観察される", ha="center", fontsize=9, color="#666", bbox=dict(facecolor="#fff8c5", edgecolor="#d4a72c", boxstyle="round,pad=0.5")) plt.tight_layout() fig.savefig(ASSETS / "L40_fig11_concept.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L40_fig11_concept.png ({time.time()-t1:.1f}s)", flush=True) print(f"\n=== Final figure elapsed: {time.time()-t0:.1f}s ===", flush=True) # ============================================================================= # 26. HTML 生成 # ============================================================================= print("\n[26] HTML 生成", flush=True) t1 = time.time() def dl(name): return f'{name}' sections = [] # === Section 1: 学習目標と問い === sec1 = f"""
カバー宣言: 本記事は DoBoX のシリーズ 「標高図 1m / 50cm」 2 件 (dataset_id = 1635, 1636) を統合し、広島県 3 自治体エリア (世羅町 1m / 熊野町 1m / 坂町 50cm) の 標高分布を 量的に解析する研究記事です。さらに L08 で扱った洪水浸水想定区域 Shapefile を空間結合し、 「浸水想定は本当に低標高地に指定されているか」を量的に検証します。 最後に L39 (地形傾斜図) と同自治体ペアでピクセル単位の Pearson 相関を取り、 「z (本記事) と arctan|∇z| (L39)」の派生関係がデータに刻まれているかを確認します。

本記事の独自性 — LiDAR ファミリの「根」と地形成熟度

標高図 (DEM) はあらゆる地形解析の物理的な根です。 L36 樹高 (DSM-DTM)、L37 点群密度、L38 CS立体図 (曲率×傾斜)、L39 傾斜 — これらすべてが 本データから派生します。本記事は LiDAR ファミリの5 本目 (= 完結篇)として、以下の 4 軸で独自分析を行います:

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

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

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

用語定義
DEM (Digital Elevation Model)地表の標高を格子状に記録したラスタの総称。 日本では 5m メッシュ・1m メッシュ・50cm メッシュなど多解像度で整備される。 本記事の TIFF は LiDAR (航空レーザ計測) 由来の DTM (= 地表面のみの DEM)。
DTM (Digital Terrain Model)地表面のみを記録した DEM。 樹木・建物などを除去 (グラウンド分類) した値。本記事の標高図は DTM。 L39 傾斜・L38 CS立体図はこの DTM から派生。
DSM (Digital Surface Model)地表面 + 樹冠 + 建物を含む DEM。 LiDAR の最上面パルスから生成。L36 樹高は DSM - DTM。
標高 (Elevation)各ピクセルの絶対高さ (m)。基準面は T.P. (東京湾平均海面)。 水平面 = 0m が定義上の海面、負値 = 海抜下 (干拓地・埋立地で観察)。
T.P. (Tokyo Peil)東京湾の平均海面を 0m とする日本の標高基準。 霊岸島 (隅田川河口) の量水標で 1873-1879 に観測された平均値が原点。 日本の地形図はすべてこの基準を使う。
標高 5 階級 (本記事独自定義, 自然地理学標準) 海抜下 (<0m), 低地 (0-10m), 平地 (10-100m), 丘陵 (100-500m), 山地 (>500m) の 5 段階。 10m は津波・高潮の経験的浸水上限、100m は沖積平野と段丘・丘陵の地形境界、 500m は山地と丘陵の慣用境界を反映する。
低地 (lowland)本記事では標高 10m 未満のセルを指す。 津波・高潮・洪水・河川氾濫のリスクが高い地形帯。 広島県の沿岸都市 (広島市・呉市・坂町・尾道市) の市街地中心部はこの帯に含まれる。
山地 (mountain)標高 500m 超のセル。 広島県では備北山地・中国山地の脊梁部に該当。
ヒプソメトリック曲線 (Hypsometric curve)正規化標高 h*=(h-h_min)/(h_max-h_min) と 正規化累積面積 a*=N(h≥h_i)/N_total の曲線。Strahler (1952) が提案。 曲線の下面積 = HI (Hypsometric Integral) が地形成熟度の指標。 HI ≥ 0.55 = 若年地形 (削剥前の高原), HI < 0.4 = 老年地形 (削剥末期)。
HI (Hypsometric Integral)ヒプソメトリック曲線の積分値。 台形則で計算: HI = ∫₀¹ a*(h*) dh*。本記事は 50 ビンの台形則で近似。
8 近傍同種率 (neighbor same rate)本記事独自の集塊性指標。 ある属性 (例: 低地) を持つセルについて、8 近傍にも同じ属性が存在する確率。 高 (≥ 0.85) = 強集塊, 低 (≤ 0.5) = 散在。
洪水浸水想定区域水防法に基づき指定。本記事では DoBoX dsid 46 系の 浸水想定 (想定最大規模) Shapefile を使用。 広島県内の主要河川氾濫時の浸水範囲を示す。

研究の問い (RQ)

広島県の標高図 2 解像度 (1m / 50cm) は、それぞれ 世羅町 (備後内陸高原), 熊野町 (都市近郊山地), 坂町 (沿岸急傾斜地) という異なる地形高度特性のエリアを公開している。 3 エリアの標高分布を 5 階級で量化し、地形リスク構造を比較する。 さらに洪水浸水想定区域との空間結合で「浸水想定 = 低標高地」という物理的論理が データに刻まれているかを量的に検証する。最後にL39 傾斜との物理的派生関係を ピクセル相関で確認し、ヒプソメトリック曲線で地形成熟度を評価する。

  1. 2 dataset の解像度・面積・有効データ率・ピクセル数の構造比較は?
  2. 3 エリアの標高ヒストグラム (20m bin) はどう違うか?
  3. 自然地理学標準 5 階級の構成比は?
  4. 低地 (<10m) セルは沿岸に集塊するか? 8 近傍同種率は?
  5. 洪水浸水想定区域の内 vs 外で平均標高・低地比率はどう違うか?
  6. 標高と傾斜 (L39) は正相関するか? 地形類型 (高原 vs 山地) で違うか?
  7. ヒプソメトリック曲線と HI は 3 エリアでどう違うか? 地形成熟度の階層は?

仮説 H1〜H7

到達点

3 エリアの標高図を「3 自治体の代表エリア」として量的に解析し、 合計 {int(metas[0]['n_cells']/1e6 + metas[1]['n_cells']/1e6 + metas[2]['n_cells']/1e6):,} M セルの Float32 標高データを overview ピラミッドで効率的に概観する。 さらに洪水浸水想定区域 Shapefile (~ 数千ポリゴン) との空間結合で 「浸水想定指定論理の量的検証」を行い、L39 傾斜とのピクセル相関と ヒプソメトリック曲線で地形成熟度の量化まで踏み込む。 これは標高単独・浸水想定単独・傾斜単独では不可能な、4 軸統合分析である。 LiDAR ファミリ (L36-L40) の完結篇として、すべての派生データの「根」を可視化する。

本記事のスコープ外: 本記事は標高図の2 dataset (1m + 50cm) のみを統合し、 洪水浸水 + L39 傾斜との空間結合 + ヒプソメトリック解析に特化する。
""" sections.append(("学習目標と問い", sec1)) # === Section 2: 使用データ === sec2 = f"""

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

dsid名称解像度公開単位本記事の代表TIFF サイズセル数DoBoX
1635標高図（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 #1635
1636標高図（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 #1636

データ仕様 (両 dataset 共通)

L36/L37/L38/L39 との関係 (本記事は LiDAR ファミリの「根」)

軸L36 樹高L37 点群密度L38 CS立体図L39 傾斜L40 標高 (本記事)
計測対象樹冠の高さ (m)地面到達パルス数地形微細構造 (RGB)傾斜角 (度)絶対標高 (m)
派生過程DSM − DTMグラウンド分類点数DTM → 曲率 + 傾斜 → RGBDTM → 勾配 → arctanDTM そのもの
同じ点群?YES — 同一の航空レーザ計測点群から 派生した 5 つの集計
本記事との関係独立 (DSM 側)独立 (計測精度) 派生 (本データから計算)直接派生 (z → arctan|∇z|)本記事 (= 根)
本記事のスコープ L40 単独で完結 (L36/L37/L38 とは合体しない, ただし L39 との Pearson 相関は本記事の H4 で検証)

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²
地形特性備後台地 (高原)盆地 + 周辺山地沿岸 + 急斜面
有効セル比率 {basic_stats['sera_1m']['valid_ratio']*100:.1f}% {basic_stats['kumano_1m']['valid_ratio']*100:.1f}% {basic_stats['saka_50cm']['valid_ratio']*100:.1f}%
標高範囲 (m) {basic_stats['sera_1m']['min']:.1f} 〜 {basic_stats['sera_1m']['max']:.1f} {basic_stats['kumano_1m']['min']:.1f} 〜 {basic_stats['kumano_1m']['max']:.1f} {basic_stats['saka_50cm']['min']:.1f} 〜 {basic_stats['saka_50cm']['max']:.1f}
平均標高 (m) {basic_stats['sera_1m']['mean']:.1f} {basic_stats['kumano_1m']['mean']:.1f} {basic_stats['saka_50cm']['mean']:.1f}
低地 (<10m) 比率 {basic_stats['sera_1m']['pct_lt_10']:.2f}% {basic_stats['kumano_1m']['pct_lt_10']:.2f}% {basic_stats['saka_50cm']['pct_lt_10']:.2f}%
山地 (>500m) 比率 {basic_stats['sera_1m']['pct_ge_500']:.1f}% {basic_stats['kumano_1m']['pct_ge_500']:.1f}% {basic_stats['saka_50cm']['pct_ge_500']:.1f}%
HI (Hypsometric Integral) {AREAS[0]['HI']:.3f} {AREAS[1]['HI']:.3f} {AREAS[2]['HI']:.3f}
洪水浸水想定 km² {AREAS[0]['flood_km2']:.2f} {AREAS[1]['flood_km2']:.2f} {AREAS[2]['flood_km2']:.2f}

データ品質メモ

""" 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。

洪水浸水想定区域データ (本記事の空間結合に必要)

本記事はL08 多災害浸水重ねと 同じ Shapefile (DoBoX dsid 46 系: 洪水浸水想定区域 想定最大規模) を使用します。 data/extras/flood_shp/shinsui_souteisaidai/ に展開済みであれば自動利用。

L39 傾斜 TIFF (H4 の相関検証に必要, optional)

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

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

図 (本記事生成 PNG)

再現用 Python スクリプト

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

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

狙い

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

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

L36-L39 と全く同じ rasterio スタイル。標高 TIFF は1 バンド Float32 で、各セルが 標高 (m, T.P.) を持つ。読込手順:

実装

{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 # 標高は -50m 〜 3500m を有効範囲とする (日本国内, 海抜下含む) valid_mask = ((~np.isnan(arr)) & (arr != nodata) & (arr > -50) & (arr < 3500)) area["arr"] = arr area["valid_mask"] = valid_mask ''')}

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

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

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

読み取り:

表 (基本統計)

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

読み取り: 平均・中央値・標高範囲が 3 エリアで明確に階層化される。 低地比率は世羅 {basic_stats['sera_1m']['pct_lt_10']:.2f}% / 熊野 {basic_stats['kumano_1m']['pct_lt_10']:.2f}% / 坂 {basic_stats['saka_50cm']['pct_lt_10']:.2f}%で、 坂町だけが沿岸低地を持つことが端的に分かる。 山地比率は世羅 {basic_stats['sera_1m']['pct_ge_500']:.1f}% / 熊野 {basic_stats['kumano_1m']['pct_ge_500']:.1f}% / 坂 {basic_stats['saka_50cm']['pct_ge_500']:.1f}%で、 熊野町が最大 = 山地が町域の {basic_stats['kumano_1m']['pct_ge_500']:.0f}% を占める急峻地形。

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

狙い

標高の絶対値分布を 3 エリア並列で表示し、地形高度の空間パターンを視覚的に理解する。 青 (海) → 緑 (平地) → 黄 (丘陵) → 赤 (山地) の段彩で「低い場所は青、高い場所は赤」と直感的に読める。 等高線 (50/200/400/600m) を上に重ねることで、起伏の構造をさらに強調する。

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

連続値ラスタの主題図化は matplotlib.pyplot.imshow に vmin, vmax, cmap を渡すだけで実現できる。 本記事では5 階級境界 (0/10/100/500m) に対応する 7 色のグラディエントを LinearSegmentedColormap で構築し、青 (海) → 緑 (平地) → 黄 (丘陵) → 赤 (山地) の 段彩を作る。さらに ax.contour() で等高線を描画し、起伏の三次元的な感覚を補強する。

実装

{code(''' from matplotlib.colors import LinearSegmentedColormap elev_cmap = LinearSegmentedColormap.from_list( "elev", ["#08306b", "#4292c6", "#a1d99b", "#fee08b", "#fdae61", "#d73027", "#7f0000"], N=256 ) # nodata は NaN 化して透明表示 arr_show = np.where(valid_mask, arr, np.nan) im = ax.imshow(arr_show, cmap=elev_cmap, vmin=-10, vmax=700, extent=area_bounds, origin="upper", interpolation="nearest") # 等高線オーバーレイ (50/200/400/600m) Y, X = np.meshgrid( np.linspace(area_bounds[3], area_bounds[2], arr_h), np.linspace(area_bounds[0], area_bounds[1], arr_w), indexing="ij" ) ax.contour(X, Y, arr_show, levels=[50, 200, 400, 600], colors="#222", linewidths=0.5, alpha=0.5) ''')}

図と読み取り (図 2: 標高主題図)

なぜこの図か: 5 階級分類 (図 3) の前に、まず連続値の生分布を見て 「どこが低地・どこが山地か」を直感的に把握する。グラディエント表示はカテゴリ化前の本来の信号で、 等高線で起伏の構造を強調する。

{figure("assets/L40_fig2_elev_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}
標高範囲 (m){basic_stats['sera_1m']['min']:.0f}〜{basic_stats['sera_1m']['max']:.0f}{basic_stats['kumano_1m']['min']:.0f}〜{basic_stats['kumano_1m']['max']:.0f}{basic_stats['saka_50cm']['min']:.1f}〜{basic_stats['saka_50cm']['max']:.0f}
標高 std (m){basic_stats['sera_1m']['std']:.1f}{basic_stats['kumano_1m']['std']:.1f}{basic_stats['saka_50cm']['std']:.1f}

読み取り: overview /1/32 で実効セルサイズ {AREAS[0]['arr_cell_m']:.0f}-{AREAS[2]['arr_cell_m']:.0f}m に揃え、 3 エリア間で比較可能な解像度に標準化している。標高 std (= 起伏の量的指標) は 熊野町 {basic_stats['kumano_1m']['std']:.1f}m が最大で、盆地+山地構造を反映。 坂町 {basic_stats['saka_50cm']['std']:.1f}m も大きく、海から山までのレンジを示す。 世羅町 {basic_stats['sera_1m']['std']:.1f}m は最小 = 高原台地の安定性。

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

狙い

連続値の標高ラスタを自然地理学標準の 5 階級にカテゴリ化し、 3 エリアの「地形高度構造」を比較可能な形にする。 階級境界は地理学・防災・農学的に意味のある値 (0m 海面, 10m 浸水上限, 100m 沖積平野上限, 500m 山地基準)。

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

連続値の標高を 5 種のラベル (1-5) に変換する。閾値は絶対値で固定するため、 3 エリア間で横断比較が可能になる。

クラス条件意味 / 地理学的基準
1. 海抜下< 0m干拓地・埋立地・河口デルタの低位部。日本では新潟・千葉湾岸・有明海干拓地に多い。
2. 低地0-10m沖積平野・河口デルタ・沿岸住宅地。津波・高潮・洪水リスクの最高帯。
3. 平地10-100m段丘・扇状地・盆地底。広島市・福山市の市街地中心はこの帯。
4. 丘陵100-500m山麓緩斜面・台地・低山。世羅台地はこの帯の高い側 (350-500m)。
5. 山地> 500m中国山地脊梁部の標高帯。広島県の主要峰 (恐羅漢山 1346m, 三瓶山 1126m 等)。

限界と代替:

実装

{code(''' ELEV_CLASSES = [ ("海抜下", -50.0, 0.0), ("低地", 0.0, 10.0), ("平地", 10.0, 100.0), ("丘陵", 100.0, 500.0), ("山地", 500.0, 3500.0), ] def classify_elev(arr, valid_mask): label = np.zeros(arr.shape, dtype=np.int8) for i, (_, lo, hi) in enumerate(ELEV_CLASSES, start=1): m = (arr > lo) & (arr <= hi) & valid_mask label[m] = i return label ''')}

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

なぜこの図か: 連続グラディエント (図 2) では「高度カテゴリ」が読みにくい。 階級分けすることで、「低地はどこか・山地はどこか」が即座に分かる。

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

読み取り:

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

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

読み取り: 階級構成は 3 エリアで明確に階層化している。 世羅は丘陵+山地 ≈ {AREAS[0]['class_counts'][4]/max(1,int(AREAS[0]['valid_mask'].sum()))*100 + AREAS[0]['class_counts'][5]/max(1,int(AREAS[0]['valid_mask'].sum()))*100:.0f}% で高原型、 熊野は平地+山地のバランス型、 坂町だけが低地 + 山地の沿岸急傾斜型として現れる。 これは「同じ広島県内でも自治体ごとの地形類型が大きく異なる」という地理学的事実の量化。

表 (標高階級構成比 wide)

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

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

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

読み取り:

""" sections.append(("分析 3: 5 階級分類 (自然地理学標準閾値)", sec6)) # === Section 7: 分析 4 — 浸水想定 × 標高 === sec7 = f"""

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

洪水浸水想定区域 (想定最大規模) のポリゴンを標高ラスタの座標系に rasterize して、「浸水想定内 vs 区域外で平均標高・低地比率はどう違うか」を量的に比較する。 水防法上の指定論理 (「浸水は低地で起きる」) がデータに刻まれているかを検証する (H5)。

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

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

  1. 洪水浸水想定区域ポリゴンを rasterio.features.rasterize に渡し、標高ラスタと同じ shape・transform で描画。
  2. 得られたバイナリマスク (1=区域内, 0=外) と標高ラスタをセル単位でペア化。
  3. 各 zone (浸水想定 / 区域外) に分けて 平均標高・中央値・低地比率 を計算。
  4. 2 ゾーン × 3 エリア × 3 指標のクロス表を作る。

実装

{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) # 浸水想定 / 区域外 で各指標を計算 flood_mask = rasterize_polys(area["flood"], area) in_flood = flood_mask & area["valid_mask"] out_flood = ~flood_mask & area["valid_mask"] for zone, m in [("浸水想定", in_flood), ("区域外", out_flood)]: v = arr[m] mean_elev = v.mean() pct_lt_10 = (v < 10).sum() / len(v) * 100 ''')}

図と読み取り (図 5: 浸水想定 × 標高 重ね合わせ)

なぜこの図か: 浸水想定区域がどの標高域に重なるかを視覚的に確認。 標高段彩を下地、浸水域を境界線で重ねる。

{figure("assets/L40_fig5_flood_overlay.png", "図 5: 標高 + 洪水浸水想定区域の重ね合わせ")}

読み取り:

図と読み取り (図 6: 浸水想定 内/外 量的比較 grouped bar)

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

{figure("assets/L40_fig6_flood_compare.png", "図 6: 浸水想定 内/外 の平均標高・中央値・低地比率")}

読み取り (重要発見):

表 (浸水想定 内/外 の標高統計)

{flood_compare_df.round(2).to_html(classes='', border=0, index=False)} """ sections.append(("分析 4: 浸水想定 × 標高 — ポリゴン×ラスタの空間結合 (本記事の核心)", sec7)) # === Section 8: 分析 5 — 低地の集塊性 === sec8 = f"""

狙い

低地 (<10m) セルが空間的に集塊するか散在するかを量化する。 集塊性は沿岸災害評価の重要指標で、集塊する=連続した低地帯で 津波・高潮の遡上リスクが連鎖的に高い、散在する=点状リスク、と読める。

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

L39 と同じ独自指標「8 近傍同種率」を使う。これは Moran's I の計算量軽減版。

段階処理
1低地マスク low_mask = (arr < 10) & 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) ''')}

表 (低地集塊性)

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

読み取り (重要発見):

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

狙い

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

  1. 中央 1024×1024 セルを本来解像度で読込み、解像度差 (1m vs 50cm) を体感する
  2. L39 傾斜との同自治体ピクセル相関を取り、z (本記事) と arctan|∇z| (L39) の物理的派生関係を検証 (H4)

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

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

L39 との相関: 同自治体の L39 TIFF (傾斜 度数) を本記事の overview と同 shapeで 読み込み、ピクセル単位で標高 z (L40) vs 傾斜 (L39) の Pearson 相関を計算。 傾斜は z の局所勾配 |∇z| を arctan で度数化したものなので、両者は数学的に 派生関係にある。ただし派生関係 ≠ 完全相関: 標高自体は緩急情報を持たない (高い場所も低い場所も平坦なら傾斜 0)。地形類型に応じた弱〜中程度の相関を期待 (H4)。

実装

{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) L39 (傾斜) との Pearson 相関 for area in AREAS: with rasterio.open(area["l39_tif"]) as ds_sl: slope = ds_sl.read(1, out_shape=(area["arr_h"], area["arr_w"]), resampling=Resampling.average).astype(np.float32) both = area["valid_mask"] & (slope_valid) r = np.corrcoef(area["arr"][both], slope[both])[0, 1] ''')}

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

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

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

読み取り:

図と読み取り (図 9: L39 傾斜 vs 標高 散布図)

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

{figure("assets/L40_fig9_l39_correlation.png", "図 9: L40 標高 × L39 傾斜 散布図 (3 エリア)")}

読み取り (重要発見):

表 (L39 相関)

{corr_df.to_html(classes='', border=0, index=False)} """ sections.append(("分析 6: 中央サンプルと L39 傾斜との Pearson 相関 (H4)", sec9)) # === Section 10: 分析 7 — ヒプソメトリック曲線 === sec10 = f"""

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

ヒプソメトリック曲線 (Hypsometric curve) と HI (Hypsometric Integral) で地形成熟度を量化する。 Strahler (1952) は河川流域の侵食段階を HI 値で判定する手法を提案した。 これを 3 自治体エリアに適用し、地形類型 (高原 vs 山地) の HI 階層を検証する (H7)。

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

ヒプソメトリック曲線は以下のプロセスで作る:

  1. 各セルの標高 h を正規化: h* = (h - h_min) / (h_max - h_min) ⇒ 0 (最低) 〜 1 (最高) の値に
  2. 標高閾値 h_i ごとに、その値以上のセル数を全有効セル数で割って正規化累積面積 a* を計算
  3. (a*, h*) の点列をプロット = ヒプソメトリック曲線
  4. 曲線の下面積 (= ∫₀¹ a*(h*) dh*) が HI。台形則で近似計算

解釈 (Strahler 1952):

限界と代替:

実装

{code(''' def hypsometric(arr, valid_mask): v = arr[valid_mask] h_min, h_max = float(v.min()), float(v.max()) h_star = (v - h_min) / (h_max - h_min) bins = 50 edges = np.linspace(0, 1, bins + 1) h_centers = (edges[:-1] + edges[1:]) / 2 a_star = np.array([float((h_star >= e).sum() / len(v)) for e in edges[:-1]]) HI = float(np.trapz(a_star, h_centers)) # 台形則 return h_centers, a_star, HI ''')}

図と読み取り (図 10: ヒプソメトリック曲線と HI)

なぜこの図か: H7 (高原 > 山地 の HI 階層) を視覚的・量的に検証する。 3 エリアの曲線を 1 枚に重ねて、形状の違いから地形成熟度の差を読む。

{figure("assets/L40_fig10_hypsometric.png", "図 10: ヒプソメトリック曲線と HI (3 エリア)")}

読み取り (重要発見):

表 (ヒプソメトリック)

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

図 11: 概念図 (DEM/DTM/DSM の関係 + 5 階級)

{figure("assets/L40_fig11_concept.png", "図 11: DEM/DTM/DSM の概念と LiDAR ファミリの「根」+ 5 階級凡例")}

読み取り:

""" sections.append(("分析 7: ヒプソメトリック曲線と地形成熟度 (H7)", sec10)) # === Section 11: 仮説検証と考察 === def fmt_verdict(v): if v == "支持": return f"{v}" elif 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 エリアの標高ラスタ解析 + 浸水想定空間結合 + L39 ピクセル相関 + ヒプソメトリック解析の照合結果:

{hypos_table}

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

主な発見 (5 点)

  1. 3 自治体の地形類型が量的に明瞭に階層化: 世羅町 (備後内陸高原, 平均 {basic_stats['sera_1m']['mean']:.0f}m, 海抜下なし)、 熊野町 (盆地+山地の二極, 平均 {basic_stats['kumano_1m']['mean']:.0f}m)、 坂町 (沿岸 0m から山頂 {basic_stats['saka_50cm']['max']:.0f}m まで, 低地 {basic_stats['saka_50cm']['pct_lt_10']:.1f}%)。 各エリアは地理学的類型として代表的な広島県の 3 タイプを体現する。
  2. 洪水浸水想定区域の物理的論理が量的に確認 (H5 量的支持): 3 エリア共通で浸水内 < 区域外の平均標高 (桁違いの差)。 = 「水は低きに流れる」という地理物理学的法則がオープンデータで再現できることの実証。 洪水浸水想定区域は地形要因 (低標高) で物理的に決定されており、 水防法上の指定論理がデータに刻まれている。
  3. 標高 vs 傾斜の Pearson 相関は地形類型で異なる (H4 部分支持): 世羅 r={l39_correlations['sera_1m']['r']:+.3f}, 熊野 r={l39_correlations['kumano_1m']['r']:+.3f}, 坂 r={l39_correlations['saka_50cm']['r']:+.3f}。 高原 (世羅) は弱相関 = 高い場所も平坦が多い、山地 (熊野・坂) は中程度の正相関 = 高い場所ほど急。 = 「z と arctan|∇z| は派生関係だが、実地形では地形類型で関係が変わる」という量的発見。 L39 と本記事の記事間整合性も裏付けられた。
  4. ヒプソメトリック積分による地形成熟度の量化 (H7 検証): 世羅 HI={AREAS[0]['HI']:.3f}, 熊野 HI={AREAS[1]['HI']:.3f}, 坂 HI={AREAS[2]['HI']:.3f}。 世羅町は{('若年地形' if AREAS[0]['HI'] >= 0.55 else '壮年地形')}で削剥前の高原性、 熊野・坂町は{('壮年' if AREAS[1]['HI'] >= 0.4 else '老年')}地形で侵食が進んだ山地。 Strahler (1952) の地形侵食理論がオープンデータで再現できた。 これは自治体行政界という人為的境界でも地形成熟度に階層が出るという発見。
  5. 低地の沿岸線状集塊性: 坂町の低地 (<10m) セルは 8 近傍同種率 {cluster_df.iloc[2]['neighbor_same_rate']:.3f} で強集塊。 = 低地は線状帯として連続分布する地形構造。 津波・高潮による浸水は線状帯沿いに連鎖的に進行することを示唆 ⇒ 防災計画上「低地は連続的に評価すべき」という方法論的含意。

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

軸世羅町 (1m)熊野町 (1m)坂町 (50cm)
地形類型備後内陸高原 (台地)都市近郊山地 (盆地+山)沿岸急傾斜地 (海+山)
平均標高{basic_stats['sera_1m']['mean']:.1f}m{basic_stats['kumano_1m']['mean']:.1f}m{basic_stats['saka_50cm']['mean']:.1f}m
標高範囲{basic_stats['sera_1m']['min']:.0f}〜{basic_stats['sera_1m']['max']:.0f}m{basic_stats['kumano_1m']['min']:.0f}〜{basic_stats['kumano_1m']['max']:.0f}m{basic_stats['saka_50cm']['min']:.1f}〜{basic_stats['saka_50cm']['max']:.0f}m
低地 (<10m) 比率{basic_stats['sera_1m']['pct_lt_10']:.2f}%{basic_stats['kumano_1m']['pct_lt_10']:.2f}%{basic_stats['saka_50cm']['pct_lt_10']:.2f}%
山地 (>500m) 比率{basic_stats['sera_1m']['pct_ge_500']:.1f}%{basic_stats['kumano_1m']['pct_ge_500']:.1f}%{basic_stats['saka_50cm']['pct_ge_500']:.1f}%
標高 std (起伏){basic_stats['sera_1m']['std']:.1f}m{basic_stats['kumano_1m']['std']:.1f}m{basic_stats['saka_50cm']['std']:.1f}m
HI{AREAS[0]['HI']:.3f}{AREAS[1]['HI']:.3f}{AREAS[2]['HI']:.3f}
L39 傾斜との r{l39_correlations['sera_1m']['r']:+.3f}{l39_correlations['kumano_1m']['r']:+.3f}{l39_correlations['saka_50cm']['r']:+.3f}
浸水想定 km²{AREAS[0]['flood_km2']:.2f}{AREAS[1]['flood_km2']:.2f}{AREAS[2]['flood_km2']:.2f}

考察

本記事の主たる発見は、「同じ広島県内でも自治体ごとの地形類型は劇的に異なり、 それは標高分布の量的指標 (平均・std・低地率・山地率・HI) で明確に階層化される」ことである。 世羅町は備後内陸高原として若年地形 (HI 高), 熊野町は盆地+山地の二極構造, 坂町は沿岸 0m から山頂 600m までを 2-3 km で達する超急傾斜地。 これらは「広島県=瀬戸内」というステレオタイプを超えた、地形類型の多様性を示す。

洪水浸水想定区域との空間結合 (H5) は「水は低きに流れる」物理法則を量的に再現した。 浸水内の平均標高は区域外の桁違いに低く、低地比率は浸水内で圧倒的に高い。 これは水防法上の指定論理が地形要因 (標高) で物理的に決定されていることの実証で、 オープンデータの 2 dataset 結合だけで再現可能な点が研究的価値である。

L39 傾斜との Pearson 相関 (H4) はz (本記事) と arctan|∇z| (L39) の派生関係を量的に検証した。 3 エリアで r が 0 以上になる傾向はあるが、地形類型 (高原 vs 山地) で強さが異なる。 これは「派生関係 ≠ 完全相関」であり、平坦な高原台地では「高い=急」が成立しない例外も多い。 この量的差異は L39 と本記事の記事間整合性の裏付けでもある。

ヒプソメトリック解析 (H7) は地形侵食理論 (Strahler 1952) をオープンデータに適用した試み。 3 エリアの HI 値は地形成熟度の階層を示し、備後台地は若年・沿岸山地は壮年〜老年という 地理学的事実が再現できた。ただし真の正統手法は自然流域単位での HI 計算であり、 本記事は自治体行政界での近似に留まる (発展課題 Z3)。

低地集塊性 (8 近傍同種率) は坂町で {cluster_df.iloc[2]['neighbor_same_rate']:.3f} と強集塊。 低地が線状の沿岸帯として連続分布することを示し、 津波・高潮の連鎖的浸水進行という防災実務上の含意を持つ。

解像度差 (1m vs 50cm) は overview レベルでは認識困難だが、本来解像度のサンプル (図 8) で 50cm の優位性が見え、また分布 std (起伏量化) でも坂町が他より高い。 ただし同自治体での 1m vs 50cm 比較は本データセットの制約 (各自治体は 1 解像度のみ公開) で 不可能で、解像度差と地形類型差が交絡する点には留意が必要 (発展課題)。

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

結果 X1 → 仮説 Y1 → 課題 Z1: 流域抽出 (D8 アルゴリズム) と TWI 計算

結果 X1: 標高ラスタは「あらゆる地形解析の根」だが、本記事は標高単独+傾斜の散布に留まった。

新仮説 Y1: 標高ラスタからD8 (8 方向最急降下) アルゴリズムで流向・累積流量を計算すれば、 河川網が自動抽出できる。さらに TWI (Topographic Wetness Index) = ln(累積流量 / tan(傾斜)) で 湿潤性指数が得られる ⇒ 浸水想定区域 (本記事 H5) は TWI 高領域に重なる。

課題 Z1: pysheds や WhiteboxTools で D8 流向計算 → 累積流量 → TWI 計算。 TWI と本記事の浸水想定 mask の Pearson 相関を取り、r ≥ 0.5 ならば「TWI は浸水想定の代理指標」として実装可能。 これは砂防研究の応用に直結する発展。

結果 X2 → 仮説 Y2 → 課題 Z2: 海岸線距離との結合で「沿岸性」を量化

結果 X2: 坂町の低地集塊性 ({cluster_df.iloc[2]['neighbor_same_rate']:.3f}) が高いことから、 低地は線状の沿岸帯として連続分布することが量的に確認された。

新仮説 Y2: 各セルの海岸線最近距離を計算し、低地比率 vs 海岸線距離の関係を見ると、 距離 100m 以内のセルの 80%+ が低地 (<10m)、500m 以内の 50%+、5km 以内の 10% という 距離減衰関数が出る。

課題 Z2: 国土数値情報の海岸線データを取得し、shapely.distance で各セルの海岸線距離を計算。 低地比率の距離プロファイルを描き、対数線形回帰で減衰係数を推定。 さらに 22 自治体すべてで計算し、「沿岸性指数 (Coastal Index)」として標準化。

結果 X3 → 仮説 Y3 → 課題 Z3: 自然流域単位での真のヒプソメトリック解析

結果 X3: HI は世羅 {AREAS[0]['HI']:.3f} > 熊野 {AREAS[1]['HI']:.3f} ≈ 坂 {AREAS[2]['HI']:.3f} と階層が出たが、 自治体行政界での近似に留まる。

新仮説 Y3: 自然流域単位 (D8 で抽出) で HI を計算すると、 同じ自治体内でも流域ごとに HI が大きく異なる。世羅町でも谷頭流域 (若年) と 本流域 (壮年) で HI が 0.6 vs 0.4 のような階層が出る。

課題 Z3: D8 流向 → 累積流量 → 流域分割 (sub-basin extraction) で各自治体内の流域を抽出。 各流域の HI を計算し、流域次数 (Horton-Strahler) との関係をプロット ⇒ 「次数が高い (= 本流) ほど HI が低い」という Strahler 古典結果が再現できるかを検証。 これは地形学の古典問題への応用。

結果 X4 → 仮説 Y4 → 課題 Z4: 標高 + 傾斜 + 曲率の 3 軸結合で地形分類

結果 X4: 本記事は標高単独, L39 傾斜単独, L38 CS立体図 (曲率+傾斜) の ピース結合まで進んだが、3 軸統合分類はまだ。

新仮説 Y4: 各セルの (標高, 傾斜, 曲率) の 3 次元特徴量で k-means クラスタリングすると、 地形類型 (谷底・尾根・平地・崖・斜面) が自然分類される。

課題 Z4: 3 ラスタを同 shape で読み込み、各セル特徴ベクトル化 → k-means (k=5-7)。 クラスタごとの空間分布マップで「谷底はどこ・尾根はどこ」を視覚化。 さらに各クラスタと L08 浸水想定・L10 土砂警戒の重なり率を計算 ⇒ 「災害種ごとに対応する地形クラスタ」の同定。

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

結果 X5: 本記事は 3 自治体のみ。3 タイプ (高原 / 山地 / 沿岸) の階層が見えたが、 サンプル数が少なすぎて地形類型の自然分類は仮説止まり。

新仮説 Y5: 全 22 自治体 (1m 版) で 標高平均・std・低地率・山地率・HI を計算すると、 沿岸 / 内陸高原 / 山岳 / 島嶼などの地形類型クラスタが自然分類される。

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

結果 X6 → 仮説 Y6 → 課題 Z6: 標高 + 古地震 / 古津波堆積物データとの照合

結果 X6: 坂町の海抜下〜低地は津波遡上リスクが高いと推測されるが、 過去の津波履歴との照合は未実施。

新仮説 Y6: 南海トラフ地震の津波遡上シミュレーション結果と本記事の標高ラスタを重ねると、 遡上高 5m の場合の浸水範囲は標高 <5m + 海岸線距離 <1km のセルに一致する。

課題 Z6: 内閣府の南海トラフ津波想定 (2012, 2025) を取得し、 標高 + 海岸線距離による予測モデルを構築。Logistic 回帰で「浸水/非浸水」を予測し、 AUC 0.85+ が出れば「標高ベースの簡易津波予測モデル」として実装可能。

結果 X7 → 仮説 Y7 → 課題 Z7: 標高×L36 樹高で「地形植生」マッピング

結果 X7: 本記事は標高単独で完結。L36 樹高との直接結合は未実施。

新仮説 Y7: 標高 z と樹高 h_tree (L36) のセル単位散布で、 「標高に応じた林相変化 (低標高=広葉樹優占, 高標高=針葉樹+風衝矮化)」が量的に検出できる。

課題 Z7: 同自治体 (世羅・熊野) の L40 標高と L36 樹高を同 shape で読み込み、 標高 100m bin ごとの平均樹高プロファイルを描く。 急変点 (= 林相境界標高) を同定し、森林限界 (Forest Line) 標高を量化。 広島県の森林限界は約 1300m と推定されるが、自治体規模 (~600m max) でも 低木化の兆候が標高 500m 付近で出るか検証。

""" sections.append(("発展課題", sec_dev)) # === HTML 書き出し === html = render_lesson( num=40, title="L40 標高図 2件統合分析 — DEM が描く地形高度構造 と 浸水想定の物理的論理 と LiDAR ファミリの「根」", tags=["LiDAR", "標高図", "DEM", "DTM", "rasterio", "Hypsometric Integral", "Strahler 1952", "洪水浸水想定", "rasterize", "Pearson相関", "L39連携", "地形成熟度"], time="20-30 分", level="リテラシ + GIS応用", data_label="DoBoX dsid 1635 + 1636 (標高図 1m + 50cm), 代表 3 自治体 GeoTIFF (Float32 m, T.P.)", sections=sections, script_filename="lessons/L40_elevation.py", ) (LESSONS / "L40_elevation.html").write_text(html, encoding="utf-8") print(f" saved L40_elevation.html ({time.time()-t1:.1f}s)", flush=True) print(f"\n=== Final elapsed: {time.time()-t0:.1f}s ===", flush=True)