Rによる樹頂点解析の試行－樹高・樹冠面積の算出やQGISへの表示
表題の備忘録。
G空間情報センター ｜ 林野庁が高精度な森林資源情報等を公開しました (geospatial.jp)
栃木県、高知県、兵庫県のレーザ計測成果が公開された。
ほかの県も進むといいな。
公開データの中にDCHMがあったので、樹頂点検出（抽出）に挑戦してみた。
※DCHM（数値樹冠高モデル）
≒DSM（数値表層モデル）-DTM（数値地形モデル）
準備
DCHMラスタの準備
DCHMラスタはG空間情報センター上で公開されている。
今回はここから栃木県のものを使わせていただくこととした。
図郭割りされているが、それでも解析に直接使うには範囲が広すぎるので、QGISに放り込んで適当に切り抜いた
「DCHM2.tif」という名前で保存した。
座標系は平面直角座標系だった。栃木県なので第9系。
（JGD2011 / Japan Plane Rectangular CS IX（EPSG:6677））
（背景地図はOpenStreetMap）
Rのインストール
略！
version 4.3.2 The Comprehensive R Archive Network (r-project.org)
Rを起動したら
「ForestTools」パッケージをインストール。
Rを起動して次のとおり入力・確定（シャープはコメント）
#ForestToolsのインストール
install.packages("ForestTools")ミラーリングサイトはとりあえずJapan Tokyoを選択した。
ForestToolsのバージョンは1.0.1だった。
作業ディレクトリはDCHMラスタを保存している場所に設定しておく
次にライブラリ読み込み
#ライブラリ「ForestTools」等を読み込む
library(ForestTools)
library(terra)
library(sf)terraは空間データ分析用パッケージ。以前はrasterパッケージだったらしい
sfはsimple featuresの略。R上で地物と属性データを扱うためのものらしい
DCHMラスタを読み込む
#DigitalCanopyHeightModel：DCHMラスタを読み込む
chm<- terra::rast("DCHM2.tif")「chm」という任意の名前を付けた
ちょっと脱線（飛ばしてOK）
ForestToolsのチュートリアルに沿いつつ、気になったことは都度調べて進めているので、樹頂点検出からやや脱線することも大いにアリ。
RでDCHMを図示してみる
#plot関数を使用してDCHMを図示する。セルの値 は、樹冠の地表からの高さと等しくなる。
#↓プロットエリアの余白設定
par(mar = rep(0.5, 4))
#図示
plot(chm, xlab = "", ylab = "", xaxt='n', yaxt = 'n')
今回のDCHMの最高値（最大樹高）は35mくらいだなーということがわかる
図示した結果をPNG形式で保存してみる
#png関数で描画デバイスを起動させる。命名したいファイル名と画像サイズを指定
png("Test.png",width=500,height=500,units="px")
#作図命令(上で作図したのと同じもの）
plot(chm, xlab = "", ylab = "", xaxt='n', yaxt = 'n')
#描画デバイスの終了命令。これで、作業ディレクトリ内に図が保存される
dev.off()
樹頂点の検出
ForestToolsでは、樹頂点を検出する手法のひとつ「Variable Window Filter （VWF）アルゴリズム」（Popescu&Wynne (2004)）を使用することができる（https://github.com/andrew-plowright/ForestTools から意訳）。
このアルゴリズムはDCHMを構成する全画素（セル）を”検索窓”でスキャニングし、そのセルの範囲内で最も高いセルを樹頂点とみなすもの。”検索窓”の大きさは、検索窓の中心セルの高さ情報によって決まる。
（https://cran.r-project.org/web/packages/ForestTools/vignettes/treetop_analysis.htmlとhttps://doi.org/10.1111/2041-210X.13860から何とか把握）
今回は、ほかの研究より日本の環境に近いと考えられるある程度複雑な構造を持つ針葉樹混交林で研究を行ったYoungら（2022）により示された次の関数を採用した。y＝0.04x
（あるセルの高さ情報がxのとき、そのセルを中心とした”検索窓”の半径はyとする）
【注意】Youngらの研究はUAV写真測量により得られたCHM等の解析であるため、今回のような航空レーザ計測成果とは異なる。
また、VWFのような間接的な検出手法には偽検出や未検出はあって当然と思うこと。いくつもの手法が検討・検証されてきた中で、より現地に近い値が得られるものとして本手法が見出されてきたと心得ること。
さて、本題。
樹頂点検出処理を行う
# まずは「検索窓」のサイズを決定する線形関数を定義
lin <- function(x){x * 0.04}
# vwf関数を用いて樹頂点（treetops）を検出
ttops <- vwf(chm, winFun = lin, minHeight = 2)vmf関数内で、minHeightに指定した値未満の箇所は樹頂点とはみなされない。
これを活用して樹頂点の偽検出を最小限にすること（所謂足切り）が可能。
今回は仮で2mとした。
検出した樹頂点を図示する
#CHMの図示（既出）
plot(chm, xlab = "", ylab = "", xaxt='n', yaxt = 'n')
# 樹頂点を書き加える
plot(ttops$geometry, col = "blue", pch = 20, cex = 0.5, add = TRUE)
けっこう偽検出あるかもしれないけど、面白い
樹高の平均、最大、最小を調べる
# 平均
mean(ttops$height)
# 最大
max(ttops$height)
# 最小
min(ttops$height)カンタン！
樹頂点解析結果をSHP形式で保存する
#上で解析した結果ttopsをベクタ形式の「s」と命名
s<-vect(ttops)
#出力時のファイル名を決める。今回は「ttop_point.shp」とした
outfile <- "ttop_point.shp"
#ベクタ形式のファイルとして作業中のディレクトリに出力。同じ名称のファイルがある場合は上書き保存
writeVector(s, outfile, overwrite=TRUE)
（右隣のwinRadiusは次回アップデートで削除予定らしい）
樹頂点のデータをcsv形式で保存する場合
（SHP形式で保存していればこれは不要かもしれない）
#作業中のディレクトリにcsv形式で樹頂点解析結果「ttops」を保存
write.csv(ttops,"test.csv")
CSVファイルから樹頂点をQGISで表示する
その前に、先ほど生成したCSVファイルを開いてみたところ、位置座標に相当するはずのセルに不要な文字が入ってしまっていることがわかった。
これを除去しないとQGISにそのまま反映させることができない。
今回はEXCELのSUBSTITUTE関数を使って手作業で加工することにした。
念のためQGISのDCHMラスタの座標を大まかに確認してみる
さて、CSVをQGISに放り込んでみる。
これまで何度も航空レーザ計測成果のテキスト形式のDEMを放り込んでいるので慣れたもの。
樹頂点検出結果の確認
ポイントデータとしてQGISに表示し、背景図と見比べてみた
本来であれば、一定区域の全林毎木調査を行い、単木レベルで解析結果と突合して偽検出・未検出・合致立木の樹高比較等を行うのが良いだろうが、今回はあくまでお試しということで簡単に航空写真と見比べてみる。
※背景図は国土地理院シームレス航空写真
DCHMデータの縁部分。
縁はこういった偽検出が出やすいのかも？
偽検出はあるが、思っていたほど悪くない結果に思える。
樹冠図の作成（ラスタ）
樹頂点検出後、各立木の樹冠範囲を推定・作図することができる。これにはmarker-controlled watershed segmentation（mcws）関数を使う。
どうやら、もともと画像解析による流域界の検出（水文解析）アルゴリズムがあって、それを生体組織や腫瘍など、形態学的な方面に応用したものを、更に樹冠の範囲推定に応用している…っぽい（多分）
# 樹冠の検出（ラスタ）
crowns_ras <- mcws(treetops = ttops, CHM = chm, minHeight = 1.5)
# 樹冠の図示（ラスタ）
plot(crowns_ras, col = sample(rainbow(50), nrow(unique(chm)), replace = TRUE), legend = FALSE, xlab = "", ylab = "", xaxt='n', yaxt = 'n')
DCHM範囲を、検出された樹頂点1本1本に対応する樹冠に分割し色分け
樹冠図の作成（ベクタ）
# 樹冠の検出（ポリゴン）
crowns_poly <- mcws(treetops = ttops, CHM = chm, format = "polygons", minHeight = 1.5)
# CHMの図示（再掲）
plot(chm, xlab = "", ylab = "", xaxt='n', yaxt = 'n')
# 樹冠輪郭を図に加える
plot(crowns_poly$geometry, border = "blue", lwd = 0.5, add = TRUE)
今度はこれをSHP形式で出力する
樹冠図をSHP形式で出力する
#樹冠ポリゴンをSHP形式で保存する。
#上で作成したcrowns_polyをベクタ形式の「s」と命名
s<-vect(crowns_poly)
#出力時のファイル名を決める。今回は「shp_test.shp」とした
outfile <- "shp_test.shp"
#ベクタ形式のファイルとして作業中のディレクトリに出力。同じ名称のファイルがある場合は上書き保存
writeVector(s, outfile, overwrite=TRUE)
QGISに表示してみよう
樹冠図（SHP形式）をQGISに表示してみる
CRSは大元のDCHMラスタと同じ（平面直角座標系）だった。
樹冠面積・樹冠直径の算出
樹冠の面積や直径を計算することができる（直径は各樹冠がほぼ円形をしていると仮定した場合、面積から機械的に算出するもの）
# 樹冠面積の算出（樹冠ポリゴンに面積データを付加）
crowns_poly[["area"]] <- st_area(crowns_poly)
#樹冠直径の算出（面積から）（樹冠ポリゴンに直径データを付加）
crowns_poly[["diameter"]] <- sqrt(crowns_poly[["area"]]/ pi) * 2
#樹冠の平均面積
mean(crowns_poly$area)
# 樹冠の平均直径
mean(crowns_poly$diameter)樹冠ポリゴンに樹高データを結合し、SHP形式で出力
#樹冠ポリゴンに"ttops"の高さデータを追加
crowns_poly[["ttopsHT"]] <- ttops$height
#-------------------------
#ポリゴンを再度SHP形式で保存する。
#crowns_polyを再度ベクタ形式の「s」として認識させる
s<-vect(crowns_poly)
#出力時のファイル名を決める。今回は「shp_test2.shp」とした
outfile <- "shp_test2.shp"
#ベクタ形式のファイルとして作業中のディレクトリに出力。同じ名称のファイルがある場合は上書き保存
writeVector(s, outfile, overwrite=TRUE)
その属性データ（右）
各樹冠面積、直径、樹高の属性が追加されている。
さいごに
樹頂点検出の一連の流れを試すことができた。
一つ一つ読み解いたり、試行錯誤するのはとても楽しい。
これを「使える」知識・技術とできるよう、努力しよう。
●本記事で以上の内容を試すにあたり、以下のページを参考にさせていただきました
ForestToolsパッケージ作者様によるチュートリアル
- Canopy Analysis in R Using ForestTools (r-project.org)
日本語で情報を残してくださった先人の方々
- RのForestToolsライブラリにより樹頂点を検出しQGIS上で表示する #R - Qiita
- Rで樹頂点を抽出 (coocan.jp)
●本敲では、栃木県森林資源データ「DCHM(数値樹冠高モデル)0.5m」（https://www.geospatial.jp/ckan/dataset/dchm_tochigi）を使用しています。