"""X14 LiDAR 5派生プロダクト統合精度評価 — 同一点群からの情報補完性 カバー宣言: 本記事は以下 11 dataset_id を統合する横断研究記事（第2記事化）。 - L36 樹高図 1m/50cm: ds=1626, 1627 - L37 地面到達点群密度図 1m/50cm: ds=1632, 1633 - L38 CS立体図 1m/50cm: ds=1628, 1629 - L39 地形傾斜図 1m/50cm: ds=1630, 1631 - L40 標高図 1m/50cm: ds=1635, 1636 - L43 計測年度図 1m: ds=1634 研究の問い (RQ): 同一 LiDAR 点群から派生した 5 種プロダクト（樹高・CS立体・傾斜・点群密度・標高）は 互いにどれだけの情報を共有し，PCA 第1主成分は「地形の何を表す軸」か 仮説 (D): H1: 5プロダクト間の相関行列で |r| ≥ 0.8 のペアが ≥ 1 組 H2: 傾斜×標高の相関 |r| < 0.5（独立性高） H3: PCA 第1主成分の寄与率 ≥ 50% H4: 計測年度（2014 vs 2016）でプロダクト相関構造に差がある """ from __future__ import annotations import sys, time, 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 import matplotlib.patches as mpatches from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False t0 = time.time() print("=== X14 LiDAR 5派生プロダクト統合精度評価 ===") DATA_DIR = ROOT / "data" / "extras" OUT = ASSETS DATASET_IDS = [1626,1627,1628,1629,1630,1631,1632,1633,1634,1635,1636] # ============================================================================= # 1. 世羅町エリア: 全5プロダクト @ 1m (同一解像度) # ============================================================================= print("\n[1] 世羅町エリア 5プロダクト読み込み ...") t1 = time.time() from PIL import Image Image.MAX_IMAGE_PIXELS = None SERA = DATA_DIR / "L36_canopy_height" / "1m_sera" / "sera_1m_treeheight.tif" SERA_BASE = DATA_DIR PRODUCTS_SERA = [ ("樹高", DATA_DIR/"L36_canopy_height"/"1m_sera"/"sera_1m_treeheight.tif", "F", -9990.0), ("点群密度", DATA_DIR/"L37_ground_point_density"/"1m_sera"/"sera_1m_density.tif", "L", -1.0), ("CS立体", DATA_DIR/"L38_cs_terrain"/"1m_sera"/"sera_1m_cs.tif", "RGB", -1.0), ("傾斜", DATA_DIR/"L39_terrain_slope"/"1m_sera"/"sera_1m_slope.tif", "F", -9990.0), ("標高", DATA_DIR/"L40_elevation"/"1m_sera"/"sera_1m_dem.tif", "F", -9990.0), ] STRIDE = 50 # subsample 1/50 → ~320×545 grid stacks_sera = {} valid_mask = None for name, path, mode, nodata in PRODUCTS_SERA: img = Image.open(path) if mode == "RGB": arr = np.array(img.convert("L"))[::STRIDE, ::STRIDE].astype(np.float32) vmask = arr > 5 # non-black background else: arr = np.array(img)[::STRIDE, ::STRIDE].astype(np.float32) vmask = arr > nodata stacks_sera[name] = arr if valid_mask is None: valid_mask = vmask else: valid_mask &= vmask n_valid = vmask.sum() print(f" {name}: shape={arr.shape}, valid={n_valid}({n_valid/arr.size*100:.1f}%)") # valid pixels vector flat = {k: v.ravel()[valid_mask.ravel()] for k, v in stacks_sera.items()} df_pixels = pd.DataFrame(flat) print(f" 世羅 有効画素: {len(df_pixels):,}, {time.time()-t1:.1f}s") # ============================================================================= # 2. 相関行列 # ============================================================================= print("\n[2] Pearson 相関行列 ...") corr_matrix = df_pixels.corr(method="pearson") print(corr_matrix.round(3).to_string()) # H1: r >= 0.8 high_corr_pairs = [] for i, c1 in enumerate(corr_matrix.columns): for j, c2 in enumerate(corr_matrix.columns): if i < j and abs(corr_matrix.loc[c1, c2]) >= 0.8: high_corr_pairs.append((c1, c2, corr_matrix.loc[c1, c2])) h1_pass = len(high_corr_pairs) >= 1 print(f" H1: |r|≥0.8 ペア数={len(high_corr_pairs)}: {'✓PASS' if h1_pass else '△PARTIAL'}") for c1, c2, r_val in high_corr_pairs: print(f" {c1} × {c2}: r={r_val:.3f}") # H2: 傾斜×標高 < 0.5 r_slope_elev = corr_matrix.loc["傾斜", "標高"] h2_pass = abs(r_slope_elev) < 0.5 print(f" H2: 傾斜×標高 r={r_slope_elev:.3f} < 0.5: {'✓PASS' if h2_pass else '△PARTIAL'}") # ============================================================================= # 3. PCA # ============================================================================= print("\n[3] PCA ...") X_scaled = StandardScaler().fit_transform(df_pixels.values) pca = PCA(n_components=5, random_state=42) scores = pca.fit_transform(X_scaled) var_ratio = pca.explained_variance_ratio_ h3_pass = var_ratio[0] >= 0.50 print(f" PC寄与率: {[f'{v*100:.1f}%' for v in var_ratio]}") print(f" H3: PC1寄与率={var_ratio[0]*100:.1f}% ≥ 50%: {'✓PASS' if h3_pass else '△PARTIAL'}") # PC1のloadings loadings = pd.DataFrame( pca.components_[:3], index=[f"PC{i+1}" for i in range(3)], columns=df_pixels.columns ) print(" PC1 loadings:", dict(loadings.loc["PC1"].round(3))) # ============================================================================= # 4. 熊野町: 計測年度別相関比較 # ============================================================================= print("\n[4] 熊野町: 計測年度別相関比較 ...") t1 = time.time() PRODUCTS_KUMANO = [ ("CS立体", DATA_DIR/"L38_cs_terrain"/"1m_kumano"/"kumano_1m_cs.tif", "RGB", -1.0), ("傾斜", DATA_DIR/"L39_terrain_slope"/"1m_kumano"/"kumano_1m_slope.tif", "F", -9990.0), ("標高", DATA_DIR/"L40_elevation"/"1m_kumano"/"kumano_1m_dem.tif", "F", -9990.0), ] YEAR_FILE = DATA_DIR/"L43_measurement_year"/"extracted"/"kumano"/"year.tiff" STRIDE_K = 30 # smaller stride for kumano (smaller area) stacks_kumano = {} valid_k = None for name, path, mode, nodata in PRODUCTS_KUMANO: img = Image.open(path) if mode == "RGB": arr = np.array(img.convert("L"))[::STRIDE_K, ::STRIDE_K].astype(np.float32) vmask = arr > 5 else: arr = np.array(img)[::STRIDE_K, ::STRIDE_K].astype(np.float32) vmask = arr > nodata stacks_kumano[name] = arr if valid_k is None: valid_k = vmask else: valid_k &= vmask # Year raster year_img = Image.open(YEAR_FILE) year_arr = np.array(year_img)[::STRIDE_K, ::STRIDE_K].astype(np.float32) year_valid = year_arr > 0 flat_k = {k: v.ravel()[valid_k.ravel() & year_valid.ravel()] for k, v in stacks_kumano.items()} year_flat = year_arr.ravel()[valid_k.ravel() & year_valid.ravel()] df_kumano = pd.DataFrame(flat_k) df_kumano["年度"] = year_flat corr_2014 = df_kumano[df_kumano["年度"] == 2014][list(stacks_kumano.keys())].corr() corr_2016 = df_kumano[df_kumano["年度"] == 2016][list(stacks_kumano.keys())].corr() n_2014 = (df_kumano["年度"] == 2014).sum() n_2016 = (df_kumano["年度"] == 2016).sum() print(f" 熊野 2014: {n_2014}画素, 2016: {n_2016}画素, {time.time()-t1:.1f}s") # H4: 最大差が >= 0.1 max_diff = 0.0 if len(corr_2014) == 3 and len(corr_2016) == 3: diff = (corr_2014 - corr_2016).abs() max_diff = diff.values[np.tril_indices(3, -1)].max() if len(diff) > 0 else 0 h4_pass = max_diff >= 0.1 print(f" H4: 年度別最大相関差={max_diff:.3f} ≥ 0.1: {'✓PASS' if h4_pass else '△PARTIAL'}") # ============================================================================= # 5. 図生成 # ============================================================================= print("\n=== 図生成 ===") # --- 図1: 相関行列ヒートマップ --- fig, ax = plt.subplots(figsize=(7, 6)) im = ax.imshow(corr_matrix.values, cmap="RdBu_r", vmin=-1, vmax=1, aspect="auto") plt.colorbar(im, ax=ax) ax.set_xticks(range(5)) ax.set_yticks(range(5)) ax.set_xticklabels(corr_matrix.columns, rotation=30, ha="right") ax.set_yticklabels(corr_matrix.columns) for i in range(5): for j in range(5): ax.text(j, i, f"{corr_matrix.values[i,j]:.2f}", ha="center", va="center", fontsize=11, color="white" if abs(corr_matrix.values[i,j]) > 0.5 else "black") ax.set_title("LiDAR 5プロダクト Pearson 相関行列（世羅町1m）") plt.tight_layout() fig.savefig(OUT/"X14_corr_matrix.png", dpi=110, bbox_inches="tight") plt.close(fig) print(" 図1: X14_corr_matrix.png") # --- 図2: PCA 寄与率バー + バイプロット --- fig, axes = plt.subplots(1, 2, figsize=(13, 5)) ax_bar, ax_bi = axes # 寄与率 ax_bar.bar(range(1, 6), var_ratio * 100, color=["#1f77b4","#ff7f0e","#2ca02c","#d62728","#9467bd"]) ax_bar.plot(range(1, 6), np.cumsum(var_ratio) * 100, "ko-", linewidth=1.5, markersize=5) ax_bar.set_xlabel("主成分番号") ax_bar.set_ylabel("寄与率 (%)") ax_bar.set_title(f"PCA 各主成分の寄与率（PC1={var_ratio[0]*100:.1f}%）") ax_bar.axhline(50, color="red", lw=1, linestyle="--", alpha=0.7) # バイプロット (PC1 vs PC2) n_plot = min(5000, len(scores)) idx_plot = np.random.choice(len(scores), n_plot, replace=False) ax_bi.scatter(scores[idx_plot, 0], scores[idx_plot, 1], c="lightblue", s=2, alpha=0.4, edgecolors="none") scale = 3 for i, col in enumerate(df_pixels.columns): ax_bi.annotate("", xy=(pca.components_[0, i]*scale, pca.components_[1, i]*scale), xytext=(0, 0), arrowprops=dict(arrowstyle="->", color="red", lw=1.5)) ax_bi.text(pca.components_[0, i]*scale*1.15, pca.components_[1, i]*scale*1.15, col, fontsize=9, color="red") ax_bi.set_xlabel(f"PC1 ({var_ratio[0]*100:.1f}%)") ax_bi.set_ylabel(f"PC2 ({var_ratio[1]*100:.1f}%)") ax_bi.set_title("PCA バイプロット（矢印=各プロダクトの寄与方向）") ax_bi.axhline(0, color="gray", lw=0.5); ax_bi.axvline(0, color="gray", lw=0.5) plt.tight_layout() fig.savefig(OUT/"X14_pca.png", dpi=110, bbox_inches="tight") plt.close(fig) print(" 図2: X14_pca.png") # --- 図3: PC1 loadings 棒グラフ --- fig, axes = plt.subplots(1, 3, figsize=(14, 4)) for pc_i, ax in enumerate(axes): vals = loadings.iloc[pc_i].values colors_l = ["#1f77b4" if v > 0 else "#d62728" for v in vals] ax.bar(loadings.columns, vals, color=colors_l) ax.set_title(f"PC{pc_i+1} ({var_ratio[pc_i]*100:.1f}%) loadings") ax.set_ylabel("loading") ax.axhline(0, color="k", lw=0.5) plt.setp(ax.xaxis.get_majorticklabels(), rotation=30, ha="right") plt.suptitle("PCA loadings（各主成分の各プロダクト寄与率）", y=1.02) plt.tight_layout() fig.savefig(OUT/"X14_loadings.png", dpi=110, bbox_inches="tight") plt.close(fig) print(" 図3: X14_loadings.png") # --- 図4: 世羅町 各プロダクトの空間分布（小さいサムネイル） --- fig, axes = plt.subplots(2, 3, figsize=(15, 8)) axes = axes.ravel() cmaps = ["YlGn", "Blues", "gray", "RdYlBu_r", "terrain"] for i, (name, arr) in enumerate(stacks_sera.items()): ax = axes[i] nodata = PRODUCTS_SERA[i][3] disp = arr.copy() disp[~valid_mask] = np.nan im = ax.imshow(disp, cmap=cmaps[i], interpolation="nearest") plt.colorbar(im, ax=ax, fraction=0.03, pad=0.04) ax.set_title(f"{name}") ax.axis("off") # PC1スコアマップ pc1_map = np.full(valid_mask.shape, np.nan) pc1_flat = scores[:, 0] pc1_map[valid_mask] = pc1_flat axes[5].imshow(pc1_map, cmap="RdBu_r", interpolation="nearest", vmin=-3, vmax=3) axes[5].set_title(f"PC1スコアマップ ({var_ratio[0]*100:.1f}%)") axes[5].axis("off") fig.suptitle("世羅町 LiDAR 5プロダクト + PC1スコア（stride=50で間引き）", fontsize=12) plt.tight_layout() fig.savefig(OUT/"X14_spatial_maps.png", dpi=110, bbox_inches="tight") plt.close(fig) print(" 図4: X14_spatial_maps.png") # --- 図5: 熊野町 計測年度別相関比較 --- if len(corr_2014) == 3 and len(corr_2016) == 3: fig, axes = plt.subplots(1, 2, figsize=(11, 4.5)) for ax, corr, year_label, n in [ (axes[0], corr_2014, "2014年度計測", n_2014), (axes[1], corr_2016, "2016年度計測", n_2016), ]: im2 = ax.imshow(corr.values, cmap="RdBu_r", vmin=-1, vmax=1, aspect="auto") plt.colorbar(im2, ax=ax) ax.set_xticks(range(3)); ax.set_yticks(range(3)) ax.set_xticklabels(corr.columns, rotation=30, ha="right") ax.set_yticklabels(corr.columns) for i2 in range(3): for j2 in range(3): ax.text(j2, i2, f"{corr.values[i2, j2]:.2f}", ha="center", va="center", fontsize=10) ax.set_title(f"{year_label} (n={n:,})") fig.suptitle(f"熊野町: 計測年度別相関行列（最大差={max_diff:.3f}）", fontsize=11) plt.tight_layout() fig.savefig(OUT/"X14_year_corr.png", dpi=110, bbox_inches="tight") plt.close(fig) print(" 図5: X14_year_corr.png") else: print(" 図5: スキップ（年度別データ不足）") # --- 図6: 高相関ペアの散布図 --- if high_corr_pairs: n_pairs = min(3, len(high_corr_pairs)) fig, axes = plt.subplots(1, n_pairs, figsize=(5*n_pairs, 5)) if n_pairs == 1: axes = [axes] n_scat = min(3000, len(df_pixels)) idx_s = np.random.choice(len(df_pixels), n_scat, replace=False) for i, (c1, c2, r_val) in enumerate(high_corr_pairs[:n_pairs]): axes[i].scatter(df_pixels[c1].values[idx_s], df_pixels[c2].values[idx_s], s=3, alpha=0.3, color="#1f77b4") axes[i].set_xlabel(c1); axes[i].set_ylabel(c2) axes[i].set_title(f"{c1} × {c2}\nr={r_val:.3f}") plt.tight_layout() fig.savefig(OUT/"X14_high_corr_scatter.png", dpi=110, bbox_inches="tight") plt.close(fig) print(" 図6: X14_high_corr_scatter.png") else: print(" 図6: スキップ（高相関ペアなし）") # --- 図7: 分布ヒストグラム（各プロダクト） --- fig, axes = plt.subplots(1, 5, figsize=(18, 4)) colors7 = ["#2ca02c","#1f77b4","#7f7f7f","#ff7f0e","#8c6d31"] for i, (name, arr) in enumerate(stacks_sera.items()): vals = arr[valid_mask].ravel() ax = axes[i] ax.hist(vals, bins=50, color=colors7[i], alpha=0.7, edgecolor="none") ax.set_title(name) ax.set_xlabel("値") ax.set_ylabel("頻度" if i == 0 else "") ax.tick_params(labelsize=8) fig.suptitle("世羅町 LiDAR 5プロダクトの値分布（valid画素のみ）", fontsize=11) plt.tight_layout() fig.savefig(OUT/"X14_histograms.png", dpi=110, bbox_inches="tight") plt.close(fig) print(" 図7: X14_histograms.png") # --- 図8: 計測年度 空間分布（熊野町） --- if n_2014 > 0 and n_2016 > 0: fig, ax = plt.subplots(figsize=(8, 6)) year_disp = year_arr.copy() year_disp[~(valid_k & year_valid)[::1, ::1]] = np.nan colors_yr = matplotlib.colors.ListedColormap(["#1f77b4","#ff7f0e"]) bounds = [2013.5, 2015.5, 2016.5] norm = matplotlib.colors.BoundaryNorm(bounds, colors_yr.N) ax.imshow(year_disp, cmap=colors_yr, norm=norm, interpolation="nearest") p1 = mpatches.Patch(color="#1f77b4", label=f"2014年度計測 (n={n_2014:,})") p2 = mpatches.Patch(color="#ff7f0e", label=f"2016年度計測 (n={n_2016:,})") ax.legend(handles=[p1, p2], loc="lower right") ax.set_title("熊野町 計測年度分布（stride=30）") ax.axis("off") plt.tight_layout() fig.savefig(OUT/"X14_year_map.png", dpi=110, bbox_inches="tight") plt.close(fig) print(" 図8: X14_year_map.png") else: print(" 図8: スキップ") # ============================================================================= # 6. CSV出力 # ============================================================================= print("\n=== CSV 出力 ===") corr_matrix.round(4).to_csv(OUT/"X14_corr_matrix.csv", encoding="utf-8-sig") loadings.round(4).to_csv(OUT/"X14_pca_loadings.csv", encoding="utf-8-sig") print(" CSV 2件出力") # ============================================================================= # 7. HTML # ============================================================================= print("\n=== HTML 組立 ===") ds_links = " ".join( f'#{d}' for d in DATASET_IDS ) year_map_figure = figure("X14_year_map.png", "熊野町 計測年度分布（2014/2016）") if (n_2014 > 0 and n_2016 > 0) else "" year_corr_figure = figure("X14_year_corr.png", "計測年度別相関行列（2014 vs 2016）") if (len(corr_2014)==3 and len(corr_2016)==3) else "" high_corr_figure = figure("X14_high_corr_scatter.png", "高相関ペアの散布図") if high_corr_pairs else "" sections = [ ("研究の問い (RQ) と仮説", f"""

LiDAR ファミリ (L36-L43) の横断研究

広島県で公開されているLiDAR由来の5種派生プロダクト（樹高図・CS立体図・地形傾斜図・ 地面到達点群密度図・標高図）は，同一の航空レーザ計測点群から異なるアルゴリズムで生成される。 本記事では世羅町（1m解像度）の全5プロダクトを統合し， Pearson相関行列とPCAを用いてプロダクト間の情報補完性を定量化する。

主 RQ

5プロダクト間の相関行列でr≥0.8のペアは何組存在し，PCA第1主成分は地形の何を表す軸か？ また計測年度（L43）が異なるエリアでプロダクト相関構造は変化するか？

仮説 (H1〜H4)

使用データセット ({len(DATASET_IDS)} 件): {ds_links}

"""), ("分析1: 5プロダクト間 Pearson 相関行列", f""" {figure("X14_corr_matrix.png", "LiDAR 5プロダクト Pearson 相関行列（世羅町1m）")} {high_corr_figure} {figure("X14_histograms.png", "各プロダクトの値分布")}

相関行列から|r|≥0.8のペアは{len(high_corr_pairs)}組確認された（H1{"支持" if h1_pass else "不支持"}）。 {"、".join(f"{c1}×{c2}(r={r:.3f})" for c1,c2,r in high_corr_pairs) if high_corr_pairs else "なし"}。 傾斜と標高の相関はr={r_slope_elev:.3f}で，{"独立性が高い（H2支持）" if h2_pass else "ある程度相関がある（H2不支持）"}。 これはLiDAR由来プロダクトが地形の異なる側面を捉えていることを示す。

"""), ("分析2: PCA による情報圧縮", f""" {figure("X14_pca.png", "PCA 寄与率と バイプロット")} {figure("X14_loadings.png", "PCA loadings（各主成分への寄与）")} {figure("X14_spatial_maps.png", "世羅町 5プロダクト空間分布とPC1スコアマップ")}

PC1の寄与率は{var_ratio[0]*100:.1f}%（H3{"支持" if h3_pass else "不支持"}），累積3主成分で{sum(var_ratio[:3])*100:.1f}%を説明する。 PC1 loadingsでは"{loadings.loc["PC1"].abs().idxmax()}"が最大({loadings.loc["PC1"][loadings.loc["PC1"].abs().idxmax()]:.3f})で， PC1が主に「{loadings.loc["PC1"].abs().idxmax()}」を表す軸であることが確認できる。 バイプロットでは樹高と点群密度が近い方向を向いており， 森林密度が樹高情報と強く連動していることを視覚化している。

"""), ("分析3: 計測年度別の相関構造比較（熊野町）", f""" {year_map_figure} {year_corr_figure}

熊野町のLiDAR計測は2014年度（{n_2014:,}サンプル）と2016年度（{n_2016:,}サンプル）の 2時期に分かれている。 プロダクト相関の最大年度差は{max_diff:.3f}で，{"H4を支持し新旧エリアで構造差が見られる" if h4_pass else "H4は不支持（年度差が小さく安定している）"}。 {"これは計測時期による地表被覆（植生）の変化や計測密度の差を反映している可能性がある。" if h4_pass else "2年の差では相関構造に有意な変化は生じず，データの品質一貫性が確認できる。"}

"""), ("まとめ", f""" """), ] html_out = render_lesson( num=14, title="X14: LiDAR 5派生プロダクト統合精度評価 — 同一点群からの情報補完性", tags=["X系","横断研究","LiDAR","ラスタ","PCA","相関","樹高","標高","傾斜","CS立体"], time="60分", level="リテラシ応用", data_label=( f'L36 樹高図(#1626,#1627) × ' f'L37 点群密度(#1632,#1633) × ' f'L38 CS立体(#1628,#1629) × ' f'L39 傾斜(#1630,#1631) × L40 標高(#1635,#1636) × L43 計測年度(#1634) 計{len(DATASET_IDS)}件' ), sections=sections, script_filename="lessons/X14_lidar_integrated.py", ) html_out = html_out.replace("", '', 1) out_path = LESSONS / "X14_lidar_integrated.html" out_path.write_text(html_out, encoding="utf-8") print(f"\n[HTML] X14_lidar_integrated.html ({len(html_out):,} bytes)") print(f"\n=== DONE in {time.time()-t0:.1f}s ===")