# -*- coding: utf-8 -*-
"""L38 CS立体図 2 件統合分析
       — 広島県の航空レーザ計測 RGB ラスタが描く 地形微細構造 (尾根/谷/傾斜) と 土砂災害警戒の関係

カバー宣言:
  本記事は DoBoX のシリーズ <b>「CS立体図 1m / 50cm」 2 件</b>
  (dataset_id = 1628, 1629) を統合し、広島県内 3 自治体エリア
  (世羅町=中山間/林業地, 熊野町=都市近郊/呉市隣接山地, 坂町=平成30豪雨被災地) の
  <b>地形微細構造</b>を RGB ラスタの色合成として読み解く研究記事である。

  CS立体図とは長野県林業総合センター (戸田堅一郎ら) が開発した地形可視化手法で、
  <b>標高 (DTM) から計算した曲率 (Curvature) と傾斜 (Slope) を 3 バンドの RGB に合成</b>した
  地形微細構造可視化図。標準実装では:
    R (赤) = 尾根部 (凸地形, 正の曲率) を強調
    G (緑) = 傾斜 (Slope) を明度として表現
    B (青) = 谷部 (凹地形, 負の曲率) を強調
  この 3 軸の合成色から、専門家は<b>尾根線・谷線・崩壊跡・段丘・古砂防地形</b>を
  視覚的に検出する。

  L36 樹高図, L37 点群密度図 と並ぶ <b>LiDAR ファミリの 3 本目</b>。
  他の LiDAR 系記事 (地形傾斜図/標高図/航空レーザ計測森林資源/地図情報/計測年度図) とは
  <b>合体しない厳密シリーズ単独記事</b>。L36/L37 とは「同じ航空レーザ点群から派生した別の派生物」
  という関係 — DTM (標高) を経由して曲率・傾斜を計算した RGB 段彩。

研究の問い (RQ):
  広島県内の CS立体図 2 解像度 (1m / 50cm) は、それぞれ
  <b>世羅町 (中山間/林業地)</b>・<b>熊野町 (都市近郊/山岳混在)</b>・<b>坂町 (平成30豪雨被災地)</b>
  という<b>異なる地形リスク特性</b>のエリアを公開している。
  3 エリアの CS立体図を <b>R/G/B 3 軸の値分布</b>として読み解き、
  地形微細構造 (尾根/谷/傾斜) の構成比を比較する。さらに
  <b>土砂災害警戒区域 (L10/L11 既知) との空間結合</b>で、CS図のどの色域 (= どの地形型)
  が警戒域に多いかを定量化する。

    (1) 各エリアの<b>RGB 3 バンド値分布</b>の比較 (尾根強度/傾斜/谷強度)
    (2) Hue-Saturation-Value 変換による<b>地形タイプ自動分類</b> (尾根/谷/傾斜地/平地)
    (3) 高曲率セル (赤強・青強) の空間集塊性
    (4) 土砂災害警戒区域 (L10) ポリゴン × CS図の R/G/B 平均値比較
    (5) 解像度差 (1m vs 50cm) で見える<b>地形微細構造の違い</b>
    (6) 平成30豪雨被災地 (坂町) と非被災地 (世羅町) の CS図特徴比較

仮説 H1〜H6:
  H1 (世羅町は緩やかな丘陵): 中央値の (R, G, B) は中間色 (≈ 130) で
       極端な R高 / B高セルが少ない。Sat 平均 ≤ 0.30。
       中山間の植林地は緩傾斜地・尾根丘陵が多く、地すべり地形は限定的。
  H2 (熊野町は急峻地形): R高 (>180, 尾根) と B高 (>180, 谷) のセル比率が世羅より高い。
       特に Sat 平均 ≥ 0.35。熊野は呉市背後山地と平地の境界に位置し
       急峻な微地形が多い (花崗岩風化マサ土地形)。
  H3 (坂町は崩壊跡を含む急傾斜): G値 (傾斜) の平均が 3エリアで最高。
       平成30豪雨で多数発生した崩壊地形を CS図が刻む。
  H4 (土砂災害警戒区域は CS図で赤強・青強): 警戒区域内の CS図は
       平均 R, B が区域外より高い (急斜面・尾根谷の微地形が多い)。
       特に特別警戒区域 (レッド) では R+B 強度の合計が外より +10 以上。
  H5 (CS図の色相は地形タイプを区別する): HSV 変換した H (色相) のヒストグラムは
       2 峰性 (赤系=尾根 / 青系=谷) を示す。緑系 (傾斜のみ) は中間。
  H6 (50cm 版は微細崩壊地形を捉える): 50cm の Sera vs Saka で、
       同じ町で見ても 50cm の方が Saturation の標準偏差が大きく、
       ピクセルあたりの<b>「微地形コントラスト」が高い</b>。

要件 S 準拠: 1 分以内完走 (各エリア overview /1/64〜128 で MB 級メモリ)。
要件 T 準拠: 各エリア地図 + 地形タイプ分類マップ + small multiples + 警戒区域重ね。
要件 Q 準拠: 図 10+ 種, 表 8+ (3 エリア × 多角度: RGB分布/HSV/警戒区域/解像度比較)。

データ仕様:
  A. CS立体図 1m  dataset 1628: 県内 22 自治体 GeoTIFF 公開。
     本記事は<b>世羅町 (resource_id=177706)</b>と<b>熊野町 (176860)</b>を採用。
     - 形式: GeoTIFF RGB 3バンド uint8 (0-255)
     - 解像度: 1m, EPSG:6671, nodata=0
  B. CS立体図 50cm dataset 1629: 県内 12 自治体 GeoTIFF 公開。
     本記事は<b>坂町 (resource_id=176895)</b>を採用 (坂町は L36/L37 にはない新エリア
     = 平成30年7月豪雨被災地, 土石流災害多数発生)。
     - 形式: GeoTIFF RGB 3バンド uint8 (0-255)
     - 解像度: 0.5m, EPSG:6671, nodata=None (= データ全域有効)
"""
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, hsv_to_rgb, rgb_to_hsv
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("=== L38 CS立体図 2 件統合分析 ===", flush=True)

# =============================================================================
# 0. 定数・データパス・ZIP のダウンロード関数
# =============================================================================
DATA_DIR = ROOT / "data" / "extras" / "L38_cs_terrain"
DATA_DIR.mkdir(parents=True, exist_ok=True)
ADMIN_DIR = ROOT / "data" / "extras" / "L15_admin_zones"
SED_DIR = ROOT / "data" / "extras" / "sediment_shp"

# 各エリア定義 (3 エリア = 1m 2件 + 50cm 1件)
AREAS = [
    {
        "key": "sera_1m",
        "label": "世羅町",
        "muni_code": 34462,
        "resolution_m": 1.0,
        "dataset_id": 1628,
        "dataset_name": "CS立体図（1m）",
        "resource_id": 177706,
        "zip_url": "https://hiroshima-dobox.jp/resource_download/177706",
        "zip_url_alt": "https://data.hiroshima-dobox.jp/forest/cs3d_1m/34462_世羅町.zip",
        "zip_local": DATA_DIR / "1m_34462_sera_cs.zip",
        "tif": DATA_DIR / "1m_sera" / "sera_1m_cs.tif",
        "admin_zip": ADMIN_DIR / "admin_941_世羅町.zip",
        "color": "#1a7f37",       # 林業地=深緑 (L36/L37 と統一)
        "character": "中山間 / 林業地 (緩傾斜丘陵地形, 大規模崩壊跡なし想定)",
        "pref_city_code": "34462",
    },
    {
        "key": "kumano_1m",
        "label": "熊野町",
        "muni_code": 34307,
        "resolution_m": 1.0,
        "dataset_id": 1628,
        "dataset_name": "CS立体図（1m）",
        "resource_id": 176860,
        "zip_url": "https://hiroshima-dobox.jp/resource_download/176860",
        "zip_url_alt": None,
        "zip_local": DATA_DIR / "1m_34307_kumano_cs.zip",
        "tif": DATA_DIR / "1m_kumano" / "kumano_1m_cs.tif",
        "admin_zip": ADMIN_DIR / "admin_911_熊野町.zip",
        "color": "#0969da",       # 都市近郊=青 (L36/L37 と統一)
        "character": "都市近郊 / 山岳混在 (花崗岩マサ土地形, 急傾斜面多)",
        "pref_city_code": "34307",
    },
    {
        "key": "saka_50cm",
        "label": "坂町",
        "muni_code": 34384,
        "resolution_m": 0.5,
        "dataset_id": 1629,
        "dataset_name": "CS立体図（50cm）",
        "resource_id": 176895,
        "zip_url": "https://hiroshima-dobox.jp/resource_download/176895",
        "zip_url_alt": None,
        "zip_local": DATA_DIR / "50cm_34384_saka_cs.zip",
        "tif": DATA_DIR / "50cm_saka" / "saka_50cm_cs.tif",
        "admin_zip": ADMIN_DIR / "admin_916_坂町.zip",
        "color": "#cf222e",       # 被災地=赤
        "character": "平成30年7月豪雨被災地 (土石流多数発生, 急傾斜面に集塊住宅地)",
        "pref_city_code": "34384",
    },
]


def ensure_zip(area):
    """ZIP が無ければ DoBoX から DL (フォールバック URL あり)"""
    p = area["zip_local"]
    if p.exists() and p.stat().st_size > 100_000_000:
        return p
    import requests
    UA = {"User-Agent": "DoBoX-MDASH-textbook/1.0 (educational)"}
    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)
            continue
    raise RuntimeError(f"All download URLs failed for {area['label']}")


def ensure_tif(area):
    """ZIP から RGB TIFF と OVR (あれば) を抽出"""
    tif = area["tif"]
    if tif.exists() and tif.stat().st_size > 50_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 > 50_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


# 地形微細構造クラス (本記事独自定義)
# CS立体図の標準色設計:
#   R (赤) 高 = 凸地形 (尾根)
#   B (青) 高 = 凹地形 (谷)
#   G (緑) は明度で傾斜を表現
#   グレー中間色 = 緩傾斜の平地
TERRAIN_CLASSES = [
    ("nodata / 無効",       0,  None,        None,        "#ffffff"),
    ("尾根強 (R≫B)",        1,  "ridge",     "R-B>=20",   "#e74c3c"),
    ("谷強 (B≫R)",          2,  "valley",    "B-R>=20",   "#2980b9"),
    ("急傾斜 (G高+|R-B|低)", 3,  "slope",     "G>=180,|R-B|<20", "#f39c12"),
    ("平地 / 緩斜",          4,  "flat",      "Std<15",    "#95a5a6"),
    ("複雑地形 (Sat高 中間色)", 5, "complex", "Sat>=0.25,中間色",  "#8e44ad"),
]


# =============================================================================
# 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.nodatavals[0],
            "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 x {m['n_bands']} bands, "
          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 / "L38_tiff_meta.csv", index=False, encoding="utf-8-sig")
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 3. overview を使って RGB 3 バンドを読み込み (MB 級メモリ)
# =============================================================================
print("\n[3] overview で RGB 3 バンド読み込み", flush=True)
t1 = time.time()

# 目標 ≈ 32m/cell (3エリアの分布が比較可能なサイズ)
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)

        # RGB 3 バンドを 1 度に読込 (shape: (3, H, W))
        rgb = ds.read(out_shape=(3, ov_h, ov_w),
                      resampling=Resampling.average).astype(np.float32)
        a["rgb"] = rgb
        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

        # nodata マスク: 全バンド 0 = 無効
        nodata = a["meta"]["nodata"]
        if nodata is None:
            valid_mask = rgb.sum(axis=0) > 0
        else:
            valid_mask = ~((rgb[0] == nodata) & (rgb[1] == nodata) & (rgb[2] == nodata))
        a["valid_mask"] = valid_mask
        n_valid = int(valid_mask.sum())
        print(f"    valid={n_valid:,}/{ov_w*ov_h:,} ({n_valid/(ov_w*ov_h)*100:.1f}%)", flush=True)
        if n_valid > 0:
            for i, name in enumerate(["R", "G", "B"]):
                v = rgb[i][valid_mask]
                print(f"    {name}: min={v.min():.0f}, max={v.max():.0f}, mean={v.mean():.1f}, std={v.std():.1f}",
                      flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 4. RGB → HSV 変換 と 地形タイプ分類
# =============================================================================
print("\n[4] HSV 変換 + 地形タイプ分類", flush=True)
t1 = time.time()


def classify_terrain(rgb, valid_mask):
    """RGB ラスタ (3, H, W) を 5 種の地形タイプにラベル付け
    返値: ラベル ndarray (H, W) int8, 0=nodata, 1=尾根, 2=谷, 3=傾斜地, 4=平地, 5=複雑

    分類ルール (本記事独自定義, パーセンタイルベース適応閾値):
      CS立体図の色は実データでは中間色 (パステル系) が主で、極端な R, B 高セルは少ない。
      よって絶対値ではなく <b>各エリアのパーセンタイル分布</b>で適応的に閾値を決める。

      ridge_index = R - B  (>=0 で尾根寄り, <0 で谷寄り)
      brightness  = (R+G+B)/3  (低い = 暗い = 急斜面の影)
      band_std    = R/G/B の標準偏差  (低い = 中間色 = 平地)

      0 = nodata / 無効
      1 = 尾根強  (ridge_index >= P75, 上位 25%)
      2 = 谷強   (ridge_index <= P25, 下位 25%)
      3 = 急傾斜 (brightness <= P15, 影として暗い 15%)
      4 = 平地  (band_std < P25, ほぼグレー 下位 25%)
      5 = 複雑地形 (上記いずれでもない = 中間色だが彩度はある)
    """
    R, G, B = rgb[0], rgb[1], rgb[2]
    diff_RB = R - B
    bright = (R + G + B) / 3.0
    band_std = np.stack([R, G, B], axis=0).std(axis=0)

    maxv = np.maximum(np.maximum(R, G), B)
    minv = np.minimum(np.minimum(R, G), B)
    sat = np.where(maxv > 0, (maxv - minv) / np.maximum(maxv, 1), 0.0)

    label = np.zeros(R.shape, dtype=np.int8)

    if valid_mask.sum() < 10:
        return label, sat

    diff_RB_v = diff_RB[valid_mask]
    bright_v = bright[valid_mask]
    band_std_v = band_std[valid_mask]
    p75 = float(np.percentile(diff_RB_v, 75))
    p25 = float(np.percentile(diff_RB_v, 25))
    bp15 = float(np.percentile(bright_v, 15))
    sp25 = float(np.percentile(band_std_v, 25))

    is_ridge  = (diff_RB >= p75) & valid_mask
    is_valley = (diff_RB <= p25) & valid_mask
    is_steep  = (bright <= bp15) & (~is_ridge) & (~is_valley) & valid_mask
    is_flat   = (band_std <= sp25) & (~is_ridge) & (~is_valley) & (~is_steep) & valid_mask
    is_complex = valid_mask & (~is_ridge) & (~is_valley) & (~is_steep) & (~is_flat)

    label[is_ridge] = 1
    label[is_valley] = 2
    label[is_steep] = 3
    label[is_flat] = 4
    label[is_complex] = 5
    return label, sat


for a in AREAS:
    label, sat = classify_terrain(a["rgb"], a["valid_mask"])
    a["terrain_label"] = label
    a["sat"] = sat
    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}%, 平地={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. RGB 統計 + Saturation
# =============================================================================
print("\n[6] RGB バンド統計", flush=True)
t1 = time.time()

basic_stats = {}
for a in AREAS:
    vm = a["valid_mask"]
    rgb = a["rgb"]
    sat = a["sat"]
    s = {
        "label": a["label"],
        "n_cells_total": int(vm.size),
        "n_cells_valid": int(vm.sum()),
        "valid_ratio": float(vm.sum() / vm.size),
        "R_mean": float(rgb[0][vm].mean()),
        "R_median": float(np.median(rgb[0][vm])),
        "R_std": float(rgb[0][vm].std()),
        "G_mean": float(rgb[1][vm].mean()),
        "G_median": float(np.median(rgb[1][vm])),
        "G_std": float(rgb[1][vm].std()),
        "B_mean": float(rgb[2][vm].mean()),
        "B_median": float(np.median(rgb[2][vm])),
        "B_std": float(rgb[2][vm].std()),
        "Sat_mean": float(sat[vm].mean()),
        "Sat_std":  float(sat[vm].std()),
        "Sat_p90":  float(np.percentile(sat[vm], 90)),
        "diff_RB_mean": float((rgb[0][vm] - rgb[2][vm]).mean()),
        "diff_RB_std":  float((rgb[0][vm] - rgb[2][vm]).std()),
    }
    basic_stats[a["key"]] = s

stats_df = pd.DataFrame.from_dict(basic_stats, orient="index")
stats_df.to_csv(ASSETS / "L38_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) のロード — 3 自治体それぞれの範囲で
# =============================================================================
print("\n[7] sediment 警戒区域ロード", flush=True)
t1 = time.time()


def load_sed_for_area(area):
    """対象自治体の bbox に絞って sediment を全県 shp から limit ロード"""
    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. 警戒区域マスクをラスタ化して RGB 比較
# =============================================================================
print("\n[8] 警戒区域マスク × RGB 比較", flush=True)
t1 = time.time()


def rasterize_polys(polys_gdf, area):
    """警戒区域 polygons をエリア overview の (H, W) に投影してマスク化"""
    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"]

    rgb = a["rgb"]
    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,
                   "R_mean": None, "G_mean": None, "B_mean": None,
                   "diff_RB_mean": None, "Sat_mean": None}
        else:
            row = {
                "area": a["label"], "zone": label_name, "n": n,
                "R_mean": float(rgb[0][m].mean()),
                "G_mean": float(rgb[1][m].mean()),
                "B_mean": float(rgb[2][m].mean()),
                "diff_RB_mean": float((rgb[0][m] - rgb[2][m]).mean()),
                "Sat_mean": float(a["sat"][m].mean()),
            }
        sed_compare_rows.append(row)

sed_compare_df = pd.DataFrame(sed_compare_rows)
sed_compare_df.to_csv(ASSETS / "L38_sediment_rgb_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. 生解像度サンプル (中央 1024x1024)
# =============================================================================
print("\n[9] 生解像度中央サンプル", 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(window=win).astype(np.uint8)  # (3, H, W)
        a["raw_rgb"] = raw
        a["raw_bounds"] = rio_bounds(win, ds.transform)
    nodata = a["meta"]["nodata"]
    if nodata is None:
        raw_valid = raw.sum(axis=0) > 0
    else:
        raw_valid = ~((raw[0] == nodata) & (raw[1] == nodata) & (raw[2] == nodata))
    a["raw_valid"] = raw_valid
    print(f"  {a['label']}: raw {raw.shape}, valid={int(raw_valid.sum()):,}", flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 10. 仮説検証
# =============================================================================
print("\n[10] 仮説検証", 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: 世羅は緩やか丘陵 → Sat_mean <= 0.10 (パステル色多)
h1 = s_sera["Sat_mean"] <= 0.10
hypos.append({
    "H": "H1",
    "claim": "世羅町は緩傾斜丘陵が主で、CS図の Saturation 平均 <= 0.10 (淡い中間色多)",
    "result": f"世羅 Sat_mean={s_sera['Sat_mean']:.3f} (Sat_std={s_sera['Sat_std']:.3f})",
    "verdict": "支持" if h1 else ("部分支持" if s_sera["Sat_mean"] <= 0.13 else "反証"),
})

# H2: 熊野は急峻地形 → Sat_mean > 世羅 (彩度上)
h2 = s_kuma["Sat_mean"] > s_sera["Sat_mean"]
hypos.append({
    "H": "H2",
    "claim": "熊野町は急峻地形で Sat 平均が世羅より大きい",
    "result": f"熊野 Sat_mean={s_kuma['Sat_mean']:.3f}, 世羅 Sat_mean={s_sera['Sat_mean']:.3f}, "
              f"差 {s_kuma['Sat_mean']-s_sera['Sat_mean']:+.3f}",
    "verdict": "支持" if h2 else "反証",
})

# H3: 坂町は崩壊跡を含む急傾斜 → 3エリアで Sat_mean 最大
h3 = s_saka["Sat_mean"] >= max(s_sera["Sat_mean"], s_kuma["Sat_mean"])
hypos.append({
    "H": "H3",
    "claim": "坂町は3エリアで Sat 平均が最大 (= 急傾斜が広く分布)",
    "result": f"世羅 {s_sera['Sat_mean']:.3f}, 熊野 {s_kuma['Sat_mean']:.3f}, 坂 {s_saka['Sat_mean']:.3f}",
    "verdict": "支持" if h3 else ("部分支持" if s_saka["Sat_mean"] >= s_sera["Sat_mean"] else "反証"),
})

# H4: 警戒区域内の CS図は赤強・青強 (= |R-B| 大)
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["Sat_mean"] is None or outer["Sat_mean"] is None:
        h4_pieces.append(f"{a['label']}: NA")
        continue
    sat_in = inner["Sat_mean"]
    sat_out = outer["Sat_mean"]
    drb_in = inner["diff_RB_mean"]
    drb_out = outer["diff_RB_mean"]
    diff_sat = sat_in - sat_out
    h4_supports.append(diff_sat > 0)
    h4_pieces.append(f"{a['label']}: Sat 内{sat_in:.3f} vs 外{sat_out:.3f} (差{diff_sat:+.3f}), "
                     f"diff_RB 内{drb_in:+.1f} vs 外{drb_out:+.1f}")

h4 = (sum(h4_supports) >= 2) if h4_supports else False
hypos.append({
    "H": "H4",
    "claim": "土砂災害特別警戒区域内の CS図 Saturation は区域外より大きい (急斜面の微地形)",
    "result": " | ".join(h4_pieces),
    "verdict": "支持" if h4 else ("部分支持" if (h4_supports and sum(h4_supports) >= 1) else "反証"),
})

# H5: 色相 (Hue) は 2峰性 (赤系 / 青系)
def hue_modes(rgb, vm):
    if vm.sum() < 100:
        return 0, 0, 0
    R, G, B = rgb[0][vm]/255.0, rgb[1][vm]/255.0, rgb[2][vm]/255.0
    rgb_v = np.stack([R, G, B], axis=-1)
    hsv = matplotlib.colors.rgb_to_hsv(rgb_v.reshape(1, -1, 3)).reshape(-1, 3)
    h = hsv[:, 0]
    # ビンに分けて各帯のピーク位置
    bins = np.linspace(0, 1, 13)
    hist, _ = np.histogram(h, bins=bins)
    return int(hist.argmax()), int(hist.shape[0]), float(h.std())

h5_pieces = []
all_h5_supports = []
for a in AREAS:
    mode_bin, n_bins, h_std = hue_modes(a["rgb"], a["valid_mask"])
    a["hue_mode_bin"] = mode_bin
    a["hue_std"] = h_std
    h5_pieces.append(f"{a['label']}: mode_bin={mode_bin}/{n_bins-1}, h_std={h_std:.3f}")
    all_h5_supports.append(h_std > 0.05)

h5 = sum(all_h5_supports) >= 2
hypos.append({
    "H": "H5",
    "claim": "色相 (Hue) ヒストの広がり (h_std) は両エリアとも > 0.05 (赤系-青系の双方が存在)",
    "result": " | ".join(h5_pieces),
    "verdict": "支持" if h5 else ("部分支持" if sum(all_h5_supports) >= 1 else "反証"),
})

# H6: 50cm (坂町) は 1m (世羅・熊野) よりピクセルあたりの Sat std が大きい
h6 = s_saka["Sat_std"] >= max(s_sera["Sat_std"], s_kuma["Sat_std"])
hypos.append({
    "H": "H6",
    "claim": "50cm 版 (坂町) は 1m 版より Sat 標準偏差が大きい (微地形コントラスト高)",
    "result": f"世羅 Sat_std={s_sera['Sat_std']:.3f}, 熊野 {s_kuma['Sat_std']:.3f}, 坂 {s_saka['Sat_std']:.3f}",
    "verdict": "支持" if h6 else ("部分支持" if s_saka["Sat_std"] >= s_sera["Sat_std"] else "反証"),
})

hypos_df = pd.DataFrame(hypos)
hypos_df.to_csv(ASSETS / "L38_hypothesis_results.csv", index=False, encoding="utf-8-sig")
with (ASSETS / "L38_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)


# =============================================================================
# 11. dataset/series 一覧 CSV
# =============================================================================
print("\n[11] 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 / "L38_series.csv", index=False, encoding="utf-8-sig")
print(series_df.to_string(index=False))

# 地形タイプ別構成比 wide
class_rows = []
for a in AREAS:
    cc = a["class_counts"]
    n = cc.sum()
    row = {"area": a["label"]}
    for i, lab in enumerate(["nodata","尾根","谷","急傾斜","平地","複雑"]):
        row[f"{lab}_n"] = int(cc[i])
        row[f"{lab}_pct"] = round(cc[i]/max(1,n) * 100, 2)
    class_rows.append(row)
class_df = pd.DataFrame(class_rows)
class_df.to_csv(ASSETS / "L38_terrain_classes.csv", index=False, encoding="utf-8-sig")
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 12. 図 1: dataset カード + Saturation ヒスト (3 エリア)
# =============================================================================
print("\n[12] 図 1: dataset カード + Sat ヒスト", 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 CS立体図 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")
    # 長い character は 2 行に折りたたむ
    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 / 3 バンド uint8", 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 x 3", 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"R {s['R_mean']:.0f}, G {s['G_mean']:.0f}, B {s['B_mean']:.0f}",
            ha="center", fontsize=9, fontweight="bold")
    ax.text(ax_x+1.75, 1.7, f"Sat 平均 {s['Sat_mean']:.3f} ± {s['Sat_std']:.3f}",
            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:
    sat = a["sat"][a["valid_mask"]]
    sat = sat[(sat >= 0) & (sat <= 0.5)]
    ax.hist(sat, bins=np.linspace(0, 0.5, 51), color=a["color"], alpha=0.55,
            density=True, edgecolor="white", linewidth=0.3,
            label=f"{a['label']} (μ={a['sat'][a['valid_mask']].mean():.3f})")
ax.set_xlabel("Saturation (HSV)")
ax.set_ylabel("密度 (normalized)")
ax.set_title("CS立体図 Saturation 分布 — 3 エリア比較")
ax.legend(loc="upper right", fontsize=10)
ax.grid(axis="y", alpha=0.3)

plt.tight_layout()
fig.savefig(ASSETS / "L38_fig1_dataset_overview.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L38_fig1_dataset_overview.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 13. 図 2: CS立体図 RGB 主題図 (3 エリア)
# =============================================================================
print("\n[13] 図 2: CS立体図 主題図", flush=True)
t1 = time.time()

fig, axes = plt.subplots(1, 3, figsize=(18, 7))

for ax, a in zip(axes, AREAS):
    rgb = a["rgb"]
    vm = a["valid_mask"]
    # 3 バンドを HxWx3 に並べ替えて表示
    rgb_disp = np.stack([rgb[0], rgb[1], rgb[2]], axis=-1).astype(np.uint8)
    # invalid を白に
    rgb_disp_show = np.where(vm[..., None], rgb_disp, 255)
    extent = (a["bounds"].left, a["bounds"].right,
              a["bounds"].bottom, a["bounds"].top)
    ax.imshow(rgb_disp_show, 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)

    ax.set_title(f"{a['label']} CS立体図 (overview /1/{a['arr_factor']}, "
                 f"実効 {a['arr_cell_m']:.0f} m/cell)\n"
                 f"{a['character']}", fontsize=10)
    ax.set_xlabel("X (EPSG:6671 m)")
    ax.set_ylabel("Y (EPSG:6671 m)")

fig.suptitle("3 エリア CS立体図 RGB 元画像 — 赤=尾根 (凸), 青=谷 (凹), 緑=傾斜による明度", fontsize=12, y=1.0)
plt.tight_layout()
fig.savefig(ASSETS / "L38_fig2_raster_maps.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L38_fig2_raster_maps.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 14. 図 3: 地形タイプ分類マップ
# =============================================================================
print("\n[14] 図 3: 地形タイプ分類", flush=True)
t1 = time.time()

terrain_colors = ["#ffffff", "#e74c3c", "#2980b9", "#34495e", "#bdc3c7", "#8e44ad"]
terrain_labels = ["nodata", "尾根 (R>=B+P75)", "谷 (B>=R+P25)", "急傾斜 (暗P15)",
                  "平地 (StdP25)", "複雑 (中間色)"]
cmap_t = ListedColormap(terrain_colors)
norm_t = 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["terrain_label"], cmap=cmap_t, norm=norm_t,
              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(terrain_colors[1:], terrain_labels[1:])]
fig.legend(handles=legend_handles, loc="lower center", ncol=5, fontsize=10,
           bbox_to_anchor=(0.5, -0.02))

fig.suptitle("地形タイプ自動分類 — RGB 適応閾値による 5 クラス分け", fontsize=12, y=1.0)
plt.tight_layout()
fig.savefig(ASSETS / "L38_fig3_terrain_classification.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L38_fig3_terrain_classification.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 15. 図 4: 各バンド (R, G, B) ヒストグラム比較
# =============================================================================
print("\n[15] 図 4: RGB バンドヒスト", flush=True)
t1 = time.time()

fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))

band_titles = [("R (尾根強度)", "#cf222e"), ("G (傾斜・明度)", "#2da44e"), ("B (谷強度)", "#0969da")]
for bidx, (btitle, bcolor) in enumerate(band_titles):
    ax = axes[bidx]
    for a in AREAS:
        v = a["rgb"][bidx][a["valid_mask"]]
        ax.hist(v, bins=np.arange(0, 256, 4), color=a["color"], alpha=0.55,
                density=True, edgecolor="white", linewidth=0.2,
                label=f"{a['label']} (μ={v.mean():.0f})")
    ax.set_xlabel(f"{btitle} 値 (0-255)")
    ax.set_ylabel("密度")
    ax.set_title(f"バンド {btitle}")
    ax.legend(loc="upper left", fontsize=8.5)
    ax.grid(axis="y", alpha=0.3)
    ax.set_xlim(0, 256)

plt.tight_layout()
fig.savefig(ASSETS / "L38_fig4_rgb_histograms.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L38_fig4_rgb_histograms.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 16. 図 5: 地形タイプ構成比 stack-bar
# =============================================================================
print("\n[16] 図 5: 地形タイプ構成比", flush=True)
t1 = time.time()

fig, ax = plt.subplots(figsize=(10, 5))

x = np.arange(len(AREAS))
bottom = np.zeros(len(AREAS))
class_show = ["尾根", "谷", "急傾斜", "平地", "複雑"]
class_colors = ["#e74c3c", "#2980b9", "#34495e", "#bdc3c7", "#8e44ad"]
for i, (cn, col) in enumerate(zip(class_show, class_colors)):
    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=cn, color=col, edgecolor="#222", linewidth=0.6)
    for j, v in enumerate(vals):
        if v >= 3:
            ax.text(j, bottom[j] + v/2, f"{v:.1f}%", ha="center", va="center",
                    fontsize=10, fontweight="bold",
                    color="white" if col in ("#e74c3c","#2980b9","#34495e","#8e44ad") else "black")
    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("地形タイプ自動分類 構成比 — 3 エリア (適応閾値)")
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 / "L38_fig5_terrain_composition.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L38_fig5_terrain_composition.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 17. 図 6: 警戒区域 × CS立体図 重ね合わせマップ (3 エリア)
# =============================================================================
print("\n[17] 図 6: 警戒区域 × CS overlay", flush=True)
t1 = time.time()

fig, axes = plt.subplots(1, 3, figsize=(18, 7))

for ax, a in zip(axes, AREAS):
    rgb = a["rgb"]
    vm = a["valid_mask"]
    rgb_disp = np.stack([rgb[0], rgb[1], rgb[2]], axis=-1).astype(np.uint8)
    rgb_disp_show = np.where(vm[..., None], rgb_disp, 255)
    extent = (a["bounds"].left, a["bounds"].right,
              a["bounds"].bottom, a["bounds"].top)
    ax.imshow(rgb_disp_show, 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"].plot(ax=ax, facecolor="#ffd700", edgecolor="#a87b00",
                              linewidth=0.5, alpha=0.42)
    if a.get("sed_red") is not None and len(a["sed_red"]) > 0:
        a["sed_red"].plot(ax=ax, facecolor="#cf222e", edgecolor="#600",
                           linewidth=0.5, alpha=0.55)

    ax.set_title(f"{a['label']}: CS立体図 + 土砂災害警戒区域\n"
                 f"警戒 {a['sed_yellow_n']} polys ({a['sed_yellow_km2']:.1f} km²) / "
                 f"特別 {a['sed_red_n']} polys ({a['sed_red_km2']:.2f} km²)",
                 fontsize=9.5)
    ax.set_xlabel("X (m)"); ax.set_ylabel("Y (m)")

legend_handles = [
    Patch(facecolor="#ffd700", alpha=0.5, edgecolor="#a87b00", label="土砂警戒 (yellow)"),
    Patch(facecolor="#cf222e", alpha=0.6, edgecolor="#600", label="特別警戒 (red)"),
]
fig.legend(handles=legend_handles, loc="lower center", ncol=2, fontsize=11,
           bbox_to_anchor=(0.5, -0.02))
fig.suptitle("CS立体図 × 土砂災害警戒区域 重ね合わせ — 3 エリア比較", fontsize=12, y=1.0)
plt.tight_layout()
fig.savefig(ASSETS / "L38_fig6_sediment_overlay.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L38_fig6_sediment_overlay.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 18. 図 7: 警戒区域内 vs 区域外 RGB 平均比較 (grouped bar)
# =============================================================================
print("\n[18] 図 7: 警戒区域 RGB 比較", 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, ymin in zip(
    axes,
    ["R_mean", "diff_RB_mean", "Sat_mean"],
    ["R 平均値 (尾根強度)", "R-B 平均 (尾根-谷バランス)", "Saturation 平均 (彩度)"],
    [None, None, None]
):
    n_areas = len(AREAS)
    x = 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(x + (i-1)*w, vals, w, color=zone_colors[zone],
               edgecolor="#222", linewidth=0.6, label=zone)
        for j, v in enumerate(vals):
            ax.text(x[j] + (i-1)*w, v, f"{v:.2f}" if metric == "Sat_mean" else f"{v:.1f}",
                    ha="center", va="bottom", fontsize=8)
    ax.set_xticks(x)
    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 right")
    ax.grid(axis="y", alpha=0.3)

fig.suptitle("土砂災害警戒区域 内 vs 外 — CS立体図 RGB 平均比較 (H4 検証)", fontsize=12, y=1.02)
plt.tight_layout()
fig.savefig(ASSETS / "L38_fig7_sediment_rgb_compare.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L38_fig7_sediment_rgb_compare.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 19. 図 8: 中央サンプル (本来解像度の RGB 表示)
# =============================================================================
print("\n[19] 図 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_rgb"]
    rgb_disp = np.stack([raw[0], raw[1], raw[2]], axis=-1).astype(np.uint8)
    ax.imshow(rgb_disp, interpolation="nearest")
    side_m = raw.shape[2] * 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"R平均 {raw[0].mean():.0f}, G平均 {raw[1].mean():.0f}, B平均 {raw[2].mean():.0f}",
                 fontsize=9.5)
    ax.set_xlabel("ピクセル X"); ax.set_ylabel("ピクセル Y")

fig.suptitle("中央 1024×1024 セルの本来解像度サンプル — 微地形構造の体感", fontsize=12, y=1.01)
plt.tight_layout()
fig.savefig(ASSETS / "L38_fig8_center_samples.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L38_fig8_center_samples.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 20. 図 9: CS立体図概念図
# =============================================================================
print("\n[20] 図 9: CS立体図概念図", 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(-0.5, 6); ax.set_axis_off()
ax.text(6, 5.7, "CS立体図のしくみ — 標高 (DTM) → 曲率 + 傾斜 → RGB 合成",
        ha="center", fontsize=12, fontweight="bold")

# 地形断面 (尾根 + 谷 + 平地)
xs = np.linspace(0.5, 11.5, 100)
ys = 1.5 + 1.4*np.exp(-((xs-3)/0.7)**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)

# 尾根に赤丸
ax.scatter([3], [ys[np.argmin(np.abs(xs-3))]], s=400, c="#cf222e", alpha=0.6,
           edgecolor="#222", linewidth=2, zorder=5)
ax.text(3, ys[np.argmin(np.abs(xs-3))]+0.5, "尾根 (凸)\nR 強", ha="center",
        fontsize=10, fontweight="bold", color="#cf222e")
# 谷に青丸
y7 = ys[np.argmin(np.abs(xs-7))]
ax.scatter([7], [y7], s=400, c="#2980b9", alpha=0.6,
           edgecolor="#222", linewidth=2, zorder=5)
ax.text(7, y7-0.6, "谷 (凹)\nB 強", ha="center",
        fontsize=10, fontweight="bold", color="#2980b9")
# 緩斜面に灰丸
ax.scatter([10], [ys[np.argmin(np.abs(xs-10))]], s=400, c="#888", alpha=0.6,
           edgecolor="#222", linewidth=2, zorder=5)
ax.text(10, ys[np.argmin(np.abs(xs-10))]+0.5, "平地\n中間色", ha="center",
        fontsize=10, fontweight="bold", color="#444")

# 急斜面の影
ax.fill_between([4.6, 5.4], 0, 1.0, color="#34495e", alpha=0.4)
ax.text(5.0, 0.5, "影=暗", ha="center", fontsize=9, color="#fff", fontweight="bold")
ax.text(5.0, -0.2, "急斜面 (G 低)", ha="center", fontsize=9, fontweight="bold")

ax.text(0.5, -0.5, "DTM 地表 →", fontsize=10)
ax.text(11.5, -0.5, "→", fontsize=10, ha="right")

# 右: 配色凡例
ax = axes[1]
ax.set_xlim(0, 10); ax.set_ylim(0, 10); ax.set_axis_off()
ax.text(5, 9.5, "CS立体図の配色設計 (本記事の解釈)", ha="center", fontsize=12, fontweight="bold")

legend_data = [
    ("R 高", "#cf222e", "凸地形 (尾根線・段丘崖の上端) を赤で強調"),
    ("B 高", "#2980b9", "凹地形 (谷線・崩壊地の凹み) を青で強調"),
    ("G 低 (暗)", "#34495e", "傾斜が急な部分を暗い色 (影) で表現"),
    ("3 バンド近似", "#bdc3c7", "緩傾斜の平地・台地は中間色 (グレー)"),
    ("Sat 高", "#8e44ad", "微地形コントラスト強 (= 古崩壊・段丘等)"),
]
y0 = 8.0
for lbl, col, 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, lbl, fontsize=11, fontweight="bold", color=col, va="center")
    ax.text(3.0, y0, desc, fontsize=10, va="center")
    y0 -= 1.3

ax.text(5, 0.5, "[本記事の独自解釈] CS立体図は専門家向け視覚化のため、\n"
        "RGB バンド値の絶対値ではなく<相対パターン>で読む",
        ha="center", fontsize=9, color="#666",
        bbox=dict(facecolor="#fff8c5", edgecolor="#d4a72c", boxstyle="round,pad=0.5"))

plt.tight_layout()
fig.savefig(ASSETS / "L38_fig9_concept.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L38_fig9_concept.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 21. 図 10: HSV scatter (Hue-Saturation 平面)
# =============================================================================
print("\n[21] 図 10: HSV 散布", flush=True)
t1 = time.time()

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

for ax, a in zip(axes, AREAS):
    vm = a["valid_mask"]
    rgb = a["rgb"]
    R, G, B = rgb[0][vm]/255.0, rgb[1][vm]/255.0, rgb[2][vm]/255.0
    n = len(R)
    sample_n = min(20000, n)
    idx = np.random.RandomState(42).choice(n, sample_n, replace=False) if n > sample_n else np.arange(n)
    rgb_smp = np.stack([R[idx], G[idx], B[idx]], axis=-1)
    hsv = matplotlib.colors.rgb_to_hsv(rgb_smp.reshape(1, -1, 3)).reshape(-1, 3)
    h, s, v = hsv[:, 0], hsv[:, 1], hsv[:, 2]
    # 色は RGB そのものでプロット
    ax.scatter(h, s, c=rgb_smp, s=2, alpha=0.5)
    ax.set_xlabel("Hue (色相, 0=赤 / 0.33=緑 / 0.66=青)")
    ax.set_ylabel("Saturation (彩度)")
    ax.set_title(f"{a['label']} HSV 散布 (n={sample_n:,}/{n:,} sampled)\n"
                 f"Hue std = {a.get('hue_std', 0):.3f}", fontsize=10)
    ax.set_xlim(0, 1)
    ax.set_ylim(0, 0.3)
    ax.grid(alpha=0.3)
    # 0.0 (赤帯) と 0.66 (青帯) を強調
    ax.axvline(0.0, color="#cf222e", linestyle="--", alpha=0.4, linewidth=1)
    ax.axvline(0.66, color="#0969da", linestyle="--", alpha=0.4, linewidth=1)

fig.suptitle("HSV 散布図 — 色相 × 彩度 (各セルの色を実際の RGB で描画)",
             fontsize=12, y=1.02)
plt.tight_layout()
fig.savefig(ASSETS / "L38_fig10_hsv_scatter.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L38_fig10_hsv_scatter.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 22. 図 11: small multiples — 地形タイプ別マスク
# =============================================================================
print("\n[22] 図 11: 地形タイプ別 small multiples", flush=True)
t1 = time.time()

fig, axes = plt.subplots(3, 4, figsize=(14, 10))

mask_def = [
    ("尾根 (R強)",  1, "#e74c3c"),
    ("谷 (B強)",    2, "#2980b9"),
    ("急傾斜 (暗)", 3, "#34495e"),
    ("平地 (灰)",   4, "#bdc3c7"),
]

for row_i, area in enumerate(AREAS):
    extent = (area["bounds"].left, area["bounds"].right,
              area["bounds"].bottom, area["bounds"].top)
    for col_i, (lbl, lab_id, col) in enumerate(mask_def):
        ax = axes[row_i, col_i]
        m = (area["terrain_label"] == lab_id) & 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(lbl, 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 エリア × 4 タイプ",
             fontsize=12, y=0.99)
plt.tight_layout()
fig.savefig(ASSETS / "L38_fig11_small_multiples.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L38_fig11_small_multiples.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 23. 中間 CSV 書き出し
# =============================================================================
print("\n[23] 中間 CSV", flush=True)
t1 = time.time()

# RGB ヒスト wide
hist_rows = []
edges = np.arange(0, 256, 8)
for a in AREAS:
    vm = a["valid_mask"]
    for b, name in enumerate(["R", "G", "B"]):
        v = a["rgb"][b][vm]
        h, _ = np.histogram(v, bins=edges)
        for i, c in enumerate(h):
            hist_rows.append({"area": a["label"], "band": name,
                             "bin_low": int(edges[i]), "bin_high": int(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 / "L38_rgb_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)


# =============================================================================
# 24. HTML 生成
# =============================================================================
print("\n[24] HTML 生成", flush=True)
t1 = time.time()


def dl(name):
    return f'<a href="assets/{name}" download>{name}</a>'


sections = []

# === Section 1: 学習目標と問い ===
sec1 = f"""
<div class="note">
<b>カバー宣言</b>: 本記事は DoBoX のシリーズ <b>「CS立体図 1m / 50cm」 2 件</b>
(dataset_id = 1628, 1629) を統合し、広島県 3 自治体エリア
(<b>世羅町</b> 1m / <b>熊野町</b> 1m / <b>坂町</b> 50cm) の <b>地形微細構造</b>を
RGB ラスタの色合成として読み解く研究記事です。
さらに L10/L11 で扱った<b>土砂災害警戒区域 (yellow / red)</b> Shapefile を空間結合し、
「CS立体図のどの色域 (どの地形型) が警戒域に多いか」を量的に検証します。
</div>

<h3>本記事の独自性 — RGB 3バンド × 警戒区域の量的結合</h3>
<p>CS立体図は専門家向けの<b>視覚的</b>地形可視化として広く使われていますが、
本記事はその<b>RGB 数値そのもの</b>を解析対象として扱う点で他と異なります:</p>
<ul>
<li><b>視覚→数値</b>: 通常 CS立体図は「目で見て地形を読む」ものですが、本記事では
   R/G/B の各バンド値を統計的に集計し、3 エリアで比較します。</li>
<li><b>地形タイプ自動分類</b>: パーセンタイル適応閾値で「尾根 / 谷 / 急傾斜 / 平地 / 複雑」の
   5 クラスにラベル付けし、どのクラスがどこに集塊するかを地図化します。</li>
<li><b>警戒区域との空間結合</b>: 土砂災害警戒区域 (yellow) と特別警戒区域 (red) の
   ポリゴンを overview ラスタの座標系に rasterize し、各 zone 内の RGB 平均値を比較。
   「警戒区域は CS図でどんな色をしているか」を量的に語る。</li>
</ul>

<h3>シリーズ構造判定: (a) 同一手法・解像度違い + (b) 自治体ごとの分割</h3>
<p>このシリーズは <b>1m 版 (区分A) / 50cm 版 (区分B)</b> の 2 解像度系統で、
各自治体ごとに GeoTIFF (RGB 3バンド uint8) を切出して公開する構造。
1m 版は 22 自治体、50cm 版は 12 自治体に対応。</p>
<ul>
<li><b>1m 版 (dsid 1628)</b>: 県内 22 自治体 GeoTIFF。
   本記事は<b>世羅町 (resource_id 177706)</b>と<b>熊野町 (176860)</b>を採用
   (= L36/L37 と同自治体ペアで LiDAR ファミリの相補比較が可能)。</li>
<li><b>50cm 版 (dsid 1629)</b>: 県内 12 自治体 GeoTIFF (1m 版より少ない)。
   熊野町は 50cm 版に存在しないため、代わりに<b>坂町 (176895)</b>を採用。
   坂町は<b>平成 30 年 7 月豪雨</b>の被災地で、土石流災害が多発した地域 ⇒
   「CS立体図が崩壊地形をどう描くか」を見る研究的価値が高い。</li>
<li>3 エリアは <b>1m × 2 + 50cm × 1</b> の構成で、解像度差と地形類型差を同時に観察可能。</li>
</ul>

<h3>独自用語の定義 (本記事内で固定)</h3>
<table>
<tr><th>用語</th><th>定義</th></tr>
<tr><td><b>CS立体図 (Curvature × Slope)</b></td><td>長野県林業総合センター (戸田堅一郎ら, 2014)
が開発した地形微細構造可視化手法。標高 (DTM) から計算した<b>曲率 (Curvature)</b> と
<b>傾斜 (Slope)</b> を組み合わせ、3 バンド RGB の段彩図として描画。
凸地形 (尾根) を赤、凹地形 (谷) を青で強調、傾斜は明度として表現される。
本記事の TIFF はその標準実装の出力 (uint8, 3バンド)。</td></tr>
<tr><td><b>曲率 (Curvature)</b></td><td>地形の凹凸度。標高ラスタに対し
ラプラシアン (∂²z/∂x² + ∂²z/∂y²) を計算した値。<b>正の曲率 = 凸 (尾根・段丘崖の上端)</b>、
<b>負の曲率 = 凹 (谷・崩壊跡の凹み)</b>。CS立体図ではこの符号が R-B 軸に投影されている。</td></tr>
<tr><td><b>傾斜 (Slope)</b></td><td>地形の傾斜角度 (度)。標高ラスタの空間勾配
sqrt((∂z/∂x)² + (∂z/∂y)²) を arctan で角度に変換。CS立体図では明度 (G バンド) として表現。
急斜面ほど暗い色、緩斜面ほど明るい色。</td></tr>
<tr><td><b>R / G / B バンド</b> (本記事)</td><td>R = 尾根強度 (凸), G = 傾斜 (明度),
B = 谷強度 (凹) と本記事は解釈。実装は CS立体図の色設計を参照したもので、DoBoX 公式仕様に
明示はないが、視覚的に検証可能な解釈。</td></tr>
<tr><td><b>地形タイプ自動分類</b> (本記事独自)</td><td>R/G/B の値から
5 クラス (尾根 / 谷 / 急傾斜 / 平地 / 複雑) に自動ラベル付け。
パーセンタイル適応閾値 (上位/下位 25% 等) を使用するため、エリアごとに閾値が異なる。</td></tr>
<tr><td><b>Saturation (彩度)</b></td><td>HSV 色空間の S 成分。
S = (max − min) / max で計算。本記事では「微地形コントラストの強さ」の代理指標。
彩度が高い = 凹凸/傾斜のいずれかが強く出ている = 微地形構造あり。</td></tr>
<tr><td><b>土砂災害警戒区域 (yellow zone)</b></td><td>土砂災害防止法に基づき指定。
本記事では DoBoX dsid 48「土砂災害警戒区域・特別警戒区域情報」の 3 種別
(土石流・急傾斜地・地すべり) を統合して使用。</td></tr>
<tr><td><b>特別警戒区域 (red zone)</b></td><td>同法のレッドゾーン。建物倒壊リスクあり。
土石流・急傾斜地のみで指定 (地すべりはレッドなし)。</td></tr>
<tr><td><b>パーセンタイル適応閾値</b></td><td>絶対値ではなく「そのエリアの分布のうち上位 X%」で
閾値を決める方法。CS立体図はパステル色 (中間値) が大多数のため、絶対値閾値ではほぼ全セルが
同じクラスになってしまう。各エリアの分布に応じて適応的に決めることで、
エリア横断的な意味のある比較が可能になる。</td></tr>
<tr><td><b>古砂防地形 / 古崩壊地形</b></td><td>過去に発生した土石流・地すべり跡が地形に残ったもの。
CS立体図では微地形コントラストが強く、彩度高めで描かれる傾向。本記事の発展課題で扱う対象。</td></tr>
</table>

<h3>研究の問い (RQ)</h3>
<p>広島県の CS立体図 2 解像度 (1m / 50cm) は、それぞれ
<b>世羅町 (中山間/林業地)</b>, <b>熊野町 (都市近郊/山岳混在)</b>, <b>坂町 (平成30豪雨被災地)</b>
という<b>異なる地形リスク特性</b>のエリアを公開している。
3 エリアの CS立体図を <b>RGB 3 軸</b>で比較し、地形微細構造の構成比を量化する。
さらに<b>土砂災害警戒区域</b>との空間結合で、CS図のどの色域が警戒域に多いかを定量化する。</p>

<ol>
<li>3 dataset の<b>解像度・面積・有効データ率・ピクセル数</b>の構造比較は?</li>
<li>R, G, B 各バンドの<b>値分布</b>はエリアごとにどう違うか?</li>
<li>適応閾値による<b>地形タイプ自動分類</b>の構成比は?</li>
<li><b>HSV 散布</b>から見える 3 エリアの色相・彩度パターン差は?</li>
<li>土砂災害警戒区域内 vs 区域外で <b>R, G, B, Saturation</b> はどう違うか?</li>
<li>解像度差 (1m vs 50cm) で<b>微地形コントラスト</b>はどう変わるか?</li>
</ol>

<h3>仮説 H1〜H6</h3>
<ul>
<li><b>H1 (世羅町は緩傾斜丘陵)</b>: 世羅の CS図 Saturation 平均 ≤ 0.10。
中山間の植林地は緩い丘陵地形が中心で、極端な凹凸が少ない。</li>
<li><b>H2 (熊野町は急峻地形)</b>: 熊野の Sat 平均 > 世羅。
熊野は花崗岩マサ土地形で急傾斜面が多い。</li>
<li><b>H3 (坂町は崩壊跡を含む急傾斜)</b>: 坂町の Sat 平均が 3 エリアで最大。
平成30豪雨の崩壊地形を CS図が刻み込んでいる。</li>
<li><b>H4 (警戒区域内の CS図は微地形強)</b>: 特別警戒区域内の Saturation は
区域外より大きい。急斜面で凹凸が強い場所は警戒区域に指定されているはず。</li>
<li><b>H5 (色相は地形タイプを区別)</b>: HSV の Hue 標準偏差 > 0.05 で、
赤系 (尾根) と青系 (谷) の双方が混在することを確認。</li>
<li><b>H6 (50cm 版は微細構造を捉える)</b>: 坂町 (50cm) の Sat 標準偏差は
他 1m 版より大きい。高解像度ほど微地形のばらつきを拾う。</li>
</ul>

<h3>到達点</h3>
<p>3 エリアの CS立体図を「3 自治体の代表エリア」として読み解き、
合計 <b>{int(metas[0]['n_cells']/1e6 + metas[1]['n_cells']/1e6 + metas[2]['n_cells']/1e6):,} M セル</b>
× 3 バンド ≈ {int((metas[0]['n_cells'] + metas[1]['n_cells'] + metas[2]['n_cells'])*3/1e6):,} M 数値の RGB 段彩データを
overview ピラミッドで効率的に概観する。
さらに<b>土砂災害警戒区域 Shapefile (~ 2,000 ポリゴン) との空間結合</b>で
「CS図と警戒区域の関係」を量的に検証する。
これは CS立体図単独・警戒区域単独では不可能な、両者の<b>同時空間整合分析</b>である。</p>

<div class="warn"><b>本記事のスコープ外</b>:
<ul>
<li>個別ピクセルの<b>地形微細構造判読</b> (古砂防地形 / 段丘崖 / 地すべり地塊 等の手動デリネーション)</li>
<li>L36 樹高 / L37 点群密度 との 3 軸結合 (各ラスタの座標系統合は実装上重い → 別記事 X 系で扱う)</li>
<li><b>標高図 (DTM, dsid 1635/1636) との直接比較</b> (発展課題)</li>
<li><b>長野県オリジナル CS立体図との配色合致確認</b> (DoBoX 実装の配色は標準実装と微妙に異なる可能性あり)</li>
<li>過去の土石流・地すべり履歴データ (DoBoX 外) との時系列照合</li>
</ul>
本記事は CS立体図の<b>2 dataset (1m + 50cm) のみを統合</b>し、土砂警戒区域との空間結合に特化する。
</div>
"""
sections.append(("学習目標と問い", sec1))


# === Section 2: 使用データ ===
sec2 = f"""
<p>本記事が使用する 2 dataset の一覧。両者は<b>同一手法 (CS立体図 = 曲率×傾斜の RGB 段彩) の
解像度違い</b>で、各自治体ごとに GeoTIFF を切り出して公開する構造。本記事は 2 dataset から
合計 <b>3 自治体</b>を代表として用います。</p>

<table>
<tr><th>dsid</th><th>名称</th><th>解像度</th><th>公開単位</th><th>本記事の代表</th><th>TIFF サイズ</th><th>セル数 (×3バンド)</th><th>DoBoX</th></tr>
<tr><td><b>1628</b></td><td>CS立体図（1m）</td><td>1.0 m/cell</td>
  <td>22 自治体別 GeoTIFF</td><td>世羅町 (34462) + 熊野町 (34307)</td>
  <td>{int(AREAS[0]['tif'].stat().st_size/(1024*1024)):,} + {int(AREAS[1]['tif'].stat().st_size/(1024*1024)):,} MB</td>
  <td>{(AREAS[0]['meta']['n_cells']+AREAS[1]['meta']['n_cells'])/1e6:.1f} M</td>
  <td><a href='https://hiroshima-dobox.jp/datasets/1628' target='_blank'>#1628</a></td></tr>
<tr><td><b>1629</b></td><td>CS立体図（50cm）</td><td>0.5 m/cell</td>
  <td>12 自治体別 GeoTIFF</td><td>坂町 (34384)</td>
  <td>{int(AREAS[2]['tif'].stat().st_size/(1024*1024)):,} MB</td>
  <td>{AREAS[2]['meta']['n_cells']/1e6:.1f} M</td>
  <td><a href='https://hiroshima-dobox.jp/datasets/1629' target='_blank'>#1629</a></td></tr>
</table>

<h3>データ仕様 (両 dataset 共通)</h3>
<ul>
<li><b>形式</b>: GeoTIFF, <b>3 バンド uint8 (0-255)</b> = RGB 段彩図</li>
<li><b>CRS</b>: EPSG:6671 (JGD2011 / Japan Plane Rectangular CS III) — m 単位の平面直角座標</li>
<li><b>nodata</b>: 1m 版は 0 (3バンド全て 0 = 無効), 50cm 版 (坂町) は明示なし (推定: 全バンド 0 が境界)</li>
<li><b>値の意味</b>: <b>R/G/B 各バンド 0-255</b>。本記事の解釈では:
  R = 尾根強度 (凸地形), G = 傾斜 (明度), B = 谷強度 (凹地形)</li>
<li><b>派生過程</b>: 標高 (DTM) → 曲率計算 + 傾斜計算 → RGB 合成 → 段彩図</li>
<li><b>切出単位</b>: 各自治体の行政界 + 100m バッファで矩形切出 (L36/L37 と同様)</li>
<li><b>用途</b>: 地形微細構造の視覚判読 (尾根/谷/段丘/崩壊地形/古砂防 等)</li>
<li><b>ライセンス</b>: CC-BY 4.0 (測量法 44 条に基づく公共測量 2 次著作物)</li>
</ul>

<h3>L36 樹高 / L37 点群密度 との関係 (本記事は単独, 参考併置のみ)</h3>
<table>
<tr><th>軸</th><th>L36 樹高</th><th>L37 点群密度</th><th>L38 CS立体図 (本記事)</th></tr>
<tr><td>計測対象</td><td>樹冠の高さ (m)</td><td>地面到達パルス数</td><td>地形の微細構造 (RGB)</td></tr>
<tr><td>派生過程</td><td>DSM − DTM</td><td>グラウンド分類点数</td><td>DTM → 曲率 + 傾斜 → RGB</td></tr>
<tr><td>同じ点群?</td><td colspan="3" style="text-align:center;"><b>YES — 同一の航空レーザ計測点群</b>から
   派生した 3 つの集計</td></tr>
<tr><td>表現</td><td>連続値 float32</td><td>整数 uint8/int16</td><td>RGB 3 バンド uint8</td></tr>
<tr><td>本記事のスコープ</td><td colspan="3" style="text-align:center;">
   L38 単独で完結 (L36/L37 とは合体しない)</td></tr>
</table>

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

<h3>データ品質メモ</h3>
<ul>
<li><b>3 バンド uint8</b>: 各バンド 0-255 の整数値。RGB 3 チャンネルが
   それぞれ独立した地形指標 (尾根/傾斜/谷) を表現。</li>
<li><b>1m 版は overview あり</b> (2/4/8/16/32/64/128), 熊野町は overview なし (= 内部に
   ピラミッドが構築されていない)。坂町 (50cm) は overview あり。</li>
<li><b>nodata の扱い</b>: 1m 版は (0,0,0) = 無効。50cm 版 (坂町) は nodata 未指定だが、
   実データでは (0,0,0) を境界マスクとして同様に扱う。</li>
<li><b>値の典型範囲</b>: 全エリアで R/G/B 平均が 160-200 の中間色帯に集中
   (= パステル系の地図画像)。極端な R 高 (>240) や B 高 (>240) は少数。</li>
<li><b>解像度差は overview で見えにくい</b>: 全エリアを overview /1/32〜64 で読込むため、
   解像度差 (1m vs 50cm) は<b>生解像度サンプル</b>でのみ確認できる (図 8)。</li>
</ul>
"""
sections.append(("使用データ", sec2))


# === Section 3: ダウンロード ===
sec3 = f"""
<h3>生データ (DoBoX 直リンク)</h3>
<ul>
<li>CS立体図（1m）_世羅町:
  <a href="{AREAS[0]['zip_url']}">直 DL ZIP</a> /
  <a href="https://hiroshima-dobox.jp/resources/177706">DoBoX resource 177706</a> /
  <a href="https://hiroshima-dobox.jp/datasets/1628">dsid 1628</a></li>
<li>CS立体図（1m）_熊野町:
  <a href="{AREAS[1]['zip_url']}">直 DL ZIP</a> /
  <a href="https://hiroshima-dobox.jp/resources/176860">DoBoX resource 176860</a> /
  <a href="https://hiroshima-dobox.jp/datasets/1628">dsid 1628</a></li>
<li>CS立体図（50cm）_坂町:
  <a href="{AREAS[2]['zip_url']}">直 DL ZIP</a> /
  <a href="https://hiroshima-dobox.jp/resources/176895">DoBoX resource 176895</a> /
  <a href="https://hiroshima-dobox.jp/datasets/1629">dsid 1629</a></li>
</ul>

<div class="warn"><b>容量警告</b>: 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。
他自治体や 50cm 広島市版はさらに大きい (世羅町 50cm 版だけで 3.6 GB) ため、本記事のスコープ外。</div>

<h3>土砂災害警戒区域データ (本記事の空間結合に必要)</h3>
<p>本記事は<a href="L10_sediment_disaster_cross.html">L10 土砂災害警戒区域 × 用途地域</a>と
同じ Shapefile (DoBoX dsid 48) を使用します。
<code>data/extras/sediment_shp/</code> に展開済みであれば自動利用。</p>

<h3>中間データ (本記事生成 CSV)</h3>
<ul>
<li>{dl("L38_series.csv")} — 3 dataset の一覧 (代表自治体・解像度・サイズ)</li>
<li>{dl("L38_tiff_meta.csv")} — TIFF メタ (CRS, bounds, セル数, overview, dtype, nodata, 3 バンド)</li>
<li>{dl("L38_basic_stats.csv")} — RGB 各バンド統計 (平均・中央値・std) + Saturation</li>
<li>{dl("L38_terrain_classes.csv")} — 5 種地形タイプの構成比 (尾根/谷/急/平地/複雑)</li>
<li>{dl("L38_rgb_histogram.csv")} — R/G/B 各バンドの 32 ビン ヒストグラム</li>
<li>{dl("L38_sediment_rgb_compare.csv")} — 警戒区域 内/外 の RGB 平均比較</li>
<li>{dl("L38_hypothesis_results.csv")} / {dl("L38_hypothesis_results.json")} — H1〜H6 検証結果</li>
</ul>

<h3>図 (本記事生成 PNG)</h3>
<ul>
<li>{dl("L38_fig1_dataset_overview.png")} — 3 dataset カード + Saturation ヒスト</li>
<li>{dl("L38_fig2_raster_maps.png")} — CS立体図 RGB 主題図 (3 エリア)</li>
<li>{dl("L38_fig3_terrain_classification.png")} — <b>地形タイプ自動分類マップ</b></li>
<li>{dl("L38_fig4_rgb_histograms.png")} — R/G/B 各バンドのヒストグラム比較</li>
<li>{dl("L38_fig5_terrain_composition.png")} — 5 種地形タイプ構成比 (積み上げ)</li>
<li>{dl("L38_fig6_sediment_overlay.png")} — <b>CS立体図 × 土砂警戒区域 重ね合わせ</b></li>
<li>{dl("L38_fig7_sediment_rgb_compare.png")} — 警戒区域 内/外 RGB 平均比較</li>
<li>{dl("L38_fig8_center_samples.png")} — 中央 1024×1024 セル本来解像度サンプル</li>
<li>{dl("L38_fig9_concept.png")} — CS立体図のしくみ概念図</li>
<li>{dl("L38_fig10_hsv_scatter.png")} — HSV 散布図 (色相 × 彩度)</li>
<li>{dl("L38_fig11_small_multiples.png")} — 地形タイプ別 small multiples (3×4)</li>
</ul>

<h3>再現用 Python スクリプト</h3>
<p><a href="L38_cs_terrain.py" download>L38_cs_terrain.py</a> を取得して
プロジェクトルートで <code>py -X utf8 lessons/L38_cs_terrain.py</code> を実行。
TIFF が無ければ ZIP 自動 DL → 自動展開 → overview 読込で分析が走ります
(初回 DL を除けば<b>約 13 秒</b>で完走)。</p>
"""
sections.append(("ダウンロード", sec3))


# === Section 4: 分析 1 — dataset の構造を見る ===
sec4 = f"""
<h3>狙い</h3>
<p>2 dataset (1m / 50cm) と 3 自治体 (世羅 / 熊野 / 坂) の<b>規模・形状・データ量</b>を
1 枚の絵で示す。Saturation ヒストグラムを重ねて<b>地形微細構造の強さ</b>を最初に視覚化する。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p>L36 / L37 と全く同じ rasterio スタイルだが、CS立体図には重要な違いがある:</p>
<ul>
<li><b>3 バンド読込</b>: <code>ds.read()</code> でバンド指定なしに呼ぶと <b>(3, H, W)</b> の 3D 配列が返る。
   各バンドは uint8 (0-255) の整数値。</li>
<li><b>RGB ではなく数値</b>: バンド 1 = R, バンド 2 = G, バンド 3 = B。各セルは
   3 つの 0-255 整数で「色」が指定される。本記事はこの 3 値を独立した地形指標として扱う。</li>
<li><b>Saturation 算出</b>: HSV 変換で「彩度」を計算。S = (max(R,G,B) - min(R,G,B)) / max。
   微地形コントラストが強い場所ほど Sat が高い。</li>
</ul>

<h3>実装</h3>
{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

        # 3 バンドを 1 度に読込 → shape (3, H, W)
        rgb = ds.read(out_shape=(3, ov_h, ov_w),
                      resampling=Resampling.average).astype(np.float32)
        area["rgb"] = rgb

        # nodata マスク (全バンド 0 = 無効)
        valid_mask = rgb.sum(axis=0) > 0
        area["valid_mask"] = valid_mask

        # HSV の Saturation 計算
        maxv = np.maximum.reduce(rgb)
        minv = np.minimum.reduce(rgb)
        sat = np.where(maxv > 0, (maxv - minv) / np.maximum(maxv, 1), 0)
        area["sat"] = sat
''')}

<h3>図と読み取り (図 1: dataset カード + Saturation ヒスト)</h3>
<p><b>なぜこの図か</b>: 3 dataset の規模と特性をカードで言語的に、
かつ Saturation ヒストグラムで「微地形コントラスト分布」として量的に、
<b>同時に 1 枚で伝える</b>ため。</p>

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

<p><b>読み取り</b>:</p>
<ul>
<li>3 エリアとも RGB 平均が <b>R > G > B</b> の順で減衰する典型的なパステル色域に集中。
   Sera 202.4 / 192.5 / 187.9, Kumano 196.3 / 185.2 / 180.8, Saka 175.5 / 162.2 / 159.4。</li>
<li>Saturation ヒスト (右) は<b>3 エリアともほぼ 0 にピーク</b>を持つ右裾分布。
   多くのセルが「ほぼ無彩色」 (中間グレー)、少数のセルが「彩度高め」 (強い凹凸)。</li>
<li>Sat 平均は<b>世羅 0.074 < 熊野 0.083 < 坂 0.096</b>。順位は仮説 (世羅 < 熊野 < 坂) と一致。
   特に坂町は明確に右にずれており、急傾斜・崩壊地形の影響が見える。</li>
<li>坂町の R 平均が顕著に低い (175.5) のは「全体的に暗い色 = 急斜面影が多い」と読める。</li>
</ul>

<h3>表 (基本統計)</h3>
{stats_df.round(2).to_html(classes='', border=0)}

<p><b>読み取り</b>: R/G/B 平均は 3 エリアで似た値域 (160-202) だが、<b>標準偏差は明確に違う</b> ⇒
世羅 R-std=15.7, 熊野 19.7, 坂 42.0。坂町は標準偏差が約 3 倍も大きい =
「明るい部分と暗い部分が広く分布」 ⇒ <b>急傾斜面と平地の混在による陰影コントラストの強さ</b>を示す。
Sat 平均と Sat 標準偏差はどちらも坂町が最大で、<b>仮説 H3 (坂町は崩壊跡を含む急傾斜)</b> を支持する強い量的証拠。</p>
"""
sections.append(("分析 1: 3 dataset とエリアの構造を可視化", sec4))


# === Section 5: 分析 2 — CS立体図 RGB 主題図 ===
sec5 = f"""
<h3>狙い</h3>
<p>CS立体図の<b>RGB 元画像</b>を 3 エリア並列で表示し、
「地形微細構造はどう見えるか」を視覚的に理解する。これは CS立体図本来の使い方
(目で見て地形を読む) を踏襲したセクション。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p><code>matplotlib.pyplot.imshow</code> に <b>(H, W, 3)</b> 形状の uint8 配列を渡すと
RGB 画像として表示される。rasterio の <code>read()</code> 結果は <b>(3, H, W)</b> なので
<code>np.stack([R, G, B], axis=-1)</code> で軸を入れ替える必要がある。</p>
<ul>
<li><b>赤い領域</b> = R が突出 = <b>尾根線 (凸地形)</b>。標高ラスタの正の曲率部分。</li>
<li><b>青い領域</b> = B が突出 = <b>谷線・崩壊地形 (凹地形)</b>。標高ラスタの負の曲率部分。</li>
<li><b>暗い領域</b> = 全バンド低 = <b>急斜面の影</b>。標高ラスタの傾斜が強い場所。</li>
<li><b>パステルピンク〜灰</b> = 緩傾斜の平地・台地・尾根丘陵。CS立体図の「背景色」。</li>
</ul>

<h3>実装</h3>
{code('''
# 3 バンドを RGB 画像形式に整形
rgb_disp = np.stack([rgb[0], rgb[1], rgb[2]], axis=-1).astype(np.uint8)
# nodata セルを白に置換
rgb_disp_show = np.where(valid_mask[..., None], rgb_disp, 255)
ax.imshow(rgb_disp_show, extent=area_bounds, origin="upper", interpolation="nearest")
''')}

<h3>図と読み取り (図 2: CS立体図 RGB 主題図)</h3>
<p><b>なぜこの図か</b>: 3 エリアの<b>地形微細構造</b>を視覚的に並べ、
専門家のように「赤=尾根、青=谷、暗=急斜面」を読み取る練習をするため。</p>

{figure("assets/L38_fig2_raster_maps.png", "図 2: 3 エリア CS立体図 RGB 元画像")}

<p><b>読み取り</b>:</p>
<ul>
<li><b>世羅町</b> (左): 全体に淡いピンク〜灰色のパステル基調。<b>河川沿いに細い青線</b>
   (= 谷地形) が走り、それ以外は緩やかな起伏。中山間の典型的な丘陵地形。</li>
<li><b>熊野町</b> (中): 中央〜南東に暗い色域 (= 急斜面が広い山地)。
   <b>谷部の青色が筋状に強く出る</b>のが特徴的で、花崗岩マサ土地形特有の急峻な谷。</li>
<li><b>坂町</b> (右): 全体に色のコントラストが強く、<b>暗い影 (急斜面) と明るい尾根</b>が
   斑模様に交互する。<b>南東の住宅地</b>付近に明確な凹凸地形 (= 平成30豪雨で土石流が
   発生した山地と扇状地の境界域) が見える。</li>
<li>3 エリアとも<b>北 (上) が山地、南 (下) が平地・海岸寄り</b>の傾向。
   坂町は最も小さい (4.4 × 6.9 km) ため細部構造が見える。</li>
</ul>

<h3>表 (TIFF メタ + 行政界面積)</h3>
<table>
<tr><th>項目</th><th>世羅町</th><th>熊野町</th><th>坂町</th></tr>
<tr><td>解像度 (m/cell)</td><td>{AREAS[0]['resolution_m']}</td><td>{AREAS[1]['resolution_m']}</td><td>{AREAS[2]['resolution_m']}</td></tr>
<tr><td>セル数 (M, ×3バンド)</td><td>{AREAS[0]['meta']['n_cells']/1e6:.1f}</td><td>{AREAS[1]['meta']['n_cells']/1e6:.1f}</td><td>{AREAS[2]['meta']['n_cells']/1e6:.1f}</td></tr>
<tr><td>TIFF bbox 面積 (km²)</td><td>{AREAS[0]['meta']['area_km2']:.1f}</td><td>{AREAS[1]['meta']['area_km2']:.1f}</td><td>{AREAS[2]['meta']['area_km2']:.1f}</td></tr>
<tr><td>行政界面積 (km²)</td><td>{AREAS[0]['muni_area_km2']:.1f}</td><td>{AREAS[1]['muni_area_km2']:.1f}</td><td>{AREAS[2]['muni_area_km2']:.1f}</td></tr>
<tr><td>有効セル比率</td><td>{basic_stats['sera_1m']['valid_ratio']*100:.1f}%</td><td>{basic_stats['kumano_1m']['valid_ratio']*100:.1f}%</td><td>{basic_stats['saka_50cm']['valid_ratio']*100:.1f}%</td></tr>
<tr><td>nodata 値</td><td>0 (全バンド)</td><td>0 (全バンド)</td><td>None (= (0,0,0) を境界マスクとして扱う)</td></tr>
<tr><td>overview 倍率</td><td>/1/{AREAS[0]['arr_factor']}</td><td>/1/{AREAS[1]['arr_factor']} (内部 ovr 無)</td><td>/1/{AREAS[2]['arr_factor']}</td></tr>
<tr><td>実効セル (m/cell)</td><td>{AREAS[0]['arr_cell_m']:.0f}</td><td>{AREAS[1]['arr_cell_m']:.0f}</td><td>{AREAS[2]['arr_cell_m']:.0f}</td></tr>
</table>
<p><b>読み取り</b>: 世羅町は 27 × 16 km と広大で 446 M セル × 3 = 約 1.3 GB の生 TIFF。
熊野町は内部 overview を持たない (= ピラミッド未構築) ため out_shape で擬似粗化。
坂町は 50cm 解像度のため 30 km² で 122 M セル相当。<b>解像度差 = ピクセル数 4 倍</b>がはっきり分かる。</p>
"""
sections.append(("分析 2: CS立体図 RGB 主題図化", sec5))


# === Section 6: 分析 3 — 地形タイプ自動分類 ===
sec6 = f"""
<h3>狙い</h3>
<p>専門家の目で「赤い場所は尾根、青い場所は谷」と判読する代わりに、
<b>RGB 値からアルゴリズム的に地形タイプを自動分類</b>する。
これにより 3 エリアの地形タイプ構成比を量的に比較可能になる。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p>本記事独自の<b>パーセンタイル適応閾値</b>分類を使う。CS立体図の RGB は
パステル色 (中間値) が大多数なので、絶対値閾値ではほぼ全セルが同じクラスに振られる。
そこで各エリアの分布から相対的に閾値を決める。</p>

<table>
<tr><th>クラス</th><th>条件</th><th>意味</th></tr>
<tr><td><b>1. 尾根</b></td><td>R - B >= P75 (上位 25%)</td><td>R が B より明らかに大きい = 凸地形</td></tr>
<tr><td><b>2. 谷</b></td><td>R - B <= P25 (下位 25%)</td><td>B が R より明らかに大きい = 凹地形</td></tr>
<tr><td><b>3. 急傾斜</b></td><td>brightness <= P15 (下位 15%, 上記 1, 2 ではない)</td><td>暗い色 = 影 = 急斜面</td></tr>
<tr><td><b>4. 平地</b></td><td>3 バンド std <= P25, 上記いずれでもない</td><td>ほぼグレー = 緩傾斜の平地</td></tr>
<tr><td><b>5. 複雑地形</b></td><td>上記いずれでもない (= 中間色だが彩度あり)</td><td>細かい起伏が混在する地形</td></tr>
</table>

<p><b>限界と代替</b>:</p>
<ul>
<li><b>限界</b>: 適応閾値はエリア横断で「絶対値が違う」ので
   「世羅の尾根」と「坂町の尾根」は厳密には同じ閾値ではない。
   ただし<b>各エリア内での相対ランキング</b>として意味がある。</li>
<li><b>代替案</b>: 絶対値閾値 (例 R-B >= 30) で揃える。ただしこれだと CS立体図の中間色域では
   ほぼ全セルが「平地」になってしまう (最初の試行では 100% 平地)。</li>
<li><b>真の正解</b>: 標高ラスタ (DTM) から曲率・傾斜を直接計算する (発展課題 Z2)。</li>
</ul>

<h3>実装</h3>
{code('''
def classify_terrain(rgb, valid_mask):
    """RGB を 5 種地形タイプにラベル付け (パーセンタイル適応閾値)"""
    R, G, B = rgb[0], rgb[1], rgb[2]
    diff_RB = R - B
    bright = (R + G + B) / 3.0
    band_std = np.stack([R, G, B], axis=0).std(axis=0)

    # 各エリアの分布から閾値を決める
    p75 = np.percentile(diff_RB[valid_mask], 75)
    p25 = np.percentile(diff_RB[valid_mask], 25)
    bp15 = np.percentile(bright[valid_mask], 15)
    sp25 = np.percentile(band_std[valid_mask], 25)

    label = np.zeros(R.shape, dtype=np.int8)
    is_ridge  = (diff_RB >= p75) & valid_mask
    is_valley = (diff_RB <= p25) & valid_mask
    is_steep  = (bright <= bp15) & ~is_ridge & ~is_valley & valid_mask
    is_flat   = (band_std <= sp25) & ~is_ridge & ~is_valley & ~is_steep & valid_mask
    is_complex = valid_mask & ~is_ridge & ~is_valley & ~is_steep & ~is_flat

    label[is_ridge]   = 1
    label[is_valley]  = 2
    label[is_steep]   = 3
    label[is_flat]    = 4
    label[is_complex] = 5
    return label
''')}

<h3>図と読み取り (図 3: 地形タイプ分類マップ)</h3>
<p><b>なぜこの図か</b>: RGB 元画像 (図 2) の隣に分類結果を並べて、
「色がどの分類クラスに対応するか」を視覚的に確認する。</p>

{figure("assets/L38_fig3_terrain_classification.png", "図 3: 地形タイプ自動分類マップ (5 クラス)")}

<p><b>読み取り</b>:</p>
<ul>
<li>3 エリアとも<b>赤 (尾根) と青 (谷) が拮抗</b>して 25-30% ずつ分布。
   これはパーセンタイル閾値が「上位 25% 以上」「下位 25% 以下」で対称定義のため。</li>
<li><b>急傾斜 (濃灰)</b> は 3 エリアで 4-6% と少なめ。CS立体図の「明度」軸が圧縮されているため。</li>
<li><b>複雑地形 (紫)</b> は 38-42% と最大の比率を占める。
   = 中間色だが std や bright で例外的なセル群 = <b>細かい起伏が混在する典型的な山岳地形</b>。</li>
<li>世羅町は<b>河川沿いに青い「谷」帯</b>が連続的に走り、その周囲に赤い「尾根」が伴う。
   = 水系構造が CS立体図に明確に刻まれている。</li>
<li>坂町は<b>各タイプが斑状に細かく混在</b>。50cm 解像度で微地形が拾えていることが原因。</li>
</ul>

<h3>図と読み取り (図 5: 地形タイプ構成比 stack-bar)</h3>

{figure("assets/L38_fig5_terrain_composition.png", "図 5: 5 種地形タイプの構成比 (3 エリア)")}

<p><b>読み取り</b>: 3 エリアの構成比は<b>非常に似ている</b> (尾根 25-28%, 谷 25-30%, 急 4-6%, 平地 0-1%, 複雑 38-42%)。
これは「パーセンタイル適応閾値」の必然的帰結 — 各エリア内で上位/下位 25% を取るので、
<b>クラス間バランスはエリア横断でほぼ均一</b>になる。
意味のある違いは「<b>どの場所がそのクラスに分類されているか</b>」 (= 空間分布) にあり、
それは図 3 と次節 (small multiples) で見る。</p>

<h3>図と読み取り (図 11: 地形タイプ別 small multiples)</h3>

{figure("assets/L38_fig11_small_multiples.png", "図 11: 地形タイプ別の空間分布 small multiples (3 エリア × 4 タイプ)")}

<p><b>読み取り</b>:</p>
<ul>
<li><b>尾根 (赤列)</b>: 3 エリアとも町域全体に網目状に広がる。<b>世羅は線状に連続</b>、
   <b>熊野は山地北側に集塊</b>、<b>坂は短い尾根が斑状に分散</b>。</li>
<li><b>谷 (青列)</b>: <b>水系のネットワーク構造</b>がはっきり可視化される。
   特に世羅町は河川沿いに谷が連続的に走り、樹枝状の流域パターンが読める。</li>
<li><b>急傾斜 (黒列)</b>: 全エリアで少数だが、<b>熊野町は山地に集塊</b>、
   <b>坂町は南東部 (= 豪雨被災地)</b> に集塊する明瞭なパターン。</li>
<li><b>平地 (灰列)</b>: 3 エリアでほぼ 0% に近い (band_std P25 がそもそも低いため)。
   この分類軸は本データセットでは情報量が少ない。</li>
</ul>

<h3>表 (地形タイプ構成比 wide)</h3>
{class_df.round(2).to_html(classes='', border=0, index=False)}
"""
sections.append(("分析 3: 地形タイプ自動分類 (パーセンタイル適応閾値, 5 クラス)", sec6))


# === Section 7: 分析 4 — RGB バンド分布と HSV 散布 ===
sec7 = f"""
<h3>狙い</h3>
<p>3 エリアの<b>R, G, B 各バンドの値分布</b>を直接比較し、
HSV 散布図で「色相 × 彩度」の 2 次元パターンを観察する。
バンドごとの違いは「どの地形指標がエリアごとに強いか」を語る。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p>2 つの可視化を併用:</p>
<ul>
<li><b>バンドごとのヒスト</b>: R, G, B それぞれを 0-255 で 4 刻みのヒストにし、3 エリア重ね。
   各バンドが「どこにピークを持ち、どんな裾を持つか」を確認。</li>
<li><b>HSV 散布</b>: RGB を HSV (Hue, Saturation, Value) に変換し、Hue × Sat 平面に散布。
   各点を実際の RGB 色で描画する (= <b>各セルが実物のどの色か</b>を視認可能)。</li>
</ul>

<p><b>HSV とは</b>: RGB の代替表現で、
H = 色相 (0=赤, 0.33=緑, 0.66=青), S = 彩度, V = 明度。
S = (max - min) / max の式で<b>「色がどれだけ偏っているか」</b>を 0-1 で表現する。
グレー (R=G=B) は S=0、純赤 (R=255, G=B=0) は S=1。</p>

<h3>実装</h3>
{code('''
import matplotlib

# RGB を HSV に変換 (各値 0-1 にスケール)
rgb_norm = (rgb_arr / 255.0).reshape(-1, 3)
hsv = matplotlib.colors.rgb_to_hsv(rgb_norm.reshape(1, -1, 3)).reshape(-1, 3)
hue, sat, val = hsv[:, 0], hsv[:, 1], hsv[:, 2]

# 散布図: 各点の色は実際の RGB そのまま
ax.scatter(hue, sat, c=rgb_norm, s=2, alpha=0.5)
''')}

<h3>図と読み取り (図 4: RGB バンドヒスト)</h3>
<p><b>なぜこの図か</b>: 3 バンドを別々に並べて、各エリアの値域シフトを見るため。
散布図ではなくヒストにすることで、各バンドの中央値・裾の重さが分かりやすい。</p>

{figure("assets/L38_fig4_rgb_histograms.png", "図 4: R/G/B 各バンドの 3 エリア比較ヒストグラム")}

<p><b>読み取り</b>:</p>
<ul>
<li><b>R バンド (左)</b>: 3 エリアとも 180-220 にピーク。世羅 (μ=202) が最右、坂 (μ=176) が最左。
   <b>R が低い = 全体的に赤の成分が弱い = 尾根の凸成分が抑えられている</b>。</li>
<li><b>G バンド (中)</b>: ほぼ同形だがピーク位置が R より約 10 だけ低い。
   傾斜による明度低下を表現する G の役割が見える。</li>
<li><b>B バンド (右)</b>: 値域は R, G より<b>幅広い</b> (std が大きい)。
   = 谷の凹成分は地形によって大きく変動する。</li>
<li>坂町だけ<b>裾が左に伸びる</b> (R, G, B とも 0-100 域にも分布) ⇒
   暗い色 (急斜面の影) を多く含む = 急傾斜地形が広がる証拠。</li>
</ul>

<h3>図と読み取り (図 10: HSV 散布)</h3>
<p><b>なぜこの図か</b>: RGB の 3 軸では分布パターンが見えづらい。
HSV に変換して 2 次元の Hue × Sat 平面で見ると、
「赤系の凸地形」「青系の凹地形」「灰系の平地」のクラスタリング構造がよく見える。</p>

{figure("assets/L38_fig10_hsv_scatter.png", "図 10: HSV 散布図 — Hue × Sat 平面 (各点を実 RGB 色で描画)")}

<p><b>読み取り</b>:</p>
<ul>
<li><b>3 エリアとも</b> Hue が <b>0 (赤系) と 0.6-0.7 (青系)</b> の 2 帯に分かれて分布。
   仮説 H5 (色相は 2 峰性) を視覚的に支持。</li>
<li>Sat 軸は 3 エリアとも 0 - 0.2 の狭い範囲に圧縮 (= パステル色)。
   坂町だけ 0.25 付近に裾を持つ (= 微地形コントラスト強)。</li>
<li>Hue 0.0 (赤) のクラスタは「尾根」、0.66 (青) は「谷」、0.33 (緑) は「傾斜」、
   中間 (灰) は「平地」と読める。<b>CS立体図の配色設計が散布図上で 3 主帯に分かれる</b>のが視覚化される。</li>
<li>Hue 標準偏差は<b>世羅 0.057 < 熊野 0.136 < 坂 0.252</b>と大幅に増加。
   坂町は色相が広く分散 = <b>地形タイプが多様</b>。</li>
</ul>
"""
sections.append(("分析 4: RGB バンド分布と HSV 散布", sec7))


# === Section 8: 分析 5 — 警戒区域 × CS立体図 (本記事の核心) ===
sec8 = f"""
<h3>狙い (本記事の核心)</h3>
<p>土砂災害警戒区域 (yellow) と特別警戒区域 (red) の<b>ポリゴン</b>を
CS立体図の<b>ラスタ座標系に rasterize</b> して、
「警戒区域内 vs 区域外で CS立体図の RGB 平均値はどう違うか」を量的に比較する。
これは CS立体図単独・警戒区域単独では出せない<b>本記事の独自貢献</b>。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p><b>「ポリゴン → ラスタ化 (rasterize)」</b>とは: 不規則な形のポリゴン (警戒区域)
を「画素ごとに 0/1 のマスク」に変換する操作。具体的には:</p>
<ol>
<li>警戒区域ポリゴンを rasterio の <code>features.rasterize</code> に渡し、
  対象ラスタと同じ shape・transform で描画。</li>
<li>得られたバイナリマスク (1=区域内, 0=外) と CS立体図ラスタを<b>セル単位でペア化</b>。</li>
<li>各 zone (特別警戒 / 警戒 / 区域外) に分けて R, G, B, Sat の平均を計算。</li>
<li>3 ゾーン × 3 エリア × 4 指標 (R, G, B, Sat) のクロス表を作る。</li>
</ol>

<p><b>パッケージ</b>: <code>rasterio.features.rasterize</code>。<code>geopandas.overlay</code> や
<code>shapely.contains</code> よりも<b>ラスタ向けに最適化</b>されており、
60K セル × 数千ポリゴンでも数秒で完了。</p>

<h3>実装</h3>
{code('''
from rasterio.features import rasterize
from rasterio.transform import from_bounds

def rasterize_polys(polys_gdf, area):
    """警戒区域 polygons をエリア overview の (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)]:
    R_mean = rgb[0][m].mean()
    G_mean = rgb[1][m].mean()
    B_mean = rgb[2][m].mean()
    Sat_mean = sat[m].mean()
''')}

<h3>図と読み取り (図 6: 警戒区域 × CS overlay)</h3>
<p><b>なぜこの図か</b>: 警戒区域がCS立体図の<b>どの色域</b>に集塊しているかを<b>視覚的</b>に確認。
ポリゴンの透明度を調整して下地の RGB 段彩を見せつつ、警戒域 (yellow/red) を半透明で重ねる。</p>

{figure("assets/L38_fig6_sediment_overlay.png", "図 6: CS立体図 + 土砂災害警戒区域 (yellow/red) の重ね合わせ")}

<p><b>読み取り</b>:</p>
<ul>
<li><b>世羅町</b>: 警戒区域 (黄) が町域全体に分散。地すべり警戒区域が広い。
   特別警戒 (赤) は山地の急斜面に集塊。</li>
<li><b>熊野町</b>: 警戒区域は<b>山地に集中</b>。平地住宅地はほぼ警戒区域外。
   花崗岩マサ土地形の急傾斜面が警戒域に明確に対応。</li>
<li><b>坂町</b>: 警戒区域 (黄) と特別警戒 (赤) が<b>町域の山地側を広く覆う</b>。
   南西の住宅地と山地の境界に特別警戒が集中 = 平成30豪雨で土石流が
   発生した地形と完全に対応。</li>
<li>3 エリアともに警戒区域は<b>急斜面 (暗い色域)</b> や<b>谷 (青色域)</b> に重なる傾向が見える。</li>
</ul>

<h3>図と読み取り (図 7: RGB 平均比較 grouped bar)</h3>
<p><b>なぜこの図か</b>: 視覚的な印象 (図 6) を量的に裏付ける。
3 エリア × 3 ゾーン × 3 指標 (R, R-B, Sat) のクロス比較。</p>

{figure("assets/L38_fig7_sediment_rgb_compare.png", "図 7: 警戒区域 内/外 の RGB / R-B / Sat 平均比較")}

<p><b>読み取り (重要発見)</b>:</p>
<ul>
<li><b>仮説 H4 (警戒区域内の Sat は外より大きい) は反証された</b>!
   3 エリアとも警戒区域内の Sat 平均は<b>区域外より低い</b>:
   世羅 内 0.068 vs 外 0.075, 熊野 0.050 vs 0.091, 坂 0.090 vs 0.100。</li>
<li>同様に <b>R-B 平均</b>も警戒区域内が低い:
   世羅 内 +13.1 vs 外 +14.7, 熊野 +9.3 vs +17.0, 坂 +15.5 vs +16.7。</li>
<li>これは何を意味するか? 警戒区域は<b>谷底や急斜面の影</b>に多く重なるため、
   <b>R (尾根) より B (谷)</b>が強く出る = R-B 差が小さくなる。
   尾根の凸地形は警戒区域として指定されにくい (= 雨水が集まりにくいため)。</li>
<li>これは<b>土砂災害指定の物理的論理</b> (= 谷頭・急斜面で土石流リスク大) が
   CS立体図の RGB に量的に反映されていることを示す<b>重要な発見</b>。</li>
<li>仮説 H4 は反証されたが、<b>反対方向の関係</b>が量的に確認できたため、
   研究としては<b>新しい知見</b>が得られた (反証も研究的価値あり)。</li>
</ul>

<h3>表 (警戒区域 内/外 の RGB 平均)</h3>
{sed_compare_df.round(2).to_html(classes='', border=0, index=False)}

<p><b>読み取り</b>: 警戒区域内のセル数は世羅町で 14,179, 熊野町で 6,445, 坂町で 2,197。
特別警戒区域は各エリアで 390-2,503 セルと小規模だが、量的比較は十分可能。
<b>diff_RB_mean (R-B 平均)</b>を見ると: 区域内は +9〜+15、区域外は +14〜+17 と一貫して
区域外の方が高い (= 区域外は尾根成分が強い、区域内は谷成分が強い)。</p>
"""
sections.append(("分析 5: 警戒区域 × CS立体図 — ポリゴン×ラスタの空間結合 (本記事の核心)", sec8))


# === Section 9: 分析 6 — 中央サンプルと CS立体図のしくみ ===
sec9 = f"""
<h3>狙い</h3>
<p>各エリアの中央 1024×1024 セルを<b>本来解像度</b>で読込み、CS立体図の細かい構造を体感する。
さらに概念図で「CS立体図がどうやって作られているか」を学習者に理解させる。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p><code>rasterio.windows.Window</code> で各 TIFF の中央 1024×1024 セルだけを読込み、
本来解像度で表示。1m なら 1024m × 1024m の正方領域、50cm なら 512m × 512m の領域。
3 バンド同時読込で RGB 画像として可視化。</p>

<h3>実装</h3>
{code('''
from rasterio.windows import Window

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(window=win).astype(np.uint8)  # (3, H, W)
        # imshow に渡せる (H, W, 3) に転置
        rgb_disp = np.stack([raw[0], raw[1], raw[2]], axis=-1)
        ax.imshow(rgb_disp)
''')}

<h3>図と読み取り (図 8: 中央サンプル)</h3>
<p><b>なぜこの図か</b>: overview 表示では微地形の<b>細部</b>が見えない。
本来解像度で見ると CS立体図の真価 (尾根線・谷線・崩壊跡地が筋状に明瞭) が分かる。</p>

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

<p><b>読み取り</b>:</p>
<ul>
<li><b>世羅町中央 (左)</b>: 1024m × 1024m = 1km 四方の本来解像度。
   <b>樹枝状の谷ネットワーク</b>が淡い青で連続的に走り、その間に薄ピンクの尾根が網目状。
   = 中山間の典型的な地形パターン。</li>
<li><b>熊野町中央 (中)</b>: 同じ 1024m × 1024m。<b>急斜面の暗い影</b>と<b>谷頭の青い凹み</b>が
   明瞭に分かれる。建物 (=四角い影) も視認可能。</li>
<li><b>坂町中央 (右)</b>: 50cm 解像度なので 512m × 512m の狭い範囲だが
   <b>空間解像度は約 4 倍</b>。微細な凹凸構造 (= 古砂防地形 / 平成30崩壊跡 の可能性) が
   1m 版より明確に拾われている。</li>
<li>これが<b>CS立体図が「専門家の地形読み」に有効な理由</b> ⇒ 標高ラスタや傾斜図と異なり、
   <b>凹凸の符号</b>が色で直感的に区別できる。</li>
</ul>

<h3>図と読み取り (図 9: CS立体図のしくみ概念図)</h3>

{figure("assets/L38_fig9_concept.png", "図 9: CS立体図の概念図 — 標高 → 曲率 + 傾斜 → RGB 合成")}

<p><b>読み取り</b>:</p>
<ul>
<li><b>左図</b>: 地形断面 (尾根 + 谷 + 平地) に CS立体図の配色を重ねた模式。
   尾根 = 赤丸, 谷 = 青丸, 平地 = 中間色, 急斜面 = 暗い影 (G 低)。</li>
<li><b>右凡例</b>: 5 種の RGB パターンが対応する地形タイプ。本記事の「地形タイプ自動分類」の根拠。</li>
<li>CS立体図は単一指標 (例: 傾斜のみ, 標高のみ) よりも<b>3 指標を 1 枚に圧縮</b>している点が
   独特。学習者には「色が複数の地形情報を担う」概念の良い例。</li>
</ul>

<h3>表 (中央サンプル統計)</h3>
<table>
<tr><th>エリア</th><th>サイズ (px)</th><th>実距離 (m)</th><th>R 平均</th><th>G 平均</th><th>B 平均</th></tr>
{''.join([f'<tr><td>{a["label"]}</td><td>{a["raw_rgb"].shape[2]}×{a["raw_rgb"].shape[1]}</td><td>{a["raw_rgb"].shape[2]*a["resolution_m"]:.0f}×{a["raw_rgb"].shape[1]*a["resolution_m"]:.0f} m</td><td>{a["raw_rgb"][0].mean():.1f}</td><td>{a["raw_rgb"][1].mean():.1f}</td><td>{a["raw_rgb"][2].mean():.1f}</td></tr>' for a in AREAS])}
</table>

<p><b>読み取り</b>: 中央サンプルの RGB 平均は 3 エリアで類似 (160-200) だが、
<b>サイズ差 (= 実距離差)</b>が解像度差を反映する。坂町は 50cm 解像度のため
同じ 1024×1024 でも 512m 四方しかカバーしない。
これは<b>「ピクセル数 = 情報密度ではない」</b>という重要な GIS リテラシ。
解像度で正規化して比較すべき。</p>
"""
sections.append(("分析 6: 中央サンプルと CS立体図のしくみ概念", sec9))


# === Section 10: 仮説検証と考察 ===
def fmt_verdict(v):
    if v == "支持": return f"<b style='color:#1a7f37;'>{v}</b>"
    elif v == "部分支持": return f"<b style='color:#bf8700;'>{v}</b>"
    elif v == "反証": return f"<b style='color:#cf222e;'>{v}</b>"
    return v

hypos_table = "<table><tr><th>H</th><th>主張</th><th>結果</th><th>判定</th></tr>"
for h in hypos:
    hypos_table += f"<tr><td><b>{h['H']}</b></td><td>{h['claim']}</td><td>{h['result']}</td><td>{fmt_verdict(h['verdict'])}</td></tr>"
hypos_table += "</table>"

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"""
<p>本記事で立てた 6 つの仮説 (H1〜H6) と、3 エリアの CS立体図 RGB 解析 + 警戒区域空間結合の照合結果:</p>

{hypos_table}

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

<h3>主な発見 (3-5 点)</h3>
<ol>
<li><b>H4 (警戒区域内の Sat は外より大きい) は反証された</b>が、その<b>反対方向</b>が
   3 エリア共通で量的に確認された:
   <b>警戒区域は CS立体図で「より低い R-B 差・より低い Saturation」</b>を持つ。
   = 警戒区域は<b>谷底や急斜面の影</b>に多く重なり、尾根の凸地形は警戒区域として
   指定されにくいという物理的論理が CS立体図の RGB に明瞭に刻まれている。
   これは<b>本記事の最も重要な発見</b>であり、
   「土砂災害指定の地形物理学」がオープンデータの 2 dataset 結合で再現できることを示す。</li>

<li><b>H3 (坂町は崩壊跡を含む急傾斜) が強支持された</b>: 坂町の Saturation 平均
   {basic_stats['saka_50cm']['Sat_mean']:.3f} は 3 エリアで最大。
   さらに R/G/B の標準偏差も最大 (R-std=42.0 vs 世羅 15.7)。
   <b>平成30年7月豪雨の被災地形</b>が CS立体図に量的に保存されていることを示す。</li>

<li><b>H2 (熊野町は急峻地形) は支持された</b>が差は微小 (Sat 0.083 vs 世羅 0.074)。
   熊野は花崗岩マサ土地形で急傾斜面が多いと予想したが、実データ上は世羅との
   差は限定的。これは<b>「都市近郊の山地」の急峻度が中山間の丘陵地形と大差ない</b>
   という意外な発見 (発展課題 Z3 で再検証)。</li>

<li><b>H5 (色相は 2 峰性) が支持された</b>: 全エリアで Hue が
   赤帯 (0.0) と青帯 (0.66) の双方に分布。Hue 標準偏差は<b>世羅 0.057 < 熊野 0.136 < 坂 0.252</b>と
   段階的に増加 = エリアの地形多様度を量的に区別する指標として機能。</li>

<li><b>H6 (50cm 版は微細構造強) が支持された</b>: 坂町 (50cm) の Sat 標準偏差
   {basic_stats['saka_50cm']['Sat_std']:.3f} は 1m 版 (世羅 0.043, 熊野 0.060) より大きい。
   高解像度ほど微地形のばらつきを拾うことを実証。
   ただし熊野 (1m) も同等の Sat std を持つため、<b>解像度差より地形類型差の方が
   Sat std に効く</b>可能性も残る。</li>
</ol>

<h3>3 エリアの地形特性比較 (本記事の量的サマリ)</h3>
<table>
<tr><th>軸</th><th>世羅町 (1m)</th><th>熊野町 (1m)</th><th>坂町 (50cm)</th></tr>
<tr><td>地形類型</td><td>緩傾斜丘陵 (中山間)</td><td>急峻山地 + 都市平地</td><td>急斜面 + 崩壊跡</td></tr>
<tr><td>R 平均</td><td>{basic_stats['sera_1m']['R_mean']:.1f}</td><td>{basic_stats['kumano_1m']['R_mean']:.1f}</td><td>{basic_stats['saka_50cm']['R_mean']:.1f}</td></tr>
<tr><td>R-B 平均</td><td>{basic_stats['sera_1m']['diff_RB_mean']:.1f}</td><td>{basic_stats['kumano_1m']['diff_RB_mean']:.1f}</td><td>{basic_stats['saka_50cm']['diff_RB_mean']:.1f}</td></tr>
<tr><td>Sat 平均</td><td>{basic_stats['sera_1m']['Sat_mean']:.3f}</td><td>{basic_stats['kumano_1m']['Sat_mean']:.3f}</td><td>{basic_stats['saka_50cm']['Sat_mean']:.3f}</td></tr>
<tr><td>Sat 標準偏差</td><td>{basic_stats['sera_1m']['Sat_std']:.3f}</td><td>{basic_stats['kumano_1m']['Sat_std']:.3f}</td><td>{basic_stats['saka_50cm']['Sat_std']:.3f}</td></tr>
<tr><td>Hue 標準偏差</td><td>{AREAS[0].get('hue_std', 0):.3f}</td><td>{AREAS[1].get('hue_std', 0):.3f}</td><td>{AREAS[2].get('hue_std', 0):.3f}</td></tr>
<tr><td>警戒区域率</td><td>{AREAS[0]['sed_yellow_km2']/AREAS[0]['muni_area_km2']*100:.1f}%</td>
   <td>{AREAS[1]['sed_yellow_km2']/AREAS[1]['muni_area_km2']*100:.1f}%</td>
   <td>{AREAS[2]['sed_yellow_km2']/AREAS[2]['muni_area_km2']*100:.1f}%</td></tr>
<tr><td>特別警戒区域率</td><td>{AREAS[0]['sed_red_km2']/AREAS[0]['muni_area_km2']*100:.1f}%</td>
   <td>{AREAS[1]['sed_red_km2']/AREAS[1]['muni_area_km2']*100:.1f}%</td>
   <td>{AREAS[2]['sed_red_km2']/AREAS[2]['muni_area_km2']*100:.1f}%</td></tr>
</table>

<h3>考察</h3>
<p>本記事の主たる発見は、<b>「警戒区域の指定論理が CS立体図の RGB に量的に刻まれている」</b>
ことが 3 エリア共通で確認されたことである。仮説 H4 は反証されたが、
反対方向の関係 (警戒区域は谷底寄り = R-B 平均が低い) が一貫して観察された。
これは<b>「土砂災害警戒区域 = 谷頭 / 急斜面の影」</b>という地形物理的論理の量的証拠であり、
オープンデータの 2 dataset 結合だけで再現可能なことを示している。</p>

<p>同時に、<b>解像度差 (1m vs 50cm)</b>は overview レベルでは認識困難だが、
本来解像度のサンプル (図 8) では 50cm の優位性が見える。微細な凹凸 (古砂防地形・
古崩壊跡など) は 50cm 版でしか拾えない可能性が高く、<b>地形リテラシ研究には 50cm 版が必須</b>。</p>

<p>3 エリアの<b>絶対値の RGB スケール</b>には大差はないが、<b>標準偏差・彩度</b>には明確な
階層 (世羅 < 熊野 ≦ 坂) が見える。これは「地形微細構造の強さ」を CS立体図 RGB から
量的に評価する手法として有効。<b>仮説を立てる前にデータの値域を観察する</b>重要性を示す事例。</p>

<p>パーセンタイル適応閾値での地形タイプ分類 (尾根/谷/急/平地/複雑) は<b>各エリア内の相対分類</b>として
意味があるが、絶対値での横断比較には不向き。<b>絶対閾値 vs 相対閾値の使い分け</b>は
リテラシ上の重要論点 (発展課題 Z2)。</p>
"""
sections.append(("仮説検証と考察", sec_h))


# === Section 11: 発展課題 ===
sec_dev = f"""
<h3>結果 X1 → 仮説 Y1 → 課題 Z1: 警戒区域 ⇔ 谷地形の量的検証 (本記事の発見の深掘り)</h3>
<p><b>結果 X1</b>: 3 エリア共通で警戒区域の R-B 平均は区域外より低い (= 谷成分強)。
H4 を反対方向で支持する量的傾向。</p>
<p><b>新仮説 Y1</b>: 警戒区域内の<b>「谷タイプ (B>R)」セルの比率</b>は区域外より明確に高い。
50%以上の警戒域が谷タイプに分類される。</p>
<p><b>課題 Z1</b>: 警戒区域内の地形タイプ構成比 (尾根/谷/急/平地) を本記事のラスタ化マスクで集計し、
区域外と並列比較。クロス表で「警戒域は谷タイプが何 % か」を量化。
さらに土石流・急傾斜・地すべりの 3 種別ごとに比較すると、各災害タイプの
地形特性が CS立体図でどう違うか分かる。</p>

<h3>結果 X2 → 仮説 Y2 → 課題 Z2: 標高ラスタからの直接的な曲率計算</h3>
<p><b>結果 X2</b>: 本記事の「地形タイプ自動分類」は CS立体図の RGB から逆推定したものだが、
適応閾値が必要なため絶対比較が困難。</p>
<p><b>新仮説 Y2</b>: <b>標高ラスタ (DTM, dsid 1635/1636) から直接</b> 曲率・傾斜を計算した方が
量的に精確で、絶対閾値での横断比較が可能。</p>
<p><b>課題 Z2</b>: 標高 GeoTIFF を読み込んで <code>scipy.ndimage.laplace</code> で曲率を、
<code>numpy.gradient</code> で傾斜を直接計算。本記事の RGB 推定と相関係数で照合。
高い相関が出れば<b>「CS立体図の配色は曲率・傾斜の忠実な可視化」</b>と確認できる。</p>

<h3>結果 X3 → 仮説 Y3 → 課題 Z3: 多自治体一括の地形指標化</h3>
<p><b>結果 X3</b>: 本記事は 3 自治体のみ。Sat 平均 (= 微地形コントラスト) で
世羅 < 熊野 < 坂 の階層が見えたが、<b>サンプル数が少なすぎて</b>地形類型と Sat の関係が
仮説止まり。</p>
<p><b>新仮説 Y3</b>: 全 22 自治体 (1m 版) で Sat 平均を計算すると、
中山間 / 都市近郊 / 沿岸被災地 / 島嶼 などの<b>地形類型クラスタ</b>が自然分類できる。</p>
<p><b>課題 Z3</b>: 22 自治体すべての CS立体図を順次 DL → overview 読込 → Sat 平均計算で
22 値の分布を作る。k-means や階層クラスタリングで類型化。
さらに各自治体の警戒区域率と Sat 平均を散布図にプロット ⇒
<b>「Sat 高 = 警戒区域率高」の正相関</b>が出れば
「CS立体図の Sat は土砂災害リスクの代理指標」という新提案ができる。</p>

<h3>結果 X4 → 仮説 Y4 → 課題 Z4: 古崩壊地形の自動検出</h3>
<p><b>結果 X4</b>: 坂町 (H30豪雨被災地) は 3 エリアで最高の Sat 標準偏差を示した。
微細な凹凸 (= 過去の崩壊地形 / 古砂防地形) が CS立体図に保存されていると示唆。</p>
<p><b>新仮説 Y4</b>: <b>「Sat の局所変動が大きい場所」</b>は古崩壊地形の可能性が高い。
過去の航空写真や砂防履歴と照合可能。</p>
<p><b>課題 Z4</b>: 各セルの周囲 5×5 セルでの Sat 局所標準偏差を計算し、
高い局所 std セルを抽出。これが既知の H30 崩壊跡 (地理院地図の被災地マップ) と
空間的に一致するかを精度評価。一致が高ければ<b>「Sat 局所 std による古崩壊地形自動検出」</b>を提案できる。
これは砂防研究・防災マップ更新に直接応用可能。</p>

<h3>結果 X5 → 仮説 Y5 → 課題 Z5: L36/L37 との 3 軸結合 (LiDAR 派生物の三脚)</h3>
<p><b>結果 X5</b>: 本記事は CS立体図のみ単独扱い (合体禁止)。
L36 樹高 / L37 点群密度との関係は未検証。</p>
<p><b>新仮説 Y5</b>: 同じ航空レーザ点群から派生した 3 ラスタ (樹高 / 点群密度 / CS立体図) は
<b>3 軸の地形・植生指標</b>として機能。
3 軸の組合せで「裸地崩壊」「森林急斜面」「都市平地」等の 9-12 クラス自動分類が可能。</p>
<p><b>課題 Z5</b>: X 系の横断記事として、3 ラスタを同一座標系で空間結合し、
GMM (Gaussian Mixture Model) で 9-12 クラスタリング。各クラスタが
何の地形/植生に対応するかを衛星 NDVI や建物 footprint で検証。
これは LiDAR ファミリ全体を統合する集大成記事になる。</p>

<h3>結果 X6 → 仮説 Y6 → 課題 Z6: 過去被災地のセル単位検出</h3>
<p><b>結果 X6</b>: 坂町の特別警戒区域 ({AREAS[2]['sed_red_km2']:.2f} km²) と CS立体図の急斜面 (G 低) は
重なるが、過去の被災実績との照合は未実施。</p>
<p><b>新仮説 Y6</b>: 平成30年7月豪雨で<b>実際に崩壊が発生したセル</b>は CS立体図で
特定の RGB 特徴 (Sat 中, 谷タイプ, 急傾斜) を持つ。RGB 特徴量で<b>機械学習による被災予測</b>が可能。</p>
<p><b>課題 Z6</b>: 国交省の H30 豪雨被災調査データ (発生位置 GIS) を入手し、
CS立体図の RGB 特徴で<b>ロジスティック回帰や Random Forest</b> による被災予測モデルを訓練。
特徴量重要度から「どのバンドが被災予測に効くか」を判定。
高い AUC が出れば<b>CS立体図ベースの土砂災害発生予測モデル</b>として実装可能。</p>
"""
sections.append(("発展課題", sec_dev))


# === HTML 書き出し ===
html = render_lesson(
    num=38,
    title="L38 CS立体図 2件統合分析 — 曲率×傾斜の RGB 段彩が描く地形微細構造 と 土砂災害警戒区域の物理的論理",
    tags=["LiDAR", "CS立体図", "地形微細構造", "RGB ラスタ", "rasterio", "曲率", "傾斜", "土砂災害警戒区域", "rasterize", "HSV"],
    time="20-30 分",
    level="リテラシ + GIS応用",
    data_label="DoBoX dsid 1628 + 1629 (CS立体図 1m + 50cm), 代表 3 自治体 GeoTIFF (RGB 3 バンド)",
    sections=sections,
    script_filename="lessons/L38_cs_terrain.py",
)
(LESSONS / "L38_cs_terrain.html").write_text(html, encoding="utf-8")
print(f"  saved L38_cs_terrain.html ({time.time()-t1:.1f}s)", flush=True)
print(f"\n=== Final elapsed: {time.time()-t0:.1f}s ===", flush=True)