Rで航空機LiDARデータ解析

はじめに

地表面の詳細な構造を把握するために、航空機によるLiDAR観測が各所で実施されており、オンライン上でも少しずつではあるがフリーで利用が可能なデータになってきた。
QGISなど、LiDARデータを解析するツールはいくつか存在するが、本記事では日本語の記事が少なく、機能が充実していそうなR言語のLidRを用いた解析についてまとめた。
本記事では処理の一工程に沿って説明しており、そこそこボリュームがあるため、目次の項目などから適宜解析に必要な項目を探して活用していただきたい。

準備・環境

環境はWindows10のWSL2（Ubuntu 20.04）。
Rのバージョンは4.0。Rは初めたばかりで試行錯誤しながら手を動かしているところ。
パッケージはいろいろと用いているが、メインはLiDARデータ処理用パッケージLidR（Github、マニュアルはこちら）。
つい先日（2020年7月）、パッケージ製作者がチュートリアルブックを作られたようなので、こちらもぜひ。

その他、本記事で必要なパッケージはすべてinstall.packages()から取得が可能。

使用データ

今回使用するデータは、静岡県ポイントクラウドデータベースよりダウンロードしたもの。静岡県内各所のLiDARデータをいくつか漁ったうち、森林域で点群密度の高いデータを選択。

とりあえず開く

点群データはlasフォーマットにて配布されている。ArcGISの説明①、②などが参考になる。

点群データのヘッダー読込

多くの場合、lasデータは大量の点群データを格納しており、数100Mbから数Gbになることもある。そのため、初めはデータ全体ではなくヘッダー情報のみを短時間で取得するのがよい。

library(lidR)
head <- readLASheader("30XXX00010005-1.las")
print(head)

出力の一部を下記に示す。一つのファイル内に300万点近くのデータを持っているということなどがわかる。
XYの座標系はJGD2000もしくはJGD2011。¹

File signature:           LASF 
File source ID:           0 
:   :   :  :
Point data record length: 34 
Num. of point records:    3055555 
Num. of points by return: 3055555 0 0 0 0 
Scale factor X Y Z:       0.001 0.001 0.001 
Offset X Y Z:             0 0 0 
min X Y Z:                -57141.35 -87000 1094.17 
max X Y Z:                -56800.01 -86683.21 1254.27 
Variable length records:  void

Google Map上で変換した緯度経度を示した図が下の通り。領域全体をほぼ同じ大きさの樹木が覆いつくしていることがわかる。

点群データの読込

点群データ本体を読み込む。lasを圧縮したフォーマットであるlazフォーマットでも全く同様に読み込むことができる。読み込む変数を制限するオプション等もある（詳しくは公式チュートリアルを参照）。

library(lidR)
las <- readLAS("30XXX00010005-1.las")

# class        : LAS (v1.2 format 3)
# memory       : 268.1 Mb 
# extent       : -57141.35, -56800.01, -87000, -86683.21 (xmin, xmax, ymin, ymax)
# coord. ref.  : NA 
# area         : 58400 units²
# points       : 3.06 million points
# density      : 52.32 points/units²

print()にて簡易的な情報が読み込める。
データを参照するには、las@dataのようにする。そのうち、特定の一項目について抽出したい場合は、las@data$項目名とするか、のちに紹介するfilter_poiを用いる。
las@dataとしたときの出力は下記の通り。今回用いたデータでは、Intensity（強度値）やReturnNumber（跳ね返り番号）には一律な定数が用いられている。また、Classification（地物分類）はすでに値が入っている。

X	Y	Z	gpstime	Intensity	ReturnNumber	NumberOfReturns	ScanDirectionFlag	EdgeOfFlightline	Classification	Synthetic_flag	Keypoint_flag	Withheld_flag	ScanAngleRank	UserData	PointSourceID	R	G	B
<dbl>	<dbl>	<dbl>	<dbl>	<int>	<int>	<int>	<int>	<int>	<int>	<lgl>	<lgl>	<lgl>	<int>	<int>	<int>	<int>	<int>	<int>
-56974.10	-86999.93	1147.04	-1	0	1	1	0	0	3	FALSE	FALSE	FALSE	0	1	1	38144	37120	32768
-56974.34	-86999.85	1147.30	-1	0	1	1	0	0	6	FALSE	FALSE	FALSE	0	1	1	44288	43008	39168
-56972.85	-86999.98	1146.44	-1	0	1	1	0	0	5	FALSE	FALSE	FALSE	0	1	1	44032	40448	37120
-56973.10	-86999.89	1146.69	-1	0	1	1	0	0	4	FALSE	FALSE	FALSE	0	1	1	43264	40704	36864

点群データの3次元描画

lidRにて、rglという3次元描画パッケージをベースにしたplot()が用意されている。²

とりあえず、デフォルト設定を用いて描画。高さ情報（Z）を用いてカラーリングがなされる。

plot(las)

# 画像の保存
library(rgl)
rgl.snapshot("ファイル名.png")

マウスを用いてぐるぐる図形を回転させることができる。どうやら矩形領域の半分近くはデータが存在しないらしい。横方向から見た画像だと、樹木で覆いつくされた層がきちんと存在していることがわかる。

今度は、lasのClassificationの項目値で色を分けて描画してみる。

plot(las, color=Classification)

下層から上層にかけてなんとなく色が異なっていることがわかる。

また、今回使用したデータにはRGB値のカラムが設定されており、それを用いた描画も可能。

plot(las, color ="RGB")

要素を用いたフィルタリング

Classificationのうち1つに絞って描画をするために、filter_poi()を用いる。
今回は、分類2（地表、参照）にフィルタリングする。

plot(filter_poi(las, Classification == 2))

最初に描画した図と比較すると点の量は減少しており、フィルタリングは行われていることがわかる。しかし、初めに入っていたClassificationの値では、きちんと地面を表せてはいないことがわかる。

点群データの領域フィルタ

もしもデータが大きすぎて処理に時間を要する場合、一度領域を切って保存してから作業を行うのがよい。領域の切り取りはclip_...()を用いて下記のように行う。

sub <- clip_rectangle(las, xmin, ymin, xmax, ymax)

点群データの保存

データ数を減らしたり、新たなカラムを追加したlasなどを保存する場合、writeLAS()を用いる。

writeLAS(sub, "ファイル名.las")

地面ポイントの抽出

地表を表すClassificationを再計算して求めることにする。
使用するのはclassify_ground()という関数。中で、pmf() (Progressive Morphological Filter)というアルゴリズムを用いているが、詳細は公式チュートリアルを参照のこと。

ws <- seq(3, 12, 3)
th <- seq(0.1, 1.5, length.out = length(ws))
las <- classify_ground(las, pmf(ws, th))

再度Classification==2の条件で描画をする。

うまく樹冠部分の点群が除かれていることがわかる。

地表マップ（標高マップ、DTM）の作成

DTM (Digital Terrain Model)は、樹冠を除いた地表の高さを表す2次元マップ。こちらも関数内でknnidw()（k近傍法+逆距離加重法）というアルゴリズムを用いている（その他の関数含め、詳細は公式チュートリアルにて）。

dtm <- grid_terrain(las, res=1, algorithm=knnidw())

ラスタ（グリッド）データの2次元描画

DTMはRasterLayerというデータ型に保持された2次元ラスタデータで、library(raster)のplot()で描画する（参考）。zlim、zlimcol（カラーバーの外れ値をNullにするフラグ）は.rasterImagePlotのオプションで、plot()を用いると内部的にこの関数のオプションが使用できるようになるらしい。³

png("ファイル名.png", width = 720, height = 720)
par(plt = c(0.1, 0.9, 0.1, 0.9), cex.axis=1.5, cex.main=1.5)
plot(dtm, main='DTM [m]', zlim = c(1090, 1250), zlimcol=FALSE)
dev.off()

内挿補完によって、点群が存在しない対角線周辺部分についても標高が計算されているのがわかる。

標高マップの3次元描画

白地図のような標高マップの描画は簡単にできる。

plot_dtm3d(dtm, bg="white")
# 画像の保存
rgl.snapshot("ファイル名.png")

適切に色を付ける方法がわからなかったため、今度はrglパッケージのrgl.surface()を用いて、先ほどの2次元マップを3次元表示する。¹
DTMはRasterLayerというデータ型で保持されており、as.array(dtm)で配列化して描画する。

x <- rev(1:dtm@nrows)
z <- 1:dtm@ncols
ylim <- range(as.array(dtm), na.rm=TRUE)
ylen <- ylim[2] - ylim[1] + 1
colorlut <- terrain.colors(ylen)
col <- colorlut[as.array(dtm) - ylim[1] + 1]
rgl.surface(x, z, as.array(dtm), color = col, back = "lines")

2次元の時と比べ、見た目の印象もなかなか異なるのがわかる。

2次元マップ（RasterLayer）の保存と読込

保存と読込はそれぞれ、library(raster)のwriteRaster()とraster()を用いる。保存形式はGeoTiffとしているが、ほかにもいくつかの形式がサポートされている。

writeRaster(dtm, 'ファイル名.tif', format = "GTiff", overwrite = TRUE)

# 読込
dtm <- raster("ファイル名.tif")

標高差をなくす

点群からDTMを減算する形で、地表面標高を均一にする。dtmがない場合でも、いくつかのアルゴリズムから均一化は行える（詳しくは公式チュートリアルにて）。

las <- normalize_height(las, dtm)

地面が平らになった。東側の木々の樹高が高いほか、西側にもぽつりと標高の高い木が認められる。

樹高マップ（CHM作成）

先ほどの点群を2次元マップ化してCHM (Canopu Height Model)を作成する。下記では、点群からグリッドがするためのアルゴリズムにdsmtin()（ドロネー三角形分割）を用いている。その他のアルゴリズムなどの詳細は公式チュートリアルにて。

chm <- grid_canopy(las, res = 1, algorithm=dsmtin())

樹頂（樹数）推定

find_trees()を用いて樹頂推定を行う。lmf() (Local Maximum Filter)はlidRの組み込み関数であり、与えるウィンドウサイズの定義によって樹数が大きく変化する（詳細は公式チュートリアルで）。

ft <- find_trees(las, lmf(5))

# 描画
x <- plot(las, color="RGB")
add_treetops3d(x, ft)

lasを狭領域にclipし、左から5m、4m、3mのウィンドウサイズで樹頂推定を行ったもの。本数としては4mのものが適切に見えるが、評価はなかなか難しい。

ベクター（座標・ポリゴン）データの保存と読込

find_trees()で得られるデータは、XY座標に樹頂高さが付与されたデータ（SpatialPointsDataFrame）であり、writeOGR()を用いてshapefileに保存できる。
この後出てくるdelineate_crowns()から得られるポリゴンデータについても同様に保存できる。

library(rgdal)
writeOGR(obj=ft, dsn="保存ディレクトリ", layer="ファイル名(.shpなし)", driver="ESRI Shapefile")

# 読込
ft <- readOGR("ファイル名.shp")

樹冠推定

点群データから単木抽出を行う。segment_trees()を実行すると点群データにtreeIDというカラムが追加され、IDごとに単木が分かれている。こちらも4種類の関数が用意されており、それらの内部でもいくつかの定数を決めてやる必要があるため、本格的に用いる場合には試行錯誤を要する（関数の説明は公式チュートリアルを参照）。

下記は、いずれもデフォルトの定数を用いた場合の実行例。

sg <- segment_trees(las, li2012()) # 左上
sg <- segment_trees(las, watershed(chm)) # 右上
sg <- segment_trees(las, dalponte2016(chm, treetops = ft)) # 左下
sg <- segment_trees(las, silva2016(chm, treetops = ft)) # 右下

平面樹冠ポリゴン作成

作成したtreeIDを2次元ポリゴンに落とし込む。

poly <- delineate_crowns(sg)

ポリゴンデータの保存については既出のfind_trees()の節を参照のこと。

ポリゴンデータの作図

各ポリゴンに樹頂高さが保持されているので、それを用いた描画を行う。

png("ファイル名.png", width = 720, height = 720)
par(plt = c(0.1, 0.9, 0.1, 0.9))
spplot(poly, "ZTOP", main=list(label="ZTOP [m]",cex=2), at = seq(5, 30, 0.5), colorkey=list(labels=list(cex=1.5)))
dev.off()

空間演算

点群全体の統計値

下記は、基本統計値のセットを事前に取り揃えた.stdmetricsを使用したもの。
公式チュートリアルでは、インベントリデータと組み合わせた統計解析を行う応用例が紹介されている。

cloud_met <- cloud_metrics(las, .stdmetrics)
print(cloud_met)
# 保存
write.csv(as.data.frame(cloud_met), "ファイル名.csv")

# $zmax
# [1] 1254.27
# $zmean
# [1] 1197.117
# $zsd
# [1] 34.80193
# :   :   :
# $itot
# [1] 0
# $imax
# [1] 0
# $imean
# [1] 0
# :   :   :
# $n
# [1] 3055555
# $area
# [1] 108133.1

水平グリッドごとの演算

今回は、自作関数を用いた例で実践してみる。今回のデータでIntensityの値が定数のため、Intensityを入力には用いていないが、Zと合わせて2変数以上の解析も可能。

myMetrics = function(z) {
    metrics = list(
        npoint  = length(z),
        zmax    = max(z),
        zmean   = mean(z)
    )
    return(metrics)
}
grid_met = grid_metrics(las, ~myMetrics(Z), res = 1)

# 保存
writeRaster(grid_met$npoint, filename="ファイル名.tif", format="GTiff", overwrite=TRUE)   

# 描画
png("ファイル名.png", width = 720, height = 720)
par(plt = c(0.1, 0.9, 0.1, 0.9), cex.axis=1.5, cex.main=1.5)
plot(grid_met$npoint, main='NPOINT', col = height.colors(255), zlim = c(0, 200), zlimcol=FALSE)
dev.off()

計算する統計値が複数の場合、出力はRasterBrickとなる。すべての変数名やデータを保持して保存する方法もあるとは思うが、調べてもわからなかったため、上の例ではLayerごとに分けて保存している。

葉面積指数（LAI）の計算

lidRにはLAD()が用意されており、各層の葉面積密度（LAD）の計算を行う。これを1mごとで計算すれば、LAIはそれを鉛直方向に和をとったものに等しい。grid_metrics()を用いて、水平ピクセルごとのLAIを計算する。

LAI <- function(z) 
{
    lad <- LAD(z, dz=1, z0=2)
    if(class(lad) == "numeric")
        return(NA_real_)
    r = sum(lad$lad, na.rm=TRUE)
    return(r)
}

lai <- grid_metrics(las, ~LAI(Z), res = 10)

# 保存
writeRaster(lai, filename="ファイル名.tif", format="GTiff", overwrite=TRUE)
# 描画
png("figs/ファイル名.png", width = 720, height = 720)
par(plt = c(0.1, 0.9, 0.1, 0.9), cex.axis=1.5, cex.main=1.5)
plot(lai, main='LAI', zlim = c(0, 6), zlimcol=FALSE)
dev.off()

3次元グリッド（voxel）ごとの演算

下記では、各グリッドごとの点の数を計算している。出力はXYZの3列と計算値（下記ではN）を含むdata.frameとなるため、保存時はcsv形式が代表例となる。

vox_met <- voxel_metrics(las, ~list(N=length(Z)), res=5)

# 保存
write.csv(vox_met, "ファイル名.csv", row.names=FALSE)
# 描画
plot(vox_met, color="N", trim=800) # trimは色の最大値

葉面積密度（LAD）の計算

LADの計算自体、鉛直方向に一度に一気に行う仕様になっているため、ボクセルごとにLADを保持するためには、下記のようにgrid_metrics()を用いる。一つの計算値に対して配列を出力できず、単独値にする必要があるため、計算値を鉛直方向1mごとに保持する形にする。

vLAD <- function(z) {
    lad <- LAD(z, z0=1.0, dz=1)
    metrics = list(
        z015 = 0,
        z025 = 0,
        （中略）
        z335 = 0,
        z345 = 0
    )
    for (i in seq(1, 34, 1)){
        metrics[i] = lad[i,"lad"]
    }
    return(metrics)
}

lad <- grid_metrics(las, ~vLAD(Z), res = 5)
# 保存
write.csv(as.data.frame(lad, xy=TRUE), "ファイル名.csv", row.names=FALSE)
# レイヤごとの描画（ファイル名にもzの値を含める例）
library(stringr)
for (i in 1:34){
    name = str_pad(i * 10 + 5, 3, pad = "0")
    v = paste("z", name, sep="")
    png(paste("ディレクトリ名/lad_z", name, ".png", sep=""), width = 720, height = 720)
    par(plt = c(0.1, 0.9, 0.1, 0.9), cex.axis=1.5, cex.main=1.5)
    plot(lad$v, main=paste("LAD (Z = ",  i, ".5 m)", sep=""), zlim = c(0, 1.5), zlimcol=FALSE)
    dev.off()
} 

その他

複数ファイルにまたがる広域lasデータセットに対しては、lasCatalogというまとめて処理する機能が備わっている。空間演算を活用すれば、より多様な解析が可能になる（屋根・電線の分類など）。

1. この記事にあるEPSG情報をもとに、この記事で紹介しているpythonの座標変換ツールを用いて、JGD2000（静岡県のEPSG:2450）から緯度経度座標に変換できることから確認した。 ↩ ↩²

2. IRkernelをインストールしてJupyter notebookにてRを用いる場合、3次元描画はうまくいかない（ウィンドウは出現するものの、インタラクティブに点群を動かすことができない）。WSL2のRStudioを用いる場合は、この記事の通りに、事前にX11設定をする必要がある。 ↩

3. <余談> zlimcolの機能を担うオプションがほとんど検索に引っかからなくて困った。公式マニュアルのplot()にも記載がないのだが、こういうのはどういう風にして調べるのがよいのだろうか。今回は、ググって何とか情報を見つけたのと、plotにわざとエラーを吐かせてtrackbackからソースコードを発見した。 ↩