More than 5 years have passed since last update.

RのForestToolsライブラリにより樹頂点を検出しQGIS上で表示する

Posted at 2019-11-26

経緯

・近年、森林・林業分野において、航空機レーザ計測データなどによる地表面のデジタルデータ化が進んでいる。
・公共測量の場合は、測量法に則った申請により計測データの取得が可能であると思われるが、データを利用価値のある形に変換する方法が必要。
・統計分析フリーソフトRにForestToolsという良さげなライブラリがあると知り、QGISに持ち込むところまでをお試ししてみた。

##今回やりたいこと
樹冠高ラスタデータ（CHM）→樹頂点の検出→QGISで表示

ForestToolsライブラリについて

・リモートセンシング技術により収集された森林情報を分析するForestToolsパッケージ（開発：Andrew Plowright氏）
・樹冠高ラスタデータ（CHM：Crown Height Model）から、主要木の樹頂点位置や個々の樹冠の輪郭の推定が可能
・指定したポリゴン内の樹高や本数情報の要約等も可能とのこと
・参考：The Comprehensive R Archive Network（"Canopy analysis in R using Forest Tools"）
　https://cran.rstudio.com/web/packages/ForestTools/vignettes/treetopAnalysis.html

手順

###１　Rを起動し、次のコマンドでライブラリをインストール
install.packages("ForestTools")
（実行するとミラーサイトの場所を聞かれるので、最も近い場所を選ぶ）

###２　ライブラリをアタッチ
library(ForestTools)
library(raster)

###３　サンプルをロード
data("kootenayCHM")

###４　樹頂点抽出処理（vwf）
####４－１　アルゴリズム概要
・vwfとは？＝variable window filterアルゴリズム（Popescu and Wynne 2004)
・CHM上を特定の大きさのwindowでスキャンし、突出して高い箇所を抽出する。
・windowの大きさは、中心となるセル値（樹冠の高さ）に依存する関数とすることで、樹木の高さよる樹冠投影面積の大きさを補完している。

####４－２　解析スクリプト
lin <- function(x){x * 0.05 + 0.6} ※後述

ttops <- vwf(CHM = kootenayCHM, winFun = lin, minHeight = 2)
引数：CHM = ラスタデータ（CHM），winFun = window sizeの関数，minHeight = 樹頂点を検出する最低値

※winFunについて
・windowの中心となるセル値を取得して、search windowの半径を返す。
・single input & single outputの式で定義（上記は最もシンプルな線形方程式）

####４－３　データの取り出し方
ttops（データフレーム）
class : SpatialPointsDataFrame
features : 1077
extent : 439689.2, 439832.2, 5526454, 5526562 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 4
names : layer, height, winRadius, treeID
min values : 1, 2.0020417869091, 0.700102089345455, 1
max values : 1, 13.4912073612213, 1.27456036806107, 1077

names(ttops)
[1] "layer" "height" "winRadius" "treeID"

ttops$x[1:3]
[1] 439689.8 439731.2 439733.2

ttops$y[1:3]
[1] 5526562 5526562 5526562

ttops$height[1:3]
[1] 3.265903 3.027651 3.463092

####４－４　解析結果を図示
プロットの枠
par(mar = rep(0.5, 4))

ラスターデータ図示
plot(kootenayCHM, xlab = "", ylab = "", xaxt='n', yaxt = 'n')

樹頂点のプロット
plot(ttops, col = "blue", pch = 20, cex = 0.5, add = TRUE)