Rで航空機LiDARデータ解析

# はじめに

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

# 準備・環境

環境はWindows10のWSL2（Ubuntu 20.04）。
Rのバージョンは4.0。Rは初めたばかりで試行錯誤しながら手を動かしているところ。
パッケージはいろいろと用いているが、メインはLiDARデータ処理用パッケージ[LidR](https://cran.r-project.org/web/packages/lidR/index.html)（[Github](https://github.com/Jean-Romain/lidR)、[マニュアル](https://cran.r-project.org/web/packages/lidR/lidR.pdf)はこちら）。
つい先日（2020年7月）、パッケージ製作者が[チュートリアルブック](https://jean-romain.github.io/lidRbook/index.html)を作られたようなので、こちらもぜひ。

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

# 使用データ

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

![図0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/2185bc28-1798-4505-5f4d-3b090063b6da.png)

![図1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/8fed231b-82d8-3c75-d42b-25361bfd7aed.png)

# とりあえず開く

点群データはlasフォーマットにて配布されている。ArcGISの説明[①](https://desktop.arcgis.com/ja/arcmap/latest/manage-data/las-dataset/storing-lidar-data.htm)、[②](https://desktop.arcgis.com/ja/arcmap/latest/manage-data/las-dataset/lidar-point-classification.htm)などが参考になる。

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

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

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

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

[^1]: [この記事](http://tmizu23.hatenablog.com/entry/20091215/1260868350)にあるEPSG情報をもとに、[この記事](https://qiita.com/HidKamiya/items/e67423a8254841ce6d46)で紹介しているpythonの座標変換ツールを用いて、JGD2000（静岡県のEPSG:2450）から緯度経度座標に変換できることから確認した。

```
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上で変換した緯度経度を示した図が下の通り。領域全体をほぼ同じ大きさの樹木が覆いつくしていることがわかる。

<img width="300" alt="図.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/df43ea95-3cec-164b-739a-6e82d6ba3e94.png">



## 点群データの読込

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

```r:
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](https://cran.r-project.org/web/packages/rgl/index.html)という3次元描画パッケージをベースにした`plot()`が用意されている。[^2]

[^2]: [IRkernel](https://cran.rstudio.com/web/packages/IRkernel/index.html)をインストールしてJupyter notebookにてRを用いる場合、3次元描画はうまくいかない（ウィンドウは出現するものの、インタラクティブに点群を動かすことができない）。WSL2のRStudioを用いる場合は、[この記事](https://qiita.com/HidKamiya/items/0f5fc5eb0bbd4d145cd5)の通りに、事前にX11設定をする必要がある。

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

```r:
plot(las)

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

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

![図10 (2).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/a41932c3-86a2-39f2-051b-24f36c2664e4.png)


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

```
plot(las, color=Classification)
```

![class.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/41a9de55-dee2-93e5-b152-3ee9cdedf084.png)

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

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

```r:
plot(las, color ="RGB")
```

![snap_rgb.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/b7dbba97-3b2c-7a54-7d76-4a9d318d2606.png)


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

Classificationのうち1つに絞って描画をするために、`filter_poi()`を用いる。
今回は、分類2（地表、[参照](https://desktop.arcgis.com/ja/arcmap/10.3/manage-data/las-dataset/lidar-point-classification.htm)）にフィルタリングする。

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

![class02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/fb156a4b-89aa-1971-3a56-61d5a97d5558.png)

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

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

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

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

## 点群データの保存

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

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

# 地面ポイントの抽出

地表を表すClassificationを再計算して求めることにする。
使用するのは`classify_ground()`という関数。中で、`pmf()` (Progressive Morphological Filter)というアルゴリズムを用いているが、詳細は[公式チュートリアル](https://jean-romain.github.io/lidRbook/gnd.html)を参照のこと。

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

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

![cg02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/ace8ca0c-322e-ad41-2f2f-6ef123f5cf6c.png)

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

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

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

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

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

DTMはRasterLayerというデータ型に保持された2次元ラスタデータで、`library(raster)`の`plot()`で描画する（[参考](http://mikemeredith.net/blog/2019/plotting_rasters.htm)）。`zlim`、`zlimcol`（カラーバーの外れ値をNullにするフラグ）は[`.rasterImagePlot`](https://rdrr.io/cran/raster/src/R/imageplot2.R)のオプションで、`plot()`を用いると内部的にこの関数のオプションが使用できるようになるらしい。[^3]

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

```r:
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()
```

![dtm2d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/d0a702da-afb1-f12b-f871-0a71c7f48ad1.png)

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

## 標高マップの3次元描画

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

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

![dtm3d10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/ba1f30ae-7469-9b8e-ecc0-770ed8948dad.png)


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

```r:
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")
```

![dtm3d03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/03be3929-8c1e-6f2f-a230-192cad321fda.png)


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

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

保存と読込はそれぞれ、library(raster)の`writeRaster()`と`raster()`を用いる。保存形式はGeoTiffとしているが、ほかにも[いくつかの形式](https://www.rdocumentation.org/packages/raster/versions/3.3-13/topics/writeRaster)がサポートされている。

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

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

# 標高差をなくす

点群からDTMを減算する形で、地表面標高を均一にする。dtmがない場合でも、いくつかのアルゴリズムから均一化は行える（詳しくは[公式チュートリアル](https://jean-romain.github.io/lidRbook/norm.html)にて）。

```r:
las <- normalize_height(las, dtm)
```

![norm.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/17f981e3-98d5-4301-0f15-2c6063e50203.png)

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

# 樹高マップ（CHM作成）

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

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

![chm3d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/5c2a738f-b4f2-3340-3432-f372e5598c94.png)


# 樹頂（樹数）推定

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

```r:
ft <- find_trees(las, lmf(5))

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

![ft_merge.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/1a222bf2-486a-75ab-2237-05bbb0f8dd8b.png)

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

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

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

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

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

# 樹冠推定

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

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

```r:
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)) # 右下
```

![sg_merge.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/9d622e3c-efd6-c396-eefb-234a224e9238.png)

# 平面樹冠ポリゴン作成

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

```r:
poly <- delineate_crowns(sg)
```

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

## ポリゴンデータの作図

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

```r:
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()
```

![de_marge.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/dcd2bac7-1843-5ecd-d518-93bb29d986b7.png)


# 空間演算

## 点群全体の統計値

下記は、基本統計値のセットを事前に取り揃えた`.stdmetrics`を使用したもの。
[公式チュートリアル](https://jean-romain.github.io/lidRbook/cba.html)では、インベントリデータと組み合わせた統計解析を行う応用例が紹介されている。

```r:
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変数以上の解析も可能。

```r:
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()
```
![grid_npoint.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/11f2df9a-9962-a9fa-7717-83afaca76636.png)


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


### 葉面積指数（LAI）の計算

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

```r:
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()
```

![grid_lai_10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/420be493-890b-9dd9-1c8f-92681d8e8377.png)


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

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

```r:
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は色の最大値
```

![vox_npoint_50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/ca166129-af6a-0a3c-b5d6-49102b8dbc28.png)


### 葉面積密度（LAD）の計算

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

```r:
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()
} 
```

![lad.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187118/ec40c782-67c6-b25f-7074-0a1312ff59ef.gif)


# その他

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