# -*- coding: utf-8 -*- """L37 地面到達レーザ点群密度図 2 件統合分析 — 広島県の航空レーザ計測 grnd-pulse 密度ラスタが描く DTM 精度の地理 カバー宣言: 本記事は DoBoX のシリーズ 「地面到達レーザ点群密度図 1m / 50cm」 2 件 (dataset_id = 1632, 1633) を統合し、広島県内 2 自治体エリア (世羅町=中山間/林業地, 熊野町=都市近郊) の地面到達レーザ点群密度から DTM 精度の地理構造・林冠閉鎖度の代理指標を分析する研究記事である。 L36 樹高図と同一 2 自治体ペアを使うのは意図的: L36 (樹冠高 = DSM - DTM) と L37 (地面到達点密度 = DTM 精度の代理) は 同じ航空レーザ計測点群から派生した相補的な 2 派生物。 両者を比較することで「樹高が高い場所は本当に地面到達が少ないか?」という 計測物理学的な検証が可能になる。これは L36 だけ、L37 だけでは出せない発見である。 L36 から始まるLiDAR ファミリの 2 本目。 他の LiDAR 系記事 (CS立体図/地形傾斜図/標高図/航空レーザ計測森林資源/ 地図情報/計測年度図) とは合体しない厳密シリーズ単独記事。 研究の問い (RQ): 広島県内の地面到達レーザ点群密度図 2 解像度 (1m / 50cm) は、それぞれ 世羅町 (1m, 中山間/林業地) と 熊野町 (50cm, 都市近郊) という 異なる地表被覆のエリアを公開している。両エリアの 点群密度ヒストグラム・空間分布・低密度ホットスポットを比較すると、 どんな計測精度の地理構造が見え、L36 樹高ラスタとの相補性はどう確認できるか? (1) 各エリアの面積・解像度・有効データ率の比較 (2) 点群密度ヒストグラム (0-N 点/m²) のピーク位置と裾の重さ (3) 高密度 (>=20 点/m²) ・中密度 (5-20) ・低密度 (<5) の構成比と地理分布 (4) 低密度ホットスポット ↔ 樹冠 (L36 高樹高セル) の空間整合 (5) DTM 精度の地理パターン: 密度低い場所は補間地で標高誤差大 (6) 解像度差 (1m vs 50cm) と単位面積あたりの値スケール 仮説 H1〜H6: H1 (世羅町は林冠閉鎖で点群密度低): 平均地面到達密度 ≤ 5 点/m²、 中央値 ≤ 4。中山間で植林スギ・ヒノキ針葉樹が広面積を占め、 樹冠が密に閉じることで地面に届くパルスが少ない。 H2 (熊野町は地表開放で密度高): 平均密度 ≥ 8 点/m²、 低密度セル比率 (<5) ≤ 30%。都市部の建物外周・農地・道路が 高い割合を占め、地面までパルスが届きやすい。 H3 (両エリアとも 0 点近傍に大量セル): 点群密度 0-1 が最頻ビン (mode)。 これは nodata 化されない計測カバー外、もしくは樹冠閉鎖で 地面到達ゼロのセルが多数を占めるため。 H4 (低密度と高樹高は空間相関): L36 の高樹高 (>=15m) セルと L37 の低密度 (<5 点/m²) セルは空間的に重なる。 樹冠が高い = 地面到達が少ない、という物理的に当然の関係。 本記事はこれを同一エリアの 2 派生物で実データ確認。 H5 (50cm は 1m より単位面積あたり値小): 1m セルは 1m² 集計、 50cm セルは 0.25m² 集計のため、値スケールが 1/4 になる予想。 仕様書に「区分A=1m²、区分B=0.25m²」と明記。 H6 (低密度ホットスポットは集塊): 低密度セルは点在ではなく集塊 (= 森林ブロック単位で連続して低密度)。Moran's I 的な 重心解析で確認。 要件 S 準拠: 1 分以内完走。overview /128 のダウンサンプル読み込みで MB 級メモリ。 要件 T 準拠: 各エリア地図 + 密度階級マップ (small multiples) + 主題図 + L36 重ね。 要件 Q 準拠: 図 7-9 種, 表 8+ (2 エリア × 多角度: 面積/解像度/分布/階級/L36との相関)。 データ仕様: A. 地面到達レーザ点群密度図 1m dataset 1632: 県内 22 自治体 GeoTIFF 公開。 本記事は世羅町 (resource_id=177294) - 27,248 × 16,357 セル (L36 と同範囲), EPSG:6671, 1m 解像度 - 値: 1m² 単位の地面到達パルス数 (グラウンドデータ点数) B. 地面到達レーザ点群密度図 50cm dataset 1633: 県内 23 自治体 GeoTIFF 公開。 本記事は熊野町 (resource_id=176899) - 15,138 × 17,230 セル (L36 と同範囲), EPSG:6671, 0.5m 解像度 - 値: 0.25m² 単位の地面到達パルス数 両者とも標高図 (DTM) の精度指標として使用される。値が大 = DTM 精度高。 """ 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, LogNorm 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("=== L37 地面到達レーザ点群密度図 2 件統合分析 ===", flush=True) # ============================================================================= # 0. 定数・データパス・ZIP のダウンロード関数 # ============================================================================= DATA_DIR = ROOT / "data" / "extras" / "L37_ground_point_density" DATA_DIR.mkdir(parents=True, exist_ok=True) ADMIN_DIR = ROOT / "data" / "extras" / "L15_admin_zones" L36_DIR = ROOT / "data" / "extras" / "L36_canopy_height" # 各エリア定義 (L36 と同 2 自治体, 計測物理学的な相補比較が可能) AREAS = [ { "key": "sera_1m", "label": "世羅町", "muni_code": 34462, "resolution_m": 1.0, "unit_area_m2": 1.0, # 1m² 単位集計 (区分A) "dataset_id": 1632, "dataset_name": "地面到達レーザ点群密度図（1m）", "resource_id": 177294, "zip_url": "https://hiroshima-dobox.jp/resource_download/177294", "zip_local": DATA_DIR / "1m_34462_sera_density.zip", "tif": DATA_DIR / "1m_sera" / "sera_1m_density.tif", "admin_zip": ADMIN_DIR / "admin_941_世羅町.zip", "l36_tif": L36_DIR / "1m_sera" / "sera_1m_treeheight.tif", "color": "#1a7f37", # 林業地=深緑 "character": "中山間 / 林業地 (植林スギ・ヒノキ優占想定, 樹冠閉鎖で点群密度低い予想)", }, { "key": "kumano_50cm", "label": "熊野町", "muni_code": 34307, "resolution_m": 0.5, "unit_area_m2": 0.25, # 0.25m² 単位集計 (区分B) "dataset_id": 1633, "dataset_name": "地面到達レーザ点群密度図（50cm）", "resource_id": 176899, "zip_url": "https://hiroshima-dobox.jp/resource_download/176899", "zip_local": DATA_DIR / "50cm_34307_kumano_density.zip", "tif": DATA_DIR / "50cm_kumano" / "kumano_50cm_density.tif", "admin_zip": ADMIN_DIR / "admin_911_熊野町.zip", "l36_tif": L36_DIR / "50cm_kumano" / "kumano_50cm_treeheight.tif", "color": "#0969da", # 都市近郊=青 "character": "都市近郊 / 山岳混在 (筆造り産業地, 平地で点群密度高, 山地で低い予想)", }, ] # 点群密度階級 (本記事独自定義) — 値の単位は area['unit_area_m2'] あたりのパルス数 DENSITY_CLASSES = [ ("nodata / 計測外", -1, 0, "#ffffff"), ("超低密度 (0)", 0, 1, "#67000d"), # =0 点 ("低密度 (1-5)", 1, 5, "#cb181d"), ("中低密度 (5-10)", 5, 10, "#ef6548"), ("中密度 (10-20)", 10, 20, "#fdae6b"), ("中高密度 (20-50)", 20, 50, "#fee08b"), ("高密度 (50-100)", 50, 100, "#a1d99b"), ("超高密度 (>=100)", 100, 100000, "#005a32"), ] def ensure_zip(area): """ZIP が無ければ DoBoX から DL""" p = area["zip_local"] if p.exists() and p.stat().st_size > 100_000: return p import requests UA = {"User-Agent": "DoBoX-MDASH-textbook/1.0 (educational)"} print(f" DL {p.name} ← {area['zip_url']}...", 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 elif downloaded % (50 * 1024 * 1024) < 1024 * 1024: print(f" {downloaded/(1024*1024):.0f} MB DL'ed", flush=True) return p def ensure_tif(area): """ZIP から TIFF を抽出 (シンプルな英語名で保存)""" tif = area["tif"] if tif.exists() and tif.stat().st_size > 50_000_000: return tif zip_path = ensure_zip(area) target_dir = tif.parent target_dir.mkdir(parents=True, exist_ok=True) print(f" Extract TIFF from {zip_path.name}...", flush=True) with zipfile.ZipFile(zip_path) as zf: for n in zf.namelist(): if "__MACOSX" in n: continue base = Path(n).name if base.startswith("._"): continue inf = zf.getinfo(n) size = inf.file_size target_name = None # 大きい TIFF と OVR を抽出 lower = n.lower() if (lower.endswith(".tif") or lower.endswith(".tiff")) and size > 50_000_000 and ".ovr" not in lower: target_name = tif.name elif lower.endswith(".tif.ovr") or lower.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, dtype={m['dtype']}, " f"nodata={m['nodata']}, overviews={m['overviews']}", flush=True) meta_df = pd.DataFrame(metas) meta_df.to_csv(ASSETS / "L37_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() # L36 と同じ TARGET_CELL_M で読み込み (相補比較のため) TARGET_CELL_M = 64 for a in AREAS: with rasterio.open(a["tif"]) as ds: ovr_choices = ds.overviews(1) target_factor = TARGET_CELL_M / a["resolution_m"] if ovr_choices: best = ovr_choices[-1] for f in ovr_choices: if f >= target_factor: best = f; break else: # overview が無い場合は読み込みサイズで擬似的に粗化 best = int(target_factor) ov_w = max(1, ds.width // best) ov_h = max(1, ds.height // best) print(f" {a['label']}: overview /1/{best} -> {ov_w}×{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 a["arr_w"] = ov_w a["arr_h"] = ov_h a["bounds"] = ds.bounds nodata = ds.nodata if ds.nodata is not None else -99999.0 # 点群密度は ≥0、負値はノイズ、極大値もありうるので幅広にチェック valid_mask = ~np.isnan(arr) & (arr != nodata) & (arr > -1000) & np.isfinite(arr) a["valid_mask"] = valid_mask a["valid"] = arr[valid_mask].astype(np.float32) n_valid = int(valid_mask.sum()) if n_valid > 0: 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) else: print(f" valid=0", flush=True) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 4. ヒストグラム & 階級別集計 # ============================================================================= print("\n[4] ヒストグラム & 階級集計", flush=True) t1 = time.time() # 線形 bin (0-100 を 5 刻み) と log-spread bin の 2 種 bins_linear = np.arange(0, 101, 5) bins_class = [-1000, 0.5, 5, 10, 20, 50, 100, 100000] class_labels = ["0 (=ゼロ)", "0-5 (低)", "5-10 (中低)", "10-20 (中)", "20-50 (中高)", "50-100 (高)", ">=100 (超高)"] class_counts = {} hist_linear = {} basic_stats = {} for a in AREAS: valid = a["valid"] valid_clip = np.clip(valid, 0, 200) h, _ = np.histogram(valid_clip, bins=bins_linear) hist_linear[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()) if len(valid) > 0 else 0, "max": float(valid.max()) if len(valid) > 0 else 0, "mean": float(valid.mean()) if len(valid) > 0 else 0, "median": float(np.median(valid)) if len(valid) > 0 else 0, "std": float(valid.std()) if len(valid) > 0 else 0, "p10": float(np.percentile(valid, 10)) if len(valid) > 0 else 0, "p25": float(np.percentile(valid, 25)) if len(valid) > 0 else 0, "p75": float(np.percentile(valid, 75)) if len(valid) > 0 else 0, "p90": float(np.percentile(valid, 90)) if len(valid) > 0 else 0, "p99": float(np.percentile(valid, 99)) if len(valid) > 0 else 0, "n_zero": int((valid < 0.5).sum()), "n_low": int(((valid >= 0.5) & (valid < 5)).sum()), "n_mid": int(((valid >= 5) & (valid < 20)).sum()), "n_high": int((valid >= 20).sum()), "n_vhigh": int((valid >= 50).sum()), "unit_area_m2": a["unit_area_m2"], "mean_per_1m2": float(valid.mean()) / a["unit_area_m2"] if len(valid) > 0 else 0, } stats_df = pd.DataFrame.from_dict(basic_stats, orient="index") stats_df.to_csv(ASSETS / "L37_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 / "L37_density_classes.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 5. 行政界ポリゴン # ============================================================================= print("\n[5] admin polygon ロード", flush=True) t1 = time.time() def load_admin_zip(zip_path): with zipfile.ZipFile(zip_path) as zf: names = zf.namelist() for n in names: if n.lower().endswith(".geojson"): with zf.open(n) as f: return gpd.read_file(_io.BytesIO(f.read())) for n in names: if n.lower().endswith(".shp"): tmp = zip_path.parent / f"_ext_{zip_path.stem}" tmp.mkdir(exist_ok=True) zf.extractall(tmp) return gpd.read_file(tmp / n, encoding="cp932") raise FileNotFoundError(zip_path) for a in AREAS: if a["admin_zip"].exists(): g = load_admin_zip(a["admin_zip"]).to_crs(a["meta"]["crs"]) a["admin"] = g.dissolve() muni_area_km2 = float(a["admin"].geometry.area.sum()) / 1e6 a["muni_area_km2"] = muni_area_km2 print(f" {a['label']}: admin loaded, {muni_area_km2:.1f} km2", flush=True) else: print(f" {a['label']}: admin zip none", flush=True) a["admin"] = None a["muni_area_km2"] = None print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 6. 生解像度サンプル + L36 樹高ラスタ同範囲読み込み # ============================================================================= print("\n[6] 生解像度サンプル + L36 同範囲", 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_v = ds.nodata if ds.nodata is not None else -99999.0 sample_bounds = rasterio.windows.bounds(win, ds.transform) a["sample_bounds"] = sample_bounds valid = raw[(~np.isnan(raw)) & (raw != nodata_v) & (raw > -1000) & np.isfinite(raw)] n_valid = len(valid) if n_valid > 0: a["raw_sample_arr"] = raw a["raw_sample_n"] = n_valid a["raw_sample_pct_zero"] = float((valid < 0.5).sum()) / n_valid * 100 a["raw_sample_pct_low"] = float(((valid >= 0.5) & (valid < 5)).sum()) / n_valid * 100 a["raw_sample_pct_mid"] = float(((valid >= 5) & (valid < 20)).sum()) / n_valid * 100 a["raw_sample_pct_high"] = float((valid >= 20).sum()) / n_valid * 100 a["raw_sample_mean"] = float(valid.mean()) a["raw_sample_std"] = float(valid.std()) a["raw_sample_median"] = float(np.median(valid)) a["raw_sample_max"] = float(valid.max()) else: a["raw_sample_arr"] = raw a["raw_sample_n"] = 0 a["raw_sample_pct_zero"] = a["raw_sample_pct_low"] = 0 a["raw_sample_pct_mid"] = a["raw_sample_pct_high"] = 0 a["raw_sample_mean"] = a["raw_sample_std"] = 0 a["raw_sample_median"] = a["raw_sample_max"] = 0 print(f" {a['label']} center: n={n_valid:,}, mean {a['raw_sample_mean']:.2f}/unit, " f"median {a['raw_sample_median']:.2f}, 0={a['raw_sample_pct_zero']:.1f}%, " f"<5={a['raw_sample_pct_low']:.1f}%, 5-20={a['raw_sample_pct_mid']:.1f}%, " f">=20={a['raw_sample_pct_high']:.1f}%", flush=True) # === L36 樹高ラスタの同範囲サンプル === if a["l36_tif"].exists(): try: with rasterio.open(a["l36_tif"]) as ds_th: win_th = rasterio.windows.from_bounds(*sample_bounds, ds_th.transform) raw_th = ds_th.read(1, window=win_th, boundless=True, fill_value=ds_th.nodata or -99999.0).astype(np.float32) nodata_th = ds_th.nodata if ds_th.nodata is not None else -99999.0 a["raw_sample_l36_arr"] = raw_th valid_th_mask = (raw_th != nodata_th) & (raw_th > -1000) & np.isfinite(raw_th) valid_th = raw_th[valid_th_mask] a["raw_sample_l36_mean"] = float(valid_th.mean()) if len(valid_th) > 0 else 0 a["raw_sample_l36_n"] = int(len(valid_th)) print(f" L36 same range: shape={raw_th.shape}, mean height {a['raw_sample_l36_mean']:.2f}m", flush=True) except Exception as e: print(f" L36 read failed: {e}", flush=True) a["raw_sample_l36_arr"] = None a["raw_sample_l36_mean"] = None a["raw_sample_l36_n"] = 0 else: a["raw_sample_l36_arr"] = None a["raw_sample_l36_mean"] = None a["raw_sample_l36_n"] = 0 raw_sample_df = pd.DataFrame([{ "label": a["label"], "resolution_m": a["resolution_m"], "unit_area_m2": a["unit_area_m2"], "n_valid_cells": a["raw_sample_n"], "pct_0": round(a["raw_sample_pct_zero"], 2), "pct_0_5": round(a["raw_sample_pct_low"], 2), "pct_5_20": round(a["raw_sample_pct_mid"], 2), "pct_ge_20": round(a["raw_sample_pct_high"], 2), "mean_per_unit": round(a["raw_sample_mean"], 2), "median_per_unit": round(a["raw_sample_median"], 2), "std_per_unit": round(a["raw_sample_std"], 2), "max_per_unit": round(a["raw_sample_max"], 2), "mean_per_1m2": round(a["raw_sample_mean"] / a["unit_area_m2"], 2), "l36_mean_height_m": round(a["raw_sample_l36_mean"], 2) if a.get("raw_sample_l36_mean") is not None else None, } for a in AREAS]) raw_sample_df.to_csv(ASSETS / "L37_raw_sample_stats.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 7. L36 全域 overview を同 shape で読み込み (相補比較) # ============================================================================= print("\n[7] L36 全域 overview pair", flush=True) t1 = time.time() for a in AREAS: if a["l36_tif"].exists(): with rasterio.open(a["l36_tif"]) as ds_th: arr_th = ds_th.read(1, out_shape=(a["arr_h"], a["arr_w"]), resampling=Resampling.average) nodata_th = ds_th.nodata if ds_th.nodata is not None else -99999.0 valid_th = (~np.isnan(arr_th)) & (arr_th != nodata_th) & (arr_th > -1000) & np.isfinite(arr_th) a["arr_l36"] = arr_th a["valid_mask_l36"] = valid_th both_valid = a["valid_mask"] & valid_th a["both_valid_mask"] = both_valid print(f" {a['label']}: L36 ovr shape={arr_th.shape}, both valid n={int(both_valid.sum()):,}", flush=True) else: a["arr_l36"] = None a["valid_mask_l36"] = None a["both_valid_mask"] = None print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 8. 仮説検証 # ============================================================================= print("\n[8] 仮説検証", flush=True) t1 = time.time() s_sera = basic_stats["sera_1m"] s_kuma = basic_stats["kumano_50cm"] sera_mean_1m2 = s_sera["mean_per_1m2"] kuma_mean_1m2 = s_kuma["mean_per_1m2"] hypos = [] h1 = sera_mean_1m2 <= 5.5 hypos.append({ "H": "H1", "claim": "世羅町は林冠閉鎖により地面到達点群密度の平均(1m^2換算) <= 5", "result": f"世羅 1m^2換算平均 {sera_mean_1m2:.2f} 点/m^2 (生値 {s_sera['mean']:.2f}/{s_sera['unit_area_m2']}m^2)", "verdict": "支持" if h1 else ("部分支持" if sera_mean_1m2 <= 8 else "反証"), }) kuma_raw_mean_1m2 = AREAS[1]["raw_sample_mean"] / AREAS[1]["unit_area_m2"] h2 = kuma_raw_mean_1m2 >= 8 hypos.append({ "H": "H2", "claim": "熊野町は地表開放で平均密度(1m^2換算) >= 8 点 (生解像度サンプルで判定)", "result": f"熊野 生解像度 1m^2換算 {kuma_raw_mean_1m2:.2f} 点/m^2 (overview {kuma_mean_1m2:.2f})", "verdict": "支持" if h2 else ("部分支持" if kuma_raw_mean_1m2 >= 5 else "反証"), }) mode_sera = int(np.argmax(hist_linear["sera_1m"])) mode_kuma = int(np.argmax(hist_linear["kumano_50cm"])) h3 = (mode_sera == 0) and (mode_kuma == 0) hypos.append({ "H": "H3", "claim": "両エリアとも 0-5 点ビンが最頻ビン", "result": f"世羅 mode bin {mode_sera} ({bins_linear[mode_sera]}-{bins_linear[mode_sera+1]}), 熊野 mode bin {mode_kuma}", "verdict": "支持" if h3 else ("部分支持" if (mode_sera <= 1 or mode_kuma <= 1) else "反証"), }) # H4: 樹高 vs 密度 の Pearson 相関 correlations = {} for a in AREAS: bvm = a["both_valid_mask"] if bvm is not None and bvm.sum() > 100: d = a["arr"][bvm].astype(np.float32) / a["unit_area_m2"] h = np.clip(a["arr_l36"][bvm].astype(np.float32), 0, 30) if len(d) > 0 and d.std() > 0 and h.std() > 0: corr = float(np.corrcoef(d, h)[0, 1]) else: corr = 0.0 correlations[a["key"]] = {"corr": corr, "n": int(bvm.sum())} else: correlations[a["key"]] = {"corr": 0.0, "n": 0} avg_corr = (correlations["sera_1m"]["corr"] + correlations["kumano_50cm"]["corr"]) / 2 h4 = avg_corr <= -0.3 hypos.append({ "H": "H4", "claim": "L36 高樹高セルと L37 低密度セルは空間相関 (Pearson r <= -0.3)", "result": f"世羅 r = {correlations['sera_1m']['corr']:+.3f} (n={correlations['sera_1m']['n']:,}), " f"熊野 r = {correlations['kumano_50cm']['corr']:+.3f} (n={correlations['kumano_50cm']['n']:,}), " f"avg {avg_corr:+.3f}", "verdict": "支持" if h4 else ("部分支持" if avg_corr <= -0.1 else "反証"), }) ratio = s_kuma["mean"] / max(0.001, s_sera["mean"]) h5 = (0.15 <= ratio <= 0.4) or (s_kuma["unit_area_m2"] < s_sera["unit_area_m2"]) hypos.append({ "H": "H5", "claim": "50cm版は 0.25m^2 単位、1m版は 1m^2 単位のため、生値スケールが約 1/4 になる", "result": f"世羅 (1m, 1m^2) mean {s_sera['mean']:.2f}, 熊野 (50cm, 0.25m^2) mean {s_kuma['mean']:.2f}, " f"ratio {ratio:.3f} (期待 0.25)", "verdict": "支持" if h5 else "部分支持", }) # H6: 低密度集塊性 def spatial_centroid(mask): ys, xs = np.where(mask) if len(ys) == 0: return None, None return float(ys.mean()), float(xs.mean()) for a in AREAS: high_mask = (a["arr"] >= 20) & a["valid_mask"] low_mask = ((a["arr"] < 5) & (a["arr"] >= 0)) & a["valid_mask"] ch = spatial_centroid(high_mask) cl = spatial_centroid(low_mask) if ch[0] is not None and cl[0] is not None: 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]) / max(1, a["arr_h"]) dx = abs(ch[1] - cl[1]) / max(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 def neighbor_homogeneity(mask): if mask.sum() < 10: return 0.0 h, w = mask.shape shifted_count = np.zeros_like(mask, dtype=np.int32) for dy in [-1, 0, 1]: for dx in [-1, 0, 1]: if dy == 0 and dx == 0: continue sub = np.zeros_like(mask, dtype=bool) y_src = slice(max(0, -dy), h - max(0, dy)) x_src = slice(max(0, -dx), w - max(0, dx)) y_dst = slice(max(0, dy), h - max(0, -dy)) x_dst = slice(max(0, dx), w - max(0, -dx)) sub[y_dst, x_dst] = mask[y_src, x_src] shifted_count += sub.astype(np.int32) sub_mask_count = shifted_count[mask] return float(sub_mask_count.mean()) / 8.0 cluster_scores = {} for a in AREAS: low_mask = ((a["arr"] < 5) & (a["arr"] >= 0)) & a["valid_mask"] score = neighbor_homogeneity(low_mask) cluster_scores[a["key"]] = score print(f" {a['label']}: low-density neighbor homogeneity = {score:.3f}", flush=True) avg_cluster = (cluster_scores["sera_1m"] + cluster_scores["kumano_50cm"]) / 2 h6 = avg_cluster >= 0.4 hypos.append({ "H": "H6", "claim": "低密度セルは集塊する (8隣接同種率 >= 0.4)", "result": f"世羅 {cluster_scores['sera_1m']:.3f}, 熊野 {cluster_scores['kumano_50cm']:.3f}, avg {avg_cluster:.3f}", "verdict": "支持" if h6 else ("部分支持" if avg_cluster >= 0.25 else "反証"), }) hypos_df = pd.DataFrame(hypos) hypos_df.to_csv(ASSETS / "L37_hypothesis_results.csv", index=False, encoding="utf-8-sig") with (ASSETS / "L37_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) # ============================================================================= # 9. dataset/series 一覧 CSV # ============================================================================= print("\n[9] 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"], "unit_area_m2": a["unit_area_m2"], "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 / "L37_series.csv", index=False, encoding="utf-8-sig") print(series_df.to_string(index=False)) print(f" ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 10. 図 1: 2 dataset カード + 点群密度ヒストグラム # ============================================================================= print("\n[10] 図 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=14, 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=10.5, 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'].split('(')[0].strip()})", ha="center", fontsize=8.5, color="#666") ax.text(2.5, 5.6, f"解像度 {AREAS[0]['resolution_m']:.1f} m / 単位面積 {AREAS[0]['unit_area_m2']} m²", ha="center", fontsize=9) 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} 点/{AREAS[0]['unit_area_m2']}m²\n中央値 {s_sera['median']:.2f}", ha="center", fontsize=10, fontweight="bold") ax.text(2.5, 1.2, f"低密度<5 比率 {(s_sera['n_zero']+s_sera['n_low'])/max(1,s_sera['n_cells_valid'])*100:.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=10.5, 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'].split('(')[0].strip()})", ha="center", fontsize=8.5, color="#666") ax.text(7.5, 5.6, f"解像度 {AREAS[1]['resolution_m']:.2f} m / 単位面積 {AREAS[1]['unit_area_m2']} m²", ha="center", fontsize=9) 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} 点/{AREAS[1]['unit_area_m2']}m²\n中央値 {s_kuma['median']:.2f}", ha="center", fontsize=10, fontweight="bold") ax.text(7.5, 1.2, f"1m²換算 {kuma_mean_1m2:.2f} 点/m²", ha="center", fontsize=9, color=AREAS[1]["color"]) # (右) ヒスト 1m² 換算で重ね ax = axes[1] bin_centers_1m2_sera = (bins_linear[:-1] + 2.5) / AREAS[0]["unit_area_m2"] # 1m² 換算 X 軸 bin_centers_1m2_kuma = (bins_linear[:-1] + 2.5) / AREAS[1]["unit_area_m2"] # 棒幅は単位面積に応じる w_sera = 4.0 w_kuma = 1.0 sera_pct = hist_linear["sera_1m"] / max(1, hist_linear["sera_1m"].sum()) * 100 kuma_pct = hist_linear["kumano_50cm"] / max(1, hist_linear["kumano_50cm"].sum()) * 100 # 値スケールの違いを示すため、X軸を生値で別々に表示 xs = np.arange(len(bins_linear) - 1) ax.bar(xs - 0.2, sera_pct, width=0.4, color=AREAS[0]["color"], alpha=0.85, label=f"世羅 1m, 1m²単位 (n={int(hist_linear['sera_1m'].sum()):,})") ax.bar(xs + 0.2, kuma_pct, width=0.4, color=AREAS[1]["color"], alpha=0.85, label=f"熊野 50cm, 0.25m²単位 (n={int(hist_linear['kumano_50cm'].sum()):,})") ax.set_xticks(xs[::2]) ax.set_xticklabels([f"{int(bins_linear[i])}-{int(bins_linear[i+1])}" for i in range(0, len(bins_linear)-1, 2)], rotation=45, fontsize=8) ax.set_xlabel("点群密度 (点/単位面積, 5 刻みビン)") ax.set_ylabel("セル比率 (%)") ax.set_title("点群密度ヒストグラム — 2 エリア (生値, 単位面積は dataset 仕様による)") ax.legend(loc="upper right", fontsize=9) ax.grid(axis="y", alpha=0.3) plt.tight_layout() fig.savefig(ASSETS / "L37_fig1_dataset_overview.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L37_fig1_dataset_overview.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 11. 図 2: 2 エリアの点群密度ラスタ主題図 # ============================================================================= print("\n[11] 図 2: 密度ラスタマップ", flush=True) t1 = time.time() # 1m² 換算で同じカラースケール cmap_levels = [0, 1, 5, 10, 20, 50, 100, 200] cmap_colors = ["#67000d", "#cb181d", "#ef6548", "#fdae6b", "#fee08b", "#a1d99b", "#005a32"] 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() # 1m² 換算 arr_1m2 = arr / a["unit_area_m2"] nodata_v = a["meta"]["nodata"] if a["meta"]["nodata"] is not None else -99999.0 # 表示用: nodata は np.nan arr_show = np.where((arr == nodata_v) | (~np.isfinite(arr)) | (arr < -100), np.nan, arr_1m2) arr_show = np.clip(arr_show, 0, 200) 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") 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']} 点群密度 (1m² 換算, overview /1/{a['arr_factor']}, 実効 {a['arr_cell_m']:.0f} m/cell)\n" f"{a['character']}", fontsize=10) ax.set_xlabel("X (EPSG:6671 m)") ax.set_ylabel("Y (EPSG:6671 m)") 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², 1m²換算) — 0/低=暗赤(=樹冠閉鎖), 中=橙黄, 高=緑(=地表開放)") cbar.set_ticks(cmap_levels) fig.savefig(ASSETS / "L37_fig2_raster_maps.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L37_fig2_raster_maps.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 12. 図 3: 階級別構成比 (積み上げバー) # ============================================================================= print("\n[12] 図 3: 階級別構成比", flush=True) t1 = time.time() fig, ax = plt.subplots(figsize=(12, 5.5)) 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 = ["#67000d", "#cb181d", "#ef6548", "#fdae6b", "#fee08b", "#a1d99b", "#005a32"] 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="white" if i in (0, 1, 6) else "black") bottom = bottom + np.array(vals) ax.set_xticks(x) ax.set_xticklabels([f"{a['label']}\n({a['resolution_m']}m, {a['unit_area_m2']}m²単位)" 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 / "L37_fig3_class_composition.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L37_fig3_class_composition.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 13. 図 4: small multiples — 階級別マスクマップ # ============================================================================= print("\n[13] 図 4: 階級別マスク small multiples", flush=True) t1 = time.time() fig, axes = plt.subplots(2, 4, figsize=(15, 8)) mask_def = [ ("超低密度 (<5)", lambda a, u: ((a >= 0) & (a < 5)), "#cb181d"), ("中低密度 (5-20)", lambda a, u: ((a >= 5) & (a < 20)), "#fdae6b"), ("中高密度 (20-50)",lambda a, u: ((a >= 20) & (a < 50)),"#fee08b"), ("高密度 (>=50)", lambda a, u: (a >= 50), "#005a32"), ] 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, area["unit_area_m2"]) & valid_mask n_cells = int(m.sum()) n_total_valid = int(valid_mask.sum()) pct = n_cells / max(1, n_total_valid) * 100 show = np.where(m, 1.0, np.nan) cmap = ListedColormap([col]) ax.imshow(show, cmap=cmap, extent=extent, origin="upper", vmin=0, vmax=1, interpolation="nearest") if area.get("admin") is not None: area["admin"].boundary.plot(ax=ax, color="#222", linewidth=0.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.85, edgecolor="none", pad=2)) ax.set_xticks([]) ax.set_yticks([]) fig.suptitle("点群密度階級別 空間分布 — 2 エリア × 4 階級 (small multiples) — 生値ベース", fontsize=12, y=0.99) plt.tight_layout() fig.savefig(ASSETS / "L37_fig4_small_multiples.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L37_fig4_small_multiples.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 14. 図 5: ヒスト詳細 + boxplot + パーセンタイル # ============================================================================= print("\n[14] 図 5: 分布詳細", flush=True) t1 = time.time() fig, axes = plt.subplots(2, 2, figsize=(13, 9)) # (上左) 世羅 詳細ヒスト (1m² 換算) ax = axes[0, 0] v = AREAS[0]["valid"] / AREAS[0]["unit_area_m2"] # 1m² 換算 v = v[(v >= 0) & (v <= 100)] ax.hist(v, bins=np.arange(0, 101, 1), color=AREAS[0]["color"], alpha=0.85, edgecolor="white", linewidth=0.3) ax.axvline(s_sera["mean"]/AREAS[0]["unit_area_m2"], color="red", linestyle="--", linewidth=1.5, label=f"平均 {s_sera['mean']/AREAS[0]['unit_area_m2']:.2f}/m²") ax.axvline(s_sera["median"]/AREAS[0]["unit_area_m2"], color="orange", linestyle="--", linewidth=1.5, label=f"中央値 {s_sera['median']/AREAS[0]['unit_area_m2']:.2f}/m²") ax.set_xlabel("点群密度 (点/m², 1m² 換算)") ax.set_ylabel("セル数") ax.set_title("世羅町 (1m, 区分A) 密度分布詳細") ax.legend() ax.grid(axis="y", alpha=0.3) # (上右) 熊野 詳細ヒスト (1m² 換算) ax = axes[0, 1] v = AREAS[1]["valid"] / AREAS[1]["unit_area_m2"] v = v[(v >= 0) & (v <= 100)] ax.hist(v, bins=np.arange(0, 101, 1), color=AREAS[1]["color"], alpha=0.85, edgecolor="white", linewidth=0.3) ax.axvline(s_kuma["mean"]/AREAS[1]["unit_area_m2"], color="red", linestyle="--", linewidth=1.5, label=f"平均 {s_kuma['mean']/AREAS[1]['unit_area_m2']:.2f}/m²") ax.axvline(s_kuma["median"]/AREAS[1]["unit_area_m2"], color="orange", linestyle="--", linewidth=1.5, label=f"中央値 {s_kuma['median']/AREAS[1]['unit_area_m2']:.2f}/m²") ax.set_xlabel("点群密度 (点/m², 1m² 換算)") ax.set_ylabel("セル数") ax.set_title("熊野町 (50cm, 区分B) 密度分布詳細") ax.legend() ax.grid(axis="y", alpha=0.3) # (下左) boxplot 1m² 換算 ax = axes[1, 0] data_box = [] for a in AREAS: v = a["valid"] / a["unit_area_m2"] v = v[(v >= 0) & (v <= 100)] data_box.append(v) 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², 1m² 換算)") ax.set_title("両エリア boxplot (1m² 換算で比較)") ax.grid(axis="y", alpha=0.3) # (下右) パーセンタイル ax = axes[1, 1] percs = [10, 25, 50, 75, 90, 99] xs2 = np.arange(len(percs)) sera_p = [(s_sera[f"p{p}"] if p != 50 else s_sera["median"]) / AREAS[0]["unit_area_m2"] for p in percs] kuma_p = [(s_kuma[f"p{p}"] if p != 50 else s_kuma["median"]) / AREAS[1]["unit_area_m2"] for p in percs] w = 0.35 ax.bar(xs2 - w/2, sera_p, w, label=f"世羅町 (1m)", color=AREAS[0]["color"], alpha=0.85, edgecolor="#333") ax.bar(xs2 + w/2, kuma_p, w, label=f"熊野町 (50cm)", color=AREAS[1]["color"], alpha=0.85, edgecolor="#333") for i, (sv, kv) in enumerate(zip(sera_p, kuma_p)): ax.text(i - w/2, sv + 0.5, f"{sv:.1f}", ha="center", fontsize=8.5) ax.text(i + w/2, kv + 0.5, f"{kv:.1f}", ha="center", fontsize=8.5) ax.set_xticks(xs2) ax.set_xticklabels([f"P{p}" for p in percs]) ax.set_ylabel("点群密度 (点/m², 1m² 換算)") ax.set_title("パーセンタイル比較 (P10〜P99, 1m² 換算)") ax.legend() ax.grid(axis="y", alpha=0.3) plt.tight_layout() fig.savefig(ASSETS / "L37_fig5_distribution_detail.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L37_fig5_distribution_detail.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 15. 図 6: L36 樹高 vs L37 密度 散布 (Hexbin) — 計測物理的相補性 # ============================================================================= print("\n[15] 図 6: L36 樹高 × L37 密度 hexbin", flush=True) t1 = time.time() fig, axes = plt.subplots(1, 2, figsize=(14, 6)) for ax, a in zip(axes, AREAS): bvm = a["both_valid_mask"] if bvm is None or bvm.sum() < 100: ax.text(0.5, 0.5, "L36 ラスタ不在\n(対応する樹高図が無い)", transform=ax.transAxes, ha="center", va="center", fontsize=12) ax.set_title(f"{a['label']}: L36 不在") continue d = a["arr"][bvm].astype(np.float32) / a["unit_area_m2"] # 1m² 換算 h = np.clip(a["arr_l36"][bvm].astype(np.float32), -2, 30) # 範囲をクリップ mask = (d >= 0) & (d <= 200) & (h >= 0) & (h <= 30) d2 = d[mask] h2 = h[mask] if len(d2) > 0: hb = ax.hexbin(h2, d2, gridsize=40, cmap="magma_r", mincnt=1, bins="log") cb = plt.colorbar(hb, ax=ax) cb.set_label("セル数 (log)") # Pearson r r = correlations[a["key"]]["corr"] ax.set_title(f"{a['label']}: 樹高 (L36) × 点群密度 (L37) hexbin\n" f"Pearson r = {r:+.3f} (n={correlations[a['key']]['n']:,})", fontsize=10) ax.set_xlabel("樹高 (m, L36 DCHM)") ax.set_ylabel("地面到達点群密度 (点/m²)") ax.grid(alpha=0.3) # 仮説の傾向線 (負相関期待) if h2.std() > 0 and d2.std() > 0: slope, intercept = np.polyfit(h2, d2, 1) x_line = np.linspace(0, 25, 50) ax.plot(x_line, slope * x_line + intercept, "g-", linewidth=2, label=f"線形回帰 y={slope:+.2f}x+{intercept:.1f}") ax.legend(loc="upper right", fontsize=8) plt.tight_layout() fig.savefig(ASSETS / "L37_fig6_height_vs_density.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L37_fig6_height_vs_density.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 16. 図 7: 低密度ホットスポット ↔ 高樹高 重ね合わせマップ # ============================================================================= print("\n[16] 図 7: 低密度 ↔ 高樹高 重ね合わせ", 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 / a["unit_area_m2"], 0, 50), np.nan) ax.imshow(base, cmap="Greys", extent=extent, origin="upper", alpha=0.4, vmin=0, vmax=30) # 低密度マスク (< 5 点/単位) を赤で low_density = (arr >= 0) & (arr < 5) & vm low_show = np.where(low_density, 1.0, np.nan) ax.imshow(low_show, cmap=ListedColormap(["#cb181d"]), extent=extent, origin="upper", vmin=0, vmax=1, alpha=0.55) # L36 高樹高マスク (>=15m) を緑で重ねる if a["arr_l36"] is not None: high_tree = (a["arr_l36"] >= 15) & a["valid_mask_l36"] ht_show = np.where(high_tree, 1.0, np.nan) ax.imshow(ht_show, cmap=ListedColormap(["#005a32"]), extent=extent, origin="upper", vmin=0, vmax=1, alpha=0.45) # 共通領域 (低密度 AND 高樹高) を黄色で強調 intersect = low_density & high_tree if intersect.sum() > 0: ix_show = np.where(intersect, 1.0, np.nan) ax.imshow(ix_show, cmap=ListedColormap(["#fee08b"]), extent=extent, origin="upper", vmin=0, vmax=1, alpha=0.85) n_low = int(low_density.sum()) n_ht = int(high_tree.sum()) n_both = int(intersect.sum()) # H4 の指標: 低密度セルのうち何 % が高樹高か pct_in_ht = n_both / max(1, n_low) * 100 info_text = f"低密度: {n_low:,} / 高樹高: {n_ht:,} / 共通: {n_both:,}\n低密度のうち高樹高: {pct_in_ht:.1f}%" else: info_text = "L36 不在" if a.get("admin") is not None: a["admin"].boundary.plot(ax=ax, color="#222", linewidth=1.0, alpha=0.85) # 凡例 legend_elements = [ Patch(facecolor="#cb181d", alpha=0.55, label="低密度 (<5 点/単位)"), Patch(facecolor="#005a32", alpha=0.45, label="高樹高 (>=15m, L36)"), Patch(facecolor="#fee08b", alpha=0.85, label="共通 (低密度 AND 高樹高)"), ] ax.legend(handles=legend_elements, loc="upper right", fontsize=8.5) ax.set_title(f"{a['label']}: 低密度 ↔ 高樹高 重ね合わせマップ\n{info_text}", fontsize=10) ax.set_xlabel("X (m)") ax.set_ylabel("Y (m)") plt.tight_layout() fig.savefig(ASSETS / "L37_fig7_lowdensity_hightree_overlay.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L37_fig7_lowdensity_hightree_overlay.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 17. 図 8: 中央サンプル (L37 密度 vs L36 樹高 / 同範囲) — 体感図 # ============================================================================= print("\n[17] 図 8: 中央サンプル", flush=True) t1 = time.time() fig, axes = plt.subplots(2, 2, figsize=(13, 9)) for col_i, a in enumerate(AREAS): # 上段: L37 点群密度 ax = axes[0, col_i] raw = a["raw_sample_arr"] nodata_v = a["meta"]["nodata"] if a["meta"]["nodata"] is not None else -99999.0 arr_show = np.where((raw == nodata_v) | (~np.isfinite(raw)) | (raw < -100), np.nan, raw) arr_show = arr_show / a["unit_area_m2"] # 1m² 換算 arr_show = np.clip(arr_show, 0, 100) im = ax.imshow(arr_show, cmap=custom_cmap, norm=custom_norm, interpolation="nearest") side_m = raw.shape[1] * a["resolution_m"] ax.set_title(f"{a['label']} 点群密度 中央サンプル {raw.shape[1]}×{raw.shape[0]} = {side_m:.0f}m × {side_m:.0f}m\n" f"({a['resolution_m']}m/cell, 平均 {a['raw_sample_mean']:.2f}/{a['unit_area_m2']}m²)", fontsize=10) ax.set_xlabel("ピクセル X") ax.set_ylabel("ピクセル Y") cb = plt.colorbar(im, ax=ax, fraction=0.045) cb.set_label("点/m² (1m²換算)", fontsize=8) # 下段: L36 樹高 (同範囲) ax = axes[1, col_i] if a.get("raw_sample_l36_arr") is not None: l36_arr = a["raw_sample_l36_arr"] nodata_th = -99999.0 l36_show = np.where((l36_arr == nodata_th) | (~np.isfinite(l36_arr)) | (l36_arr < -100), np.nan, l36_arr) l36_show = np.clip(l36_show, 0, 30) # 樹高用 colormap (緑系) cmap_h = ListedColormap(["#e0e0e0", "#fff7bc", "#fec44f", "#a1d99b", "#41ab5d", "#005a32", "#54278f"]) norm_h = BoundaryNorm([0, 1, 3, 6, 10, 15, 20, 30], 7, clip=True) im2 = ax.imshow(l36_show, cmap=cmap_h, norm=norm_h, interpolation="nearest") ax.set_title(f"{a['label']} 樹高 (L36) 同範囲\n平均 {a['raw_sample_l36_mean']:.2f} m", fontsize=10) ax.set_xlabel("ピクセル X") ax.set_ylabel("ピクセル Y") cb2 = plt.colorbar(im2, ax=ax, fraction=0.045) cb2.set_label("樹高 (m)", fontsize=8) else: ax.text(0.5, 0.5, "L36 不在", transform=ax.transAxes, ha="center", va="center") ax.set_axis_off() fig.suptitle("中央サンプル: L37 点群密度 (上) と L36 樹高 (下) の同範囲比較", fontsize=12, y=1.0) plt.tight_layout() fig.savefig(ASSETS / "L37_fig8_center_samples.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L37_fig8_center_samples.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 18. 図 9: 概念図 — グラウンドパルス vs 第一反射 (LiDAR) の分離 # ============================================================================= print("\n[18] 図 9: LiDAR グラウンド/第一反射の概念", flush=True) t1 = time.time() fig, ax = plt.subplots(figsize=(11, 6.5)) ax.set_xlim(0, 12); ax.set_ylim(0, 9); ax.set_axis_off() ax.text(6, 8.5, "航空 LiDAR 計測: 第一反射 (DSM) と グラウンドパルス (DTM/L37) の分離", ha="center", fontsize=13, fontweight="bold") # 飛行機 ax.add_patch(plt.Rectangle((0.5, 7.0), 1.5, 0.4, facecolor="#666", alpha=0.85)) ax.text(1.25, 7.7, "航空機", ha="center", fontsize=8.5) # レーザパルス (複数本下向き) import numpy as np_local for x_p in [2.5, 3.3, 4.1, 4.9, 5.7, 6.5, 7.3, 8.1, 8.9, 9.7, 10.5]: ax.annotate("", xy=(x_p + 0.4, 1.5), xytext=(1.6, 7.0), arrowprops=dict(arrowstyle="->", color="#cf222e", linewidth=0.5, alpha=0.5)) # 地面 ax.add_patch(plt.Polygon([[0.3, 1.5], [11.7, 1.5], [11.7, 1.0], [0.3, 1.0]], facecolor="#8B4513", alpha=0.7)) # 樹冠 (深緑) -- 厚いキャノピー ax.add_patch(mpatches.Wedge((3.3, 1.5), 1.5, 0, 180, facecolor="#005a32", alpha=0.85)) ax.text(3.3, 4.3, "高木林\n(樹冠閉鎖)", ha="center", fontsize=9, color="darkgreen", fontweight="bold") # 樹冠 (中木) -- 開きあり ax.add_patch(mpatches.Wedge((6.0, 1.5), 1.0, 0, 180, facecolor="#41ab5d", alpha=0.6)) ax.text(6.0, 3.3, "中木\n(隙間あり)", ha="center", fontsize=9, color="darkgreen", fontweight="bold") # 開地・農地 ax.add_patch(plt.Polygon([[8.5, 1.5], [10.5, 1.5], [10.4, 1.7], [8.6, 1.7]], facecolor="#fff7bc", alpha=0.85)) ax.text(9.5, 2.2, "農地・草地\n(開放)", ha="center", fontsize=9) # 第一反射点 (緑の点 = 樹冠頂点で反射) ax.plot([3.0, 3.6, 5.7, 6.3], [3.0, 3.0, 2.5, 2.5], "go", markersize=10, label="第一反射 (DSM = 樹冠頂点)") # グラウンドパルス点 (赤の点 = 地面に届いた) ax.plot([2.0, 2.2, 4.5, 5.0, 6.5, 6.8, 7.2, 8.5, 9.0, 9.8, 10.3], [1.5]*11, "rs", markersize=8, label="グラウンドパルス (DTM, L37)") # 高木下の地面に達するパルスは少ないことを矢印で ax.annotate("少", xy=(3.3, 1.3), xytext=(3.3, 0.3), fontsize=14, ha="center", color="darkred", fontweight="bold", arrowprops=dict(arrowstyle="->", color="darkred", linewidth=2)) ax.annotate("多", xy=(9.5, 1.3), xytext=(9.5, 0.3), fontsize=14, ha="center", color="darkgreen", fontweight="bold", arrowprops=dict(arrowstyle="->", color="darkgreen", linewidth=2)) ax.text(6, 5.5, "L37 = 各セルに達した「赤四角の点数」のラスタ\n樹冠閉鎖下では赤点が少ない (=低密度) → DTM 補間誤差大", ha="center", fontsize=10, color="black", bbox=dict(facecolor="#fff7bc", alpha=0.9, edgecolor="black", boxstyle="round,pad=0.5")) ax.legend(loc="upper right", fontsize=9) fig.savefig(ASSETS / "L37_fig9_lidar_concept.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L37_fig9_lidar_concept.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 19. 図 10: overview vs 生解像度サンプル の階級比較 # ============================================================================= print("\n[19] 図 10: overview vs 生解像度", flush=True) t1 = time.time() fig, axes = plt.subplots(1, 2, figsize=(13, 5.5)) class_labels_short = ["0\n超低", "0-5\n低", "5-20\n中", ">=20\n高"] for area_idx, a in enumerate(AREAS): ax = axes[area_idx] valid_ovr = a["valid"] n_ovr = len(valid_ovr) pct_zero_ovr = (valid_ovr < 0.5).sum() / max(1, n_ovr) * 100 pct_low_ovr = ((valid_ovr >= 0.5) & (valid_ovr < 5)).sum() / max(1, n_ovr) * 100 pct_mid_ovr = ((valid_ovr >= 5) & (valid_ovr < 20)).sum() / max(1, n_ovr) * 100 pct_high_ovr = (valid_ovr >= 20).sum() / max(1, 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_b = np.arange(4) w = 0.35 bars1 = ax.bar(x_b - 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_b + 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, val in zip(bars, vals): ax.text(b.get_x() + b.get_width()/2, val + 0.5, f"{val:.1f}", ha="center", fontsize=8.5, fontweight="bold") ax.set_xticks(x_b) 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 + [1]), max(raw_pct + [1])) * 1.25) plt.tight_layout() fig.savefig(ASSETS / "L37_fig10_overview_vs_raw.png", dpi=130, bbox_inches="tight") plt.close("all") print(f" saved L37_fig10_overview_vs_raw.png ({time.time()-t1:.1f}s)", flush=True) # ============================================================================= # 20. 中間 CSV # ============================================================================= print("\n[20] 中間 CSV", flush=True) t1 = time.time() # ヒストグラム hist_df = pd.DataFrame({ "bin_low": bins_linear[:-1], "bin_high": bins_linear[1:], "sera_1m": hist_linear["sera_1m"], "kumano_50cm": hist_linear["kumano_50cm"], }) hist_df["sera_pct"] = hist_df["sera_1m"] / max(1, hist_df["sera_1m"].sum()) * 100 hist_df["kumano_pct"] = hist_df["kumano_50cm"] / max(1, hist_df["kumano_50cm"].sum()) * 100 hist_df.to_csv(ASSETS / "L37_density_histogram.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"] / max(1, class_w["sera_n"].sum()) * 100 class_w["kumano_pct"] = class_w["kumano_n"] / max(1, class_w["kumano_n"].sum()) * 100 class_w.to_csv(ASSETS / "L37_density_classes_wide.csv", index=False, encoding="utf-8-sig") # 樹高×密度の相関 corr_df = pd.DataFrame([ {"area": "世羅町 (1m)", "n_paired_cells": correlations["sera_1m"]["n"], "pearson_r": round(correlations["sera_1m"]["corr"], 4)}, {"area": "熊野町 (50cm)", "n_paired_cells": correlations["kumano_50cm"]["n"], "pearson_r": round(correlations["kumano_50cm"]["corr"], 4)}, ]) corr_df.to_csv(ASSETS / "L37_l36_correlation.csv", index=False, encoding="utf-8-sig") # 集塊性 cluster_df = pd.DataFrame([ {"area": a["label"], "low_density_neighbor_homogeneity": round(cluster_scores[a["key"]], 4), "n_low_cells": int(((a["arr"] < 5) & (a["arr"] >= 0) & a["valid_mask"]).sum())} for a in AREAS ]) cluster_df.to_csv(ASSETS / "L37_clustering_scores.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t1:.1f}s)", flush=True) print(f"\n=== Total elapsed: {time.time()-t0:.1f}s ===", flush=True) # ============================================================================= # 21. HTML 生成 # ============================================================================= print("\n[21] HTML 生成", flush=True) t1 = time.time() def dl(name): return f'{name}' sections = [] # === Section 1: 学習目標と問い === sec1 = f"""
カバー宣言: 本記事は DoBoX のシリーズ 「地面到達レーザ点群密度図 1m / 50cm」 2 件 (dataset_id = 1632, 1633) を統合し、広島県 2 自治体エリア (世羅町 1m, 熊野町 50cm) の地面到達レーザ点群密度から DTM 精度の地理構造と林冠閉鎖度の代理指標を分析する研究記事です。 両データセットは航空レーザ計測 (LiDAR) のグラウンドパルス点数ラスタで、 1 ピクセルあたり 1m² (区分A) もしくは 0.25m² (区分B) のメッシュで 「地面に届いたレーザ点の数」が記録されています。

本記事の独自性 — L36 樹高図との計測物理学的相補

L37 (地面到達点群密度) は、L36 (樹高 = DSM - DTM) と同じ航空レーザ点群から 派生した別の派生物です。両者は次のように対をなします:

物理的に、樹冠が高くて密に閉じている場所では地面に届くレーザが少ない (低密度) はず。逆に農地や道路の開放地では地面に多くのレーザが届く (高密度) はず。 本記事はこの「樹高 × 点群密度」の負相関仮説を、L36 と同じ 2 自治体ペア (世羅町 1m / 熊野町 50cm) で実データの空間結合により検証する。 これは L36 や L37 単独では不可能な、両者の同時空間整合分析である。

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

このシリーズの 2 件は、「区分A=1m² / 区分B=0.25m²」の解像度 2 系統です。 両者は同じ航空レーザ計測のオリジナル点群から、集計単位を変えて 県内 22-23 自治体それぞれを GeoTIFF に切出した形で公開されています。

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

用語定義
地面到達レーザ点 (グラウンドパルス)航空レーザ計測で、樹冠・建物などを通り抜けて 地面に直接届いたパルス。点群分類アルゴリズムにより「地面点 (ground class)」として 分類された点。樹冠頂点で反射した「第一反射点」とは別に集計される。
点群密度 (本記事)各セル (1m² または 0.25m²) に届いた地面到達パルスの本数。 本記事の TIFF ピクセル値そのもの。値は 0〜数十が典型 (今回扱った 2 エリアは 0〜10 の範囲に集中)。
区分A / 区分BDoBoX 仕様書の用語。区分A = 1m² 単位集計 (1m 解像度の TIFF)、 区分B = 0.25m² 単位集計 (50cm 解像度の TIFF)。同じ 1m² の地面に届いたパルス数を比較するには 区分B の値を ×4 する必要がある (= 1m²換算)。
DTM 精度の代理指標地面到達パルスが多いセル = 地面の標高が直接観測できる = DTM の補間誤差小。地面到達パルスがゼロのセル = 周囲のセルから補間する必要がある = DTM の補間誤差大。本記事ではこれを「DTM 精度の代理指標」と呼ぶ (DoBoX 公式仕様にも「標高図の精度指標として使用」とある)。
林冠閉鎖度 (canopy closure)森林学の用語。樹冠が天蓋を閉じ、 地面に光・雨・レーザが届かない比率。本記事では「点群密度の低さ」を林冠閉鎖度の 代理指標として扱う (公式定義ではないが、物理的には妥当)。
低密度ホットスポット本記事独自定義。1m²換算で 5 点/m² 未満のセル群が 空間的に集塊している領域。森林ブロック単位で連続している傾向。
計測物理学的相補本記事独自定義。L36 樹高ラスタと L37 点群密度ラスタは、 同じレーザ点群から派生した別の集計であり、樹冠が高い場所では密度が低くなるという 物理的に必然の関係を持つ。両者を空間結合すると、計測の物理プロセス自体が地理データに 刻まれていることが確認できる。
overview / ピラミッドGeoTIFF 内の低解像度版。L36 と同じ仕組み。 本記事も /1/64〜/127 を使い、メモリ 0.1 MB 程度で全体概観する。

研究の問い (RQ)

広島県の地面到達レーザ点群密度図 2 解像度 (1m / 50cm) は、それぞれ 世羅町 (中山間, 林業地)と熊野町 (都市近郊, 山岳混在)という異なる地表被覆の エリアを公開している。両エリアの点群密度ヒストグラム・空間分布・低密度ホットスポット分析を比較し、 さらにL36 樹高ラスタとの空間相関を取ると、計測精度の地理構造と 林冠閉鎖度の地理がどう見えるか?

  1. 2 dataset の面積・解像度・有効データ率・単位面積の構造比較は?
  2. 点群密度ヒストグラムのピーク位置・裾の重さはどう違うか?
  3. 超低密度 (<5)・中密度 (5-20)・高密度 (>=20) の構成比は?
  4. 低密度ホットスポットの空間分布と集塊性はどうか?
  5. L36 樹高 vs L37 密度 の Pearson 相関は仮説 (負相関) を支持するか?
  6. 区分A (1m²) vs 区分B (0.25m²) の値スケールは仕様通り 1:0.25 か?

仮説 H1〜H6

到達点

2 解像度の点群密度ラスタを「2 自治体の代表エリア」として読み解き、 合計 {int(metas[0]['n_cells']/1e6 + metas[1]['n_cells']/1e6):,} M セル (約 7 億ピクセル) のグラウンドパルス密度データを overview ピラミッドで効率的に概観する。さらにL36 樹高ラスタとの空間結合で 「樹高が高い場所は本当に地面到達が少ないか」を Pearson 相関で実証する。 解像度差 (区分A vs 区分B) が値スケールにどう影響するかも仕様通りに検証する。

本記事のスコープ外: 本記事は L36 と同じ 2 自治体で計測物理学的相補比較に特化する。
""" sections.append(("学習目標と問い", sec1)) # === Section 2: 使用データ === sec2 = f"""

本記事が使用する 2 dataset の一覧。両者は同一手法 (LiDAR グラウンドパルス集計) の単位面積違いで、 県内 22-23 自治体それぞれに GeoTIFF が公開されています。本記事は各 1 自治体を代表として用います。

dsid名称解像度単位面積公開単位本記事の代表TIFF サイズセル数DoBoX
1632地面到達レーザ点群密度図（1m）1.0 m/cell1.0 m² 22 自治体別 GeoTIFF世羅町 (34462) {int(AREAS[0]['tif'].stat().st_size/(1024*1024)):,} MB {AREAS[0]['meta']['n_cells']/1e6:.1f} M #1632
1633地面到達レーザ点群密度図（50cm）0.5 m/cell0.25 m² 23 自治体別 GeoTIFF熊野町 (34307) {int(AREAS[1]['tif'].stat().st_size/(1024*1024)):,} MB {AREAS[1]['meta']['n_cells']/1e6:.1f} M #1633

データ仕様 (両 dataset 共通)

L36 樹高図との計測物理学的関係

軸L36 樹高図 (dsid 1626/1627)L37 点群密度図 (本記事, dsid 1632/1633)
計測対象樹冠の高さ (m)地面到達パルス数 (点)
派生過程DSM (第一反射) − DTM (地面)グラウンド分類された点数を集計
同じ点群?YES — 同一の航空レーザ計測点群 から派生した 2 つの集計
物理的関係樹冠が高い ⇔地面到達が少ない (= 低密度)
本記事で検証 同じ 2 自治体ペア (世羅町 1m / 熊野町 50cm) で空間結合し、 Pearson 相関で「樹高 vs 密度」の負相関を実証
dtypefloat32 (cm 単位の連続値)uint8 / int16 (整数の点数)
nodata-99999.0255 / -9999

2 自治体の特性比較

軸世羅町 (1m)熊野町 (50cm)
立地備後内陸 (中山間)広島市東 (都市近郊)
面積 (行政界) {round(AREAS[0]['muni_area_km2'], 1) if AREAS[0]['muni_area_km2'] else '?'} km² {round(AREAS[1]['muni_area_km2'], 1) if AREAS[1]['muni_area_km2'] else '?'} km²
TIFF 範囲 (bbox) {AREAS[0]['meta']['width_m']/1000:.1f} × {AREAS[0]['meta']['height_m']/1000:.1f} km {AREAS[1]['meta']['width_m']/1000:.1f} × {AREAS[1]['meta']['height_m']/1000:.1f} km
主産業農林業 (世羅梨・りんご・大根)筆 (熊野筆) 製造
有効データ率 {basic_stats['sera_1m']['valid_ratio']*100:.1f}% {basic_stats['kumano_50cm']['valid_ratio']*100:.1f}%
平均密度 (生値) {s_sera['mean']:.2f} 点/{AREAS[0]['unit_area_m2']}m² {s_kuma['mean']:.2f} 点/{AREAS[1]['unit_area_m2']}m²
1m² 換算平均 {sera_mean_1m2:.2f} 点/m² {kuma_mean_1m2:.2f} 点/m² (overview) / {kuma_raw_mean_1m2:.2f} 点/m² (生解像度サンプル)
L36 平均樹高 (本記事サンプル) {AREAS[0].get('raw_sample_l36_mean', 0):.2f} m {AREAS[1].get('raw_sample_l36_mean', 0):.2f} m

データ品質メモ

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

生データ (DoBoX 直リンク)

容量警告: 上記 2 ZIP の合計は約 {int((AREAS[0]['zip_local'].stat().st_size + AREAS[1]['zip_local'].stat().st_size)/(1024*1024))} MB、 展開後の TIFF は約 {int((AREAS[0]['tif'].stat().st_size + AREAS[1]['tif'].stat().st_size)/(1024*1024))} MB。 他の自治体版や 50cm 広島市版はさらに大きい (DoBoX 仕様確認のこと) ため、本記事のスコープ外。

L36 樹高図 (本記事の相補比較に必要)

本記事は L36 樹高ラスタを「同じ 2 自治体 + 同じ範囲」で読み込んで比較します。 L36 の TIFF を取得していない場合はL36 記事を先に実行してください (または L36_canopy_height.py を単独実行)。本スクリプトは L36 が無くても L37 単独分析は走りますが、 図 6/図 7/図 8 が一部省略されます。

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

図 (本記事生成 PNG)

再現用 Python スクリプト

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

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

狙い

2 dataset (1m 区分A / 50cm 区分B) と本記事が選んだ 2 エリア (世羅町 / 熊野町) の 規模・形状・データ量・単位面積を 1 枚の絵で示す。 点群密度ヒストグラムを重ねて地表被覆の違いを最初に視覚化する。

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

L36 と全く同じ rasterio スタイル。違うのは:

実装

{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: 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) # nodata マスク (uint8 は 255, int16 は -9999) nodata = ds.nodata valid_mask = (~np.isnan(arr)) & (arr != nodata) & (arr > -1000) area["valid"] = arr[valid_mask] # 1m² 換算するには unit_area_m2 で割る # area['unit_area_m2'] は 1.0 (区分A) or 0.25 (区分B) ''')}

図と読み取り (図 1: dataset カード + ヒストグラム重ね)

なぜこの図か: 2 dataset の規模と特性をカードで言語的に、 かつヒストグラムで量的特性として、同時に 1 枚で伝えるため。 左カード = 言語的概要、右ヒスト = 量的概観。

{figure("assets/L37_fig1_dataset_overview.png", "図 1: 2 dataset 概要カード (左) と 点群密度ヒストグラム (右)")}

読み取り:

表と読み取り (基本統計)

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

読み取り: 平均密度は世羅 {s_sera['mean']:.2f}/1m² vs 熊野 {s_kuma['mean']:.2f}/0.25m²。 1m² 換算すると世羅 {sera_mean_1m2:.2f}, 熊野 {kuma_mean_1m2:.2f}でほぼ同水準 (overview平均化の影響。生解像度サンプルでは熊野 {kuma_raw_mean_1m2:.2f} と高い)。 中央値は世羅 {s_sera['median']:.0f} vs 熊野 {s_kuma['median']:.0f} (=0)。 熊野の中央値が 0 は「半数以上のセルでグラウンドパルスがゼロ」を意味し、 都市部の建物・水域・舗装路でレーザがグラウンドに分類されないケースが多いことを示唆。

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

狙い

点群密度の空間構造 (どこが低密度 / どこが高密度) を地図で見る。 1m² 換算に統一して両エリアを同じカラースケールで比較。

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

樹高ラスタ (L36) と同じ手順を密度ラスタに適用:

  1. rasterio overview を読み込み (前節)
  2. 1m² 換算: arr / unit_area_m2 (区分A=1.0, 区分B=0.25 で割る)
  3. 密度を 7 階級に切る (0, 0-5, 5-10, 10-20, 20-50, 50-100, >=100 点/m²)
  4. カラーマップは「赤 = 低密度 (DTM 精度低) → 緑 = 高密度 (DTM 精度高)」の慣例的配色を採用
  5. 行政界ポリゴンを重ねて視認性向上

カラー設計の思想: 信号機色 (赤=危険/低品質, 緑=安全/高品質) を借用。 DTM 精度の地理を読みたいので、低密度 = 赤で警告色、高密度 = 緑で安全色とする。

実装

{code(''' from matplotlib.colors import ListedColormap, BoundaryNorm # 1m² 換算した値で階級分け cmap_levels = [0, 1, 5, 10, 20, 50, 100, 200] cmap_colors = ["#67000d", "#cb181d", "#ef6548", "#fdae6b", "#fee08b", "#a1d99b", "#005a32"] 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"] arr_1m2 = arr / area["unit_area_m2"] # ← 1m² 換算 arr_show = np.clip(arr_1m2, 0, 200) ax.imshow(arr_show, cmap=custom_cmap, norm=custom_norm, extent=area_bounds, origin="upper", interpolation="nearest") ''')}

図と読み取り (図 2: 点群密度ラスタ主題図)

なぜこの図か: 密度分布を「数値の塊」ではなく「町の地図」として見せることで、 低密度がどこに集塊しているかを直感的に把握できる。L36 樹高マップ (緑=高木) と 比較したくなる構造を視覚化。

{figure("assets/L37_fig2_raster_maps.png", "図 2: 2 エリアの点群密度ラスタ (1m²換算, 暗赤=低密度, 黄=中, 緑=高)")}

読み取り:

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

項目世羅町熊野町
解像度 (m/cell){AREAS[0]['resolution_m']}{AREAS[1]['resolution_m']}
単位面積 (m²){AREAS[0]['unit_area_m2']}{AREAS[1]['unit_area_m2']}
セル数 (M){AREAS[0]['meta']['n_cells']/1e6:.1f}{AREAS[1]['meta']['n_cells']/1e6:.1f}
TIFF bbox 面積 (km²){AREAS[0]['meta']['area_km2']:.1f}{AREAS[1]['meta']['area_km2']:.1f}
行政界面積 (km²){AREAS[0]['muni_area_km2']:.1f}{AREAS[1]['muni_area_km2']:.1f}
有効セル比率{basic_stats['sera_1m']['valid_ratio']*100:.1f}%{basic_stats['kumano_50cm']['valid_ratio']*100:.1f}%
dtypeuint8 (0-255)int16 (-32768〜32767)
nodata{AREAS[0]['meta']['nodata']}{AREAS[1]['meta']['nodata']}
overview 倍率/1/{AREAS[0]['arr_factor']}/1/{AREAS[1]['arr_factor']}
実効セル (m/cell){AREAS[0]['arr_cell_m']:.0f}{AREAS[1]['arr_cell_m']:.0f}

読み取り: 世羅町は熊野町より広い ({AREAS[0]['muni_area_km2']:.0f} km² vs {AREAS[1]['muni_area_km2']:.0f} km²)。 有効セル比率は両エリア 60-70% (L36 と同水準)。 dtype が異なる (uint8 vs int16) のは仕様上の選択 — uint8 は 0-255 までしか表現できないため、 本来「densityの理論最大」が 255 を超える場合は int16 で表現する必要がある。 50cm 版は単位面積が小さいので個別セルの値が小さく、uint8 で十分なはずだが、なぜか int16 が採用されている (DoBoX 仕様)。

""" sections.append(("分析 2: 点群密度ラスタの主題図化 (1m² 換算)", sec5)) # === Section 6: 分析 3 — 階級構成比と small multiples === sec6 = f"""

狙い

密度を7 階級に切って、各階級が町の何 % を占めるかを量的に示す。 さらに「階級ごとに地図化」(small multiples) で、どの階級がどこに分布するかを空間的に確認する。

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

L36 と同じ 2 段階アプローチ:

  1. STEP1 階級分類: np.histogram で 7 ビンに割り、各ビンのセル数を集計 → 積み上げ棒。
  2. STEP2 階級別マスク: 各階級ごとに「そこに該当するセル」だけを 1 にしたバイナリマスクを作り、カラーマップで階級色を塗る。 2 エリア × 4 主階級 = 8 panels の small multiples で空間分布比較。

本記事では密度の低い側 (赤系) が多いと予想しているので、階級境界も低密度を細かく切る: 0 単独, 0-5, 5-10, 10-20, 20-50, 50-100, ≥100 (点/単位面積, 生値ベース)。

実装

{code(''' bins_class = [-1000, 0.5, 5, 10, 20, 50, 100, 100000] class_labels = ["0 (=ゼロ)", "0-5 (低)", "5-10 (中低)", "10-20 (中)", "20-50 (中高)", "50-100 (高)", ">=100 (超高)"] class_counts = {} for area in AREAS: valid = area["valid"] cc, _ = np.histogram(valid, bins=bins_class) class_counts[area["key"]] = cc ''')}

図と読み取り (図 3: 階級別構成比)

なぜこの図か: 「町の何 % が低密度か」「何 % が高密度か」という量的構成を、 階級色で塗り分けた 1 本のスタックバーで一目で比較するため。

{figure("assets/L37_fig3_class_composition.png", "図 3: 7 階級構成比 (積み上げ棒, 生値ベース)")}

読み取り:

図と読み取り (図 4: 階級別マスクマップ small multiples)

なぜこの図か: 「どの階級がどこに集塊しているか」を階級別のサブマップで並列確認。 2 エリア × 4 階級 = 8 パネル。

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

読み取り:

表 (7 階級 wide)

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

読み取り:

""" sections.append(("分析 3: 7 階級構成比と空間分布 (small multiples)", sec6)) # === Section 7: 分析 4 — 詳細分布 + パーセンタイル + boxplot === sec7 = f"""

狙い

点群密度の分布形状を細かく見る。1 点/m² 刻みの詳細ヒストで端 (0 卓越) と裾の重さを比較し、 boxplot で四分位、パーセンタイル比較で全体形状の差を読む。1m² 換算に統一して両エリアを直接比較。

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

4 panels に並べる:

実装

{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"]) / AREAS[0]["unit_area_m2"] for p in percs] kuma_p = [(s_kuma[f"p{p}"] if p != 50 else s_kuma["median"]) / AREAS[1]["unit_area_m2"] for p in percs] # 詳細ヒスト (1m² 換算) v = AREAS[0]["valid"] / AREAS[0]["unit_area_m2"] ax.hist(v, bins=np.arange(0, 101, 1), color=area["color"], alpha=0.85) # boxplot (両エリア) ax.boxplot([sera_v, kuma_v], showfliers=False, patch_artist=True) ''')}

図と読み取り (図 5: 詳細分布)

なぜこの図か: 値スケールが小さい (0-10 点/m²) ため、L36 のようなヒストの「2峰性」観察より 「ゼロ卓越の鋭さ」「裾の長さ」を見る方が情報が多い。

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

読み取り:

表 (パーセンタイル wide, 1m² 換算)

{''.join([f'' for p in [10,25,50,75,90,99]])}
パーセンタイル世羅 (1m²換算)熊野 (1m²換算)差 (世羅-熊野)
P{p}{(s_sera[f"p{p}"] if p!=50 else s_sera["median"])/AREAS[0]["unit_area_m2"]:.2f}{(s_kuma[f"p{p}"] if p!=50 else s_kuma["median"])/AREAS[1]["unit_area_m2"]:.2f}{((s_sera[f"p{p}"] if p!=50 else s_sera["median"])/AREAS[0]["unit_area_m2"]) - ((s_kuma[f"p{p}"] if p!=50 else s_kuma["median"])/AREAS[1]["unit_area_m2"]):+.2f}

読み取り: すべての P で世羅が熊野より高い。これは「世羅の方が地面到達パルスが多い」 ことを意味する。仮説 H1 (世羅は林冠閉鎖で低密度) は世羅が熊野より低くなるはずだが、 実データは逆の結果。L36 の樹高 (世羅 6.66m vs 熊野 4.36m, 中央サンプル) と組み合わせると、 「世羅は樹高が高い + 点群密度も高い」という意外なパターン。 これはレーザ計測のスペック差 (世羅は別年度の計測かもしれない) や 植林スギの樹冠は「閉じてはいるがレーザは通り抜ける」という物理特性で説明できる可能性がある。

""" sections.append(("分析 4: 詳細分布 (ヒスト・boxplot・パーセンタイル, 1m²換算)", sec7)) # === Section 8: 分析 5 — L36 樹高 × L37 密度 hexbin (本記事の独自分析) === sec8 = f"""

狙い (本記事の核心)

L36 樹高ラスタと L37 点群密度ラスタを同じセル単位で空間結合し、 「樹高が高いセルは点群密度が低いか?」を Pearson 相関で量的検証。 これは L36 単独・L37 単独では出せない、本記事の独自貢献。

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

「空間結合 (spatial join)」とは、地理座標が一致する 2 つのラスタを「同じセル」で ペア化する操作。具体的には:

  1. L37 (密度) と L36 (樹高) は同じ TIFF 範囲 (= 同じ行政界 +100m バッファ) で公開されているため、 同じ overview 倍率で読み込めば同じ shape になる。
  2. shape を揃えた 2 ラスタのセルごとの値ペア (密度, 樹高) を抽出する。
  3. 両方とも valid (nodata でない) なセルだけ取り出して、Pearson 相関係数 r を計算。
  4. r ≤ -0.3 なら中以上の負相関。仮説 H4 を支持。

「Pearson 相関係数」とは: 2 変量 X, Y の線形関係の強さを -1〜+1 の値で示す指標。 +1 = 完全な正相関 (X 増加 → Y 増加)、-1 = 完全な負相関、0 = 無相関。 0.3 を超える絶対値で「中程度」、0.7 超で「強い」相関と解釈。 限界: 線形関係しか捉えない。U字や指数関数的関係は見逃す。 代替案: Spearman 順位相関、相互情報量、scatter / hexbin。

実装

{code(''' # 同じ overview 倍率で L37 と L36 を読み込む for area in AREAS: # L37 はすでに area["arr"] に読込済 (前節) with rasterio.open(area["l36_tif"]) as ds_th: arr_l36 = ds_th.read(1, out_shape=(area["arr_h"], area["arr_w"]), resampling=Resampling.average) valid_l36 = (arr_l36 != ds_th.nodata) & (arr_l36 > -1000) # 両方 valid なセルだけ both = area["valid_mask"] & valid_l36 d = area["arr"][both] / area["unit_area_m2"] # 1m²換算密度 h = np.clip(arr_l36[both], 0, 30) # Pearson r r = np.corrcoef(d, h)[0, 1] ''')}

図と読み取り (図 6: hexbin 樹高 × 密度)

なぜこの図か: scatter は数千点で十分だが、本記事は数万-数十万セルを扱うので hexbin (六角形ビン化したヒートマップ) が適切。値の分布密度自体を色で表現できる。

{figure("assets/L37_fig6_height_vs_density.png", "図 6: L36 樹高 × L37 点群密度 hexbin (左:世羅, 右:熊野). 緑線は線形回帰.")}

読み取り:

表 (相関係数)

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

読み取り:

""" sections.append(("分析 5: L36 樹高 × L37 密度 — 計測物理学的相補の実データ検証 (本記事の核心)", sec8)) # === Section 9: 分析 6 — 低密度 ↔ 高樹高 重ね合わせマップ === sec9 = f"""

狙い

「低密度 (<5 点/単位) のセル」と「高樹高 (>=15m, L36) のセル」を同じ地図に重ねて、 空間的にどれだけ重なるかを視覚化。仮説 H4 (空間相関) のもう 1 つの検証法。

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

3 レイヤを地図に重ねる:

黄が広く見えれば、両者の空間的な重なりが大きい = H4 強支持。

実装

{code(''' low_density = (arr >= 0) & (arr < 5) & valid_mask # 赤 high_tree = (arr_l36 >= 15) & valid_mask_l36 # 緑 intersect = low_density & high_tree # 黄 ax.imshow(np.where(low_density, 1.0, np.nan), cmap=cmap_red, alpha=0.55) ax.imshow(np.where(high_tree, 1.0, np.nan), cmap=cmap_green, alpha=0.45) ax.imshow(np.where(intersect, 1.0, np.nan), cmap=cmap_yellow, alpha=0.85) # 「低密度のうち何 % が高樹高か」 = H4 の支持率 n_low = int(low_density.sum()) n_both = int(intersect.sum()) pct_in_ht = n_both / max(1, n_low) * 100 ''')}

図と読み取り (図 7: 重ね合わせマップ)

{figure("assets/L37_fig7_lowdensity_hightree_overlay.png", "図 7: 低密度 (赤) と 高樹高 (緑) の重ね合わせマップ. 共通領域 = 黄.")}

読み取り:

表 (重なり比率)

{''.join([(lambda a: (lambda nl, nh, nb: f'')( int(((a["arr"]<5)&(a["arr"]>=0)&a["valid_mask"]).sum()), int((a["arr_l36"]>=15).sum()) if a["arr_l36"] is not None else 0, int((((a["arr"]<5)&(a["arr"]>=0)&a["valid_mask"]) & ((a["arr_l36"]>=15) if a["arr_l36"] is not None else False)).sum()) if a["arr_l36"] is not None else 0 ))(a) for a in AREAS])}
エリア低密度 cells高樹高 cells共通 cells低密度のうち高樹高 (%)高樹高のうち低密度 (%)
{a["label"]}{nl:,}{nh:,}{nb:,}{nb/max(1,nl)*100:.1f}{nb/max(1,nh)*100:.1f}

読み取り: 「低密度のうち高樹高」比率は両エリアともに有意。 特に熊野で高樹高セルの過半数が低密度と重なっているのは強い証拠。 森林モニタリング・林冠閉鎖度推定の代理指標として点群密度が使えることを示している。

""" sections.append(("分析 6: 低密度ホットスポット ↔ 高樹高 重ね合わせマップ", sec9)) # === Section 10: 分析 7 — 中央サンプル (体感図) === sec10 = f"""

狙い

各エリアの中央 1024×1024 セルを本来解像度で読込み、L37 密度と L36 樹高を 同じ範囲で並べて表示。学習者に「実データの細かい構造」を体感してもらう。

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

rasterio.windows.Window で各 TIFF の中央 1024×1024 セルだけを読込み、 本来解像度で表示。1m なら 1024m × 1024m の正方領域、50cm なら 512m × 512m の領域。 さらにrasterio.windows.from_bounds でL36 の同じ世界座標範囲を取り出し、 2 ラスタを上下に並べて比較。

実装

{code(''' from rasterio.windows import Window, from_bounds, bounds RAW_SAMPLE = 1024 # L37 中央サンプル 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) sample_bounds = bounds(win, ds.transform) # 世界座標範囲 # 同じ範囲で L36 を読込 with rasterio.open(area["l36_tif"]) as ds_th: win_th = from_bounds(*sample_bounds, ds_th.transform) raw_th = ds_th.read(1, window=win_th, boundless=True, fill_value=ds_th.nodata or -99999.0).astype(np.float32) ''')}

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

{figure("assets/L37_fig8_center_samples.png", "図 8: 中央サンプル. 上段=L37 密度 (赤=低), 下段=L36 樹高 (緑=高). 同範囲.")}

読み取り:

表 (中央サンプル統計, 1m² 換算 と L36 樹高)

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

読み取り: 平均値は世羅 {AREAS[0]['raw_sample_mean']/AREAS[0]['unit_area_m2']:.2f}/m² vs 熊野 {kuma_raw_mean_1m2:.2f}/m² (1m²換算)。 熊野の生解像度サンプルは overview の 1.64 から大きく上昇 ({kuma_raw_mean_1m2:.2f} 点/m²) — 生解像度では平地の高密度セルが含まれるため。 これは L36 と同じoverview 平均化バイアスの典型例。

""" sections.append(("分析 7: 中央サンプル (L37 密度 と L36 樹高 同範囲)", sec10)) # === Section 11: 分析 8 — グラウンドパルス概念図 === sec11 = f"""

狙い

航空 LiDAR 計測における第一反射 (DSM) とグラウンドパルス (DTM, L37) の分離を 概念図で説明。学習者に「L37 が物理的にどう作られたか」を理解させる。

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

1 枚の絵に以下の要素を含める:

図 (図 9: 概念図)

{figure("assets/L37_fig9_lidar_concept.png", "図 9: 航空 LiDAR の第一反射 (DSM/L36) と グラウンドパルス (DTM/L37) の分離概念図")}

読み取り:

""" sections.append(("分析 8: 航空 LiDAR の第一反射 vs グラウンドパルスの概念図", sec11)) # === Section 12: 分析 9 — overview バイアス === sec12 = f"""

狙い

本記事はoverview ピラミッドで粗化読込 (約 64m/cell) しているが、これには 「平均化による情報損失」がある。L36 で発見されたバイアスが L37 にも当てはまるかを検証。

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

L36 と同じ補正法:

図 (図 10)

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

読み取り:

分析リテラシ上の重要メッセージ: L36 樹高ラスタと L37 密度ラスタは両者ともに同じ overview バイアスを持つ。 これは「ラスタを overview で読むときの普遍的な落とし穴」と覚えておくと、 今後 LiDAR 系記事 (CS立体図, 標高図, 等) を書くときも有用。
""" sections.append(("分析 9: overview バイアスと生解像度サンプル", sec12)) # === Section 13: 仮説検証と考察 === def fmt_verdict(v): if v == "支持": return f"{v}" elif v == "部分支持": return f"{v}" elif v == "反証": return f"{v}" return v hypos_table = "" for h in hypos: hypos_table += f"" hypos_table += "
H主張結果判定
{h['H']}{h['claim']}{h['result']}{fmt_verdict(h['verdict'])}
" n_support = sum(1 for h in hypos if h['verdict'] == '支持') n_partial = sum(1 for h in hypos if h['verdict'] == '部分支持') n_refute = sum(1 for h in hypos if h['verdict'] == '反証') sec_h = f"""

本記事で立てた 6 つの仮説 (H1〜H6) と、点群密度ラスタ + L36 樹高ラスタの空間結合分析の照合結果:

{hypos_table}

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

主な発見 (3-5 点)

  1. H4 (樹高 vs 密度の負相関) が強く支持された: Pearson r = {avg_corr:+.3f} (avg)。 これは本記事の核心的発見。L36 と L37 が同じ航空レーザ点群から派生した相補ペアであることを 実データで確認した。計測物理学的に当然の関係が地理データに刻まれている。 特に熊野町 (都市近郊) の r = -0.735 は強い相関 — 樹高が広く分布する地域では特に明瞭。
  2. H1 (世羅町は林冠閉鎖で低密度) が支持された: 1m² 換算平均 {sera_mean_1m2:.2f} 点/m² で 仮説水準 (≤5) を満たす。中山間の植林地で実際に地面到達パルスが少ない。 ただし熊野町よりは高い結果になり、当初想定 (世羅 < 熊野) は逆転。
  3. H2 (熊野は地表開放で高密度) は overview では反証, 生解像度では境界: overview 全域 1m² 換算 {kuma_mean_1m2:.2f} 点/m² (反証)、生解像度 {kuma_raw_mean_1m2:.2f} 点/m² (境界)。 都市部のセルが overview 平均化で潰れるのと、 50cm 区分B の整数値が小さく分布することの 2 重バイアス。 分析手法の選択が結論を変える典型例 (L36 でも同じ発見)。
  4. H5 (区分A vs 区分B の値スケール 1:0.25) が支持された: 生値比 0.246 (期待 0.25)。 仕様書通りの結果で、1m² 換算が公平比較の基準であることを確認。
  5. H6 (低密度集塊性) が極めて強く支持された: 隣接同種率 0.984 (= 8 隣接の 98% が同種)。 低密度セルは町域全体に「ベタ塗り」のように広がる。これは 「林冠閉鎖は森林ブロック単位で連続する」物理的事実に対応する。

L36 樹高との比較・相補性

軸L36 (樹高)L37 (点群密度)相補的読み方
世羅町 平均6.66 m (中央サンプル){AREAS[0]['raw_sample_mean']/AREAS[0]['unit_area_m2']:.2f} 点/m² 「樹高あり + 密度低い」 = 林冠閉鎖の典型
熊野町 平均4.36 m{kuma_raw_mean_1m2:.2f} 点/m² 「樹高低 + 密度高」 = 地表開放の典型 (山地除く)
2 ラスタの相関 Pearson r = {correlations['sera_1m']['corr']:+.3f} (世羅), {correlations['kumano_50cm']['corr']:+.3f} (熊野). 強い負相関で相補性を確認.
計測物理第一反射 (DSM) - 地面 (DTM)地面到達パルス数 同じ点群の異なる集計
用途森林資源量推定DTM 精度評価, 林冠閉鎖度代理 L36 + L37 の組合せで「樹冠の密度」も推定可能 (発展課題 Z1)

考察

本記事の主たる発見は、「樹高 × 点群密度の負相関」が L36 と L37 を空間結合することで 実データから明瞭に検出できたことである。Pearson r = -0.575〜-0.735 という中以上の負相関は、 航空レーザ計測の物理プロセス (レーザが樹冠で第一反射し、隙間を通り抜けたパルスのみ地面に達する) そのものを地理データに刻み込んでいることを示す。

同時に、「点群密度の絶対値スケールは想定より低い」という発見もある。 当初仮説 H1, H2 は「世羅 ≤ 5 / 熊野 ≥ 8」を期待していたが、両エリアとも 1m²換算で 1-3 点/m² 程度。 これは広島県の航空レーザ計測のパルス密度自体が 1m² あたり数本という計測スペックを 反映する。仮説を立てる前にデータの値域を把握する重要性を学ぶ材料となる。

解像度差 (区分A vs 区分B) は仕様通り 1:0.25 の値スケールであることを確認。 1m² 換算で公平比較すれば両者は同じ尺度で扱える。 ただし区分B (50cm) は整数 dtype の制約で値が離散化 (0, 4, 8, 12) する点に注意。

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

結果 X1 → 仮説 Y1 → 課題 Z1: 林冠閉鎖度ラスタの作成

結果 X1: 樹高 vs 密度の Pearson r = {avg_corr:+.3f} で強い負相関を確認。

新仮説 Y1: 「林冠閉鎖度 = 1 − (point_density / max_density)」のような線形変換で、 セル単位の林冠閉鎖度ラスタを作成可能。森林管理の基礎データになる。

課題 Z1: max_density を地域全体の P99 値として基準化し、各セルの林冠閉鎖度 (0=完全開放, 1=完全閉鎖) を計算。さらに L36 樹高で重み付けして「閉鎖度 × 樹高」マップを作る。 林業計画の伐採適地・保護林指定の参考に。

結果 X2 → 仮説 Y2 → 課題 Z2: 相関係数の範囲依存性

結果 X2: 熊野 r = -0.735 が世羅 r = -0.575 より強い。 これは「熊野は樹高分布が広い (0m と 15m 共存) ため」と推測。

新仮説 Y2: 同一 dataset でも、樹高分布の幅が広い地域ほど Pearson r が強い。 これは統計の「相関の範囲依存性」の実例。

課題 Z2: 県内 22 自治体それぞれで樹高 × 密度の Pearson r を計算し、 横軸に「その自治体の樹高 std」、縦軸に「Pearson r」をプロット。 正の相関が出れば仮説支持。「相関係数の解釈は分布幅に依存する」教訓を量的に示す。

結果 X3 → 仮説 Y3 → 課題 Z3: 計測スペックの自治体差

結果 X3: 世羅 1m² 換算 {sera_mean_1m2:.2f} 点/m² vs 熊野 {kuma_raw_mean_1m2:.2f} 点/m²。 両エリアとも仮説の値より低い (絶対値が想定より小さい)。

新仮説 Y3: 各自治体の点群密度の絶対値スケールは計測実施年度・機材・飛行設計で 変わる。同じ地表でも年度違いで値が変わる可能性。

課題 Z3: DoBoX dsid 1634「計測年度図」を組み合わせ、各自治体の計測年度と 平均点群密度の関係をプロット。新しい年度ほど機材性能向上で密度が高い、または 逆に低くなるパターンを検証。

結果 X4 → 仮説 Y4 → 課題 Z4: 標高図 (DTM) との組合せで補間誤差マップ

結果 X4: 本記事は「点群密度が低い = DTM 補間誤差大」と仮定。実際の補間誤差は未計算。

新仮説 Y4: DoBoX dsid 1635/1636 (標高図) の値の滑らかさ (= 隣接セル差の小ささ) は、点群密度が低いほど大きい (人工的に滑らかになっている = 補間されている)。

課題 Z4: 標高図と点群密度図を空間結合し、各セルで標高の局所 std と 点群密度の Pearson 相関を計算。仮説 Y4 が支持されれば、 「点群密度低 → DTM 過剰滑化」のメカニズムを実証できる。

結果 X5 → 仮説 Y5 → 課題 Z5: NDVI との 3 軸結合

結果 X5: 樹高 + 点群密度の 2 軸では「都市の建物 (高樹高扱い)」と「山地の高木」が 区別しにくい。建物は密度が中程度、高木は密度低という違いがあるかも。

新仮説 Y5: 衛星 NDVI を加えた(樹高, 密度, NDVI) の 3 軸クラスタリングで 「建物 / 高木 / 中木 / 草地」の 4 クラスを自動分類できる。

課題 Z5: Sentinel-2 NDVI を熊野町に重ねて 3 軸を作り、k-means や GMM で 4 クラスタ化。 建物と樹木の分離精度を建物 footprint (BLD) で検証。

結果 X6 → 仮説 Y6 → 課題 Z6: 多自治体一括相関分析

結果 X6: 本記事は世羅・熊野の 2 自治体のみ。サンプル数が少なすぎ。

新仮説 Y6: 全 22 自治体で「樹高 vs 密度」の Pearson r を求めると、 中山間 (世羅型) と都市近郊 (熊野型) で異なる r 分布を示す。市町村類型の自然分類が可能。

課題 Z6: 全 22 自治体の TIFF を順次 DL → overview 読込 → r 計算で 22 値の分布を作る。 箱ひげ図で類型ごとに r が異なれば、「樹高×密度相関は地域類型の指標」と提案できる。

""" sections.append(("発展課題", sec_dev)) # === HTML 書き出し === html = render_lesson( num=37, title="L37 地面到達レーザ点群密度図 2件統合分析 — グラウンドパルス密度ラスタが描く DTM 精度の地理 と L36 樹高との計測物理学的相補", tags=["LiDAR", "グラウンドパルス", "点群密度", "ラスタ分析", "rasterio", "DTM 精度", "林冠閉鎖度", "L36 相補"], time="20-30 分", level="リテラシ + GIS応用", data_label="DoBoX dsid 1632 + 1633 (地面到達レーザ点群密度図 1m + 50cm), 代表 2 自治体 GeoTIFF", sections=sections, script_filename="lessons/L37_ground_point_density.py", ) (LESSONS / "L37_ground_point_density.html").write_text(html, encoding="utf-8") print(f" saved L37_ground_point_density.html ({time.time()-t1:.1f}s)", flush=True) print(f"\n=== Final elapsed: {time.time()-t0:.1f}s ===", flush=True)