L41 航空レーザ森林資源 3件統合分析 — 樹種ポリゴン・単木ポイント・20m mesh の地理単位差を読む

データ取得手順

✅ このスクリプトは初回実行時にデータを自動取得します（DoBoX からの直接ダウンロード）。

実行コマンド:

cd "2026 DoBoX 教材"
python -X utf8 lessons/L41_forest_resources.py

DoBoX のオープンデータは申請不要・商用/非商用とも利用可。 data/extras/ は .gitignore 対象（約 57 GB のキャッシュ）。 スクリプト実行で自動再生成されます。

学習目標と問い

研究の問い (RQ)

広島県内の航空レーザ計測に基づく森林資源 3 dataset:

• (a) 樹種ポリゴン (dsid 1519): 市町村単位の polygon で「どの面が何の樹種か」を示す
• (b) 単木ポイント (dsid 1520): 市町村単位の point で「個々のスギ・ヒノキ単木の位置・樹高・直径」を示す
• (c) 森林資源量集計メッシュ (dsid 1521): 国土基本図図郭の 20m mesh で「林分指標 (材積/密度/形状比など 10+ 列)」を示す

これら 3 dataset は同じ航空レーザ計測 (2018-2019) から派生したが、 地理単位が異なる (市町別 vs 図郭別) ため、単純に空間結合できない。 本記事は「3 dataset を統合すると、広島県の森林タイプ・林相成熟度・林業生産性の地理がどう描けるか」を量的に検証する。

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

立てた仮説 (H1〜H7)

到達点

• 3 dataset の地理編成方式の違いを量的・視覚的に把握する
• 樹種構成 (件数 vs 面積) と単木形状 (樹高/直径) のアロメトリを実データで読む
• 20m mesh の10 種以上の林分指標を樹種別に統計化し、収穫期・健全度を判定
• ha材積の地理集塊性を量化 (8 近傍同種率)
• 3 dataset の統合上の制約 (= 地理単位の差) を量的に確認し、方法論的提言を示す

使用データ

本記事では DoBoX の「航空レーザ計測に基づく森林資源データ」シリーズ 3 件を統合する。 すべて広島県林業課が 森林管理基盤情報 (公共測量 2 次著作物, 測量法 44 条) として公開した GeoPackage 形式のベクタデータ。CRS はEPSG:6671 (JGD2011 平面直角座標系 III 系)。

3 dataset の地理編成の違い (本記事の核心)

• 1519 樹種ポリゴン: 市町村単位 polygon。10 市町公開。 1 市町ファイル ≈ 数千〜2 万 polygon。属性は樹種・面積・計測年・計測法。
• 1520 単木ポイント: 市町村単位 point。10 市町公開。 1 市町ファイル ≈ 数十万〜数百万 point (世羅町は 3,397,186 point)。 属性は樹種・樹高・胸高直径・単木材積・形状比・樹冠長率・樹冠面積。 注: スギ・ヒノキと判読された区域内のみ抽出。
• 1521 森林資源量集計メッシュ: 国土基本図図郭単位の 20m mesh polygon。10 図郭公開。 1 図郭ファイル = 750,000 mesh (= 20km × 15km ÷ 20m × 20m)。 うち valued (= 林分指標が計測された) は 約 13 万件。 属性は樹種・面積・立木本数・立木密度・平均樹高・平均直径・合計材積・ha材積・ 収量比数・相対幹距比・形状比・樹冠長率 (12 列)。

重要な構造的差異: 1519 と 1520 は同じ市町村で公開されているため 直接統合可能 (本記事は世羅町ペアを採用)。一方 1521 は図郭単位なので、 市町村と必ずしも一致しない。本記事の 1521 は 03NF2 図郭 (庄原市域中心) を採用したため、 1520 (世羅町) とは地理が重ならない。これが本データセットの最大の方法論的制約であり、 本記事 H7 で量的に確認する。

共通メタデータ

• 計測期間: 2018-2019 (広島県・国土交通省・林野庁の合同観測)
• 森林計測法: 1 (= 第一手法, 単一バージョン)
• ライセンス: CC-BY 4.0 (測量法 44 条 2 次著作物として広島県林業課が取りまとめ)
• 使用承認: 広島県砂防課 (令和5年砂防第136号), 林野庁治山課 (令和5年林整治第1041号), 国土交通省中国地方整備局 (令和4年国中整水予第37号)

L36-L40 との関係 (LiDAR ファミリ)

本記事の 3 dataset は L36-L40 のラスタ系 5 dataset と同じ航空レーザ計測から派生する:

• L36 樹高図 (DCHM, ラスタ) → 樹頂点抽出 → 本記事 1520 単木ポイント
• L40 標高図 (DTM, ラスタ) → DCHM 計算 → 樹高 → 本記事 1520, 1521
• L37 点群密度図 (ラスタ) → 計測精度の指標 → 本記事 1521 形状比・収量比数の信頼性
• L38 CS立体図 / L39 傾斜 (ラスタ) → 微地形・斜面方位 → 本記事の樹種立地条件

つまり LiDAR ファミリは「点群 → ラスタ (L36-L40) → 林業ベクタ (本記事)」という派生階層を成す。 本記事の 3 dataset は林業実務に直結する応用層であり、森林経営計画 (森林法第 5 条) の基礎資料となる。

ダウンロード（再現用データ・中間データ・図）

本記事の再現に必要なデータ・中間データ・図はすべて以下から直 DL 可能。

原データ (DoBoX 直リンク)

• 1519 樹種ポリゴン (世羅町): 直DL (45 MB ZIP → 91 MB GeoPackage)
• 1520 単木ポイント (世羅町): 直DL (107 MB ZIP → 552 MB GeoPackage)
• 1521 森林資源量メッシュ (03NF2): 直DL (20 MB ZIP → 159 MB GeoPackage)

中間データ (本スクリプトが生成)

• L41_dataset_meta.csv — 3 dataset のメタ情報比較
• L41_species_summary.csv — 1519 樹種別件数・面積
• L41_tree_point_by_species.csv — 1520 樹種別単木統計
• L41_allometry.csv — 1520 樹高 vs 胸高直径 アロメトリ係数
• L41_mesh_by_species.csv — 1521 樹種別 mesh 集計
• L41_volume_class.csv — 1521 ha材積 5 階級構成
• L41_high_volume_pct.csv — 1521 高材積 (≥500) 比率
• L41_harvest_ready.csv — 1521 収穫期セルの集計
• L41_clustering.csv — 1521 樹種別 8 近傍同種率

図 (PNG)

• 図1: 3 dataset 編成構造
• 図2: 樹種ポリゴン主題図
• 図3: 樹種構成 棒グラフ
• 図4: 単木 樹高・直径 ヒスト
• 図5: 樹高 vs 直径 アロメトリ
• 図6: 単木サンプル位置マップ
• 図7: ha材積 choropleth
• 図8: 林分指標 box plot
• 図9: ha材積 階級分類
• 図10: 収量比数 vs 平均樹高
• 図11: 形状比 ヒスト
• 図12: 樹種別 small multiples

再現スクリプト

L41_forest_resources.py をダウンロードし、data/extras/L41_forest_resources/ にデータがなければ自動取得から始まります。
実行:

cd "2026 DoBoX 教材"
py -X utf8 lessons\L41_forest_resources.py

分析 1: 3 dataset の構造比較 (本記事の核心)

狙い

3 dataset の地理編成・属性・スケールを 1 枚の表にまとめ、視覚的にも対比する。 本記事の核心命題「同じ航空レーザから派生したが地理単位が異なる」を量的に確認する。

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

pyogrio.read_info() で各 GeoPackage のメタ (要素数・列数・bbox) を取得する。 これは GIS データ俯瞰の標準手法で、ファイル全体を読まずにスキーマだけを軽量に確認できる。

• 入力: GeoPackage パス + レイヤ名
• 出力: dict (features, fields, dtypes, crs, total_bounds)
• 限界: feature 値の分布や geometry 詳細は読めない (それは read_dataframe で別途)

実装

↑ L41_forest_resources.py 行 1294–1309

1294
1295
1296
1297
1298

import pyogrio
for k, d in DATASETS.items():
    info = pyogrio.read_info(d["gpkg"], layer=d["layer"])
    print(f"{d['label']}: features={info['features']:,}, "
          f"fields={len(info['fields'])}, bbox_w={info['total_bounds'][2]-info['total_bounds'][0]:.0f}m")

図と読み取り (図 1: 3 dataset 編成構造概念図)

なぜこの図か: 3 dataset の bbox 大きさ・要素数・属性数を 1 枚で比較。 矩形の面積で地理範囲、注釈テキストで要素数・列数を読み取る。

図 1: 3 dataset の地理編成 — 同じ航空レーザから派生したが異なる地理単位

読み取り (重要発見):

• 1519 樹種ポリゴン (世羅町): bbox 27.0×16.3km, polygon 13,025, 列 7。 = 行政界スケール、属性は素朴 (樹種・面積)。
• 1520 単木ポイント (世羅町): bbox 26.7×16.0km, point 3,397,186, 列 12。 1519 とほぼ同じ bbox だが、要素数は 261 倍 (≈ 各ポリゴン内の単木抽出のため)。
• 1521 mesh (03NF2 図郭): bbox 20×15km, mesh 750,000, 列 17。 属性数が圧倒的に多い (12 vs 7) = 林分指標が豊富。
• 3 dataset の地理単位は重ならない: 1519/1520 は世羅町, 1521 は 03NF2 図郭 = 庄原市域。 = 「同じ広島県林業データだが、3 dataset を 1 自治体で全部見ることはできない」。
• これは森林データの設計思想として「市町村単位 (地方自治) vs 国土基本図単位 (測量行政)」の 2 系統が並立していることの量的反映。

表 (3 dataset メタ情報比較)

読み取り: 列 geographic_unit の差を読むと、本記事の核心命題が量的に確認できる。 1519/1520 は「市町村」(自治単位), 1521 は「国土基本図図郭」(測量単位) で公開される。 1521 のメッシュは EPSG:6671 (JGD2011 III 系) の原点から 20m 等間隔で規則格子を切るため、 mesh の position だけで bbox が決まり (ファイル名で識別可能)、行政界とは無関係。 これは GIS の「主題ベース vs 規則格子ベース」という古典的二元論の例である。

分析 2: 樹種ポリゴン (1519) の構成 (H1)

狙い

1519 樹種ポリゴンを読み込み、世羅町の森林の樹種構成を量化する。 件数ベースと面積ベースの 2 角度で見て、ヒノキ/スギ/その他 の地理的優占を確認 (H1 検証)。

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

geopandas で GeoPackage を全件読込 → groupby で樹種別の count と area を合計。 さらに件数ベースと面積ベースの構成比を別々に計算 = 「件数の多い樹種」 ≠ 「面積の大きい樹種」 である可能性を見る。

• 入力: GeoPackage (Polygon, 13,025 行 × 7 列 + geometry)
• 出力: 樹種別の polygon 件数, 総面積_ha, 構成比 (count_pct, area_pct)
• 限界: 「その他」カテゴリ内訳 (広葉樹各種) は本データには無い

実装

↑ L41_forest_resources.py 行 1351–1378

 1
 2
 3
 4
 5
 6
 7
 8
 9
1360
1361

import geopandas as gpd
gdf_ts = gpd.read_file("tree_species_34462.gpkg", layer="tree_species_34462")
species_count = gdf_ts["樹種"].value_counts()
species_area = gdf_ts.groupby("樹種")["面積_ha"].sum().sort_values(ascending=False)
species_summary = pd.concat([
    species_count.rename("polygon_count"),
    species_area.rename("total_ha"),
], axis=1).fillna(0)
species_summary["count_pct"] = species_summary["polygon_count"] / species_summary["polygon_count"].sum() * 100
species_summary["area_pct"]  = species_summary["total_ha"]      / species_summary["total_ha"].sum()      * 100
species_summary["mean_ha"]   = species_summary["total_ha"]      / species_summary["polygon_count"]

図と読み取り (図 2: 樹種ポリゴン主題図)

なぜこの図か: 樹種構成の地理分布を見る。緑=ヒノキ、青=スギ、灰=その他 で塗り分け、 世羅町行政界を黒線で囲む。集落・地形との関係も読める。

図 2: 樹種ポリゴン主題図 (1519 世羅町)

読み取り:

• 世羅町全域に緑 (ヒノキ) と灰 (その他) が斑模様に分布。スギ (青) は局所集塊で散在。
• 町の周辺部 (山地) に針葉樹植林が多く、町中央 (台地・農地) は「その他」(広葉樹混交林・里山) が多い傾向。
• polygon のサイズ分布が広い: 小さな点状から町を横断する大きな polygon まで。 = 「林分」が均一ではなく、管理単位や微地形に応じて細分化されている。

図と読み取り (図 3: 樹種構成 棒グラフ)

なぜこの図か: 件数と面積の 2 角度で樹種の優占を確認。 件数 vs 面積の差が「polygon サイズの樹種偏り」を示す。

図 3: 樹種構成 — 件数 vs 面積 (1519 世羅町)

読み取り (重要発見):

• 件数ベース: ヒノキ 50.8% > その他 23.5% > スギ 19.8%。 ヒノキが約半数を占める優占樹種。
• 面積ベース: ヒノキ 15.6% > その他 74.2% > スギ 2.5%。 件数比 ≈ 面積比で、樹種ごとの平均 polygon サイズはほぼ同等。
• 各樹種の平均 polygon 面積: ヒノキ 0.46ha, スギ 0.19ha, その他 4.72ha。 = 「その他」polygon は10.3倍大きく、 広葉樹混交林は連続性が高い (大きな塊で分布する) ことを示唆。
• H1 (ヒノキ ≥ 40% かつ スギ ≤ 30% かつ その他 ≤ 30%) の判定: 支持。 ヒノキ 50.8% と「その他」23.5% の関係次第で「ヒノキ優占」の主張強度が変わる。
• 世羅町の森林は「針葉樹植林 (ヒノキ+スギ) ≈ 71% + 広葉樹等 ≈ 24%」と読める。 これは中山間林業地として典型的な構成。

表 (1519 樹種別集計)

分析 3: 単木ポイント (1520) の樹高・直径分布 (H2)

狙い

1520 単木ポイント (3,397,186 点) を全件属性読込し、 スギとヒノキの樹高・胸高直径・単木材積の分布を比較。 樹種ごとの形状特性と生長到達点の差を確認 (H2 検証)。

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

pyogrio.read_dataframe(read_geometry=False) で属性のみ全件読込。 geometry を読まないので 552 MB GPKG でも~10 秒で全件取得可能。 これは GIS データ「巨大なら属性だけ読む」という標準的高速戦略。

• 入力: GeoPackage (Point, 3.4M 行 × 12 列 + geometry)
• 出力 (属性のみ): pandas DataFrame (3.4M 行 × 12 列)
• 限界: 空間結合や地図表示には別途 geom 読込が必要 (本記事 図 6 のサブサンプル)

実装

↑ L41_forest_resources.py 行 1429–1462

 1
 2
 3
 4
 5
 6
 7
 8
 9
1438
1439
1440
1441
1442

import pyogrio
# geometry なしで全件属性読込 (高速)
tp_attr = pyogrio.read_dataframe("tree_point_34462.gpkg",
                                 layer="tree_point_34462",
                                 read_geometry=False)
# 樹種別集計
tp_by_species = tp_attr.groupby("樹種").agg(
    n=("樹高", "count"),
    height_mean=("樹高", "mean"),
    height_median=("樹高", "median"),
    dbh_mean=("胸高直径", "mean"),
    vol_mean=("単木材積", "mean"),
    vol_total=("単木材積", "sum"),
)

図と読み取り (図 4: 単木 樹高・胸高直径ヒスト)

なぜこの図か: 樹種ごとの分布形状を比較する。 ヒストグラムのピーク位置・裾の長さ・複峰性が樹種特性を示す。

図 4: 単木 樹高・胸高直径分布 (1520 世羅町, 樹種別)

読み取り (重要発見):

• 樹高分布: スギ平均 21.0m (中央 21.6m, std 5.7m), ヒノキ平均 14.1m (中央 14.1m, std 5.0m)。 差 = +6.8mでスギが高い。
• 胸高直径分布: スギ平均 28.1cm, ヒノキ平均 20.2cm。 スギは胸高直径も+7.9cm大きい。
• 単木材積: スギ平均 0.711 m³ (中央 0.574), ヒノキ平均 0.311 m³ (中央 0.210)。 倍率 = ×2.3でスギが大きい。
• 分布形状: 両樹種ともピーク約 8-12m 付近の右裾長分布。 = 多くは中木 (10-20m), 少数の高木 (30m+) を持つ。
• H2 (スギ +4m 以上, 材積 ×2 以上) の判定: 支持。 樹高条件は満たす (差 6.8m), 材積条件は満たす (倍 2.3)。
• これは林学の常識 (スギは早生・大径) を 3.4M 単木で量的に裏付ける結果。 かつサンプル数が桁違い (従来の毎木調査は数百本程度) なので統計的検定力が圧倒的。

表 (1520 樹種別単木統計)

分析 4: 樹高 × 胸高直径 アロメトリ (H3, 本記事の生物物理学的核心)

狙い (本記事の生物物理学的核心)

樹木は「樹高 h ∝ 胸高直径 d^β」というべき乗法則 (アロメトリ) に従う。 1520 単木 3.4M 点で h-d 関係の Pearson 相関と回帰係数を計算し、 スギとヒノキで係数 c, β がどう違うかを検証 (H3)。

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

アロメトリ (allometry) = 「異速成長」。樹木では:

• 線形モデル: h = a + b · d (簡易だが大径木で乖離)
• べき乗モデル: h = c · d^β (生物学的に正当)
• 両辺対数を取れば log(h) = log(c) + β · log(d) という線形回帰

係数 β は樹種固有で、典型的に 0.5-1.0。β < 1 は「太くなる速度が高くなる速度より早い」(寸胴), β = 1 は等比 (細長), β > 1 は「高さが速い」(超細長)。

• 入力: 樹種ごとの (h, d) ペア (スギ 33万件, ヒノキ 307万件)
• 出力: Pearson r, 線形係数 (a, b), 対数べき乗係数 (c, β)
• 限界: アロメトリ係数は地域・管理状態で変動、本研究の値は世羅町固有

実装

↑ L41_forest_resources.py 行 1504–1554

 1
 2
 3
 4
 5
 6
 7
 8
 9
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525

import numpy as np
allometry_rows = []
for sp in ["スギ", "ヒノキ"]:
    sub = tp_clean[tp_clean["樹種"] == sp]
    h = sub["樹高"].values.astype(np.float64)
    d = sub["胸高直径"].values.astype(np.float64)
    # 線形回帰 h = a + b * d
    X = np.vstack([np.ones(len(d)), d]).T
    coef, *_ = np.linalg.lstsq(X, h, rcond=None)
    a, b = coef
    # 対数べき乗モデル h = c * d^β
    log_h, log_d = np.log(h), np.log(d)
    Xl = np.vstack([np.ones(len(log_d)), log_d]).T
    coef2, *_ = np.linalg.lstsq(Xl, log_h, rcond=None)
    log_c, exp_d = coef2
    allometry_rows.append({
        "species": sp,
        "n": len(sub),
        "r_pearson": float(np.corrcoef(h, d)[0, 1]),
        "lin_intercept": float(a), "lin_slope": float(b),
        "log_c": float(np.exp(log_c)), "log_exp": float(exp_d),
    })

図と読み取り (図 5: 樹高 × 胸高直径 散布 + 対数回帰線)

なぜこの図か: アロメトリ関係を視覚的に確認する。 散布点で雲の形を、回帰線でべき乗モデルの中心線を見る。 30k 点の subsample 表示と回帰線で「3.4M 全体の傾向」を読む。

図 5: 樹高 × 胸高直径 アロメトリ (1520 世羅町, スギ/ヒノキ)

読み取り (重要発見):

• スギ: n=329,746, Pearson r = 0.829, 線形 h = 6.72 + 0.507 d, べき乗 h = 1.63·d^0.77。
• ヒノキ: n=3,067,440, Pearson r = 0.805, 線形 h = 3.73 + 0.515 d, べき乗 h = 1.08·d^0.85。
• 両樹種とも r ≥ 0.81 と強正相関。 = 樹高と胸高直径は不可分にリンクする (大きい木は太い木)。
• べき乗指数 β: スギ=0.77, ヒノキ=0.85。 β < 1 = 樹高の生長速度より直径の生長速度が早い ("寸胴"傾向)。 β > 1 = 細長傾向。
• スギの β > ヒノキの β: スギは細長く育ちやすい。 これは林学常識 (スギは樹高優先、ヒノキは直径太り) と整合 (またはずれを発見)。
• H3 (r ≥ 0.7) の判定: 支持。
• これは「アロメトリ法則 (= 生物物理学)」がオープンデータの 3.4M 単木で量的に再現できることの実証。 研究的価値: 従来は数百本の毎木調査でしか得られなかった係数を、 航空レーザの単木抽出でデータ単位を 4 桁拡大。 将来的に「自治体ごと・斜面方位ごと・標高ごと」のアロメトリ係数が計算可能。

表 (アロメトリ係数)

図と読み取り (図 6: 単木サンプル位置マップ)

なぜこの図か: 単木 3.4M 点の地理分布を視覚的に把握する。 点を樹高で色分け (黄=低木 → 緑=高木) して、世羅町内の樹高地理パターンを読む。

図 6: 単木サンプル {len(gdf_tp):,} 点 (1520 世羅町, 色=樹高 m)

読み取り:

• 50,000 点の subsample (= 全 3.4M の 1.5%) でも、 世羅町の森林の輪郭が明瞭に描かれる。 = 単木点群は準連続的に分布し、空間的サンプリングが密。
• 樹高色 (緑=高木) は町の周辺部 (山地) に集中、 町中央 (= 高原台地・農地) には森林がない (空白)。
• これは「世羅町は高原台地で、周辺山地に植林帯」という地理学的事実の量的視覚化。
• L36 樹高ラスタ (DCHM) と比較すると、1520 単木は L36 の樹頂点抽出結果であることが視覚的に裏付けられる。 = 派生関係 (L36 → 1520) の量的実証。

分析 5: mesh (1521) ha材積分布と林業生産性 (H4)

狙い (本記事の林業生産性核心)

1521 mesh 750k セル (うち numeric 40,361 セル) で、 ha材積 (m³/ha) の地理分布と樹種別構成を分析する。 「どこに高材積エリア (= 林業価値の高い森林) があるか」を量化する (H4 集塊性)。

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

20m mesh の各セルには林分指標が 12 列付与されている。 geopandas で全件読込し、ha材積 で choropleth マップを描く。 さらにha材積 5 階級 (0-100, 100-250, 250-500, 500-1000, >1000 m³/ha) で量的比較する。

• 入力: GeoPackage 750k mesh, 列 ha材積
• 出力: 樹種別 ha材積 統計, 5 階級構成, 高材積比率
• 限界: mesh は 20m 規則格子のため境界アーティファクトあり (林分の自然境界とは無関係)

実装

↑ L41_forest_resources.py 行 1601–1636

 1
 2
 3
 4
 5
 6
 7
 8
 9
1610
1611
1612
1613
1614
1615

import geopandas as gpd
gdf_mesh_full = gpd.read_file("fr_mesh20m_03NF2.gpkg", layer="fr_mesh20m_03NF2")
mesh_valued = gdf_mesh_full[gdf_mesh_full["樹種"].notna()].copy()
mesh_numeric = mesh_valued.dropna(subset=["ha材積", "立木密度", "平均樹高"]).copy()

# 5 階級分類
VOLUME_CLASSES = [
    ("極低", 0, 100), ("低", 100, 250), ("中", 250, 500),
    ("高", 500, 1000), ("極高", 1000, 5000)
]
def classify_volume(v):
    for i, (_, lo, hi) in enumerate(VOLUME_CLASSES, start=1):
        if lo <= v < hi: return i
    return 0
mesh_numeric["volume_class"] = mesh_numeric["ha材積"].apply(classify_volume)

図と読み取り (図 7: ha材積 choropleth)

なぜこの図か: 林業価値の地理分布を視覚化。色の濃い領域 = 高材積。 庄原市行政界を黒線で囲み、図郭外との関係も読む。

図 7: ha材積 choropleth (1521 03NF2, 庄原市域)

読み取り:

• 03NF2 図郭 (20×15km) 内に40,361 valued meshが分布。 ha材積の地理パターンが明瞭。
• 濃赤の領域 (高材積 ≥ 500 m³/ha) は線状・斑状に集塊。 = 同一林班・同一樹種・同一管理単位を反映。
• 図郭の 4 隅は valued mesh が少ない (= 山頂・川・農地) ことが多く、 森林分布の外縁を反映する。
• 03NF2 図郭は庄原市を中心に含む内陸山地で、針葉樹植林が連続的に広がる。

図と読み取り (図 9: ha材積 5 階級分類)

なぜこの図か: 連続値 ha材積 を慣用 5 階級にカテゴリ化し、 階級マップ + 樹種別構成比 stack-bar の 2 角度で見る。

図 9: ha材積 5 階級 — マップ × 樹種別構成

読み取り (重要発見):

• (a) 階級マップ: 高 (赤) と極高 (暗赤) が領域内に集塊。 中 (橙) と低 (薄赤) が周辺に広がる。
• (b) 樹種別構成: スギは極高比率が 1.1%, ヒノキは 2.7%。 = スギ林分が高材積に偏る傾向 (スギは早生・大径のため)。
• 重要な構造的観察 (本記事独自発見): 1521 mesh で 「その他」 (= 広葉樹混交林) は樹種判読されているが林分指標 (ha材積等) は欠損している (89,696 mesh, 全 valued の 68.3%)。 つまり 1521 の数値指標は管理対象林 (スギ・ヒノキ) のみ計算される設計。 = 広葉樹混交林は林業生産性指標としては未集計 (生物多様性価値は別軸として別途測定要)。 この発見は本データセットの仕様上の重要点であり、利用者は「mesh で広葉樹資源量も同時に得られる」と 誤解しないよう注意が必要。
• 高材積セル (≥500 m³/ha) の比率: スギ=43.9%, ヒノキ=11.7%。
• 8 近傍同種率 (高材積 mask): 0.564。 中程度の集塊。 H4 の判定: 部分支持。
• これは「高材積エリアは点在するのではなく、林班単位の塊で連続分布する」ことの量的実証。 林業実務上は「収穫期到来エリアの一括皆伐計画」の根拠データになる。

表 (1521 樹種別 ha材積 5 階級構成)

表 (1521 樹種別 高材積比率)

図と読み取り (図 12: 樹種別 ha材積 small multiples)

なぜこの図か: 樹種ごとの ha材積地理分布を並列で比較。 スギとヒノキの林分立地の差を視覚的に確認。

図 12: 樹種別 ha材積マップ (1521 small multiples)

読み取り:

• スギ: 高材積 (赤) が斑状に集塊。比較的「川沿い・谷頭」の好条件立地に分布。
• ヒノキ: 中-高材積が広範に分布。「斜面・尾根」など多様な立地に対応。
• その他 (広葉樹等): 林分指標が欠損のためマップ上は空白 (= 樹種は判読されているが ha材積等の数値計算は未実施)。 = 林業生産対象外の「広葉樹混交林・里山」は本データセットでは空間的位置のみ与えられ、量化は別途必要。
• スギ vs ヒノキの立地選好 (スギ=湿潤好・尾根回避, ヒノキ=斜面適応・分布広い) が視覚化される。

分析 6: 林分指標 box plot と健全度 (H5, H6)

狙い

1521 mesh の10+ 列の林分指標を樹種別に箱ひげ図で並べて、 スギ・ヒノキ・その他の形状特性 (高さ・直径)と 健全度 (収量比数・形状比・樹冠長率) を 1 枚で比較する。

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

箱ひげ図 (box plot) は分布の中央値・四分位範囲・外れ値を 1 つの図で示す。 ヒストグラムより圧縮的で、複数群の比較に向く。

• 箱の中線 = 中央値、箱の上下 = 第 1/第 3 四分位 (= 中央 50%), ヒゲ = 1.5 × IQR 内
• 外れ値は本記事では非表示 (showfliers=False) で本体の差に集中
• 限界: ヒストグラムの形状情報は失う

実装

↑ L41_forest_resources.py 行 815–848

 1
 2
 3
 4
 5
 6
 7
 8
 9
824
825
826
827
828

fig, axes = plt.subplots(2, 3, figsize=(15, 9))
metrics = [
    ("ha材積", "ha材積 (m³/ha)", (0, 1500)),
    ("立木密度", "立木密度 (本/ha)", (0, 2500)),
    ("平均樹高", "平均樹高 (m)", (0, 35)),
    ("平均直径", "平均直径 (cm)", (0, 50)),
    ("収量比数", "収量比数", (0, 1.2)),
    ("形状比", "形状比", (30, 110)),
]
for ax, (col, title, ylim) in zip(axes.flat, metrics):
    data = [mesh_numeric[mesh_numeric["樹種"]==sp][col].dropna().values
            for sp in ["スギ","ヒノキ","その他"]]
    ax.boxplot(data, labels=["スギ","ヒノキ","その他"],
               patch_artist=True, showfliers=False)

図と読み取り (図 8: 林分指標 box plot)

なぜこの図か: 6 つの林分指標を樹種別に並列比較し、 「樹種ごとの林分特性プロファイル」を 1 枚で読む。

図 8: 林分指標 6 種 × 樹種 (1521 mesh)

読み取り (重要観察):

• 本図は スギ・ヒノキ のみ表示: 1521 mesh で「その他」(= 広葉樹混交林) は 樹種だけ判読されて林分指標 (材積・密度・樹高...) は欠損している (89,696 mesh, 全 valued の 68.3%)。 = 1521 の数値指標は管理対象林 (スギ・ヒノキ) のみ計算される設計。
• ha材積: スギ > ヒノキ の階層が明確 (中央値 スギ 468, ヒノキ 278 m³/ha)。 スギの方が単位面積あたり林業価値が高い。
• 立木密度: ヒノキが平均的に高い (中央値 ヒノキ 950 本/ha, スギ 725 本/ha)。 ヒノキは小径木が密植される慣行を反映。
• 平均樹高: スギが圧倒的に高い (中央値 22m vs ヒノキ 15m)。 1520 単木と整合 (図 4 の樹高ヒスト)。
• 平均直径: スギ 28cm vs ヒノキ 22cm。 こちらもスギが太い。
• 収量比数: スギ・ヒノキとも中央 0.6 付近で似ているが、 スギの方が分布が右寄り (= 収穫期到来割合が多い)。
• 形状比: スギ中央値 78, ヒノキ 69。 スギはやや細長 (高木で gravity-driven)。 100 超の不健全林分は両樹種とも少数 = 広島県の管理林分は概ね健全。

図と読み取り (図 10: 収量比数 vs 平均樹高 散布)

なぜこの図か: 収穫期判定の閾値 (収量比数 0.7) と樹高の関係を視覚化。 赤線 = 収量比数 0.7, 橙線 = 平均樹高 18m。両線の右上 = 「収穫期到来かつ大樹」エリア。

図 10: 収量比数 × 平均樹高 (1521 mesh)

読み取り (重要発見):

• 収量比数 0.7 以上のセルは合計 10,220 (25.3% of 40,361)。 これらの平均樹高は 19.86m。
• H5 (収量比数≥0.7 で平均樹高≥18m) の判定: 支持。
• 樹種別 収穫期セル比率: スギ=29.2%, ヒノキ=23.2%。
• 収量比数 ≥ 0.7 のセルは右上の散布雲 (高材積・高木) に集塊することが視覚的に確認できる。 = 「収穫期到来は樹高でも検出可能」(代理指標として機能)。
• これは林業実務上「収量比数を持たないデータ (例: L36 樹高ラスタ単独) でも、 平均樹高で代替的に収穫期エリアを推定できる」可能性を示唆する。

表 (1521 樹種別 mesh 集計)

表 (1521 収穫期セル集計)

図と読み取り (図 11: 形状比ヒスト — 林相健全度)

なぜこの図か: 形状比 70-90 が健全、>100 が不健全という慣用閾値を視覚化。 mesh (集計値) と単木 (個別値) を並べて、スケール依存性を見る。

図 11: 形状比分布 — mesh (集計) vs 単木 (個別)

読み取り:

• (a) mesh 形状比: 中央値 72, 70-90 範囲内 55.5%, 100 超 0.07%。 = 大部分の林分は健全範囲に収まる。
• (b) 単木 形状比: 中央値 71, 分布の裾が mesh より広い。 = 個別樹はばらつきが大きい (集計でならされる)。
• H6 (70-90 が健全, >100 が過密で不健全) の判定: 支持。 健全範囲の比率が想定通り
• 過密林分 (形状比 >100) は少数であり、 広島県の管理林分は概ね適正範囲。 ただし不健全セルの地理集塊性は別途要検証 (= 特定地区での放置林リスク)。

分析 7: 3 dataset 統合の制約と方法論的提言 (H7, 本記事の方法論的核心)

狙い (本記事の方法論的核心)

本記事の最大の発見は「3 dataset は同じ航空レーザから派生したが、地理単位が異なるため直接統合できない」 ことの量的確認である。具体的には:

• 1519 樹種ポリゴン: 世羅町 (市町村code 34462)
• 1520 単木ポイント: 世羅町 (市町村code 34462) → 1519 と同自治体
• 1521 森林資源量メッシュ: 03NF2 図郭 (市町村code 34210 = 庄原市) → 1519/1520 と異なる地理

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

本セクションは新規計算ではなく、これまでの分析結果から整合性を量的に確認する。 特に「同じ自治体ペア」と「異なる地理単位」での指標差を整理。

• 1519 + 1520 (同自治体): 樹種別の polygon 件数 (1519) と単木件数 (1520) の比率を確認
• 1521 (別地理): bbox 重なりが世羅町とゼロであることを地理座標で確認
• 研究的含意: 「県全体の森林資源図を作るには 10 市町 × 10 図郭 = 100 ファイル必要」

1519 と 1520 の整合性 (同自治体ペア)

世羅町の 1519 polygon (n=13,025) と 1520 point (n=3,397,186) は同一自治体のため、 直接重ね合わせ可能。1520 はスギ・ヒノキ抽出のため、1519 のスギ・ヒノキ polygon 内にのみ存在する:

読み取り (重要発見):

• 世羅町スギ密度: 1520 単木 ÷ 1519 polygon 面積 ≈ 669本/ha。 林業の慣用密度 (1500-2500本/ha 植栽 → 間伐後 600-1200本/ha) と比較すると、現在の世羅町スギ林は中程度の管理状態。
• 世羅町ヒノキ密度: 1008本/ha。 ヒノキはスギより低密度植栽が一般的 (1500本/ha 程度)。
• これは1519 + 1520 の組合せでしか得られない派生指標。 1521 mesh 単独では (世羅町とは別地理のため) この値は出せない。

1521 と 1519/1520 の地理重なり (なし)

読み取り (本記事の方法論的核心):

• X 範囲は部分的に重なる、 Y 範囲は完全に分離。
• 市町村code が異なる (世羅 34462 vs 庄原 34210) ため、 1519/1520 と 1521 を 1 枚の図で重ねることは物理的に不可能。
• H7 (1520 と 1521 は地理単位が異なるため直接統合不可) の判定: 支持。
• これは本データセットの最大の方法論的制約であり、 利用者は「県全体の森林資源図を作るには 10 市町 × 10 図郭 = ファイル整備の手間が大きい」ことに留意すべき。
• 提言: 「広島県は 1521 を市町村単位でも公開する」か、 「3 dataset を統合した県全域 GeoPackage を 1 ファイルで配布する」のが 利用者にとって最も望ましい (発展課題 Z3)。

仮説検証と考察

本記事で立てた 7 つの仮説 (H1〜H7) と、3 dataset 統合分析の照合結果:

判定サマリ: 支持 6 / 部分支持 1 / 反証 0 (全 7 仮説中)

主な発見 (5 点)

1. 3 dataset の地理単位差が量的に確認 (H7 支持): 1519 (市町村 polygon) + 1520 (市町村 point) は同自治体ペアとして直接統合可能、 1521 (図郭 mesh) は独立した地理単位で 1519/1520 とは bbox が重ならない。 これは「同じ航空レーザ計測から派生したが、設計思想 (主題ベース vs 規則格子ベース) の違いで地理単位が分離」 という GIS 古典問題の例。
2. 世羅町は針葉樹植林 + 広葉樹混交林の中山間林業地 (H1 支持): 1519 で ヒノキ 51% + スギ 20% = 針葉樹植林 71%, その他 (広葉樹等) 24%。 = 「中山間林業地として典型的な構成」が量的に確認できた。
3. スギはヒノキより高木・大径・大材積 (H2 支持): 1520 単木 3.4M で スギ平均樹高 21.0m, ヒノキ 14.1m (差 +6.8m)。 単木材積はスギ ×2.3。 林学常識をオープンデータで桁違いの標本数 (従来の毎木調査 ~10² → 単木抽出 ~10⁶) で再現。
4. 樹高 vs 胸高直径のアロメトリ法則がオープンデータで再現 (H3 支持): スギ r=0.83, β=0.77, ヒノキ r=0.81, β=0.85。 = 樹木のべき乗法則 (allometry) が 3.4M 単木で量的に成立。 研究的意義: 自治体別・斜面方位別・標高別のアロメトリ係数算出への展開可能性。
5. 高材積エリアは集塊し、収穫期セルは平均樹高で代替推定可能 (H4 部分支持, H5 支持): 1521 mesh で ha材積 ≥ 500 m³/ha のセルは 8 近傍同種率 0.564。 = 林班単位の塊で連続分布。 収量比数 ≥ 0.7 のセル平均樹高 19.9m。 = 「収穫期は樹高でも代理判定可能」という方法論的応用可能性。

世羅町・庄原市域の森林資源量サマリ (本記事の量的サマリ)

考察

本記事の主たる発見は、「広島県の航空レーザ森林資源 3 dataset は同じ計測 (2018-2019) から派生したが、 地理編成 (市町村別 vs 図郭別) と属性深度 (基本 vs 詳細) が大きく異なるため、 利用者は dataset ごとに異なるアプローチが必要」であることの量的実証である。

1519 樹種ポリゴンは森林の輪郭マップとして最も基本的、属性は素朴 (樹種・面積) だが 全市町村で公開済みなので県全域の樹種分布把握に最適。 1520 単木ポイントは桁違いのデータ量 (~10⁶ point/市町村) でアロメトリ等の生物物理学的研究に向く。 1521 mesh は10+ 列の林分指標で林業実務 (収穫期判定・健全度評価) に直結する。

樹種構成 (H1) は世羅町ではヒノキ優占 (51%), スギ 20%, その他 24%。 針葉樹植林 (ヒノキ+スギ) が約半分を占める中山間林業地として典型的。 これは森林経営計画の基礎情報であり、林野庁 GIS でも国土数値情報「森林地域」と整合する。

単木スケールでのスギ vs ヒノキ (H2) は、3.4M 点という統計的検定力の極めて高い標本で、 スギの優位性 (高木・大径・大材積) を明瞭に示した。 特にスギ単木材積はヒノキの×2.3であり、 林業経済的にスギ林分の方が単位面積あたり価値が高い (= 高林齢で高材積) ことの量的根拠。

アロメトリ (H3) は樹木の生物物理学法則を 3.4M 単木で実証した重要な発見である。 スギ h=1.63·d^0.77, ヒノキ h=1.08·d^0.85 というべき乗関係は、従来は数百本の毎木調査でしか得られなかった係数を 4 桁拡大したサンプルで安定推定。 研究的意義として、地域別・斜面方位別アロメトリ係数の算出が現実的になった。

1521 mesh の林分指標 (H4-H6) は林業実務への直接的応用を示す。 高材積エリアの集塊性 (8 近傍 0.564) は林班単位の連続性を反映し、 収量比数 ≥ 0.7 セルの平均樹高 19.9m は収穫期の樹高代理可能性を示唆。 形状比 70-90 健全範囲が 56% と多数派であり、 広島県の管理林分の概ね健全性が量的に確認できた。

本記事最大の方法論的貢献は3 dataset の地理単位差 (H7) の量的確認である。 1519/1520 (市町村) と 1521 (図郭) は構造的に分離されており、 県全域の森林資源 1 枚画像を作るには 10 市町 × 10 図郭 = 大量のファイル統合作業が必要。 オープンデータ提供者 (広島県林業課) への提言として、「市町村単位 1521 mesh の追加公開」または 「県全域統合 GeoPackage の単一ファイル配布」が利用者の研究効率を大幅に向上させる (発展課題 Z3 で扱う)。

発展課題

結果 X1 → 仮説 Y1 → 課題 Z1: L36 樹高ラスタ (DCHM) と 1520 単木ポイントの直接対応関係検証

結果 X1: 1520 単木ポイントは「樹種ポリゴン (1519) でスギ・ヒノキと判読された区域内において DCHM データ等を利用して樹頂点を抽出した点」と公式説明されている (DoBoX dataset 1520 説明文)。 本記事 図 6 で世羅町内に単木点が密に分布することを確認したが、 L36 樹高ラスタ (DCHM) との 1:1 対応は未検証。

新仮説 Y1: 1520 単木ポイントの座標 (x, y) における L36 DCHM ラスタのピクセル値は、 1520 の樹高 (m) 列とほぼ一致するはず (= L36 ラスタ値が 1520 単木抽出の元データ)。 ただし overview 解像度差や時期ずれ (DCHM 計測時期 vs 1520 計測 2019-01-19) で±1m 程度の誤差あり。

課題 Z1: L36 樹高ラスタ (世羅町 1m, 446M セル) を rasterio で開き、 1520 の各単木座標で rasterio.sample によりピクセル値を取得。 単木の樹高列 vs DCHM 値の Pearson 相関を計算し、r ≥ 0.95 で 1:1 対応が量的に裏付けられる。 さらに残差 (DCHM - 1520樹高) のヒストグラムで抽出アルゴリズムのバイアス (= 樹頂点検出の系統誤差) を同定可能。 これは LiDAR 林業応用層と LiDAR ラスタ層の派生関係の量的実証であり、 研究的意義は L36 と本記事の記事間整合性の裏付け。

結果 X2 → 仮説 Y2 → 課題 Z2: 1521 mesh の収量比数を地形要因で予測

結果 X2: 1521 mesh の収量比数は 0-1.08 の範囲で分布、≥0.7 セルが 25% を占める。 ただし収量比数の地理パターンは本記事では未分析。

新仮説 Y2: 収量比数は地形要因 (標高・傾斜・斜面方位) で予測できる。 すなわち 「南向き・緩斜面・標高 200-500m」のセルは収量比数 ≥ 0.7 の確率が高い。 これは林業経営学の常識 (= スギ・ヒノキの好立地条件) を量的に検証することに相当。

課題 Z2: 1521 mesh と L40 標高ラスタ・L39 傾斜ラスタをセル単位で空間結合。 ただし mesh は 20m, ラスタは 1m なので解像度ダウンサンプルが必要。 標高 + 傾斜 + 斜面方位 (= 北/南/東/西の 4 区分) を特徴量として、収量比数を目的変数に Random Forest 回帰。 重要度ランキングで「林業生産性に最も効く地形要因」を同定。 さらに 03NF2 だけでなく 10 図郭で繰り返して、地域横断的な収量比数モデルを構築。

結果 X3 → 仮説 Y3 → 課題 Z3: 県全域の森林資源 1 ファイル統合 (オープンデータ提供者への提言)

結果 X3: 1519/1520 (市町村単位 10 ファイル) と 1521 (図郭単位 10 ファイル) は地理編成が異なり、 県全域図を作るには 20+ ファイルの個別 DL → 結合処理が必要。 これは利用者の研究効率を大きく阻害する。

新仮説 Y3: 県全域統合 GeoPackage (= 全市町村 + 全図郭をマージしたシングルファイル) を提供すれば、 利用者の前処理コストが劇的に下がり、研究の生産性が向上する。 GeoPackage は 1 ファイル内に複数レイヤ (tree_species_pref, tree_point_pref, fr_mesh_pref) を保持できるので、技術的には可能。

課題 Z3: 本記事のスクリプトを 10 市町 × 3 dataset で繰り返して実行し、 gpd.GeoDataFrame.to_file で県全域マージファイルを生成。 ファイルサイズの試算 (10 市町 × 552 MB ≈ 5 GB の単木) で実用上の制約を確認し、 「Cloud Optimized GeoPackage (COG-style chunked) 配布」または 「市町村別 ZIP のままで提供 + 統合スクリプトを公開」のどちらが現実的か検討。 広島県林業課・DoBoX 運営への政策提言として論文化可能。

結果 X4 → 仮説 Y4 → 課題 Z4: アロメトリ係数の自治体横断比較

結果 X4: 世羅町でスギ β=0.77, ヒノキ β=0.85。 これらは世羅町固有の値であり、他自治体での再現性は未検証。

新仮説 Y4: アロメトリ係数 β は気候帯・標高帯・管理状態で変動する。 広島県の 10 自治体で計算すると、沿岸 (温暖) vs 内陸 (寒冷) で β が ±0.05 程度異なる。 特に内陸高標高 (= 庄原・北広島町) では樹高生長が抑制され、β が小さくなる。

課題 Z4: 1520 単木データを 10 自治体すべてで取得 (= 計 ~30M 単木) し、 自治体ごとに β を計算。気温平年値 (気象庁データ) ・標高 (L40 平均) との回帰で 「地域要因 → β」のモデルを構築。アロメトリ係数の地域分布図を作り、 広島県の森林管理ゾーニングへの活用 (低 β 地域は早期間伐推奨等) を提案。

結果 X5 → 仮説 Y5 → 課題 Z5: ha材積の時系列トレンド (再計測時の変化検出)

結果 X5: 本記事の計測年は 2019-01-19 単一時点。森林の材積は 1 年で 5-10 m³/ha 増加するため、 再計測 (例: 2024 年) で5 年間の材積増加量が観測できる。

新仮説 Y5: 2019 → 2024 の ha材積差分は、 10-50 m³/ha の正の増加 (= 通常成長) を示すセルが大半。 ただし皆伐・間伐・台風被害のセルは負の差分 (- 100 m³/ha+) を示し、空間的に集塊する。

課題 Z5: 広島県が 2024-2025 に再計測を予定 (DoBoX dataset 1527, 1434 等参照)。 新計測の 1521 mesh と 2019 版を空間結合し、ha材積差分マップを作成。 正の差分 (成長) vs 負の差分 (伐採・被害) の地理パターンを分析。 これは森林炭素吸収量 (CO2) 推定の基礎データになり、 脱炭素政策・J-クレジット制度への直接応用可能。

結果 X6 → 仮説 Y6 → 課題 Z6: 単木 形状比の不健全林分マッピング

結果 X6: 単木 1520 で形状比 >100 (= 過密で不健全) のセルは 6.1%。 これらの地理集塊性は未分析。

新仮説 Y6: 不健全林分 (形状比 >100) は、「植栽後 30-50 年で間伐が遅れた」放置林に集塊する。 すなわち管理単位 (林班) 単位で連続して不健全。

課題 Z6: 1520 単木の形状比 >100 を mask とし、 本記事と同じ 8 近傍同種率を計算。集塊性 ≥ 0.5 ならば「林班単位の放置林」が量的に同定。 さらに森林経営計画 GIS (国土数値情報「森林地域」 + 広島県林業 GIS) と空間結合し、 不健全林分の所有者・管理者を特定 (政策提言: 間伐補助金重点投下対象の地理ターゲティング)。

結果 X7 → 仮説 Y7 → 課題 Z7: 樹冠面積データの欠損分析

結果 X7: 単木 1520 の 樹冠面積 列は世羅町データで全件 NaN。 他自治体や mesh 1521 の面積_ha 列とは別概念で、未充填の理由は未明。

新仮説 Y7: 樹冠面積は計測初期段階では未集計で、後の自治体 (= 2024 以降公開分) では充填されている可能性が高い。 あるいは樹頂点抽出アルゴリズムの仕様上、樹冠半径推定が信頼できないため未公開。

課題 Z7: 10 自治体の 1520 で樹冠面積列の充填率を確認。 未充填なら「DCHM ラスタから個別樹冠を抽出する Watershed セグメンテーション」を独自実装。 これは林業計測学の標準手法 (Popescu 2007 等) で、L36 樹高ラスタ + L37 点群密度 から派生可能。 新計算した樹冠面積 vs 樹冠長率 vs 単木材積で樹冠形態と材積の関係を量的に検証。