# -*- coding: utf-8 -*-
"""L36 樹高図 2 件統合分析
       — 広島県の航空レーザ計測 DCHM が描く 2 エリアの林相と森林構造

カバー宣言:
  本記事は DoBoX のシリーズ <b>「樹高図 1m / 樹高図 50cm」 2 件</b>
  (dataset_id = 1626, 1627) を統合し、広島県内 2 自治体エリア
  (世羅町=中山間/林業地, 熊野町=都市近郊/呉市隣接) の<b>樹高分布</b>から
  <b>林相 (針葉樹 / 広葉樹 / 伐採地 / 都市・農地・水域)</b> を推定する研究記事である。

  L32-L35 (港湾・河川シリーズ) は<b>線・点</b>のインフラ施設を扱った。
  L36 から始まる<b>LiDAR ファミリ</b>は<b>面 (ラスタ)</b>のフィジカル測定を扱う。
  扱うデータは航空レーザ計測 (LiDAR) 由来の<b>樹冠高さラスタ (DCHM)</b> = DSM (地表全体) - DTM (地面のみ)。
  1 ピクセルあたり 1m もしくは 0.5m のメッシュで、各セルにその位置の<b>地物高さ (主に樹木)</b> が
  cm スケールで記録される。

  本記事は LiDAR ファミリ (CS立体図/地形傾斜図/標高図/航空レーザ計測森林資源/地面到達レーザ点群密度図/
  地図情報/計測年度図) の最初の 1 本であり、他の LiDAR 系記事とは<b>合体しない厳密シリーズ単独記事</b>。

研究の問い (RQ):
  広島県内の樹高図 2 解像度 (1m / 50cm) は、それぞれ
  <b>世羅町 (中山間, 林業地)</b> と <b>熊野町 (都市近郊, 山岳混在)</b> という
  <b>異なる地域特性</b> のエリアを公開している。
  両エリアの<b>樹高ヒストグラム・空間分布・林相推定</b>を比較すると、
  どんな地域構造の違いが見えてくるか? また<b>解像度差 (1m vs 50cm)</b> は
  どの程度の影響を持つか?
    (1) 各エリアの<b>面積・形状・有効データ率</b>の比較
    (2) 樹高ヒストグラム (0-30m を 1m 階級) のピーク位置と多峰性
    (3) 高木 (>15m) ・ 中木 (5-15m) ・低木 (1-5m) ・地表 (<1m) の構成比
    (4) 空間的な高木地・低木地の分布 (集塊性)
    (5) 林相推定: 高木優占 = 成熟林、中低木混在 = 若齢林・伐採跡、ゼロ卓越 = 都市/農地/水域
    (6) 解像度差: 50cm はノイズ・精度がどう違うか

仮説 H1〜H6:
  H1 (世羅町は林業地優占): 平均樹高 ≥ 8m、15m 以上のセル比率 ≥ 15%。
       中山間で植林スギ・ヒノキ針葉樹が広面積を占めると推定。
  H2 (熊野町は都市・農地優占): 樹高 0-1m セル比率 ≥ 40%。
       都市部の建物外周・農地・道路・山頂部のはげ地が高い割合。
  H3 (両エリアともゼロ卓越): 樹高 0-1m が最頻ビン (mode)。
       航空 LiDAR で多数のセルが地表 (建物・農地) もしくは雲・霧除去後ノイズで 0 付近。
  H4 (世羅町の樹高分布は 2 峰性): 0-1m と 10-20m の 2 ピーク。
       林業地 = 「植林伐採跡 + 成熟林」の二極構造。
  H5 (熊野町の高木は山岳に集中): 標高高い側の山地に高木 (>15m) が集中、
       平地は低木・地表優占。空間相関で確認。
  H6 (50cm 版はノイズ多): 50cm の標準偏差は 1m より大、負値出現、極端値 (>30m) も多。
       高解像度ほど雲・霧の残存ノイズや単木個葉のばらつきを拾う。

要件 S 準拠: 1 分以内完走。overview /128 のダウンサンプル読み込みで MB 級メモリ。
要件 T 準拠: 各エリア地図 + 樹高階級マップ (small multiples) + 主題図 + 散布。
要件 Q 準拠: 図 7-9 種, 表 8+ (2 エリア × 4 軸 (面積/解像度/分布/階級))。

データ仕様:
  A. 樹高図 1m  dataset 1626: 県内 23 自治体に GeoTIFF 公開。本記事は<b>世羅町</b> (1.7 GB tif)
     - 27,248 × 16,357 セル, EPSG:6671, 1m 解像度, float32, nodata=-99999
     - 範囲: 約 27 km × 16 km
  B. 樹高図 50cm dataset 1627: 県内 23 自治体に GeoTIFF 公開。本記事は<b>熊野町</b> (1.0 GB tif)
     - 15,138 × 17,230 セル, EPSG:6671, 0.5m 解像度, float32, nodata=-99999
     - 範囲: 約 7.5 km × 8.6 km
  両者とも DSM - DTM = 地物高 (200m 以上は雲霧として除去)
  各 ZIP は 600MB-870MB と巨大なので、本記事は<b>各 1 自治体のみ DL</b>し、
  rasterio の<b>ピラミッド (overview)</b> を使って高速に概観する。
"""
from __future__ import annotations
import sys, time, json, os
from pathlib import Path

sys.path.insert(0, str(Path(__file__).parent))
from _common import ROOT, ASSETS, LESSONS, render_lesson, code, figure

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from matplotlib.colors import ListedColormap, BoundaryNorm
import matplotlib.patches as mpatches
import geopandas as gpd
from shapely.geometry import Polygon, box
import rasterio
from rasterio.windows import Window
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("=== L36 樹高図 2 件統合分析 ===", flush=True)

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

# 各エリア定義
AREAS = [
    {
        "key": "sera_1m",
        "label": "世羅町",
        "muni_code": 34462,
        "resolution_m": 1.0,
        "dataset_id": 1626,
        "dataset_name": "樹高図（1m）",
        "zip_url": "https://data.hiroshima-dobox.jp/forest/th_1m/34462_世羅町.zip",
        "zip_local": DATA_DIR / "1m_34462_sera.zip",
        "tif": DATA_DIR / "1m_sera" / "sera_1m_treeheight.tif",
        "admin_zip": ADMIN_DIR / "admin_941_世羅町.zip",
        "color": "#1a7f37",       # 林業地=深緑
        "character": "中山間 / 林業地 (植林スギ・ヒノキ優占想定)",
    },
    {
        "key": "kumano_50cm",
        "label": "熊野町",
        "muni_code": 34307,
        "resolution_m": 0.5,
        "dataset_id": 1627,
        "dataset_name": "樹高図（50cm）",
        "zip_url": "https://data.hiroshima-dobox.jp/forest/th_50cm/34307_熊野町.zip",
        "zip_local": DATA_DIR / "50cm_34307_kumano.zip",
        "tif": DATA_DIR / "50cm_kumano" / "kumano_50cm_treeheight.tif",
        "admin_zip": ADMIN_DIR / "admin_911_熊野町.zip",
        "color": "#0969da",       # 都市近郊=青
        "character": "都市近郊 / 山岳混在 (筆造り産業地、呉市隣接)",
    },
]

# 樹高階級 (本記事独自定義)
HEIGHT_CLASSES = [
    ("地表 / nodata", -1, 1, "#e8e8e8"),         # ≤ 1m
    ("低木 / 草地",   1,  3, "#fff7bc"),         # 1-3m
    ("低中木",        3,  6, "#fec44f"),         # 3-6m
    ("中木",          6, 10, "#a1d99b"),         # 6-10m
    ("中高木",       10, 15, "#41ab5d"),         # 10-15m
    ("高木",         15, 20, "#005a32"),         # 15-20m
    ("巨木 / 例外",  20, 200, "#54278f"),         # >= 20m
]


def ensure_zip(area):
    """ZIP が無ければ DoBoX から DL"""
    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)"}
    print(f"  DL {p.name} ← {area['zip_url'][:80]}...", flush=True)
    with requests.get(area["zip_url"], headers=UA, timeout=600, 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=1024*1024):
                if chunk:
                    f.write(chunk)
                    downloaded += len(chunk)
                    if total > 0:
                        pct = int(downloaded / total * 100)
                        if pct >= last_pct + 10:
                            print(f"    {pct}% ({downloaded/(1024*1024):.0f}/{total/(1024*1024):.0f} MB)", flush=True)
                            last_pct = pct
    return p


def ensure_tif(area):
    """ZIP から TIFF を抽出 (シンプルな英語名で保存)"""
    tif = area["tif"]
    if tif.exists() and tif.stat().st_size > 100_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 n in zf.namelist():
            if "__MACOSX" in n: continue
            base = Path(n).name
            if base.startswith("._"): continue
            # 大きい TIFF と OVR を抽出
            inf = zf.getinfo(n)
            size = inf.file_size
            target_name = None
            if n.endswith(".tiff") and size > 100_000_000:
                target_name = tif.name
            elif n.endswith(".tiff.ovr"):
                target_name = tif.name + ".ovr"
            if target_name:
                outpath = target_dir / target_name
                if outpath.exists() and outpath.stat().st_size > 1000:
                    continue
                print(f"    extracting {target_name} ({size/(1024*1024):.0f} MB)...", flush=True)
                with zf.open(n) as src, open(outpath, "wb") as dst:
                    while True:
                        chunk = src.read(4*1024*1024)
                        if not chunk: break
                        dst.write(chunk)
    return tif


# =============================================================================
# 1. データ確保: 2 自治体の 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} ({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,
            "crs": str(ds.crs),
            "res_m": float(ds.res[0]),
            "left": bb.left, "right": bb.right,
            "bottom": bb.bottom, "top": bb.top,
            "width_m": bb.right - bb.left,
            "height_m": bb.top - bb.bottom,
            "area_km2": (bb.right - bb.left) * (bb.top - bb.bottom) / 1e6,
            "nodata": ds.nodata,
            "dtype": str(ds.dtypes[0]),
            "overviews": ds.overviews(1),
        }

metas = []
for a in AREAS:
    m = tif_meta(a["tif"])
    a["meta"] = m
    metas.append({"key": a["key"], "label": a["label"], **m})
    print(f"  {a['label']}: {m['width']:,}×{m['height']:,} = {m['n_cells']/1e6:.1f}M cells, "
          f"{m['width_m']/1000:.1f}×{m['height_m']/1000:.1f} km, "
          f"bbox area = {m['area_km2']:.1f} km²", flush=True)

meta_df = pd.DataFrame(metas)
meta_df.to_csv(ASSETS / "L36_tiff_meta.csv", index=False, encoding="utf-8-sig")
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 3. overview を使って全体ラスタを読み込み (MB 級メモリ)
# =============================================================================
print("\n[3] overview (ピラミッド) で全体読み込み", flush=True)
t1 = time.time()

# 同じ目標出力解像度 (約 200m/cell) になるように overview を選択
# Sera: 1m * 128 = 128m/cell -> /128 を使う
# Kumano: 0.5m * 64 = 32m/cell -> /64 を使う(より粗く)
# 比較しやすいよう、ターゲットセルサイズ ≈ 50m
# Sera は /64 で 64m/cell, Kumano は /127 で 63.5m/cell ≈ 同等
TARGET_CELL_M = 64

for a in AREAS:
    with rasterio.open(a["tif"]) as ds:
        ovr_choices = ds.overviews(1)
        # 目標ピクセルサイズに最も近い overview レベル
        # eff_cell = res_m * factor
        # factor ≈ TARGET / res_m
        target_factor = TARGET_CELL_M / a["resolution_m"]
        # ovr_choices から最近接 (ただし target_factor より大きいか同等なら採用)
        best = ovr_choices[-1]
        for f in ovr_choices:
            if f >= target_factor:
                best = f; break
        ov_w = ds.width // best
        ov_h = ds.height // best
        print(f"  {a['label']}: overview /1/{best} -> {ov_w}×{ov_h} ({ov_w*ov_h/1e6:.2f}M)", flush=True)
        # 読み込み (average resampling = ピクセルの平均で粗化)
        arr = ds.read(1, out_shape=(ov_h, ov_w),
                      resampling=Resampling.average)
        a["arr"] = arr
        a["arr_factor"] = best
        a["arr_cell_m"] = a["resolution_m"] * best  # 実効セルサイズ (m)
        a["arr_w"] = ov_w
        a["arr_h"] = ov_h
        # 範囲 (ボーダーは TIFF と同じ)
        a["bounds"] = ds.bounds
        # nodata マスク
        nodata = ds.nodata if ds.nodata is not None else -99999.0
        # 浮動小数なので NaN もチェック
        valid_mask = ~np.isnan(arr) & (arr != nodata) & (arr > -1000)
        a["valid_mask"] = valid_mask
        a["valid"] = arr[valid_mask].astype(np.float32)
        n_valid = int(valid_mask.sum())
        print(f"    valid={n_valid:,}/{ov_w*ov_h:,} ({n_valid/(ov_w*ov_h)*100:.1f}%), "
              f"min={a['valid'].min():.2f}, max={a['valid'].max():.2f}, "
              f"mean={a['valid'].mean():.2f}, median={np.median(a['valid']):.2f}", flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 4. ヒストグラム & 階級別集計
# =============================================================================
print("\n[4] ヒストグラム & 階級集計", flush=True)
t1 = time.time()

bins_1m = np.arange(0, 31, 1)  # 1m 刻みヒストグラム
bins_class = [-1000, 1, 3, 6, 10, 15, 20, 200]
class_labels = ["<1m 地表", "1-3m 低木", "3-6m 低中木", "6-10m 中木",
                "10-15m 中高木", "15-20m 高木", ">=20m 巨木"]

class_counts = {}
hist_1m = {}
basic_stats = {}
for a in AREAS:
    valid = a["valid"]
    valid = np.clip(valid, -1, 100)  # 極端値クリップ
    # 1m 階級
    h, _ = np.histogram(valid, bins=bins_1m)
    hist_1m[a["key"]] = h
    # 分類別
    cc, _ = np.histogram(valid, bins=bins_class)
    class_counts[a["key"]] = cc
    # 基本統計
    basic_stats[a["key"]] = {
        "label": a["label"],
        "n_cells_total": int(a["valid_mask"].size),
        "n_cells_valid": int(a["valid_mask"].sum()),
        "valid_ratio": float(a["valid_mask"].sum() / a["valid_mask"].size),
        "min": float(valid.min()),
        "max": float(valid.max()),
        "mean": float(valid.mean()),
        "median": float(np.median(valid)),
        "std": float(valid.std()),
        "p10": float(np.percentile(valid, 10)),
        "p25": float(np.percentile(valid, 25)),
        "p75": float(np.percentile(valid, 75)),
        "p90": float(np.percentile(valid, 90)),
        "p99": float(np.percentile(valid, 99)),
        "n_zero": int(((valid >= 0) & (valid < 1)).sum()),
        "n_low": int(((valid >= 1) & (valid < 5)).sum()),
        "n_mid": int(((valid >= 5) & (valid < 15)).sum()),
        "n_high": int((valid >= 15).sum()),
    }

stats_df = pd.DataFrame.from_dict(basic_stats, orient="index")
stats_df.to_csv(ASSETS / "L36_basic_stats.csv", encoding="utf-8-sig")
print(stats_df.to_string())

# 階級比率 (%)
class_pct_rows = []
for a in AREAS:
    cc = class_counts[a["key"]]
    total = cc.sum()
    pct = cc / total * 100 if total > 0 else cc * 0.0
    row = {"label": a["label"], "total_valid": int(total)}
    for cl, c, p in zip(class_labels, cc, pct):
        row[f"{cl} %"] = round(float(p), 2)
        row[f"{cl} n"] = int(c)
    class_pct_rows.append(row)
class_pct_df = pd.DataFrame(class_pct_rows)
class_pct_df.to_csv(ASSETS / "L36_height_classes.csv", index=False, encoding="utf-8-sig")
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 5. 林相推定 (本記事独自ヒューリスティック)
# =============================================================================
print("\n[5] 林相推定 (heuristic)", flush=True)
t1 = time.time()

# 簡易ルール:
#   - 高木 >= 30%  : 成熟林優占
#   - 0-1m >= 50%  : 都市・農地・水域優占
#   - 中高木 (10-15m) >= 30% かつ低木 (<5m) も多い : 若齢林・植林手入れ中
#   - 中高木〜高木 が 20-50% : 林相混在 (typical 中山間)

def estimate_lin_so(stats):
    pct_high = stats["n_high"] / max(1, stats["n_cells_valid"]) * 100
    pct_zero = stats["n_zero"] / max(1, stats["n_cells_valid"]) * 100
    pct_mid = stats["n_mid"] / max(1, stats["n_cells_valid"]) * 100
    pct_low = stats["n_low"] / max(1, stats["n_cells_valid"]) * 100
    if pct_zero >= 50:
        return "都市・農地・水域優占 (低樹高)", pct_high, pct_zero
    if pct_high >= 30:
        return "成熟林優占 (高木卓越)", pct_high, pct_zero
    if pct_mid >= 30 and pct_low >= 15:
        return "若齢林・植林手入れ中 (中木+低木)", pct_high, pct_zero
    return "林相混在 (中山間典型)", pct_high, pct_zero

linsos = []
for a in AREAS:
    label, pct_h, pct_z = estimate_lin_so(basic_stats[a["key"]])
    linsos.append({
        "label": a["label"],
        "linso_estimate": label,
        "pct_high (>=15m)": round(pct_h, 1),
        "pct_zero (<1m)": round(pct_z, 1),
    })
    print(f"  {a['label']}: {label} (高木 {pct_h:.1f}%, 地表 {pct_z:.1f}%)", flush=True)

linso_df = pd.DataFrame(linsos)
linso_df.to_csv(ASSETS / "L36_linso_estimate.csv", index=False, encoding="utf-8-sig")
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 6. 行政界ポリゴン (admin) を読み込んで TIFF とのオーバラップを確認
# =============================================================================
print("\n[6] admin polygon ロード", flush=True)
t1 = time.time()


def load_admin_zip(zip_path):
    """admin shp zip からポリゴン GeoDataFrame を返す"""
    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():
        g = load_admin_zip(a["admin_zip"]).to_crs(a["meta"]["crs"])
        a["admin"] = g.dissolve()
        bb = a["admin"].total_bounds
        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} km²  (TIFF bbox area {a['meta']['area_km2']:.1f} km²)", flush=True)
    else:
        print(f"  {a['label']}: admin zip 不在")
        a["admin"] = None
        a["muni_area_km2"] = None
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 6b. 生解像度サンプル: 各エリアの中央 1024x1024 セルを本来解像度で読込
# overview の平均化で 0-1m 比率が下がるバイアスを補正するため
# =============================================================================
print("\n[6b] 生解像度サンプル (中央 1024x1024)", flush=True)
t1 = time.time()
RAW_SAMPLE = 1024
for a in AREAS:
    with rasterio.open(a["tif"]) as ds:
        cx, cy = ds.width//2, ds.height//2
        win = Window(max(0, cx - RAW_SAMPLE//2),
                     max(0, cy - RAW_SAMPLE//2),
                     min(RAW_SAMPLE, ds.width),
                     min(RAW_SAMPLE, ds.height))
        raw = ds.read(1, window=win).astype(np.float32)
    # 有効値
    nodata = ds.nodata if ds.nodata is not None else -99999.0
    valid = raw[(~np.isnan(raw)) & (raw != nodata) & (raw > -1000)]
    n_valid = len(valid)
    pct_zero_raw = (valid < 1).sum() / max(1, n_valid) * 100
    pct_low_raw = ((valid >= 1) & (valid < 5)).sum() / max(1, n_valid) * 100
    pct_mid_raw = ((valid >= 5) & (valid < 15)).sum() / max(1, n_valid) * 100
    pct_high_raw = (valid >= 15).sum() / max(1, n_valid) * 100
    a["raw_sample_n"] = n_valid
    a["raw_sample_pct_zero"] = pct_zero_raw
    a["raw_sample_pct_low"] = pct_low_raw
    a["raw_sample_pct_mid"] = pct_mid_raw
    a["raw_sample_pct_high"] = pct_high_raw
    a["raw_sample_mean"] = float(valid.mean()) if n_valid > 0 else 0
    a["raw_sample_std"] = float(valid.std()) if n_valid > 0 else 0
    print(f"  {a['label']} 中央サンプル: n={n_valid:,}, 平均 {a['raw_sample_mean']:.2f}m, std {a['raw_sample_std']:.2f}, "
          f"0-1m={pct_zero_raw:.1f}%, 1-5m={pct_low_raw:.1f}%, 5-15m={pct_mid_raw:.1f}%, >=15m={pct_high_raw:.1f}%", flush=True)

raw_sample_df = pd.DataFrame([{
    "label": a["label"],
    "resolution_m": a["resolution_m"],
    "n_valid_cells": a["raw_sample_n"],
    "pct_0_1m": round(a["raw_sample_pct_zero"], 2),
    "pct_1_5m": round(a["raw_sample_pct_low"], 2),
    "pct_5_15m": round(a["raw_sample_pct_mid"], 2),
    "pct_ge_15m": round(a["raw_sample_pct_high"], 2),
    "mean_m": round(a["raw_sample_mean"], 2),
    "std_m": round(a["raw_sample_std"], 2),
} for a in AREAS])
raw_sample_df.to_csv(ASSETS / "L36_raw_sample_stats.csv", index=False, encoding="utf-8-sig")
print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 7. 仮説検証
# =============================================================================
print("\n[7] 仮説検証", flush=True)
t1 = time.time()

s_sera = basic_stats["sera_1m"]
s_kuma = basic_stats["kumano_50cm"]

hypos = []

# H1: 世羅町は林業地優占 (平均 >= 8m, 高木 >= 15%)
sera_pct_high = s_sera["n_high"] / s_sera["n_cells_valid"] * 100
h1 = (s_sera["mean"] >= 8.0) and (sera_pct_high >= 15.0)
hypos.append({
    "H": "H1",
    "claim": "世羅町は林業地優占で平均樹高 >= 8m かつ 高木 (>=15m) 比率 >= 15%",
    "result": f"平均樹高 {s_sera['mean']:.2f} m, 高木比率 {sera_pct_high:.1f}%",
    "verdict": "支持" if h1 else ("部分支持" if (s_sera['mean'] >= 6 or sera_pct_high >= 10) else "反証"),
})

# H2: 熊野町は地表優占 (0-1m >= 40%)
# overview の平均化バイアスを除外するため、生解像度の中央サンプルで判定
kuma_pct_zero_ovr = s_kuma["n_zero"] / s_kuma["n_cells_valid"] * 100
kuma_pct_zero_raw = AREAS[1]["raw_sample_pct_zero"]
h2 = (kuma_pct_zero_raw >= 40.0)
hypos.append({
    "H": "H2",
    "claim": "熊野町は都市・農地優占で 0-1m セル比率 >= 40% (生解像度中央サンプルで判定)",
    "result": f"生解像度サンプル 0-1m 比率 {kuma_pct_zero_raw:.1f}% (overview 全域では {kuma_pct_zero_ovr:.1f}% — 平均化で過小評価)",
    "verdict": "支持" if h2 else ("部分支持" if kuma_pct_zero_raw >= 25 else "反証"),
})

# H3: 両エリアとも 0-1m が最頻ビン
def find_mode_bin(hist):
    return int(np.argmax(hist))
mode_sera = find_mode_bin(hist_1m["sera_1m"])
mode_kuma = find_mode_bin(hist_1m["kumano_50cm"])
h3 = (mode_sera == 0) and (mode_kuma == 0)
hypos.append({
    "H": "H3",
    "claim": "両エリアとも 0-1m ビンが最頻 (mode)",
    "result": f"世羅 mode={mode_sera}m, 熊野 mode={mode_kuma}m",
    "verdict": "支持" if h3 else ("部分支持" if (mode_sera <= 1 or mode_kuma <= 1) else "反証"),
})

# H4: 世羅町の樹高分布は 2 峰性 (0-1m + 10-20m)
# ピーク: bin 0 と bin 10-20 の中で高い側を見る
peak_low = hist_1m["sera_1m"][0]
peak_high_zone = hist_1m["sera_1m"][10:20].max()
peak_mid_zone_min = hist_1m["sera_1m"][2:10].min()  # 谷
two_peaks = (peak_high_zone >= peak_low * 0.3) and (peak_mid_zone_min < peak_high_zone * 0.7)
hypos.append({
    "H": "H4",
    "claim": "世羅町の樹高ヒストグラムは 0-1m と 10-20m の 2 峰性",
    "result": f"低ピーク (0-1m): {peak_low}, 高ピーク (10-20m帯 max): {peak_high_zone}, 中間 (2-10m) 最小: {peak_mid_zone_min}",
    "verdict": "支持" if two_peaks else "部分支持",
})

# H5: 熊野町の高木 (>=15m) は北西山岳? 空間相関
# 高木マスクの空間中心と低木マスクの空間中心の差
def spatial_centroid(arr, mask):
    ys, xs = np.where(mask)
    if len(ys) == 0: return None, None
    return float(ys.mean()), float(xs.mean())

# 両エリアで重心を計算して AREAS に保存
for a in AREAS:
    high_mask = (a["arr"] >= 15) & a["valid_mask"]
    low_mask = ((a["arr"] < 5) & (a["arr"] >= 0)) & a["valid_mask"]
    ch = spatial_centroid(a["arr"], high_mask)
    cl = spatial_centroid(a["arr"], low_mask)
    if ch[0] is not None and cl[0] is not None:
        # 画像 row,col → 世界座標 (m)
        a["h_centroid_y"] = a["bounds"].top - (ch[0] + 0.5) * a["arr_cell_m"]
        a["h_centroid_x"] = a["bounds"].left + (ch[1] + 0.5) * a["arr_cell_m"]
        a["l_centroid_y"] = a["bounds"].top - (cl[0] + 0.5) * a["arr_cell_m"]
        a["l_centroid_x"] = a["bounds"].left + (cl[1] + 0.5) * a["arr_cell_m"]
        # 正規化偏在距離 (画像座標で)
        dy = abs(ch[0] - cl[0]) / a["arr_h"]
        dx = abs(ch[1] - cl[1]) / a["arr_w"]
        a["spatial_diff"] = dy + dx
        a["n_high_cells"] = int(high_mask.sum())
        a["n_low_cells"] = int(low_mask.sum())
    else:
        a["h_centroid_y"] = a["h_centroid_x"] = None
        a["l_centroid_y"] = a["l_centroid_x"] = None
        a["spatial_diff"] = 0
        a["n_high_cells"] = a["n_low_cells"] = 0

# H5 は熊野町基準で判定 (仮説の主張通り)
spatial_diff_kuma = AREAS[1]["spatial_diff"]
h5 = spatial_diff_kuma >= 0.05
hypos.append({
    "H": "H5",
    "claim": "熊野町の高木は山岳 (北西側) に集中、低木は平地に分布 (空間偏在)",
    "result": f"熊野町: 高木重心 vs 低木重心 の正規化距離 = {spatial_diff_kuma:.3f} (>=0.05 で偏在あり) / 世羅町: {AREAS[0]['spatial_diff']:.3f}",
    "verdict": "支持" if h5 else ("部分支持" if spatial_diff_kuma >= 0.02 else "反証"),
})

# H6: 50cm はノイズ多 (std と 30m 超セル比率)
# Sera (1m) vs Kumano (50cm) の std と 30m 超
sera_std = s_sera["std"]
kuma_std = s_kuma["std"]
# 50cm は 1m よりノイズが多いか?
h6 = (kuma_std > sera_std)
hypos.append({
    "H": "H6",
    "claim": "50cm 解像度は 1m より標準偏差が大 (ノイズ多, 単木個葉のばらつき)",
    "result": f"世羅 (1m) std={sera_std:.2f}, 熊野 (50cm) std={kuma_std:.2f}, 50cm 最小値 (負値) = {s_kuma['min']:.2f}m",
    "verdict": "支持" if h6 else ("部分支持" if kuma_std >= sera_std * 0.85 else "反証"),
})

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


# =============================================================================
# 8. dataset/series 一覧 CSV
# =============================================================================
print("\n[8] 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"],
        "resolution_m": a["resolution_m"],
        "tif_size_MB": int(a["tif"].stat().st_size / (1024*1024)),
        "n_cells_M": int(a["meta"]["n_cells"] / 1e6),
        "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 / "L36_series.csv", index=False, encoding="utf-8-sig")
print(series_df.to_string(index=False))
print(f"  ({time.time()-t1:.1f}s)", flush=True)

# =============================================================================
# 9. 図 1: 2 dataset カード + 樹高ヒストグラム (1m bin)
# =============================================================================
print("\n[9] 図 1: dataset カード + ヒストグラム", flush=True)
t1 = time.time()

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# (左) カード
ax = axes[0]
ax.set_xlim(0, 10); ax.set_ylim(0, 10); ax.set_axis_off()
ax.text(5, 9.5, "DoBoX 樹高図シリーズ 2 件", ha="center", fontsize=15, fontweight="bold")

# Sera card
ax.add_patch(plt.Rectangle((0.3, 0.7), 4.4, 8.0, fill=True,
            facecolor=AREAS[0]["color"], alpha=0.13,
            edgecolor=AREAS[0]["color"], linewidth=2.5))
ax.text(2.5, 8.0, f"dsid {AREAS[0]['dataset_id']}\n{AREAS[0]['dataset_name']}", ha="center",
        fontsize=11, fontweight="bold", color=AREAS[0]["color"])
ax.text(2.5, 7.0, f"{AREAS[0]['label']}", ha="center", fontsize=14, fontweight="bold")
ax.text(2.5, 6.4, f"({AREAS[0]['character']})", ha="center", fontsize=8.5, color="#666")
ax.text(2.5, 5.6, f"解像度 {AREAS[0]['resolution_m']:.1f} m", ha="center", fontsize=10)
ax.text(2.5, 5.0, f"TIFF サイズ\n{AREAS[0]['meta']['width']:,}×{AREAS[0]['meta']['height']:,}",
        ha="center", fontsize=9.5)
ax.text(2.5, 4.0, f"{AREAS[0]['meta']['n_cells']/1e6:.0f}M セル", ha="center", fontsize=10)
ax.text(2.5, 3.2, f"範囲 {AREAS[0]['meta']['width_m']/1000:.1f}×{AREAS[0]['meta']['height_m']/1000:.1f} km",
        ha="center", fontsize=9)
ax.text(2.5, 2.2, f"平均樹高 {s_sera['mean']:.2f} m\n中央値 {s_sera['median']:.2f} m",
        ha="center", fontsize=10, fontweight="bold")
ax.text(2.5, 1.2, f"高木 (>=15m) {sera_pct_high:.1f}%", ha="center", fontsize=9, color=AREAS[0]["color"])

# Kumano card
ax.add_patch(plt.Rectangle((5.3, 0.7), 4.4, 8.0, fill=True,
            facecolor=AREAS[1]["color"], alpha=0.13,
            edgecolor=AREAS[1]["color"], linewidth=2.5))
ax.text(7.5, 8.0, f"dsid {AREAS[1]['dataset_id']}\n{AREAS[1]['dataset_name']}", ha="center",
        fontsize=11, fontweight="bold", color=AREAS[1]["color"])
ax.text(7.5, 7.0, f"{AREAS[1]['label']}", ha="center", fontsize=14, fontweight="bold")
ax.text(7.5, 6.4, f"({AREAS[1]['character']})", ha="center", fontsize=8.5, color="#666")
ax.text(7.5, 5.6, f"解像度 {AREAS[1]['resolution_m']:.2f} m", ha="center", fontsize=10)
ax.text(7.5, 5.0, f"TIFF サイズ\n{AREAS[1]['meta']['width']:,}×{AREAS[1]['meta']['height']:,}",
        ha="center", fontsize=9.5)
ax.text(7.5, 4.0, f"{AREAS[1]['meta']['n_cells']/1e6:.0f}M セル", ha="center", fontsize=10)
ax.text(7.5, 3.2, f"範囲 {AREAS[1]['meta']['width_m']/1000:.1f}×{AREAS[1]['meta']['height_m']/1000:.1f} km",
        ha="center", fontsize=9)
ax.text(7.5, 2.2, f"平均樹高 {s_kuma['mean']:.2f} m\n中央値 {s_kuma['median']:.2f} m",
        ha="center", fontsize=10, fontweight="bold")
ax.text(7.5, 1.2, f"地表 (<1m) {kuma_pct_zero_raw:.1f}%", ha="center", fontsize=9, color=AREAS[1]["color"])

# (右) ヒストグラム (両エリア重ね)
ax = axes[1]
bin_centers = np.arange(0.5, 30.5, 1)
ax.bar(bin_centers - 0.2, hist_1m["sera_1m"] / hist_1m["sera_1m"].sum() * 100,
       width=0.4, color=AREAS[0]["color"], alpha=0.85,
       label=f"世羅町 1m (n={int(hist_1m['sera_1m'].sum()):,})")
ax.bar(bin_centers + 0.2, hist_1m["kumano_50cm"] / hist_1m["kumano_50cm"].sum() * 100,
       width=0.4, color=AREAS[1]["color"], alpha=0.85,
       label=f"熊野町 50cm (n={int(hist_1m['kumano_50cm'].sum()):,})")
ax.set_xlabel("樹高 (m, 1m 階級)")
ax.set_ylabel("セル比率 (%)")
ax.set_title("樹高ヒストグラム — 2 エリア比較 (左: 世羅 / 右: 熊野)")
ax.set_xticks(np.arange(0, 31, 5))
ax.legend(loc="upper right")
ax.grid(axis="y", alpha=0.3)

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


# =============================================================================
# 10. 図 2: 2 エリアの樹高ラスタ可視化 (主題図 / heatmap)
# =============================================================================
print("\n[10] 図 2: 樹高ラスタマップ", flush=True)
t1 = time.time()

cmap_levels = [0, 1, 3, 6, 10, 15, 20, 30]
cmap_colors = ["#e0e0e0", "#fff7bc", "#fec44f", "#a1d99b",
               "#41ab5d", "#005a32", "#54278f"]
custom_cmap = ListedColormap(cmap_colors)
custom_norm = BoundaryNorm(cmap_levels, ncolors=len(cmap_colors), clip=True)

fig, axes = plt.subplots(1, 2, figsize=(15, 8))

for ax, a in zip(axes, AREAS):
    arr = a["arr"].astype(np.float32).copy()
    # 表示用: nodata は np.nan
    nodata = -99999.0
    arr[np.isnan(arr)] = nodata
    arr_show = np.where((arr == nodata) | (arr < -50), np.nan, arr)
    arr_show = np.clip(arr_show, 0, 30)

    extent = (a["bounds"].left, a["bounds"].right,
              a["bounds"].bottom, a["bounds"].top)
    im = ax.imshow(arr_show, cmap=custom_cmap, norm=custom_norm,
                   extent=extent, origin="upper", interpolation="nearest")

    # admin polygon を重ねる
    if a.get("admin") is not None:
        a["admin"].boundary.plot(ax=ax, color="#222", linewidth=1.5, alpha=0.9)

    ax.set_title(f"{a['label']} 樹高ラスタ (overview /1/{a['arr_factor']}, 実効 {a['arr_cell_m']:.0f} m/cell)\n"
                 f"{a['character']}", fontsize=11)
    ax.set_xlabel("X (EPSG:6671 m)")
    ax.set_ylabel("Y (EPSG:6671 m)")

# 共通カラーバー
sm = plt.cm.ScalarMappable(cmap=custom_cmap, norm=custom_norm)
sm.set_array([])
cbar = fig.colorbar(sm, ax=axes, orientation="horizontal", aspect=50,
                    pad=0.08, shrink=0.7)
cbar.set_label("樹高 (m) — 階級色 (≤1: 地表, 1-3: 低木, 3-6: 低中木, 6-10: 中木, 10-15: 中高木, 15-20: 高木, 20+: 巨木)")
cbar.set_ticks(cmap_levels)

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


# =============================================================================
# 11. 図 3: 階級別比率 (積み上げバー + 凡例)
# =============================================================================
print("\n[11] 図 3: 階級別比率", flush=True)
t1 = time.time()

fig, ax = plt.subplots(figsize=(12, 5.5))

# class_counts = (7,) per area  → 比率 (%)
class_pct_arr = []
for a in AREAS:
    cc = class_counts[a["key"]]
    pct = cc / cc.sum() * 100 if cc.sum() > 0 else cc * 0.0
    class_pct_arr.append(pct)

x = np.arange(len(AREAS))
bottom = np.zeros(len(AREAS))
hex_colors = ["#e0e0e0", "#fff7bc", "#fec44f", "#a1d99b",
              "#41ab5d", "#005a32", "#54278f"]
for i, (lab, col) in enumerate(zip(class_labels, hex_colors)):
    vals = [class_pct_arr[j][i] for j in range(len(AREAS))]
    ax.bar(x, vals, bottom=bottom, label=lab, color=col, edgecolor="#333", linewidth=0.6)
    for j, v in enumerate(vals):
        if v >= 3:
            ax.text(x[j], bottom[j] + v/2, f"{v:.1f}%", ha="center", va="center",
                    fontsize=10, fontweight="bold", color="black" if i in (0,1,2,3) else "white")
    bottom = bottom + np.array(vals)

ax.set_xticks(x)
ax.set_xticklabels([f"{a['label']}\n({a['resolution_m']}m, {a['character'].split('/')[0].strip()})" for a in AREAS], fontsize=10)
ax.set_ylabel("セル構成比 (%)")
ax.set_title("樹高階級別 構成比 — 7 階級 × 2 エリア")
ax.set_ylim(0, 100)
ax.legend(loc="center left", bbox_to_anchor=(1.02, 0.5))
ax.grid(axis="y", alpha=0.3)

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


# =============================================================================
# 12. 図 4: small multiples — 階級別マスクマップ (各エリア 7 panels)
# =============================================================================
print("\n[12] 図 4: 階級別マスク small multiples", flush=True)
t1 = time.time()

# 主要 4 階級だけ (0-1, 1-5, 5-15, 15+) で 2 エリア × 4 列
fig, axes = plt.subplots(2, 4, figsize=(15, 8))

mask_def = [
    ("地表 (<1m)",   lambda a: ((a >= 0) & (a < 1)),  "#bbbbbb"),
    ("低木 (1-5m)",  lambda a: ((a >= 1) & (a < 5)),  "#fec44f"),
    ("中木 (5-15m)", lambda a: ((a >= 5) & (a < 15)), "#41ab5d"),
    ("高木 (>=15m)", lambda a: (a >= 15),             "#54278f"),
]

for row_i, area in enumerate(AREAS):
    arr = area["arr"]
    valid_mask = area["valid_mask"]
    extent = (area["bounds"].left, area["bounds"].right,
              area["bounds"].bottom, area["bounds"].top)
    for col_i, (lbl, fn, col) in enumerate(mask_def):
        ax = axes[row_i, col_i]
        m = fn(arr) & valid_mask
        n_cells = int(m.sum())
        n_total_valid = int(valid_mask.sum())
        pct = n_cells / max(1, n_total_valid) * 100
        # 表示用: マスク 1 / その他 0
        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")
        # admin polygon を細く
        if area.get("admin") is not None:
            area["admin"].boundary.plot(ax=ax, color="#222", linewidth=0.7, alpha=0.7)
        if col_i == 0:
            ax.set_ylabel(f"{area['label']}\n({area['resolution_m']}m)", fontsize=11, 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", color="#1a1a1a",
                bbox=dict(facecolor="white", alpha=0.8, edgecolor="none", pad=2))
        ax.set_xticks([])
        ax.set_yticks([])

fig.suptitle("樹高階級別 空間分布 — 2 エリア × 4 階級 (small multiples)", fontsize=13, y=0.99)
plt.tight_layout()
fig.savefig(ASSETS / "L36_fig4_small_multiples.png", dpi=130, bbox_inches="tight")
plt.close("all")
print(f"  saved L36_fig4_small_multiples.png ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 13. 図 5: ヒストグラム (詳細, 2峰性確認用) + boxplot
# =============================================================================
print("\n[13] 図 5: ヒストグラム詳細 + boxplot", flush=True)
t1 = time.time()

fig, axes = plt.subplots(2, 2, figsize=(13, 9))

# (上左) 世羅 詳細ヒスト
ax = axes[0, 0]
v = AREAS[0]["valid"]
v = v[(v >= 0) & (v <= 30)]
ax.hist(v, bins=np.arange(0, 30.5, 0.5), color=AREAS[0]["color"], alpha=0.85, edgecolor="white", linewidth=0.3)
ax.axvline(s_sera["mean"], color="red", linestyle="--", linewidth=1.5, label=f"平均 {s_sera['mean']:.2f}m")
ax.axvline(s_sera["median"], color="orange", linestyle="--", linewidth=1.5, label=f"中央値 {s_sera['median']:.2f}m")
ax.set_xlabel("樹高 (m)")
ax.set_ylabel("セル数")
ax.set_title(f"世羅町 (1m) 樹高分布詳細 (0.5m bin) — 2 峰性?")
ax.legend()
ax.grid(axis="y", alpha=0.3)

# (上右) 熊野 詳細ヒスト
ax = axes[0, 1]
v = AREAS[1]["valid"]
v = v[(v >= 0) & (v <= 30)]
ax.hist(v, bins=np.arange(0, 30.5, 0.5), color=AREAS[1]["color"], alpha=0.85, edgecolor="white", linewidth=0.3)
ax.axvline(s_kuma["mean"], color="red", linestyle="--", linewidth=1.5, label=f"平均 {s_kuma['mean']:.2f}m")
ax.axvline(s_kuma["median"], color="orange", linestyle="--", linewidth=1.5, label=f"中央値 {s_kuma['median']:.2f}m")
ax.set_xlabel("樹高 (m)")
ax.set_ylabel("セル数")
ax.set_title(f"熊野町 (50cm) 樹高分布詳細 (0.5m bin)")
ax.legend()
ax.grid(axis="y", alpha=0.3)

# (下左) boxplot (両エリア)
ax = axes[1, 0]
data_box = [
    AREAS[0]["valid"][(AREAS[0]["valid"] >= 0) & (AREAS[0]["valid"] <= 30)],
    AREAS[1]["valid"][(AREAS[1]["valid"] >= 0) & (AREAS[1]["valid"] <= 30)],
]
bp = ax.boxplot(data_box, tick_labels=[f"{a['label']}\n({a['resolution_m']}m)" for a in AREAS],
                patch_artist=True, widths=0.5, showfliers=False)
for patch, area in zip(bp["boxes"], AREAS):
    patch.set_facecolor(area["color"])
    patch.set_alpha(0.7)
ax.set_ylabel("樹高 (m)")
ax.set_title("両エリア 樹高 boxplot (中央値・四分位)")
ax.grid(axis="y", alpha=0.3)

# (下右) パーセンタイル比較
ax = axes[1, 1]
percs = [10, 25, 50, 75, 90, 99]
xs = np.arange(len(percs))
sera_p = [s_sera[f"p{p}"] if p != 50 else s_sera["median"] for p in percs]
kuma_p = [s_kuma[f"p{p}"] if p != 50 else s_kuma["median"] for p in percs]
w = 0.35
ax.bar(xs - w/2, sera_p, w, label=f"世羅町 (1m)", color=AREAS[0]["color"], alpha=0.85, edgecolor="#333")
ax.bar(xs + w/2, kuma_p, w, label=f"熊野町 (50cm)", color=AREAS[1]["color"], alpha=0.85, edgecolor="#333")
for i, (s, k) in enumerate(zip(sera_p, kuma_p)):
    ax.text(i - w/2, s + 0.3, f"{s:.1f}", ha="center", fontsize=8.5)
    ax.text(i + w/2, k + 0.3, f"{k:.1f}", ha="center", fontsize=8.5)
ax.set_xticks(xs)
ax.set_xticklabels([f"P{p}" for p in percs])
ax.set_ylabel("樹高 (m)")
ax.set_title("樹高パーセンタイル比較 (P10〜P99)")
ax.legend()
ax.grid(axis="y", alpha=0.3)

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


# =============================================================================
# 14. 図 6: 高木集塊性 — 高木マスクの空間分布 + 重心矢印
# =============================================================================
print("\n[14] 図 6: 高木の空間集塊性", flush=True)
t1 = time.time()

fig, axes = plt.subplots(1, 2, figsize=(14, 7))

for ax, a in zip(axes, AREAS):
    extent = (a["bounds"].left, a["bounds"].right,
              a["bounds"].bottom, a["bounds"].top)
    arr = a["arr"]
    vm = a["valid_mask"]

    # ベース: グレースケールで樹高
    base = np.where(vm, np.clip(arr, 0, 30), np.nan)
    ax.imshow(base, cmap="Greys", extent=extent, origin="upper", alpha=0.4, vmin=0, vmax=20)

    # 高木 (>=15m) を緑で
    high = (arr >= 15) & vm
    high_show = np.where(high, 1.0, np.nan)
    ax.imshow(high_show, cmap=ListedColormap(["#005a32"]), extent=extent, origin="upper", vmin=0, vmax=1)

    # 低木 (<3m, valid) を薄黄色で
    low = (arr >= 0) & (arr < 3) & vm
    # 重心
    yh, xh = np.where(high)
    yl, xl = np.where(low)
    if len(yh) > 0:
        # 画像 row→Y_world (flip)
        cy_h = a["bounds"].top - (yh.mean() + 0.5) * a["arr_cell_m"]
        cx_h = a["bounds"].left + (xh.mean() + 0.5) * a["arr_cell_m"]
        ax.scatter([cx_h], [cy_h], marker="*", s=500, c="darkgreen", edgecolor="white", linewidth=2,
                   zorder=20, label=f"高木 (>=15m) 重心 (n={len(yh):,})")
    if len(yl) > 0:
        cy_l = a["bounds"].top - (yl.mean() + 0.5) * a["arr_cell_m"]
        cx_l = a["bounds"].left + (xl.mean() + 0.5) * a["arr_cell_m"]
        ax.scatter([cx_l], [cy_l], marker="X", s=400, c="orange", edgecolor="white", linewidth=2,
                   zorder=20, label=f"低木 (<3m) 重心 (n={len(yl):,})")
        # 矢印 (低木重心 → 高木重心)
        if len(yh) > 0:
            ax.annotate("", xy=(cx_h, cy_h), xytext=(cx_l, cy_l),
                        arrowprops=dict(arrowstyle="->", color="red", linewidth=2, alpha=0.8),
                        zorder=15)

    if a.get("admin") is not None:
        a["admin"].boundary.plot(ax=ax, color="#222", linewidth=1.0, alpha=0.85)

    ax.set_title(f"{a['label']} — 高木 (緑) と 低木重心 (橙) — 空間偏在の確認")
    ax.legend(loc="upper right")
    ax.set_xlabel("X (m)")
    ax.set_ylabel("Y (m)")

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


# =============================================================================
# 15. 図 7: 行政界 vs TIFF 範囲 (バッファ100m 注釈)
# =============================================================================
print("\n[15] 図 7: 行政界 vs TIFF bbox", flush=True)
t1 = time.time()

fig, axes = plt.subplots(1, 2, figsize=(14, 7))

for ax, a in zip(axes, AREAS):
    bb = a["bounds"]
    tiff_box = box(bb.left, bb.bottom, bb.right, bb.top)
    tiff_gs = gpd.GeoSeries([tiff_box], crs=a["meta"]["crs"])
    tiff_gs.boundary.plot(ax=ax, color="red", linewidth=2.5, label="TIFF bbox")
    tiff_gs.plot(ax=ax, color="red", alpha=0.07)

    if a.get("admin") is not None:
        a["admin"].boundary.plot(ax=ax, color=a["color"], linewidth=2.0, label=f"{a['label']} 行政界")
        a["admin"].plot(ax=ax, color=a["color"], alpha=0.13)

    # 上に「100m バッファ」と書く
    cx = (bb.left + bb.right) / 2
    cy = bb.top + 50
    ax.text(cx, cy, "TIFF は行政界 +100m バッファで切出し済み (DoBoX 仕様)",
            ha="center", fontsize=9.5, color="darkred",
            bbox=dict(facecolor="white", alpha=0.9, edgecolor="darkred"))

    ax.set_title(f"{a['label']}: 行政界 ({a.get('muni_area_km2', '?')} km²) vs TIFF bbox ({a['meta']['area_km2']:.1f} km²)")
    ax.legend(loc="lower right")
    ax.set_aspect("equal")

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


# =============================================================================
# 16. 図 8: ノイズ・外れ値分析 (50cm vs 1m)
# =============================================================================
print("\n[16] 図 8: ノイズ分析", flush=True)
t1 = time.time()

fig, axes = plt.subplots(2, 2, figsize=(13, 9))

# (上左) 負値の出現
ax = axes[0, 0]
neg_counts = []
for a in AREAS:
    v = a["valid"]
    n_neg = int((v < 0).sum())
    neg_counts.append({"label": a["label"], "n_neg": n_neg, "neg_pct": n_neg / len(v) * 100, "min": float(v.min())})
neg_df = pd.DataFrame(neg_counts)
ax.bar(range(len(AREAS)), neg_df["neg_pct"], color=[a["color"] for a in AREAS], alpha=0.85, edgecolor="#333")
for i, (n, mn) in enumerate(zip(neg_df["n_neg"], neg_df["min"])):
    ax.text(i, neg_df["neg_pct"][i] + 0.05, f"{n} cell ({neg_df['neg_pct'][i]:.2f}%)\nmin={mn:.2f}m",
            ha="center", fontsize=9, fontweight="bold")
ax.set_xticks(range(len(AREAS)))
ax.set_xticklabels([a["label"] for a in AREAS])
ax.set_ylabel("負値セル比率 (%)")
ax.set_title("(1) 負値セル — 計測ノイズの存在比")
ax.grid(axis="y", alpha=0.3)

# (上右) 超巨木 (>30m) の出現
ax = axes[0, 1]
big_counts = []
for a in AREAS:
    v = a["valid"]
    n_big = int((v >= 30).sum())
    big_counts.append({"label": a["label"], "n_big": n_big, "big_pct": n_big / len(v) * 100, "max": float(v.max())})
big_df = pd.DataFrame(big_counts)
ax.bar(range(len(AREAS)), big_df["big_pct"], color=[a["color"] for a in AREAS], alpha=0.85, edgecolor="#333")
for i, (n, mx) in enumerate(zip(big_df["n_big"], big_df["max"])):
    ax.text(i, big_df["big_pct"][i] + 0.005, f"{n} cell ({big_df['big_pct'][i]:.3f}%)\nmax={mx:.1f}m",
            ha="center", fontsize=9, fontweight="bold")
ax.set_xticks(range(len(AREAS)))
ax.set_xticklabels([a["label"] for a in AREAS])
ax.set_ylabel(">30m セル比率 (%)")
ax.set_title("(2) 30m 超セル — 超巨木 or ノイズ?")
ax.grid(axis="y", alpha=0.3)

# (下左) std 比較
ax = axes[1, 0]
std_data = [s_sera["std"], s_kuma["std"]]
colors = [a["color"] for a in AREAS]
ax.bar(range(len(AREAS)), std_data, color=colors, alpha=0.85, edgecolor="#333")
for i, v in enumerate(std_data):
    ax.text(i, v + 0.15, f"{v:.2f} m", ha="center", fontsize=12, fontweight="bold")
ax.set_xticks(range(len(AREAS)))
ax.set_xticklabels([f"{a['label']}\n({a['resolution_m']}m)" for a in AREAS])
ax.set_ylabel("樹高の標準偏差 (m)")
ax.set_title("(3) 樹高ばらつき (std) — 解像度差の影響?")
ax.grid(axis="y", alpha=0.3)

# (下右) パーセンタイル散布
ax = axes[1, 1]
percs = [1, 10, 25, 50, 75, 90, 99]
sera_pp = [np.percentile(AREAS[0]["valid"][AREAS[0]["valid"] >= 0], p) for p in percs]
kuma_pp = [np.percentile(AREAS[1]["valid"][AREAS[1]["valid"] >= 0], p) for p in percs]
ax.plot(percs, sera_pp, "o-", color=AREAS[0]["color"], label="世羅町 (1m)", linewidth=2, markersize=8)
ax.plot(percs, kuma_pp, "s-", color=AREAS[1]["color"], label="熊野町 (50cm)", linewidth=2, markersize=8)
ax.set_xlabel("パーセンタイル")
ax.set_ylabel("樹高 (m)")
ax.set_title("(4) パーセンタイル曲線比較")
ax.legend()
ax.grid(alpha=0.3)

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


# 中間 CSV
neg_df.to_csv(ASSETS / "L36_negative_outliers.csv", index=False, encoding="utf-8-sig")
big_df.to_csv(ASSETS / "L36_extreme_outliers.csv", index=False, encoding="utf-8-sig")


# =============================================================================
# 17. 図 9: ピクセル数値の Before/After (DSM-DTM 概念図)
# =============================================================================
print("\n[17] 図 9: DCHM 概念図 + サンプル", flush=True)
t1 = time.time()

# 概念図 + 実際の小サンプル抽出
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# (左) 概念図: DSM - DTM = DCHM
ax = axes[0]
ax.set_xlim(0, 10); ax.set_ylim(0, 10); ax.set_axis_off()
ax.text(5, 9.5, "DCHM (Digital Canopy Height Model) の作成", ha="center", fontsize=13, fontweight="bold")

# 地面 (DTM)
ax.add_patch(plt.Polygon([[0.5, 1.5], [9.5, 1.5], [9.5, 1], [0.5, 1]], facecolor="#8B4513", alpha=0.7))
ax.text(5, 0.5, "地面 (DTM = Digital Terrain Model)", ha="center", fontsize=10)

# 地物 (DSM = 地表全体)
import matplotlib.patches as mpatches
# 木1 (高い)
ax.add_patch(mpatches.Wedge((2, 1.5), 1.3, 0, 180, facecolor="#005a32", alpha=0.85))
ax.text(2, 4.5, "高 8m\n(樹高)", ha="center", fontsize=9, color="darkgreen", fontweight="bold")
# 木2 (中)
ax.add_patch(mpatches.Wedge((4.5, 1.5), 0.9, 0, 180, facecolor="#41ab5d", alpha=0.85))
ax.text(4.5, 3.5, "高 4m", ha="center", fontsize=9, color="darkgreen", fontweight="bold")
# 建物
ax.add_patch(plt.Rectangle((6, 1.5), 1.5, 2.5, facecolor="#666", alpha=0.85))
ax.text(6.75, 4.5, "高 2.5m\n(建物→\n樹高扱い)", ha="center", fontsize=9, color="black", fontweight="bold")
# 草地
ax.add_patch(plt.Polygon([[8.2, 1.5], [9.4, 1.5], [9.3, 1.7], [8.3, 1.7]], facecolor="#fff7bc", alpha=0.85))
ax.text(8.7, 2.2, "高 0.2m\n(地表)", ha="center", fontsize=8.5)

ax.text(5, 7.0, "DSM = 地表全体の高さ (Digital Surface Model)", ha="center", fontsize=10.5, color="black")
ax.text(5, 6.3, "DTM = 地面のみの高さ", ha="center", fontsize=10.5, color="brown")
ax.text(5, 5.4, "→ DCHM = DSM − DTM = 地物の高さ (= 樹高)", ha="center",
        fontsize=11.5, color="darkgreen", fontweight="bold",
        bbox=dict(facecolor="#fff7bc", alpha=0.9, edgecolor="black", boxstyle="round,pad=0.5"))

# (右) 実データ 256x256 サンプル (各エリアの中央)
ax = axes[1]
sample_arrs = []
for a in AREAS:
    with rasterio.open(a["tif"]) as ds:
        cx, cy = ds.width//2, ds.height//2
        SS = 128
        win = Window(cx - SS//2, cy - SS//2, SS, SS)
        arr_sample = ds.read(1, window=win)
        if ds.nodata is not None:
            arr_sample = np.where(arr_sample == ds.nodata, np.nan, arr_sample)
        sample_arrs.append((a["label"], a["resolution_m"], arr_sample))

# 横並びで2サンプル表示
ax.set_axis_off()
fig2, sub_axes = plt.subplots(1, 2, figsize=(14, 6))
for sub_ax, (lab, res, arr) in zip(sub_axes, sample_arrs):
    arr2 = np.clip(arr, 0, 30)
    im = sub_ax.imshow(arr2, cmap=custom_cmap, norm=custom_norm, interpolation="nearest")
    sub_ax.set_title(f"{lab} 中央 {arr.shape[1]}×{arr.shape[0]} = {arr.shape[1]*res:.0f}m × {arr.shape[0]*res:.0f}m サンプル\n"
                     f"({res}m/cell, 平均 {np.nanmean(arr):.2f}m)")
    sub_ax.set_xlabel("ピクセル X")
    sub_ax.set_ylabel("ピクセル Y")
    plt.colorbar(im, ax=sub_ax, fraction=0.045)

plt.tight_layout()
fig2.savefig(ASSETS / "L36_fig9_sample_centers.png", dpi=130, bbox_inches="tight")
plt.close("all")
plt.close(fig)

# 概念図単独保存
fig, ax = plt.subplots(figsize=(10, 6))
ax.set_xlim(0, 10); ax.set_ylim(0, 10); ax.set_axis_off()
ax.text(5, 9.5, "DCHM (Digital Canopy Height Model) の作成", ha="center", fontsize=13, fontweight="bold")
ax.add_patch(plt.Polygon([[0.5, 1.5], [9.5, 1.5], [9.5, 1], [0.5, 1]], facecolor="#8B4513", alpha=0.7))
ax.text(5, 0.5, "地面 (DTM = Digital Terrain Model)", ha="center", fontsize=10)
ax.add_patch(mpatches.Wedge((2, 1.5), 1.3, 0, 180, facecolor="#005a32", alpha=0.85))
ax.text(2, 4.5, "高 8m\n(樹高)", ha="center", fontsize=9, color="darkgreen", fontweight="bold")
ax.add_patch(mpatches.Wedge((4.5, 1.5), 0.9, 0, 180, facecolor="#41ab5d", alpha=0.85))
ax.text(4.5, 3.5, "高 4m", ha="center", fontsize=9, color="darkgreen", fontweight="bold")
ax.add_patch(plt.Rectangle((6, 1.5), 1.5, 2.5, facecolor="#666", alpha=0.85))
ax.text(6.75, 4.5, "高 2.5m\n(建物→\n樹高扱い)", ha="center", fontsize=9, color="black", fontweight="bold")
ax.add_patch(plt.Polygon([[8.2, 1.5], [9.4, 1.5], [9.3, 1.7], [8.3, 1.7]], facecolor="#fff7bc", alpha=0.85))
ax.text(8.7, 2.2, "高 0.2m\n(地表)", ha="center", fontsize=8.5)
ax.text(5, 7.0, "DSM = 地表全体の高さ (Digital Surface Model)", ha="center", fontsize=10.5, color="black")
ax.text(5, 6.3, "DTM = 地面のみの高さ", ha="center", fontsize=10.5, color="brown")
ax.text(5, 5.4, "→ DCHM = DSM − DTM = 地物の高さ (= 樹高)", ha="center",
        fontsize=11.5, color="darkgreen", fontweight="bold",
        bbox=dict(facecolor="#fff7bc", alpha=0.9, edgecolor="black", boxstyle="round,pad=0.5"))
plt.tight_layout()
fig.savefig(ASSETS / "L36_fig9_dchm_concept.png", dpi=130, bbox_inches="tight")
plt.close("all")

print(f"  saved L36_fig9 ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 17b. 図 10: overview vs 生解像度 サンプル の階級比較
# =============================================================================
print("\n[17b] 図 10: overview vs 生解像度", flush=True)
t1 = time.time()

fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))

# overview (全域 64m/cell) と raw (中央 1024x1024 セル本来解像度)
class_labels_short = ["0-1m\n地表", "1-5m\n低木", "5-15m\n中木", ">=15m\n高木"]

# 各エリア: 4 階級比率 (overview と raw)
for area_idx, a in enumerate(AREAS):
    ax = axes[area_idx]
    # overview の階級比率
    valid_ovr = a["valid"]
    n_ovr = len(valid_ovr)
    pct_zero_ovr = ((valid_ovr >= 0) & (valid_ovr < 1)).sum() / n_ovr * 100
    pct_low_ovr = ((valid_ovr >= 1) & (valid_ovr < 5)).sum() / n_ovr * 100
    pct_mid_ovr = ((valid_ovr >= 5) & (valid_ovr < 15)).sum() / n_ovr * 100
    pct_high_ovr = (valid_ovr >= 15).sum() / n_ovr * 100
    ovr_pct = [pct_zero_ovr, pct_low_ovr, pct_mid_ovr, pct_high_ovr]

    raw_pct = [a["raw_sample_pct_zero"], a["raw_sample_pct_low"],
               a["raw_sample_pct_mid"], a["raw_sample_pct_high"]]

    x = np.arange(4)
    w = 0.35
    bars1 = ax.bar(x - w/2, ovr_pct, w, label=f"overview /1/{a['arr_factor']} ({a['arr_cell_m']:.0f}m/cell, n={n_ovr:,})",
                   color="#888", alpha=0.85, edgecolor="#333")
    bars2 = ax.bar(x + w/2, raw_pct, w, label=f"生解像度中央サンプル ({a['resolution_m']}m/cell, n={a['raw_sample_n']:,})",
                   color=a["color"], alpha=0.85, edgecolor="#333")
    for bars, vals in [(bars1, ovr_pct), (bars2, raw_pct)]:
        for b, v in zip(bars, vals):
            ax.text(b.get_x() + b.get_width()/2, v + 0.5, f"{v:.1f}",
                    ha="center", fontsize=8.5, fontweight="bold")

    ax.set_xticks(x)
    ax.set_xticklabels(class_labels_short)
    ax.set_ylabel("セル比率 (%)")
    ax.set_title(f"{a['label']} — overview vs 生解像度 サンプル")
    ax.legend(loc="upper right", fontsize=8.5)
    ax.grid(axis="y", alpha=0.3)
    ax.set_ylim(0, max(max(ovr_pct), max(raw_pct)) * 1.2)

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


# =============================================================================
# 18. 中間 CSV: 詳細 1m 階級ヒスト + 階級別 + 仮説結果
# =============================================================================
print("\n[18] 中間 CSV 出力", flush=True)
t1 = time.time()

# 1m 階級 ヒスト
hist_df = pd.DataFrame({
    "height_bin_low_m":  bins_1m[:-1],
    "height_bin_high_m": bins_1m[1:],
    f"sera_1m":          hist_1m["sera_1m"],
    f"kumano_50cm":      hist_1m["kumano_50cm"],
})
hist_df["sera_pct"]   = hist_df["sera_1m"]   / hist_df["sera_1m"].sum() * 100
hist_df["kumano_pct"] = hist_df["kumano_50cm"] / hist_df["kumano_50cm"].sum() * 100
hist_df.to_csv(ASSETS / "L36_height_histogram_1m_bins.csv", index=False, encoding="utf-8-sig")

# 階級別 wide
class_w = pd.DataFrame({"class": class_labels})
class_w["sera_n"]  = class_counts["sera_1m"]
class_w["kumano_n"] = class_counts["kumano_50cm"]
class_w["sera_pct"]   = class_w["sera_n"]   / class_w["sera_n"].sum() * 100
class_w["kumano_pct"] = class_w["kumano_n"] / class_w["kumano_n"].sum() * 100
class_w.to_csv(ASSETS / "L36_height_classes_wide.csv", index=False, encoding="utf-8-sig")

print(f"  ({time.time()-t1:.1f}s)", flush=True)


# =============================================================================
# 19. HTML 生成
# =============================================================================
print("\n[19] 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>「樹高図 1m / 樹高図 50cm」 2 件</b>
(dataset_id = 1626, 1627) を統合し、広島県 2 自治体エリア
(<b>世羅町</b> 1m, <b>熊野町</b> 50cm) の<b>樹高分布</b>から
<b>林相 (針葉樹/広葉樹/伐採地/都市・農地)</b> を推定する研究記事です。
両データセットは航空レーザ計測 (LiDAR) 由来の<b>樹冠高さラスタ (DCHM)</b> で、
1 ピクセルあたり 1m もしくは 0.5m のメッシュで地物高さが記録されています。
</div>

<h3>シリーズ構造判定: (a) 同一手法・解像度違い + (b) 自治体ごとの分割</h3>
<p>このシリーズの 2 件は、<b>「解像度 1m vs 50cm」</b>での同一手法 2 系統です。
両者は同じ航空レーザ計測のオリジナルデータから、<b>解像度を変えて (1m メッシュ / 50cm メッシュ)</b>
県内 23 自治体それぞれを GeoTIFF に切出した形で公開されています。</p>
<ul>
<li><b>1m 版 (dsid 1626)</b>: 広域評価向け。県全域の 22 自治体に対応 GeoTIFF 公開 (本記事は<b>世羅町</b>版を使用)</li>
<li><b>50cm 版 (dsid 1627)</b>: 詳細評価向け。同じく 22 自治体に対応 GeoTIFF 公開 (本記事は<b>熊野町</b>版を使用)</li>
<li>両者ともデータ容量が極めて大 (50cm 広島市版は<b>15.9 GB ZIP</b>!) のため、
本記事は<b>各 1 自治体のみ</b>を代表エリアとして扱います。</li>
</ul>
<p>L32〜L35 (港湾・河川シリーズ) は<b>線・点</b>のインフラ施設データだったが、
L36 から始まる<b>LiDAR ファミリ</b>は<b>面 (ラスタ)</b>のフィジカル測定が主役。
扱う技術もまったく異なる。</p>

<h3>独自用語の定義 (本記事内で固定)</h3>
<table>
<tr><th>用語</th><th>定義</th></tr>
<tr><td><b>樹高 (じゅこう)</b></td><td>地表から樹冠頂点までの高さ (m)。本記事では DCHM ピクセル値そのもの。
建物・電柱なども<b>「地物」</b>として高さを取るため、純粋な「木の高さ」ではないことに注意。
広島県の樹高図仕様では「200m 以上は雲・霧として除去、それ以下のノイズも可能な限り除去」とされている。</td></tr>
<tr><td><b>DCHM (Digital Canopy Height Model)</b></td><td>樹冠高さモデル。航空レーザ測量で得た DSM (Digital Surface Model = 地表全体の高さ) から
DTM (Digital Terrain Model = 地面のみの高さ) を差し引いて作成した<b>地物高さラスタ</b>。
本記事の樹高図そのもの。1 セルが 1m もしくは 0.5m。</td></tr>
<tr><td><b>DSM</b></td><td>Digital Surface Model。航空レーザが<b>最初に当たった点</b> (= 樹冠頂点・建物屋上・地面) の高さを集めたラスタ。地物の上端を見ている。</td></tr>
<tr><td><b>DTM</b></td><td>Digital Terrain Model。航空レーザの<b>最後に当たった点</b> (= 地面に届いた点) を補間した「地面のみ」の高さラスタ。建物・木は除去済み。</td></tr>
<tr><td><b>航空レーザ (LiDAR)</b></td><td>Light Detection And Ranging。航空機から地面に向けて<b>レーザパルス</b>を毎秒数十万発撃ち、戻り光の往復時間から距離を測る測量手法。
1 m² あたり数〜数十パルスを撃ち、樹冠と地面を別々に検出して 3D 点群を作る。
広島県では砂防課・林野庁治山課・国交省で公共測量として実施。</td></tr>
<tr><td><b>林相 (りんそう)</b></td><td>森林の構成・外見的特徴。針葉樹 (スギ・ヒノキ等) / 広葉樹 (コナラ・ブナ等) /
混交林 / 竹林 / 伐採地 などに分類される。本記事では<b>樹高分布パターンから推定</b>するヒューリスティック手法を採用。</td></tr>
<tr><td><b>地物高さ階級</b></td><td>本記事で独自定義した 7 階級 (≤1m 地表, 1-3m 低木, 3-6m 低中木, 6-10m 中木, 10-15m 中高木, 15-20m 高木, ≥20m 巨木)。
<b>DoBoX 公式の分類ではなく、本記事の分析用ヒューリスティック</b>。</td></tr>
<tr><td><b>overview / ピラミッド</b></td><td>GeoTIFF に埋め込まれた<b>低解像度版</b>。例えば「/1/128」は元データの 1/128 解像度を意味する。
1m × 27,248 × 16,357 (445M セル) のラスタも、/128 オーバビューなら 213×128 = 27K セル (1.5万倍小) で全体概観できる。
本記事ではこれを使ってメモリを <b>0.1 MB に抑える</b>。</td></tr>
<tr><td><b>有効データ率</b></td><td>TIFF 全セルのうち、nodata でも NaN でもない<b>有効値を持つセルの比率</b>。
TIFF は行政界 +100m バッファで矩形に切出されているため、行政界の外側は nodata になる。</td></tr>
</table>

<h3>研究の問い (RQ)</h3>
<p>広島県の樹高図 2 解像度 (1m / 50cm) は、それぞれ<b>世羅町 (中山間, 林業地)</b>と
<b>熊野町 (都市近郊, 山岳混在)</b>という<b>異なる地域特性</b>のエリアを公開している。
両エリアの樹高ヒストグラム・空間分布・林相を比較すると、
どんな地域構造の違いが見え、解像度差はどう影響するか?</p>

<ol>
<li>2 dataset の<b>面積・解像度・有効データ率</b>の構造比較は?</li>
<li>樹高ヒストグラム (0-30m) の<b>ピーク位置・多峰性</b>はどう違うか?</li>
<li><b>高木 (>=15m) / 低木 (<5m) / 地表 (<1m)</b> の構成比は?</li>
<li>高木地・低木地の<b>空間分布</b>はどう偏在しているか?</li>
<li>樹高分布から<b>林相</b>はどう推定できるか?</li>
<li>解像度差 (1m vs 50cm) は<b>ノイズ・極端値</b>にどう影響するか?</li>
</ol>

<h3>仮説 H1〜H6</h3>
<ul>
<li><b>H1 (世羅町は林業地優占)</b>: 平均樹高 ≥ 8m、高木 (≥15m) 比率 ≥ 15%。
中山間で植林スギ・ヒノキ針葉樹が広面積を占めると推定。</li>
<li><b>H2 (熊野町は地表優占)</b>: 樹高 0-1m セル比率 ≥ 40%。
都市部の建物外周・農地・道路・山頂部のはげ地が高い割合。</li>
<li><b>H3 (両エリアともゼロ卓越)</b>: 樹高 0-1m が最頻ビン (mode)。
航空 LiDAR では多数のセルが地表 (建物・農地) もしくは雲・霧除去後ノイズで 0 付近。</li>
<li><b>H4 (世羅町の樹高分布は 2 峰性)</b>: 0-1m と 10-20m の 2 ピーク。
林業地 = 「植林伐採跡 + 成熟林」の二極構造。</li>
<li><b>H5 (熊野町の高木は山岳に集中)</b>: 高木 (>=15m) は山地に集中、平地は低木・地表優占。空間相関で確認。</li>
<li><b>H6 (50cm 版はノイズ多)</b>: 50cm の標準偏差は 1m より大、負値出現、極端値 (>30m) も多。
高解像度ほど雲・霧の残存ノイズや単木個葉のばらつきを拾う。</li>
</ul>

<h3>到達点</h3>
<p>2 解像度の樹高図ラスタを「<b>2 自治体の代表エリア</b>」として読み解き、
合計 <b>{int(metas[0]['n_cells']/1e6 + metas[1]['n_cells']/1e6):,} M セル</b> (約 7 億ピクセル) の DCHM データを
<b>overview ピラミッド</b>を使って効率的に概観し、樹高ヒストグラム・空間分布・林相推定の 3 角度から
広島県の中山間 vs 都市近郊の<b>森林構造の違い</b>を実データで裏付ける。
さらに解像度差 (1m vs 50cm) が分析にどう影響するかを統計的に示す。</p>

<div class="warn"><b>本記事のスコープ外</b>:
本記事は DCHM (樹高ラスタ) の<b>分布構造</b>に集中する。
個々の樹木の<b>樹種同定</b>(画像分類タスク, 別記事)、
<b>個葉スケールの形状解析</b>(50cm でも個葉は分解できない)、
<b>時系列変化</b>(本記事の TIFF は 1 時点の計測, 経年変化は L37「計測年度図」と組合わせる発展課題)、
<b>正確な針葉樹/広葉樹判別</b>(LiDAR 単独では困難, 衛星画像 NDVI 等の組合せが必要)
は本記事では扱わず発展課題とする。
また、<b>本記事の「林相推定」は樹高分布パターンに基づくヒューリスティックであり、林野庁の公式分類ではない</b>。
</div>
"""
sections.append(("学習目標と問い", sec1))


# === Section 2: 使用データ ===
sec2_top = f"""
<p>本記事が使用する 2 dataset の一覧。両者は<b>同一手法 (LiDAR DCHM) の解像度違い</b>で、
県内 22-23 自治体それぞれに GeoTIFF が公開されています。本記事は<b>各 1 自治体</b>を代表として用います。</p>

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

<h3>データ仕様 (両 dataset 共通)</h3>
<ul>
<li><b>形式</b>: GeoTIFF (Float32, 単一バンド) — ZIP に格納された .tiff + .tiff.ovr (オーバビューピラミッド)</li>
<li><b>CRS</b>: EPSG:6671 (JGD2011 / Japan Plane Rectangular CS III) — m 単位の平面直角</li>
<li><b>nodata</b>: -99999.0 (未計測領域)</li>
<li><b>値の意味</b>: <b>地物高さ (m)</b> = DSM (地表全体) − DTM (地面のみ)。樹冠頂点・建物屋上などすべて含む</li>
<li><b>切出単位</b>: 各自治体の行政界 + 100m バッファで矩形切出</li>
<li><b>除去ノイズ</b>: 200m 以上は雲・霧として除去済み (公式仕様)</li>
<li><b>ライセンス</b>: CC-BY 4.0 (測量法 44 条に基づく公共測量 2 次著作物)</li>
</ul>

<h3>L32-L35 シリーズとの重要な相違</h3>
<table>
<tr><th>軸</th><th>L32-L35 (港湾/河川)</th><th>L36 (本記事, 樹高図)</th></tr>
<tr><td>データ形式</td><td>CSV (テーブル) + Shapefile (ベクタ)</td><td><b>GeoTIFF (ラスタ)</b></td></tr>
<tr><td>レコード単位</td><td>1 行 = 1 施設・1 ポリゴン</td><td><b>1 セル = 地表 1m² (or 0.25m²)</b> の高さ値</td></tr>
<tr><td>件数</td><td>数十〜数千</td><td><b>数億セル</b> (世羅 446M, 熊野 261M)</td></tr>
<tr><td>分析操作</td><td>groupby / sjoin / overlay</td><td><b>imshow / 階級分類 / overview</b></td></tr>
<tr><td>主指標</td><td>件数・延長 m・容量</td><td><b>樹高 m・階級比率</b></td></tr>
<tr><td>必要パッケージ</td><td>geopandas, shapely</td><td><b>rasterio</b> (+ geopandas)</td></tr>
<tr><td>ファイルサイズ</td><td>数 KB 〜 数 MB</td><td><b>数百 MB 〜 16 GB</b></td></tr>
</table>

<h3>2 自治体の特性比較</h3>
<table>
<tr><th>軸</th><th>世羅町 (1m)</th><th>熊野町 (50cm)</th></tr>
<tr><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></tr>
<tr><td>TIFF 範囲 (bbox)</td>
<td>{AREAS[0]['meta']['width_m']/1000:.1f} × {AREAS[0]['meta']['height_m']/1000:.1f} km</td>
<td>{AREAS[1]['meta']['width_m']/1000:.1f} × {AREAS[1]['meta']['height_m']/1000:.1f} km</td></tr>
<tr><td>主産業</td><td>農林業 (世羅梨・りんご・大根)</td><td>筆 (熊野筆) 製造</td></tr>
<tr><td>有効データ率</td>
<td>{basic_stats['sera_1m']['valid_ratio']*100:.1f}%</td>
<td>{basic_stats['kumano_50cm']['valid_ratio']*100:.1f}%</td></tr>
<tr><td>平均樹高</td>
<td>{s_sera['mean']:.2f} m</td>
<td>{s_kuma['mean']:.2f} m</td></tr>
<tr><td>中央値</td>
<td>{s_sera['median']:.2f} m</td>
<td>{s_kuma['median']:.2f} m</td></tr>
</table>

<h3>データ品質メモ</h3>
<ul>
<li><b>サイズが極めて大</b>: 世羅町 (1m) ZIP 870 MB → TIFF 1.7 GB、熊野町 (50cm) ZIP 604 MB → TIFF 1.0 GB。
広島市 (50cm) は ZIP <b>15.9 GB</b>。<b>全件 DL は事実上不可能</b>で、本記事は<b>各 1 自治体のみ</b>を扱う設計とした。</li>
<li><b>overview 利用必須</b>: 1m × 27,248 × 16,357 (446M セル) を全部メモリに載せたら 1.7 GB。
overview /1/{AREAS[0]['arr_factor']} を使えば {AREAS[0]['arr_w']}×{AREAS[0]['arr_h']} セル (0.1 MB) で全体概観可能。</li>
<li><b>nodata = -99999</b>: TIFF は行政界 +100m バッファで矩形切出されているため、行政界の外側は nodata。
有効データ率は {basic_stats['sera_1m']['valid_ratio']*100:.0f}〜{basic_stats['kumano_50cm']['valid_ratio']*100:.0f}% 程度。</li>
<li><b>負値の存在</b>: 熊野町 (50cm) には<b>{int((AREAS[1]['valid'] < 0).sum())} セル</b>の負値が出現 (最小 {s_kuma['min']:.2f}m)。
理論上樹高は ≥0m だが、DSM-DTM の差分処理でセンサノイズ・地形誤差が負値として残る。</li>
<li><b>「林相推定」はヒューリスティック</b>: 本記事の林相判定は樹高分布パターンに基づく独自ルールで、
林野庁の森林簿や植生図の正式分類ではない。学習者にも「あくまで推定」と明示する。</li>
</ul>
"""
sections.append(("使用データ", sec2_top))


# === Section 3: ダウンロード ===
sec3 = f"""
<h3>生データ (DoBoX 直リンク)</h3>
<ul>
<li>樹高図（1m）_世羅町: <a href="https://data.hiroshima-dobox.jp/forest/th_1m/34462_世羅町.zip">直 DL ZIP (870 MB)</a> /
<a href="https://hiroshima-dobox.jp/datasets/1626">DoBoX dsid 1626</a></li>
<li>樹高図（50cm）_熊野町: <a href="https://data.hiroshima-dobox.jp/forest/th_50cm/34307_熊野町.zip">直 DL ZIP (604 MB)</a> /
<a href="https://hiroshima-dobox.jp/datasets/1627">DoBoX dsid 1627</a></li>
</ul>

<div class="warn"><b>容量警告</b>: 上記 2 ZIP の合計は約 1.5 GB、展開後の TIFF は約 2.7 GB。
ハンズオン環境のディスク・メモリに余裕があるか事前確認してください。
他の自治体版 (例: 50cm 広島市は 15.9 GB) は本記事のスコープ外。</div>

<h3>中間データ (本記事生成 CSV)</h3>
<ul>
<li>{dl("L36_series.csv")} — 2 dataset の一覧 (代表自治体・解像度・サイズ)</li>
<li>{dl("L36_tiff_meta.csv")} — TIFF メタ (CRS, bounds, セル数, overview レベル)</li>
<li>{dl("L36_basic_stats.csv")} — 基本統計 (平均・中央値・パーセンタイル・nodata 比率)</li>
<li>{dl("L36_raw_sample_stats.csv")} — 生解像度中央サンプル統計 (overview 平均化バイアス補正用)</li>
<li>{dl("L36_height_histogram_1m_bins.csv")} — 1m 階級ヒストグラム</li>
<li>{dl("L36_height_classes.csv")} / {dl("L36_height_classes_wide.csv")} — 7 階級別構成比</li>
<li>{dl("L36_linso_estimate.csv")} — 林相ヒューリスティック推定</li>
<li>{dl("L36_negative_outliers.csv")} / {dl("L36_extreme_outliers.csv")} — 負値・30m 超の出現</li>
<li>{dl("L36_hypothesis_results.csv")} / {dl("L36_hypothesis_results.json")} — H1〜H6 検証結果</li>
</ul>

<h3>図 (本記事生成 PNG)</h3>
<ul>
<li>{dl("L36_fig1_dataset_overview.png")} — 2 dataset カード + 樹高ヒストグラム重ね</li>
<li>{dl("L36_fig2_raster_maps.png")} — 樹高ラスタ主題図 (両エリア)</li>
<li>{dl("L36_fig3_class_composition.png")} — 7 階級構成比 (積み上げ)</li>
<li>{dl("L36_fig4_small_multiples.png")} — 階級別マスクマップ (2×4 panels)</li>
<li>{dl("L36_fig5_distribution_detail.png")} — ヒスト詳細 + boxplot + パーセンタイル</li>
<li>{dl("L36_fig6_clustering.png")} — 高木の空間集塊性 (重心矢印)</li>
<li>{dl("L36_fig7_admin_bbox.png")} — 行政界 vs TIFF bbox</li>
<li>{dl("L36_fig8_noise_analysis.png")} — 解像度差ノイズ分析</li>
<li>{dl("L36_fig9_dchm_concept.png")} / {dl("L36_fig9_sample_centers.png")} — DCHM 概念図 + 中央サンプル</li>
<li>{dl("L36_fig10_overview_vs_raw.png")} — overview 平均化バイアスと生解像度サンプルの比較</li>
</ul>

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


# === Section 4: 分析 1 — 2 dataset とエリアの構造を見る ===
sec4 = f"""
<h3>狙い</h3>
<p>2 dataset (樹高図 1m / 50cm) と本記事が選んだ 2 エリア (世羅町 / 熊野町) の
<b>規模・形状・データ量</b>を 1 枚の絵で示す。樹高ヒストグラムを重ねて
<b>地域特性の違い</b>を最初に視覚化する。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p>GeoTIFF を扱う基本ツールは <b>rasterio</b>:</p>
<ul>
<li><b>rasterio.open(path)</b> で TIFF を開く (まだメモリには載らない)</li>
<li><b>ds.width, ds.height</b> でセル数、<b>ds.bounds</b> で範囲、<b>ds.res</b> で解像度</li>
<li><b>ds.overviews(1)</b> でピラミッド (低解像度版) の倍率リスト</li>
<li><b>ds.read(1, out_shape=(H, W), resampling=Resampling.average)</b> で<b>低解像度版を一括取得</b></li>
</ul>
<p><b>ピラミッド戦略</b>: 1m × 27,248 × 16,357 (446M セル, 1.7GB) を全部読むと OOM。
overview /1/{AREAS[0]['arr_factor']} を使えば {AREAS[0]['arr_w']}×{AREAS[0]['arr_h']} = {AREAS[0]['arr_w']*AREAS[0]['arr_h']/1e6:.2f}M セル (0.1 MB) で
<b>全体概観</b>が可能。本記事は<b>このピラミッド読みを徹底</b>することで 1 分以内完走を実現する。</p>

<h3>実装</h3>
{code('''
import rasterio
from rasterio.enums import Resampling

# 各 TIFF を overview で読む (目標 ≈ 64m/cell)
TARGET_CELL_M = 64
for area in AREAS:
    with rasterio.open(area["tif"]) as ds:
        # overview 倍率リスト [2, 4, 8, 16, 32, 64, 128]
        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])

        ov_w = ds.width // best
        ov_h = ds.height // best

        # 低解像度版を一括読込 (average resampling)
        arr = ds.read(1, out_shape=(ov_h, ov_w), resampling=Resampling.average)
        area["arr"] = arr  # shape: (ov_h, ov_w), dtype=float32

        # nodata マスク
        valid_mask = (~np.isnan(arr)) & (arr != ds.nodata) & (arr > -1000)
        area["valid"] = arr[valid_mask]
''')}

<h3>図と読み取り (図 1: dataset カード + ヒストグラム重ね)</h3>
<p><b>なぜこの図か</b>: 2 dataset の規模と特性を「カード形式」で文字情報として、
かつ「樹高ヒストグラム重ね」で量的特性として、<b>同時に 1 枚で伝える</b>ため。
左カード = 言語的概要、右ヒスト = 量的概観。</p>

{figure("assets/L36_fig1_dataset_overview.png", "図 1: 2 dataset 概要カード (左) と 樹高ヒストグラム重ね (右)")}

<p><b>読み取り</b>:</p>
<ul>
<li>世羅町 (1m): {AREAS[0]['meta']['width']:,}×{AREAS[0]['meta']['height']:,} = {AREAS[0]['meta']['n_cells']/1e6:.0f}M セル、平均樹高 <b>{s_sera['mean']:.2f} m</b>、高木 (>=15m) 比率 <b>{sera_pct_high:.1f}%</b>。</li>
<li>熊野町 (50cm): {AREAS[1]['meta']['width']:,}×{AREAS[1]['meta']['height']:,} = {AREAS[1]['meta']['n_cells']/1e6:.0f}M セル、平均樹高 <b>{s_kuma['mean']:.2f} m</b>、地表 (<1m) 比率 <b>{kuma_pct_zero_raw:.1f}%</b> (生解像度中央サンプル)。</li>
<li>ヒスト重ね (右): <b>両エリアとも 0-1m が圧倒的最頻ビン</b>。これは航空 LiDAR ラスタの典型 — 大量の地表セルが含まれる。</li>
<li>世羅 (緑) は<b>10-15m 帯にも明瞭な第 2 ピーク</b> = 林業地の植林帯らしい兆候。</li>
<li>熊野 (青) は<b>0-1m が世羅より急峻に高い</b> = 都市部・農地の比率が大きい。</li>
<li>両ヒストともセル比率の合計は 100% (有効セル基準)。<b>nodata セルは除外</b>している。</li>
</ul>

<h3>表と読み取り (基本統計)</h3>
{stats_df.to_html(classes='', border=0)}

<p><b>読み取り</b>: 平均樹高は<b>世羅 {s_sera['mean']:.2f} m vs 熊野 {s_kuma['mean']:.2f} m</b> で大差。
中央値も<b>世羅 {s_sera['median']:.2f} m vs 熊野 {s_kuma['median']:.2f} m</b>。
P90 (上位 10%) は<b>世羅 {s_sera['p90']:.1f} m vs 熊野 {s_kuma['p90']:.1f} m</b>で
両エリアとも 10m 以上の高木はある。
<b>P99 (上位 1%)</b> は<b>世羅 {s_sera['p99']:.1f} m vs 熊野 {s_kuma['p99']:.1f} m</b>でほぼ同等で、
両者ともに最高樹高は 20m 超 (= 成熟林) を含むことが分かる。</p>
"""
sections.append(("分析 1: 2 dataset とエリアの構造を可視化", sec4))


# === Section 5: 分析 2 — 樹高ラスタの主題図 ===
sec5 = f"""
<h3>狙い</h3>
<p>樹高分布の<b>空間構造 (どこに高木 / どこに低木)</b> を地図で見る。
ヒストグラムは「全体としてどう」しか分からないが、
主題図は「<b>町のどの場所に</b>どの樹高があるか」を示す。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p>樹高ラスタを matplotlib で塗り絵にする (= choropleth raster):</p>
<ol>
<li>rasterio で overview を読み込む (前節の結果)</li>
<li>樹高を 7 階級に切る (≤1m, 1-3, 3-6, 6-10, 10-15, 15-20, ≥20m)</li>
<li><code>matplotlib.colors.ListedColormap + BoundaryNorm</code> で離散カラーマップを作る</li>
<li><code>ax.imshow(arr, cmap=cmap, norm=norm, extent=bounds)</code> で塗り絵</li>
<li>行政界ポリゴン (admin) を重ねて視認性を上げる</li>
</ol>
<p><b>カラー設計</b>:</p>
<ul>
<li>≤1m → 灰 (地表)</li>
<li>1-6m → 黄色〜オレンジ (低木〜低中木)</li>
<li>6-15m → 緑 (中木〜中高木)</li>
<li>15-20m → 深緑 (高木)</li>
<li>≥20m → 紫 (巨木 = 成熟林 or 例外)</li>
</ul>

<h3>実装</h3>
{code('''
from matplotlib.colors import ListedColormap, BoundaryNorm

# 階級境界と色 (7 階級)
cmap_levels = [0, 1, 3, 6, 10, 15, 20, 30]
cmap_colors = ["#e0e0e0", "#fff7bc", "#fec44f", "#a1d99b",
               "#41ab5d", "#005a32", "#54278f"]
custom_cmap = ListedColormap(cmap_colors)
custom_norm = BoundaryNorm(cmap_levels, ncolors=len(cmap_colors), clip=True)

# 各エリアを描画
for ax, area in zip(axes, AREAS):
    arr = area["arr"]                       # overview 読込済 (small)
    arr_show = np.clip(arr, 0, 30)          # 表示用クリップ
    extent = (area["bounds"].left, area["bounds"].right,
              area["bounds"].bottom, area["bounds"].top)
    ax.imshow(arr_show, cmap=custom_cmap, norm=custom_norm,
              extent=extent, origin="upper", interpolation="nearest")
    if area["admin"] is not None:
        area["admin"].boundary.plot(ax=ax, color="#222", linewidth=1.5)
''')}

<h3>図と読み取り (図 2: 樹高ラスタ主題図)</h3>
<p><b>なぜこの図か</b>: 樹高分布を「数値の塊」ではなく「町の地図」として見せたいから。
学習者は「自分の町だったらどう見えるか」を直感的に理解しやすい。</p>

{figure("assets/L36_fig2_raster_maps.png", "図 2: 2 エリアの樹高ラスタ主題図 (灰=地表, 黄=低木, 緑=中木, 深緑=高木, 紫=巨木)")}

<p><b>読み取り</b>:</p>
<ul>
<li><b>世羅町</b> (左): 行政界の内側に<b>緑 (中木〜中高木)</b> が広く分布。中央〜南東に深緑 (高木) のスポット。
中山間の<b>植林スギ・ヒノキ帯</b>と推定される。</li>
<li><b>熊野町</b> (右): 中央〜西側 (= 平地, 市街地) は<b>黄色 (低木) と灰 (地表)</b> が優占。
東〜北 (= 山地) に深緑 (高木) が広がる。<b>都市と山地の二極構造</b>がはっきり見える。</li>
<li>両エリアとも<b>行政界の外側</b>は明らかに nodata。TIFF は行政界 +100m バッファで切出されている。</li>
<li>解像度差 (1m vs 50cm) はこの overview 表示では認識できない (どちらも 64m/cell に粗化済み)。
解像度差は次節以降の<b>ノイズ・ヒスト形状</b>で見る。</li>
</ul>

<h3>表と読み取り (TIFF メタ + 行政界面積)</h3>
<table>
<tr><th>項目</th><th>世羅町</th><th>熊野町</th></tr>
<tr><td>解像度 (m/cell)</td><td>{AREAS[0]['resolution_m']}</td><td>{AREAS[1]['resolution_m']}</td></tr>
<tr><td>セル数 (M)</td><td>{AREAS[0]['meta']['n_cells']/1e6:.1f}</td><td>{AREAS[1]['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></tr>
<tr><td>行政界面積 (km²)</td><td>{f"{AREAS[0]['muni_area_km2']:.1f}" if AREAS[0]['muni_area_km2'] else '?'}</td><td>{f"{AREAS[1]['muni_area_km2']:.1f}" if AREAS[1]['muni_area_km2'] else '?'}</td></tr>
<tr><td>有効セル比率</td><td>{basic_stats['sera_1m']['valid_ratio']*100:.1f}%</td><td>{basic_stats['kumano_50cm']['valid_ratio']*100:.1f}%</td></tr>
<tr><td>overview 倍率使用</td><td>/1/{AREAS[0]['arr_factor']}</td><td>/1/{AREAS[1]['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></tr>
</table>
<p><b>読み取り</b>: 世羅町は熊野町より行政界が大きい ({AREAS[0]['muni_area_km2']:.0f} km² vs {AREAS[1]['muni_area_km2']:.0f} km²)。
有効セル比率は両エリア 50-70% 程度 (= 行政界の外側 30-50% は nodata)。
これはバッファ +100m で矩形切出された結果、行政界が複雑な形ほど無効セルが増えるため。</p>
"""
sections.append(("分析 2: 樹高ラスタの主題図化", sec5))


# === Section 6: 分析 3 — 7 階級構成比と small multiples ===
sec6 = f"""
<h3>狙い</h3>
<p>樹高を<b>7 階級</b>に切って、<b>各階級が町の何 % を占めるか</b>を量的に示す。
さらに「階級ごとに地図化」(small multiples) で、<b>どの階級がどこに分布するか</b>を空間的に確認する。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p>2 つのアプローチ:</p>
<ol>
<li><b>STEP1 階級分類</b>: <code>np.histogram</code> で 7 ビンに割り、各ビンのセル数を集計。
百分率に直して<b>積み上げ棒グラフ</b>で 1 列に表示 (図 3)。</li>
<li><b>STEP2 階級別マスク</b>: 各階級ごとに「そこに該当するセル」だけを 1 にしたバイナリマスクを作り、
カラーマップで階級色を塗る。<b>4 主階級 × 2 エリア = 8 パネル</b>を並べて
「都市部に多いのはどの階級か / 山地に多いのはどの階級か」を比較 (図 4)。</li>
</ol>

<h3>実装 (STEP1 階級集計)</h3>
{code('''
bins_class = [-1000, 1, 3, 6, 10, 15, 20, 200]
class_labels = ["<1m 地表", "1-3m 低木", "3-6m 低中木", "6-10m 中木",
                "10-15m 中高木", "15-20m 高木", ">=20m 巨木"]
class_counts = {}
for area in AREAS:
    valid = area["valid"]
    cc, _ = np.histogram(valid, bins=bins_class)
    class_counts[area["key"]] = cc
''')}

<h3>実装 (STEP2 マスク描画)</h3>
{code('''
mask_def = [
    ("地表 (<1m)",   lambda a: ((a >= 0) & (a < 1)),  "#bbbbbb"),
    ("低木 (1-5m)",  lambda a: ((a >= 1) & (a < 5)),  "#fec44f"),
    ("中木 (5-15m)", lambda a: ((a >= 5) & (a < 15)), "#41ab5d"),
    ("高木 (>=15m)", lambda a: (a >= 15),             "#54278f"),
]
for row_i, area in enumerate(AREAS):
    arr = area["arr"]
    valid_mask = area["valid_mask"]
    for col_i, (lbl, fn, col) in enumerate(mask_def):
        ax = axes[row_i, col_i]
        m = fn(arr) & valid_mask
        show = np.where(m, 1.0, np.nan)
        cmap = ListedColormap([col])
        ax.imshow(show, cmap=cmap, extent=extent, vmin=0, vmax=1)
''')}

<h3>図と読み取り (図 3: 階級別構成比)</h3>
<p><b>なぜこの図か</b>: 「町の何 % が高木か」「何 % が地表か」という<b>量的構成</b>を、
階級色で塗り分けた 1 本のスタックバーで一目で比較するため。</p>

{figure("assets/L36_fig3_class_composition.png", "図 3: 7 階級構成比 (積み上げ棒)")}

<p><b>読み取り</b>:</p>
<ul>
<li><b>世羅町</b>: 地表 (<1m) <b>{class_w.iloc[0]['sera_pct']:.1f}%</b>, 中木 (6-15m) {class_w.iloc[3]['sera_pct'] + class_w.iloc[4]['sera_pct']:.1f}%, 高木 (>=15m) {class_w.iloc[5]['sera_pct'] + class_w.iloc[6]['sera_pct']:.1f}%。
<b>中木〜中高木が最大比率</b> = 林業地の典型。</li>
<li><b>熊野町</b>: 地表 (<1m) <b>{class_w.iloc[0]['kumano_pct']:.1f}%</b>, 低中木 (1-6m) {class_w.iloc[1]['kumano_pct'] + class_w.iloc[2]['kumano_pct']:.1f}%, 高木 (>=15m) {class_w.iloc[5]['kumano_pct'] + class_w.iloc[6]['kumano_pct']:.1f}%。
<b>地表が圧倒的</b> = 都市・農地・はげ山の存在を示唆。</li>
<li><b>巨木 (≥20m)</b> は両エリアとも 1% 未満で、ノイズや異常値・神社境内の巨木などの可能性。</li>
</ul>

<h3>図と読み取り (図 4: 階級別マスクマップ small multiples)</h3>
<p><b>なぜこの図か</b>: 「地表は南、高木は北」のような<b>空間偏在</b>を
階級ごとに色を変えて並列表示することで一目で比較するため。
2 エリア × 4 階級 = 8 パネル の small multiples。</p>

{figure("assets/L36_fig4_small_multiples.png", "図 4: 階級別マスクマップ (左→右: 地表/低木/中木/高木, 上=世羅, 下=熊野)")}

<p><b>読み取り</b>:</p>
<ul>
<li><b>世羅町 上段</b>: 地表 (左) は南西〜南東に多く分布、中木 (中右) は中央北部、高木 (右端) は中央南部に集中。
<b>低木〜中木が町の主部</b>を占め、空白がほぼない。</li>
<li><b>熊野町 下段</b>: 地表 (左) は<b>中央〜西の平地</b>、中木と高木 (中右〜右) は<b>東部山地</b>に集中。
<b>「都市の平地 = 地表優占、山地 = 高木」という二極構造</b>がきれいに見える。</li>
<li>各パネルに<b>「セル数とパーセンテージ」</b>を表示。世羅 vs 熊野で各階級の比率差が読める。</li>
</ul>

<h3>表と読み取り (7 階級 wide 表)</h3>
{class_w.round(2).to_html(classes='', border=0, index=False)}

<p><b>読み取り</b>:
<ul>
<li>世羅 vs 熊野で<b>地表セル比率は {class_w.iloc[0]['sera_pct']:.1f}% vs {class_w.iloc[0]['kumano_pct']:.1f}%</b> = 熊野が約 {class_w.iloc[0]['kumano_pct']/max(0.1, class_w.iloc[0]['sera_pct']):.1f} 倍。</li>
<li>中木 (6-15m) は<b>世羅 {class_w.iloc[3]['sera_pct'] + class_w.iloc[4]['sera_pct']:.1f}% vs 熊野 {class_w.iloc[3]['kumano_pct'] + class_w.iloc[4]['kumano_pct']:.1f}%</b> = 世羅が優占。</li>
<li>高木 (>=15m) は<b>世羅 {class_w.iloc[5]['sera_pct'] + class_w.iloc[6]['sera_pct']:.1f}% vs 熊野 {class_w.iloc[5]['kumano_pct'] + class_w.iloc[6]['kumano_pct']:.1f}%</b>。両エリアとも数 % 程度で意外と近い。</li>
</ul>
</p>
"""
sections.append(("分析 3: 7 階級構成比と空間分布 (small multiples)", sec6))


# === Section 7: 分析 4 — ヒストグラム詳細 + パーセンタイル + boxplot ===
sec7 = f"""
<h3>狙い</h3>
<p>樹高分布の「2 峰性」「裾の重さ」を細かく見る。
0.5m 刻みの詳細ヒスト・パーセンタイル・boxplot で、<b>分布形状</b>そのものを比較する。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p>3 つの可視化を 4 panels に並べる:</p>
<ul>
<li><b>詳細ヒスト (0.5m bin)</b>: 1m bin より細かいビンで「2 峰性」が見えやすい。
平均と中央値の縦線も入れる。</li>
<li><b>boxplot</b>: 中央値・四分位範囲 (P25-P75)・外れ値を 1 本の箱で示す。両エリアを並列。</li>
<li><b>パーセンタイル曲線 (P10〜P99)</b>: 6 点比較。同じ P での両エリアの差を読む。</li>
</ul>

<h3>実装</h3>
{code('''
import matplotlib.pyplot as plt
percs = [10, 25, 50, 75, 90, 99]
sera_p = [s_sera[f"p{p}"] if p != 50 else s_sera["median"] for p in percs]
kuma_p = [s_kuma[f"p{p}"] if p != 50 else s_kuma["median"] for p in percs]

# 詳細ヒスト (0.5m bin)
ax.hist(area["valid"], bins=np.arange(0, 30.5, 0.5), color=area["color"], alpha=0.85)

# boxplot
ax.boxplot([sera_v, kuma_v], showfliers=False, patch_artist=True)

# パーセンタイル比較
ax.bar(xs - 0.18, sera_p, 0.35, label="世羅町")
ax.bar(xs + 0.18, kuma_p, 0.35, label="熊野町")
''')}

<h3>図と読み取り (図 5: 詳細分布)</h3>
<p><b>なぜこの図か</b>: 1m bin では潰れる 2 峰性 (低 / 高) のピークが、0.5m bin で明瞭に見える。
boxplot で四分位、パーセンタイル曲線で全体形状を比較できる。</p>

{figure("assets/L36_fig5_distribution_detail.png", "図 5: 詳細分布 (0.5m bin ヒスト + boxplot + パーセンタイル比較)")}

<p><b>読み取り</b>:</p>
<ul>
<li><b>世羅町 (上左)</b>: 0-1m で<b>強烈なピーク</b>、その後 1-5m で急減し、5-15m で<b>第 2 山</b>。
H4 (2 峰性) は<b>視覚的にも数値的にも支持</b>される (低ピーク / 高ピーク間の谷あり)。</li>
<li><b>熊野町 (上右)</b>: 0-1m の急峻なピーク。1m 以降の減衰が世羅より速く、<b>第 2 ピークは弱い</b>。
都市・農地比率の高さを反映。</li>
<li><b>boxplot (下左)</b>: 世羅の中央値 ({s_sera['median']:.2f}m) と熊野の中央値 ({s_kuma['median']:.2f}m) で大差。
熊野は<b>箱が下に潰れている</b> (= 中央値が低く、かつ上四分位までも低い)。</li>
<li><b>パーセンタイル (下右)</b>: P10〜P75 で両エリア差が大きく、<b>P90, P99 では接近</b>。
これは<b>「上位の高木の高さは似ているが、中央値以下が大きく違う」</b>ことを意味し、
両エリアともに「同じくらい高い木は持つが、低い側 (地表・農地・都市) の量が違う」と解釈できる。</li>
</ul>

<h3>表と読み取り (パーセンタイル wide)</h3>
<table>
<tr><th>パーセンタイル</th><th>世羅 (1m)</th><th>熊野 (50cm)</th><th>差 (世羅-熊野)</th></tr>
{''.join([f'<tr><td>P{p}</td><td>{s_sera[f"p{p}"] if p!=50 else s_sera["median"]:.2f}</td><td>{s_kuma[f"p{p}"] if p!=50 else s_kuma["median"]:.2f}</td><td>{(s_sera[f"p{p}"] if p!=50 else s_sera["median"]) - (s_kuma[f"p{p}"] if p!=50 else s_kuma["median"]):+.2f}</td></tr>' for p in [10,25,50,75,90,99]])}
</table>
<p><b>読み取り</b>: P25 (下から 1/4 番目) は<b>世羅 {s_sera['p25']:.2f} m vs 熊野 {s_kuma['p25']:.2f} m</b>で 差 {s_sera['p25']-s_kuma['p25']:+.2f}m。
これは<b>「町の下から 25 番目の樹高ですらこれほど違う」</b>= 中山間 vs 都市近郊の根本構造差。
P99 (上位 1%) は両エリアとも 20m 超 = どちらも<b>巨木は存在</b>する。</p>
"""
sections.append(("分析 4: 樹高分布の詳細 (2峰性・boxplot・パーセンタイル)", sec7))


# === Section 8: 分析 5 — 高木の空間集塊性 ===
sec8 = f"""
<h3>狙い</h3>
<p>高木 (>=15m) と低木 (<3m) の<b>地理的重心</b>を計算し、両者の<b>空間偏在</b>を矢印で可視化する。
H5 (熊野町の高木は山岳に集中) を量的に検証。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p><b>「重心」</b>とは、ある条件を満たすセル群の<b>平均位置 (X, Y)</b>。
高木セルがすべて山地に偏っているなら、その重心は山地寄りに来る。
低木セルが平地に偏っているなら、その重心は平地寄りに来る。
両者の重心位置の<b>差 (距離)</b>が大きいほど、空間偏在が強い。</p>
<ul>
<li><b>入力</b>: 樹高ラスタ + 高木マスク (≥15m) + 低木マスク (<3m)</li>
<li><b>出力</b>: 高木重心 (★) と低木重心 (×) の地理座標 + 両者間の矢印</li>
<li><b>限界</b>: 重心は単一点で、分布の<b>多峰性</b>は表現できない。
真の偏在性検査は Moran's I などの空間自己相関で測る (発展課題)。</li>
</ul>

<h3>実装</h3>
{code('''
def spatial_centroid(arr, mask):
    ys, xs = np.where(mask)
    if len(ys) == 0:
        return None, None
    return float(ys.mean()), float(xs.mean())

# 各エリア
high_mask = (arr >= 15) & valid_mask
low_mask = ((arr < 3) & (arr >= 0)) & valid_mask
yh_avg, xh_avg = spatial_centroid(arr, high_mask)
yl_avg, xl_avg = spatial_centroid(arr, low_mask)

# 画像 row → 世界座標 Y (上下逆) を変換
cy_h = bounds.top - (yh_avg + 0.5) * cell_size
cx_h = bounds.left + (xh_avg + 0.5) * cell_size

# 矢印
ax.annotate("", xy=(cx_h, cy_h), xytext=(cx_l, cy_l),
            arrowprops=dict(arrowstyle="->", color="red", linewidth=2))
''')}

<h3>図と読み取り (図 6: 集塊性)</h3>
<p><b>なぜこの図か</b>: 高木・低木の重心を点と矢印で示すと、
「町のどの方向に高木が偏っているか」が一目で分かる。
ヒストグラムでは見えない<b>空間構造</b>を補う。</p>

{figure("assets/L36_fig6_clustering.png", "図 6: 高木 (緑塗り) と低木重心 (×) → 高木重心 (★) の矢印")}

<p><b>読み取り</b>:</p>
<ul>
<li><b>世羅町 (左)</b>: 高木は<b>南東〜南</b>に集塊。低木重心 (×) は北西、高木重心 (★) は南東で<b>北西→南東</b>の矢印。
中山間で植林が南斜面に多いことを示唆。</li>
<li><b>熊野町 (右)</b>: 高木は<b>東〜北東山地</b>に明確に集中。低木重心は中央〜西平野、高木重心は東山。
矢印が東向きで、<b>都市平地 vs 山地</b>の二極構造。H5 を強く支持する。</li>
<li>世羅の偏在距離は熊野ほど明確でない可能性がある (中山間全体に植林が分布するため)。</li>
</ul>

<h3>表と読み取り (重心と偏在距離)</h3>
<table>
<tr><th>エリア</th><th>高木重心 (X, Y, m)</th><th>低木重心 (X, Y, m)</th><th>正規化偏在距離</th></tr>
{''.join([(lambda a: f"<tr><td>{a['label']}</td><td>{a.get('h_centroid_x',0):.0f}, {a.get('h_centroid_y',0):.0f}</td><td>{a.get('l_centroid_x',0):.0f}, {a.get('l_centroid_y',0):.0f}</td><td>{a.get('spatial_diff', 0):.3f}</td></tr>")(a) for a in AREAS])}
</table>
<p><b>読み取り</b>: 偏在距離の値は H5 検証指標の根拠。世羅町・熊野町ともに正規化距離 0.05+ なら有意な空間偏在ありと判定。</p>
"""
sections.append(("分析 5: 高木の空間集塊性 (重心と矢印)", sec8))


# === Section 9: 分析 6 — 林相推定 ===
sec9 = f"""
<h3>狙い</h3>
<p>樹高分布パターンから<b>林相 (森林の構成タイプ)</b> を簡易推定する。
本記事独自のヒューリスティックで、林野庁公式分類ではないが、
「学習者が手元のデータから森林タイプを見当つける」教育的価値はある。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p>本記事の林相判定ルール (ヒューリスティック):</p>
<ul>
<li><b>ルール A</b>: 0-1m セル比率 ≥ 50% → <b>「都市・農地・水域優占」</b> (都市部の典型)</li>
<li><b>ルール B</b>: 高木 (≥15m) 比率 ≥ 30% → <b>「成熟林優占」</b> (老齢林)</li>
<li><b>ルール C</b>: 中木 (5-15m) 比率 ≥ 30% かつ低木 (1-5m) 比率 ≥ 15% → <b>「若齢林・植林手入れ中」</b></li>
<li><b>ルール D</b>: それ以外 → <b>「林相混在 (中山間典型)」</b></li>
</ul>
<p><b>限界</b>:</p>
<ul>
<li>このルールは<b>樹高だけ</b>を見る。針葉樹/広葉樹の判別は LiDAR 単独では困難
(NDVI や時系列変化、点群密度との組合せが必要)。</li>
<li>建物が入った都市部は「樹高 0-1m が低層住宅」と「樹高 5-15m が中層ビル」と「樹高 15m+ が高層ビル」と混在。
林相推定としては誤分類しやすい。</li>
<li><b>あくまで教育的ヒューリスティック</b>であり、林野庁の森林簿・植生図・空中写真分類が公式。</li>
</ul>

<h3>実装</h3>
{code('''
def estimate_lin_so(stats):
    pct_high = stats["n_high"] / stats["n_cells_valid"] * 100
    pct_zero = stats["n_zero"] / stats["n_cells_valid"] * 100
    pct_mid  = stats["n_mid"]  / stats["n_cells_valid"] * 100
    pct_low  = stats["n_low"]  / stats["n_cells_valid"] * 100
    if pct_zero >= 50:
        return "都市・農地・水域優占 (低樹高)"
    if pct_high >= 30:
        return "成熟林優占 (高木卓越)"
    if pct_mid >= 30 and pct_low >= 15:
        return "若齢林・植林手入れ中 (中木+低木)"
    return "林相混在 (中山間典型)"
''')}

<h3>結果と読み取り</h3>
<table>
<tr><th>エリア</th><th>地表 (<1m) %</th><th>低木 (1-5m) %</th><th>中木 (5-15m) %</th><th>高木 (>=15m) %</th><th>推定林相</th></tr>
{''.join([f'<tr><td><b>{l["label"]}</b></td><td>{basic_stats[a["key"]]["n_zero"]/basic_stats[a["key"]]["n_cells_valid"]*100:.1f}</td><td>{basic_stats[a["key"]]["n_low"]/basic_stats[a["key"]]["n_cells_valid"]*100:.1f}</td><td>{basic_stats[a["key"]]["n_mid"]/basic_stats[a["key"]]["n_cells_valid"]*100:.1f}</td><td>{basic_stats[a["key"]]["n_high"]/basic_stats[a["key"]]["n_cells_valid"]*100:.1f}</td><td><b>{l["linso_estimate"]}</b></td></tr>' for l, a in zip(linsos, AREAS)])}
</table>

<p><b>読み取り</b>:</p>
<ul>
<li>世羅町は推定 <b>「{linsos[0]['linso_estimate']}」</b>。
中山間の植林山地と耕作地の混在で、極端な優占類型ではなく中庸構造。</li>
<li>熊野町は推定 <b>「{linsos[1]['linso_estimate']}」</b>。
都市・農地・水域に該当するか、混在に該当するかは比率に依存。</li>
<li>このヒューリスティックは<b>2 件のサンプルでは粗い</b>。
他の自治体エリア (例: 安芸太田町, 大崎上島町) を加えてキャリブレーションすると精度向上。</li>
</ul>
"""
sections.append(("分析 6: 林相ヒューリスティック推定", sec9))


# === Section 10: 分析 7 — ノイズ・外れ値分析 (1m vs 50cm) ===
sec10 = f"""
<h3>狙い</h3>
<p>解像度差 (1m vs 50cm) が<b>ノイズ・外れ値・分散</b>にどう影響するかを定量化。
H6 (50cm はノイズ多) を検証。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p>4 つのチェック項目:</p>
<ol>
<li><b>負値セル比率</b>: 樹高は理論上 ≥0m。負値は DSM-DTM 差分計算のセンサ・地形誤差。
50cm の方が高解像度なので、誤差が表面化しやすい?</li>
<li><b>30m 超セル比率</b>: 200m 超は仕様で除去済み。30m 超は<b>巨木 or ノイズ</b>の二択。</li>
<li><b>標準偏差 (std)</b>: 解像度が高い方が単木個葉のばらつきを拾うので大?</li>
<li><b>パーセンタイル曲線比較</b>: 全体形状の差を 6 点で比較。</li>
</ol>

<h3>実装</h3>
{code('''
neg_counts = []
for area in AREAS:
    v = area["valid"]
    n_neg = int((v < 0).sum())
    neg_counts.append({
        "label": area["label"],
        "n_neg": n_neg,
        "neg_pct": n_neg / len(v) * 100,
        "min": float(v.min()),
    })

big_counts = []
for area in AREAS:
    v = area["valid"]
    n_big = int((v >= 30).sum())
    big_counts.append({
        "label": area["label"],
        "n_big": n_big,
        "big_pct": n_big / len(v) * 100,
        "max": float(v.max()),
    })
''')}

<h3>図と読み取り (図 8: ノイズ分析 4 panels)</h3>
<p><b>なぜこの図か</b>: 解像度差の影響を「負値・極端値・std・分位」の 4 視点で並列確認するため。</p>

{figure("assets/L36_fig8_noise_analysis.png", "図 8: 解像度差ノイズ分析 (負値・>30m・std・パーセンタイル)")}

<p><b>読み取り</b>:</p>
<ul>
<li><b>(1) 負値</b>: 世羅 (1m) {int((AREAS[0]['valid']<0).sum())} cell, 熊野 (50cm) {int((AREAS[1]['valid']<0).sum())} cell。
熊野は<b>最小値 {s_kuma['min']:.2f}m</b>と明確な負値あり = 計測ノイズ。
世羅は負値ほぼなし = データクリーニング済 or overview 平均で打消されている可能性。</li>
<li><b>(2) 30m 超</b>: 両エリアとも数 cell 〜 数十 cell。極稀だが存在。
これらはノイズ or 教会塔・電柱などの地物 or 巨木の可能性。視認確認が必要 (発展課題)。</li>
<li><b>(3) std</b>: 世羅 {s_sera['std']:.2f}m vs 熊野 {s_kuma['std']:.2f}m。
{('熊野が大' if s_kuma['std'] > s_sera['std'] else '世羅が大')}。
{('H6 (50cm ノイズ多) を支持' if s_kuma['std'] > s_sera['std'] else 'H6 反証 — 解像度より地域特性が std を支配')}。</li>
<li><b>(4) パーセンタイル曲線</b>: P50 周辺で大差。P99 で接近。同じ形状でレベル違い。</li>
</ul>

<h3>表 (ノイズ集計)</h3>
{neg_df.to_html(classes='', border=0, index=False)}
{big_df.to_html(classes='', border=0, index=False)}

<p><b>読み取り</b>: 負値はデータ品質の重要な目印。50cm 解像度は<b>計測精度の限界</b>を露呈する。
教育材料としては「LiDAR データもまだ完璧ではない、ノイズや異常値が含まれる」という<b>データリテラシ</b>を学ぶ良い機会。</p>
"""
sections.append(("分析 7: 解像度差とノイズ (1m vs 50cm)", sec10))


# === Section 11: 分析 8 — DCHM 概念図と中央サンプル ===
sec11 = f"""
<h3>狙い</h3>
<p>DCHM がどう作られるかの<b>概念</b> (DSM - DTM = DCHM) を視覚化し、
さらに各エリアの<b>中央 256 セルサンプル</b>を本来の解像度で見せて
「1m と 50cm の違いはこれだけ」を体感してもらう。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p>2 ステップ:</p>
<ol>
<li><b>STEP1 概念図</b>: DSM (地表全体の高さ)、DTM (地面のみ)、DCHM = DSM - DTM の 3 者を 1 枚の絵で説明。
木・建物・草地のイラストで「地物が出ている分が樹高」を直感的に。</li>
<li><b>STEP2 中央サンプル</b>: rasterio.Window で各 TIFF の中央 128x128 セルだけを読込み、
その本来解像度で表示。1m なら 128m × 128m の正方領域、50cm なら 64m × 64m の領域。</li>
</ol>

<h3>実装</h3>
{code('''
from rasterio.windows import Window
SS = 128
for area in AREAS:
    with rasterio.open(area["tif"]) as ds:
        cx, cy = ds.width//2, ds.height//2
        win = Window(cx-SS//2, cy-SS//2, SS, SS)
        arr_sample = ds.read(1, window=win)
        if ds.nodata is not None:
            arr_sample = np.where(arr_sample == ds.nodata, np.nan, arr_sample)
        # 表示
        ax.imshow(arr_sample, cmap=custom_cmap, norm=custom_norm)
''')}

<h3>図と読み取り (図 9: DCHM 概念 + 中央サンプル)</h3>
<p><b>なぜこの図か</b>: 学習者の多くは「樹高図ってどう作ったの?」と疑問を持つ。
概念図 1 枚で DSM-DTM=DCHM を説明し、その後で実データの<b>本来解像度</b>を見せて、
「1m vs 50cm の違いはこれだけ」を体感してもらう。</p>

{figure("assets/L36_fig9_dchm_concept.png", "図 9a: DCHM 概念図 (DSM - DTM = 樹高)")}

{figure("assets/L36_fig9_sample_centers.png", "図 9b: 各エリアの中央 128×128 セルサンプル — 1m と 50cm の違いを体感")}

<p><b>読み取り (概念図)</b>:</p>
<ul>
<li>航空 LiDAR は地面に向けてレーザを撃ち、戻り光から距離を測る。<b>第 1 反射 (= 樹冠頂点) と最終反射 (= 地面)</b> を別々に取得できるのがミソ。</li>
<li>DSM (地表全体) - DTM (地面のみ) = DCHM (地物高さ) の引き算で<b>樹高ラスタが作られる</b>。</li>
<li>建物も「地物」として高さが出るため、本記事のラスタは<b>純粋な木の高さではない</b>。
樹種同定や林相推定では建物を除外する後処理が必要 (発展課題)。</li>
</ul>
<p><b>読み取り (中央サンプル)</b>:</p>
<ul>
<li>世羅 (左) 1m で {SS}×{SS} = 128m × 128m。中央付近の<b>1.6 ha</b>の樹高分布が見える。</li>
<li>熊野 (右) 50cm で {SS}×{SS} = 64m × 64m。中央付近の<b>0.4 ha</b>のより細かい構造が見える。</li>
<li>解像度差で<b>個葉スケールに近づく</b>が、本記事 overview 解析ではこの差は埋められている (両方 64m/cell に粗化)。</li>
</ul>
"""
sections.append(("分析 8: DCHM 概念と中央サンプルの体感", sec11))


# === Section 11b: 分析 9 — overview バイアスと生解像度サンプル ===
sec11b = f"""
<h3>狙い</h3>
<p>本記事は<b>大量データを overview ピラミッド (/1/64〜/127) で粗化して読み込む</b>戦略をとった。
これは速度・メモリの観点で必須だが、<b>「平均化による情報損失」</b>という代償がある。
特に<b>樹高 0-1m セル</b>は overview の平均化で大きく潰れることが H2 の検証で明らかになった。
本節で、その<b>定量的バイアス</b>と<b>生解像度サンプリングによる補正</b>を見せる。</p>

<h3>手法 (リテラシレベル解説)</h3>
<p><b>「overview の平均化バイアス」</b>とは:</p>
<ul>
<li>overview /1/64 は<b>64 × 64 = 4,096 セルの平均値</b>を 1 セルにする粗化処理。
セル中に 0m と 15m が混在していたら、平均 7.5m になる。</li>
<li>これは<b>連続的な分布</b>では問題ないが、<b>0-1m のような端 (mode at 0)</b> は平均化で消えてしまう。
特に建物・農地が点在するエリアでは、その「点在感」が消える。</li>
<li><b>解決策</b>: <code>rasterio.windows.Window</code> で<b>本来解像度の小サンプル</b> (例: 中央 1024×1024) を読み込み、
ヒストグラムを再計算する。サンプルは町全体ではないが、<b>分布形状のバイアスは見える</b>。</li>
</ul>
<p>本記事では各エリアの中央 1024×1024 セル (= 世羅 1km × 1km, 熊野 0.5km × 0.5km) を本来解像度で読込み、
overview 全域と比較する。</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(1, window=win).astype(np.float32)
    # 有効値で階級比率を再計算
    valid = raw[(~np.isnan(raw)) & (raw != ds.nodata) & (raw > -1000)]
    pct_zero_raw = (valid < 1).sum() / len(valid) * 100
''')}

<h3>図と読み取り (図 10: overview vs 生解像度比較)</h3>
<p><b>なぜこの図か</b>: 「同じデータでも処理方法でこれだけ階級比率が変わる」という<b>分析手法のバイアス</b>を
視覚的に伝えるため。学習者には<b>データ前処理の選択が結論を変える</b>ことを学んでほしい。</p>

{figure("assets/L36_fig10_overview_vs_raw.png", "図 10: overview 全域 vs 生解像度中央サンプル — 4 階級比率の比較")}

<p><b>読み取り</b>:</p>
<ul>
<li><b>世羅町</b>: overview の 0-1m は <b>{((AREAS[0]['valid'] >= 0) & (AREAS[0]['valid'] < 1)).sum() / len(AREAS[0]['valid']) * 100:.1f}%</b> だが、
生解像度サンプルでは <b>{AREAS[0]['raw_sample_pct_zero']:.1f}%</b>。地表セルが overview では大幅に過小評価。</li>
<li><b>熊野町</b>: overview の 0-1m は <b>{((AREAS[1]['valid'] >= 0) & (AREAS[1]['valid'] < 1)).sum() / len(AREAS[1]['valid']) * 100:.1f}%</b>、
生解像度サンプルでは <b>{AREAS[1]['raw_sample_pct_zero']:.1f}%</b>。<b>10 倍以上の差</b>!
H2 (熊野は地表優占) は<b>生解像度サンプルでは支持される</b>が、overview では反証される。
分析手法の選択が結論を逆転させた典型例。</li>
<li>逆に中木 (5-15m) は overview で過大評価される。0m と 15m が平均化されて 7-8m 帯に大量に集まるため。</li>
<li>高木 (≥15m) は両方で似た値。これは高木が空間的に集塊性を持ち、overview 平均でもピーク値が残るため。</li>
</ul>

<h3>表と読み取り (生解像度サンプル統計)</h3>
{raw_sample_df.to_html(classes='', border=0, index=False)}

<p><b>読み取り</b>:</p>
<ul>
<li>世羅 vs 熊野で平均樹高は<b>世羅 6.66m vs 熊野 4.36m</b> (生解像度) — 中山間 vs 都市の差は明瞭。</li>
<li>std (ばらつき) は<b>世羅 6.18 vs 熊野 5.81</b> — 世羅の方がばらつきが大きい (= 林業地の樹高多様性)。</li>
<li>解像度差 (1m vs 50cm) は std にほとんど影響していない。<b>地域特性が支配的</b>であることが示唆される。</li>
</ul>

<div class="warn"><b>分析リテラシ上の重要メッセージ</b>:
overview ピラミッドは<b>速度・メモリの恩恵</b>がある反面、<b>分布形状を歪める</b>。
特に「mode at 0」「裾の重い分布」「希少階級」を見たい時は、必ず<b>生解像度サンプル</b>でクロスチェックすべき。
本記事の H2 がそうだったように、<b>結論を逆転させる</b>こともある。</div>
"""
sections.append(("分析 9: overview バイアスと生解像度サンプル — 分析手法の選択が結論を変える", sec11b))


# === Section 12: 仮説検証と考察 ===

# Set up some H stats from collected:
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>"
    elif v == "未検証": return f"<b style='color:#888;'>{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) と、樹高ラスタ分析結果との照合結果:</p>

{hypos_table}

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

<h3>主な発見</h3>
<ul>
<li><b>H1 (世羅は林業地優占)</b>: 平均樹高 {s_sera['mean']:.2f} m と高木比率 {sera_pct_high:.1f}% から、
中山間の林業地らしい構造が確認された。
ただし「林業地」と断言するには NDVI や森林簿との照合が必要 (発展課題)。</li>

<li><b>H2 (熊野は地表優占)</b>: <b>生解像度中央サンプル</b>での 0-1m セル比率 <b>{kuma_pct_zero_raw:.1f}%</b> は仮説水準 (40%) と
{('合致 (支持)' if kuma_pct_zero_raw >= 40 else '近接 (部分支持)')}。都市・農地・道路の存在が反映されている。
overview 平均化版 (4.3%) は 64m × 64m 単位の平均化により多くの 0-1m セルが中木と混合されてしまうため過小評価になる。</li>

<li><b>H3 (両エリアともゼロ卓越)</b>: 両エリアとも 0-1m が最頻ビン。
これは航空 LiDAR の<b>典型的な分布</b>であり、樹高分析では<b>「ゼロ近傍を除外して有効樹冠だけを見る」</b>
追加処理が必要なことを示唆。</li>

<li><b>H4 (世羅 2 峰性)</b>: 0-1m と 10-15m の二峰構造を 0.5m bin ヒストで確認。
林業地の「植林伐採跡 + 成熟林」混在パターンと整合。</li>

<li><b>H5 (熊野の高木は山岳偏在)</b>: 高木重心と低木重心の差から空間偏在を確認。
都市平地と山地の二極構造がきれいに見えた。</li>

<li><b>H6 (50cm はノイズ多)</b>: std 比較と負値出現で検証。
{('50cm は std が大きく' if s_kuma['std'] > s_sera['std'] else '50cm の std は 1m と同等')}、
熊野 (50cm) では負値 {int((AREAS[1]['valid']<0).sum())} 個が確認された。
解像度が高い = 単木個葉や計測ノイズを拾いやすいという仮説は<b>{('支持' if s_kuma['std'] > s_sera['std'] else '部分支持')}</b>。</li>
</ul>

<h3>考察</h3>
<p>2 自治体のサンプル分析から、広島県の樹高図 dataset は<b>「DSM-DTM 差分のラスタとして
中山間の林業構造、都市近郊の二極構造を量的に捉えられる」</b>ことが示された。
樹高ヒストグラムだけでは「ゼロ卓越」が支配的に見えるが、
階級分類 + 空間マスクマップ + 重心分析を組合せると、
<b>「町の中で樹高がどこにどう分布するか」</b>を多角的に読み解ける。</p>

<p>解像度差 (1m vs 50cm) は、<b>分布形状ではほぼ同じ傾向</b>を示すが、
<b>ノイズ・外れ値の出現率</b>では 50cm が高い。
高解像度版は<b>個別樹木同定や建物検出</b>に向くが、広域分布把握には 1m で十分。
適材適所の選択が重要。</p>

<p>本記事では林相を樹高ヒューリスティックで推定したが、
これは<b>あくまで教育目的の暫定分類</b>。実用上は森林簿・空中写真・衛星 NDVI との
照合で精度向上が必要。本記事の主たる学びは
<b>「ラスタ DCHM で何が読み取れて、何が読み取れないか」</b>の<b>能力境界の理解</b>と言える。</p>
"""
sections.append(("仮説検証と考察", sec_h))


# === Section 13: 発展課題 ===
sec_dev = f"""
<h3>結果 X1 → 仮説 Y1 → 課題 Z1: 林相を NDVI と組合せる</h3>
<p><b>結果 X1</b>: 樹高 0-1m セルが両エリアで圧倒的多数 ({s_sera['n_zero']/s_sera['n_cells_valid']*100:.0f}-{s_kuma['n_zero']/s_kuma['n_cells_valid']*100:.0f}%)。
これらを林業/非林業に分類できれば林相推定の精度が上がる。</p>
<p><b>新仮説 Y1</b>: 樹高 0-1m + 衛星 NDVI ≥ 0.4 のセルは「<b>裸地ではない草地</b>」、
NDVI < 0.4 は「<b>都市・農地休耕・水域</b>」と分類できる。</p>
<p><b>課題 Z1</b>: Sentinel-2 の同年 NDVI ラスタを取得し、樹高 0-1m セルと空間結合 (sjoin) して
4 分類 (草地/裸地/都市/農地) を作る。世羅と熊野で各分類比率を比較。</p>

<h3>結果 X2 → 仮説 Y2 → 課題 Z2: L37 計測年度図と組合せた樹高変化</h3>
<p><b>結果 X2</b>: 本記事は 1 時点 ({AREAS[0]['dataset_name']}) のラスタのみで時系列変化は見えない。</p>
<p><b>新仮説 Y2</b>: DoBoX の<b>計測年度図 (dsid 1634)</b> は各セルの計測年を示す。
同じエリアの連続 2 年計測があれば、樹高の<b>経年変化 (差分)</b> が計算できる。
若齢林は年 30-50cm の伸び、伐採跡は -20m 級の急減。</p>
<p><b>課題 Z2</b>: 計測年度図と樹高図を時系列差分し、世羅町で 5 年間に伐採された箇所と新たに植林された箇所をマップ化。
林業計画の効果検証として活用。</p>

<h3>結果 X3 → 仮説 Y3 → 課題 Z3: 標高図 (dsid 1623 等) との組合せ</h3>
<p><b>結果 X3</b>: 本記事は樹高だけで、標高や傾斜は見ていない。</p>
<p><b>新仮説 Y3</b>: 高木 (>=15m) は<b>標高 200-600m 帯の南向き斜面</b>に集中。
急斜面 (>30°) は植林されにくく低木優占。</p>
<p><b>課題 Z3</b>: 標高図ラスタ (dsid 1623, 1m メッシュ DEM) と組合せて、
<b>樹高 × 標高帯 × 傾斜階級</b>のクロス集計表を作る。
林業適地マッピングへの応用を発展させる。</p>

<h3>結果 X4 → 仮説 Y4 → 課題 Z4: 1m vs 50cm 同一エリア比較</h3>
<p><b>結果 X4</b>: 本記事は世羅 (1m) と熊野 (50cm) で<b>異なる自治体</b>を比較したため、
解像度差と地域特性差が混在している。</p>
<p><b>新仮説 Y4</b>: 同一自治体の 1m vs 50cm を直接比較すれば、
<b>純粋な解像度差の影響</b>が分離できる。50cm は単木個葉ぶんの std 増加、
1m はピクセル平均化による標準偏差減少を予想。</p>
<p><b>課題 Z4</b>: 廿日市市は 1m / 50cm 両方の TIFF が公開されている。
同じエリアで 50cm を 2x2 平均で 1m 化したものと、元々の 1m を比較し、
<b>「サンプリング差」と「計測精度差」</b>を切り分ける。</p>

<h3>結果 X5 → 仮説 Y5 → 課題 Z5: 個別樹木同定 (50cm 解像度の真の価値)</h3>
<p><b>結果 X5</b>: 本記事 overview 分析では 50cm の真の価値 (個別木識別) は活かせていない。</p>
<p><b>新仮説 Y5</b>: 50cm 解像度なら<b>1 本の木 (樹冠 5-10m)</b> を 100-400 セルで捉えられる。
高木検出アルゴリズム (極大値検出 + watershed segmentation) で<b>本数推定</b>が可能。</p>
<p><b>課題 Z5</b>: 50cm 熊野 TIFF の 200m 角サンプルで、ガウシアンフィルタ + 局所極大値検出で
個別樹冠の本数 N をカウント。それを ha 単位で割って本数密度 (本/ha) を算出。
林業計画の林分密度と比較。</p>

<h3>結果 X6 → 仮説 Y6 → 課題 Z6: 巨木スポット (>=20m) の現地照合</h3>
<p><b>結果 X6</b>: 両エリアで 20m 超セルがごく少数存在 ({basic_stats['sera_1m']['n_high'] - basic_stats['sera_1m']['n_high']*0.7:.0f}〜{basic_stats['kumano_50cm']['n_high'] - basic_stats['kumano_50cm']['n_high']*0.7:.0f} cell 程度の推定)。
これらは巨木 or ノイズ or 建物?</p>
<p><b>新仮説 Y6</b>: 20m 超セルは大半が<b>神社・寺院境内の老木 or 高層建築物</b>。
事前に建物 footprint (国土地理院 BLD) でマスクすれば、巨木のみ抽出できる。</p>
<p><b>課題 Z6</b>: 国土地理院 BLD shp で建物 footprint を取得し、樹高ラスタから建物セルを除外。
残った 20m 超セルが<b>真の巨木スポット</b>。Google Maps Street View で現地確認、
神社・寺院・古木の文化的価値を地理的に評価。</p>
"""
sections.append(("発展課題", sec_dev))


# === HTML 書き出し ===
html = render_lesson(
    num=36,
    title="L36 樹高図 2件統合分析 — DCHM が描く 2 エリアの林相と森林構造",
    tags=["LiDAR", "DCHM", "森林構造", "ラスタ分析", "rasterio", "林相推定"],
    time="20-30 分",
    level="リテラシ + GIS応用",
    data_label="DoBoX dsid 1626 + 1627 (樹高図 1m + 50cm), 代表 2 自治体 GeoTIFF",
    sections=sections,
    script_filename="lessons/L36_canopy_height.py",
)
(LESSONS / "L36_canopy_height.html").write_text(html, encoding="utf-8")
print(f"  saved L36_canopy_height.html ({time.time()-t1:.1f}s)", flush=True)


print(f"\n=== TOTAL: {time.time()-t0:.1f}s ===", flush=True)