三次元点群データをQGISでひり出してみた

![](https://i.imgur.com/uF4vvTN.gif)

```
0.4秒くらいでこの記事を読み終えたいかた向けにオチを書いておきますと
【三次元点群データ処理…それQGISでもできますよ】
です(約2万字・読了15分)。
```

## はじめに

先日、自由でオープンなジオツールやデータを東京と札幌でいじる[GeoSaturday](https://hackmd.io/@nyampire/SyM94GE5S/)というワイガヤ勉強会を開催し、東京札幌で50名を超える参加がありました。皆様ありがとうございます。

そのなかで、林業クラスタの方が、レーザー航空測量の3次元データ（LiDARデータ）がテキストデータであるのだけど、オープンソースの地図ソフト[QGIS](https://qgis.org/ja/site/)で表示したり変換するにはどうすれば…という課題をお持ちでした。

海外ではオープンソースソフトとLiDARデータをビシバシいじって可視化したり、予想される木材の量などの算出がされているようです。しかしこのような手法の日本語情報は少なく、可視化するだけでも業者さんに頼んだり高いソフトを購入したりという状況のようです。

また、行政機関などによるオープンデータの推進が進む中、[ロックな静岡県庁さん](https://pointcloud.pref.shizuoka.jp/)のように今後も3Dデータの公開はガンガン進んでいくとみられます。

3Dデータに関係ありそうな業界（農林・建設・自動運転・MaaS…etc）の方はいまのうちにおうちのQGISで3Dデータのハンドリングのコツに触れておくのも良いことです。

ということで、贖罪と義憤()と未来への期待から、QGISでLiDARデータをいじる方法を調べ倒してまとめてみました。

なお、あれこれ試しながらのトライアルだったので、実は玄米からピザを作ろうとしているような回り道をしているかもしれませんのでご容赦ください🙇。
### こんな方に

本記事は

- オラの山のLiDARデータがあるのに扱えるツールがない
- ツールがないのでコンサルに可視化作業見積とったら股間が震え上がった
- もっとデータとQGISを駆使して楽に儲けたい
- LiDARデータなんかオープンにして意味あんの?と予算査定で酷評された
- 定番宴会ネタ「QGISでできること100連発」に新しい項目を足したい

といったバックグラウンドをもったかたを想像して書いています。


### やること

ここでは

- QGISでのLiDAR関係ツールの設定
- 3次元テキストからのデータ処理
- 処理結果から樹冠と樹高を算出

…を紹介します。なお、

- LiDARデータの美しいクリーニング
- 樹冠構造から樹種を推定する
- たくましい木こりを伴侶にする秘訣
- 亡くなった愛犬を団子鼻の模様まで3Dプリンタで復元

などは今回扱いません、すみません。


### 3次元テキストデータとは

さて、今回相談で見せていただいたレーザー測量による３次元テキストデータこんな感じでした

```
（ヘッダなし）
962129.89,3973446.45,977.09,2,160640.130496,10,73,1,1,21,0,1,32,0,0,0
962129.92,3973445.91,977.00,2,160640.130506,10,62,1,1,21,0,1,32,0,0,0
962129.94,3973445.36,976.97,2,160640.130516,10,67,1,1,21,0,1,32,0,0,0
962129.97,3973444.81,976.94,1,160640.130526,10,62,1,1,21,0,1,32,0,0,0
…
```
データの仕様書がないので概ね推測ですが、たいていこのようなデータは位置のXYZとその他諸々の情報（色やオブジェクト分類フラグ）などが入っており、たいてい1ファイル数百M〜数GBほどの大きさになります。

これらは[LiDAR](https://ja.wikipedia.org/wiki/LIDAR)というリモートセンシング手法で測定された3次元データです。最近ではドローンにも機器が搭載されたり、自動運転技術にも応用されつつあります。[国土地理院さんの解説もわかりやすいです。](https://www.gsi.go.jp/kankyochiri/Laser_senmon.html)

<img src="https://i.imgur.com/6RQlyFJ.png" width=80%>!|
|:-:|
出典：[航空レーザ測量の仕組み（国土地理院）](https://www.gsi.go.jp/kankyochiri/Laser_senmon.html)

これらのデータは近年の高精度化によりさらにファイルサイズが増える傾向にあるため、デスクトップのQGISで何でもこれらのデータを処理しまくるのはママチャリでツールドフランスに出場を目指すようなものかもしれません。おそらく今後は[クラウドでの処理がベースになってくるでしょう。](https://hackmd.io/@colspan/one-stop-3dtiles-on-GoogleColab#/)

とはいえ、現場の方が既存の地図やオープンデータなどの位置データとさっと重ねて確認したり、すぐ作業に使うためには使い慣れたQGISで活用したほうが良い場合もあります。



ε｀ ).oO(15年ほど前ですが、レーザー計測結果のtxtデータがCD-ROM()20枚くらい届き、どうにかしてラスタにした思い出がフラッシュバックしました。たしか業者さんに頼んだか有償のソフトでむりくりSHPやラスタにしたものの当時のPCスペックでは開くことすらままならず…結局どうやったのか思い出せないです。


### QGISにおけるLiDARツールの現状

QGISでこの手の処理を行うには[LAStools](http://lastools.org/)という外部ツールがよく使われています。

LAStoolsはノースカロライナ大学チャペルヒル校で開発されたプログラムでLiDARデータの処理に特化したツールが数多く登録されています。

LAStoolsに関しては[このwiki](http://wiki.awf.forst.uni-goettingen.de/wiki/index.php/Canopy_Height_Model_based_on_Airborne_Laserscanning_using_LAStools#Create_a_DSM_with_abolute_heights:_las2dem)にいくつか機能が紹介されているのでけっこう役に立ちました。応用的な利用をしたい方にはおすすめです。英語が苦手な方に[Google翻訳版](http://translate.google.com/translate?langpair=en%7Cja&hl=en&u=http://wiki.awf.forst.uni-goettingen.de/wiki/index.php/Canopy_Height_Model_based_on_Airborne_Laserscanning_using_LAStools#Create_a_DSM_with_abolute_heights:_las2dem)も貼っておきます。

また、[QGISのオフィシャル教本](https://docs.qgis.org/3.4/ja/docs/training_manual/forestry/basic_lidar.html)でもLAStoolsの使い方は紹介されているのですが、スクショが古いままだったり~~かなり難解~~少々ややこしいのでそこから意訳(超訳)して、本記事の参考にしました。

ここまで紹介しといてなんですが、LAStoolsはいろいろ使えて大変便利そうな反面、（仕方ないのですが）一部のツールの[ライセンスは有償](http://lastools.org/LICENSE.txt)とされています。

![](https://i.imgur.com/hiGZaAQ.png)

たとえば、この記事の後半で紹介するLASデータからDEMラスタに変換する処理を無料ライセンスで行うと、このような焼印的斜め線検閲が"意図的に"ガッツリ入れられてしまうことがあります。
(フィールド調査でよく見た[トドマツに残されたヒグマのツメ痕](https://www.google.com/search?sxsrf=ACYBGNTFJzRz67FNjoCzEK5LaJrpMekZLg:1577011041882&q=%E3%83%92%E3%82%B0%E3%83%9E+%E7%88%AA%E7%97%95+%E3%83%88%E3%83%89%E3%83%9E%E3%83%84&tbm=isch&source=univ&sa=X&ved=2ahUKEwjx9NT3h8nmAhVD7GEKHY8zBbgQsAR6BAgJEAE&biw=1422&bih=711)を思い出しました)

回避策はいくつかあるようで、

- お金出しましょう(すみません)
- QGISにバンドルされてるSAGAとかGRASS使え（ですよね）
- QGISに[FUSION](http://forsys.cfr.washington.edu/fusion/fusionlatest.html)とか[pktools](http://pktools.nongnu.org/html/index.html)入れるとよい

…などといった情報がありました。
ということで、今回はLAStoolsのほかにLiDARデータ処理ツールも色々入っているFUSIONもQGISにいれてしまいましょう。

※SAGAやGRASSはマシンやバージョン条件によってけっこう不安定な挙動をするので今回は割愛しました。

参考：
https://gis.stackovernet.com/ja/q/49582
https://gis.stackovernet.com/ja/q/41447
https://gis.stackexchange.com/questions/181712/creating-dem-from-las-file-without-using-lastools/206332#206332


### 作業フロー

作業ゴールとしては、木のてっぺんである樹冠とその高さをオープンデータとオープンソース（QGIS+その他）で抽出することを目的としました。

これらのデータがあると、空からの情報で作業区域内の樹種や材積などが推定でき、効率的な施業にも貢献するされています。

樹冠の算出に必要なデータとしては、DSM/DTM/CHMの３項目になります。

急に謎の3文字が飛び交うようになってしまい、読み手を[22世紀の傑作なぞなぞと称される「アルファがベータをカッパらったらイプシロンした。なぜだろう」を聞いたのび太](https://www.google.com/search?q=%E3%82%A2%E3%83%AB%E3%83%95%E3%82%A1%E3%81%8C%E3%83%99%E3%83%BC%E3%82%BF%E3%82%92%E3%82%AB%E3%83%83%E3%83%91%E3%82%89%E3%81%A3%E3%81%9F%E3%82%89%E3%82%A4%E3%83%97%E3%82%B7%E3%83%AD%E3%83%B3%E3%81%97%E3%81%9F%E3%80%82%E3%81%AA%E3%81%9C%E3%81%A0%E3%82%8D%E3%81%86%E3%80%82&sxsrf=ACYBGNSgRsrHPpgRIvTMbIXNQ8X_suigmA:1577367305257&source=lnms&tbm=isch&sa=X&ved=2ahUKEwiO0qOPt9PmAhWXF4gKHanlDW0Q_AUoAXoECAwQAw&biw=1288&bih=884)状態にさせてしまい恐縮ですが、このDSM/DTM/CHMの関係はこのようなイメージです。

<img src="https://i.imgur.com/CH9m9LF.png" width=60%>
via https://www.earthdatascience.org/courses/earth-analytics/lidar-raster-data-r/lidar-chm-dem-dsm/

つまり、樹冠算出の元になるCHMを出すために、LiDARから算出したDSMとDTMが必要ということになります。

今回行う処理フローはこんなかんじになります

![](https://i.imgur.com/YykyLpp.png)|
|:-:|

※Silva, Carlos Alberto, et al. "LiDAR remote sensing to individual tree processing: A comparison between high and low pulse density in Florida, United States of America." XVI Brazilian Remote Sensing Symposium, Foz do Iguaçú, Brazil. 2013.
の[図2を参考にしました。](https://www.researchgate.net/figure/Steps-to-individual-tree-processing-DTM-Digital-Terrain-model-DSM-Digital-Surface_fig2_236620625)


本来ですともうすこし手間を減らしたりLiDARデータのクリーニングなどをしないといけないようなのですが（例えば佐々木(2012)）、わたしはLiDAR専門家でもないので、本記事では荒くても樹冠算出のトライをQGISでやってみました。

フローが多段階なのでうまくいったらQGISのプロセッシングツールのグラフィカルモデラーなどで一括でできるのではないかなと思います。ただしLiDARのデータフォーマットがけっこう[わんぱくなので](https://hackmd.io/@colspan/one-stop-3dtiles-on-GoogleColab#/11)パラメータ設定などは柔軟な設定にしておいたほうがよいかもしれません。

## 環境構築

ではまずQGISの環境を整えましょう、ここからの作業は

- Windws10-64bit
- QGIS 3.10.11

で行います。安定版(LTR)のQGIS3.4でもOKだと思います。

なお最初わたしのMacでやろうとして全然動かなくて3日ほどほったらかしたのですが、よく見るとツール内で*.exeを多用してました…。
どうやらMacOSはLAStools/FUSIONが正式に対応してないようなので今回はMacでのやり方は割愛します。 

※[wineとか使うとMacでも使える](https://gis.ubc.ca/2018/02/installing-lastools-for-lidar-data-in-qgis-for-mac/)ような情報もありますが、今回は試してないですごめんなさい。
※[wineのセットアップはこのあたりを御覧ください](https://weblabo.oscasierra.net/wine4-macos/)

### LAStoolsのセットアップ

では環境構築をおこなっていきます。
まずQGISのプラグインリポジトリを見ると、[LAStools](https://plugins.qgis.org/plugins/LAStools/)がすでにあります

|![](https://i.imgur.com/lZu4jV9.png) |
|:-:|

まずこれをQGISにインストールしましょう。プラグインの追加方法がわからない方は[このへん](https://gis-oer.github.io/gitbook/book/materials/QGIS/QGIS.html#%E3%83%97%E3%83%A9%E3%82%B0%E3%82%A4%E3%83%B3)をどうぞ。

これでQGISのプロセッシングツールボックスにLAStoolsのグループが追加されます。

|![](https://i.imgur.com/IkYezda.png)|
|:-:|

お！これですぐ使える！…と思いがちなのですが、このプラグインは「ガワ」のみです。いわば「中身のない水餃子」状態です。つまり単にwindows内のLAStoolsプログラムにパラメータを投げてくれる役目だけでして、プラグインを起動させるとウイザードが出ますが、処理をかけようとするとエラーが出て計算してくれません。ということでLAStools本体をwindowsにインストールしましょう。
LASToolsの本体は[こちらからDLできます]( http://rapidlasso.com/LAStools/)。zipを解凍してできた `LAStools`フォルダをCドライブ直下に置きます。

![](https://i.imgur.com/eFsE5cb.png)|
|:-:|

QGISにもどり、プロセッシングツールボックスのスパナアイコン（オプション）をクリックします。

![](https://i.imgur.com/m0DQVfP.png)

「プロバイダ」からLAStoolsをクリックしてActivateに✔をいれて、LASTools folder行の横をクリックして`C:\LAStools`を入力します。

![](https://i.imgur.com/Gfe31VH.png)|
|:-:|

ウインドウ右下のOKを押すとQGISからLAStoolsでの処理が可能になります。（OK押して適用忘れをしやすいのでご注意を）

.oO(今気づいたのですがwineの設定の行もあるのでMacな方はここにwineのパス入れると起動するかも…)

### LAStoolsの起動確認

ほんとにLAStoolsがQGISから起動するか確認するために、[ロックな静岡県庁さん](https://pointcloud.pref.shizuoka.jp/)がオープンデータとして提供されているLASデータを読んでみましょう（Sさんいただきます！）。

![](https://i.imgur.com/KssS0me.png)|
|:-:|


ここからたどって、たとえば

という[`平成29年度建設工事技術管理事業UAV計測及び活用に係る説明資料作成業務（浜松市発注）（浜松城_2）`](https://pointcloud.pref.shizuoka.jp/lasmap/ankendetail?ankenno=29XXX00010002)というガチ現場オープンデータから29XXX00010002-3.lasをDLしてきます（200Mほどあるのでテザリングやテレホーダイタイムにご注意を）。


テストということでこのデータをQGISからLAStoolsを操って可視化してみましょう。QGISのプロセッシングツールボックスからLAStools→lasviewを選択します。

![](https://i.imgur.com/GbrUdjF.png)|
|:-:|

input LAS/LAZ fileには
さっきDLしてきた`29XXX00010002-3.las`を指定して、とりあえずその他設定はデフォルトのままにして

![](https://i.imgur.com/9Ax5kkR.png)

OKを押すと別ウインドウで何か立ち上がり…

|![](https://i.imgur.com/BUkPbzn.png)|
|:-:|

*出典：「静岡県ポイントクラウドデータベース」under CC BY*

ドラッグして視点を変えると浜松城が見えるはずです。ドラッグで位置も変更でき、スペースキーで画面の動かし方も変えれます。
また右クリックで背景を黒にするなどのオプションも選べます。これで、QGISでのLAStoolsのセットアップが完了しました。


![](https://i.imgur.com/YLWDI6i.gif)|
|:-:|
*出典：「静岡県ポイントクラウドデータベース」under CC BY*

### FUSIONのセットアップ

次にLAStoolと同様に様々な処理ができるFUSIONもQGISに入れてみます。[FUSION](http://forsys.cfr.washington.edu/FUSION/fusion_overview.html)はアメリカ農務省の国有林管理部門太平洋岸北西部事務所で開発されたLiDAR解析ソフトウェアで、アメリカ合衆国著作権法(17 U.S.C. §105)に則って[パブリックドメインとされています](https://ja.wikipedia.org/wiki/%E3%82%A2%E3%83%A1%E3%83%AA%E3%82%AB%E5%90%88%E8%A1%86%E5%9B%BD%E6%94%BF%E5%BA%9C%E3%81%AE%E8%91%97%E4%BD%9C%E7%89%A9)。

QGISのプラグインリポジトリには[FUSION for processing](https://plugins.qgis.org/plugins/processing_fusion/)があり、まずこれをインストールします。
![](https://i.imgur.com/ROQ9NM0.png)


これでプロセッシングツールボックスにFUSIONのグループが追加されます。

![](https://i.imgur.com/iqedJer.png)|
|:-:|

ご多分にもれず、LAStoolsと同様にこれもいわば「ガワ(中身なし水餃子)」状態ですので、FUSION本体をwindowsにインストールします。

FUSION本体はちょっと90'sな感じのする[こちらのサイト](http://forsys.cfr.washington.edu/fusion/fusionlatest.html)の`Install file` からインストーラーがダウンロードできます。

![](https://i.imgur.com/47Dlezs.png)|
|:-:|


インストーラーの`fusionlatest.exe`を管理者権限で起動し、すべてOKを押すと`C:\FUSION\`にFUSIONのプログラム一式がインストールされます。

また、LASToolsと同様に
QGISにもどり、プロセッシングツールボックスのスパナアイコン（オプション）をクリックします。

![](https://i.imgur.com/m0DQVfP.png)


Activateに✔をいれて、FUSION folder行の横をクリックして`C:\FUSION`を入力します。
![](https://i.imgur.com/enQgpUy.png)

ウインドウ右下のOKを押すとQGISからのFUSIONを使った処理が可能になります（これも右下の`OK`を押すのを忘れやすいのでご注意を）

### FUSIONの起動確認

起動確認のために、さっきと同様にロックな静岡県庁さんによる浜松城を見てみましょう。
FUSIONの場合は、一番下の
Visualizations→Open viewer
から読み込みます。 

![](https://i.imgur.com/vwLwOqY.png)

PDQというビューワーが起動します

![](https://i.imgur.com/qFa1ZtG.png)

メニューのFile→Open…
で先程の`29XXX00010002-3.las`を読み込んでみます。LAStoolsより表示までの時間がかかりますが、こんな絵が出るはずです。点の色情報を読んでないのでちょっとそっけないですが、その分表示が早いような気がします。

|![](https://i.imgur.com/XtXDrrm.png)|
|:-:|

これで、QGISでのFUSIONのセットアップも完了しました。


## 処理してみる
### 3Dテキストファイルの準備

やっとツールの準備が整いましたので、では材料となるサンプルの3Dテキストファイルを準備します。

世界各地の地形に関するオープンデータがダウンロードできる[OpenTopography](https://opentopography.org/)のサイトから[熊本のレーザー測量テキストファイル](http://opentopo.sdsc.edu/lidarDataset?opentopoID=OTLAS.052018.2444.2
)を頂戴してきます。データ諸元によると平面の投影法はEPSG:2444です。

|![](https://i.imgur.com/hmaWBzC.png)|
|:-:|



データの出典はこちらになります、OpenScienceですね、ありがとうございます。

>Scott, C. P., Arrowsmith, J. R., Nissen, E., Lajoie, L., Maruyama, T., & Chiba, T. (2018). The M7 2016 Kumamoto, Japan, Earthquake: 3-D Deformation Along the Fault and Within the Damage Zone Constrained From Differential Lidar Topography. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2018JB015581

>Scott, C., Champenois, J., Klinger, Y., Nissen, E., Maruyama, T., Chiba, T., & Arrowsmith, R. (2019). 2016 M7 Kumamoto, Japan, Earthquake Slip Field Derived From a Joint Inversion of Differential Lidar Topography, Optical Correlation, and InSAR Surface Displacements. Geophysical Research Letters. https://doi.org/10.1029/2019GL082202
>

この[OpenTopography](https://opentopography.org/)サイト、測定範囲一括DLだとデータ量が膨大になるので矩形選択で一部を切り取ってDLすることもできるので、今回のようなちょっとした練習にはありがたいです。

データを切り取る矩形選択はちょっと分かりづらいですが地図の左上のボタンになります。今回は北東の朝来山らへんを選択しました（山あいの任意場所でOKです）

<img src="https://i.imgur.com/6wJ8TTx.png)" width=60%>

ファイル形式はLASやTIFFをDLしちゃうと練習にならないｗので、ASCIIフォーマット（テキストファイル）を選択します

![](https://i.imgur.com/5D80L19.png)

生成処理が終わったら教えてくれるメールアドレスなどを入力し、`SUBMIT`を押します（メールの宛先を実家のオカンや気になるあの子にしてはいけませんよ）。

![](https://i.imgur.com/4MB3l0I.png)|
|:-:|


生成中はこのような画面が表示されます

![](https://i.imgur.com/iT8pCM7.png)|
|:-:|



メールを確認するまでもなく生成がすぐ終わるとおもいますので
Download point cloud data in ASCII format   **points.tar.gz** (xx MB)

をクリックしてダウンロードを行います。

![](https://i.imgur.com/cOON9TE.png)|
|:-:|

このファイルはtar.gz圧縮されていますので、[7-zip](https://sevenzip.osdn.jp/)などで解凍するとこのようなpoints.txtのテキストができます。
（これをメモ帳で開くとハングアップしたりしますので、わたしは巨大テキストファイルでもさくっと開いてくれる[VSCode](https://azure.microsoft.com/ja-jp/products/visual-studio-code/)で確認してみました、まあエディタはお好きなやつで…）

![](https://i.imgur.com/eRNhGyQ.png)

今回のテキストファイルはファイルサイズ120MB、約313万行でした（VSCodeでは表示まで５秒ほどでした）。

### *txt→ *.lasに変換

では、このテキストを汎用的なlas形式にするために、QGISのLAStoolsで変換してみます。QGISのプロセッシングツールボックスにある`LAStools→file-conversion→txt2las`を選択します。ここではASCIIファイル→テキストファイル、という意味です。あとででてくるASCIIラスタとはちょっと異なりますのでご注意を。

![](https://i.imgur.com/evwcoqs.png)|
|:-:|

埋めるべき情報は

- run new 64 bit executable(処理の高速化)…`✔`
- Input ASCII file(材料テキストファイル)…お手持ちのtxtなど
- parse line as(列の並び)
…たいていは左からXYZ列なので`xyz`で良いのですが、ときどき~~おてんば~~独自のフォーマットで左端にIDとかついている場合は**S**kipをいれて
`sxyz`※１と記述します。

このほかの設定は

- skip the first n lines(スキップする行数)
…今回は列名ヘッダがついてたので１を入力
- projection
…今回はEPSGのデータを使うので`epsg`を選択する
- EPSG code
…ファイル諸元にあったようにEPSGコード`2444`を入力
- Output LAS/lAZ file
…出力先を選択してファイル名を入力。
(ファイル名だけ入れるとエラーが出る場合、パスとファイル名のあとに手打ちで.lasの拡張子まで入れてあげると通ります)


<img src="https://i.imgur.com/9hpIlMx.png" width=60%>

「実行」を押すと、すこし計算後変換されたlasファイルが生成されます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/100567/21356599-068b-488b-9b49-3e2ad02b133c.png)

ためしに先述のlasviewで見てみると、おお、できてますね！
<img src="https://i.imgur.com/2sTPh9d.png" width=80%>

※１
カラムごとの独自情報をtxtから取り込みたい場合は下記のような記述をします

    * i - intensity
    * n - number of returns of given pulse
    * r - number of return
    * c - classification
    * u - user data
    * p - point source ID
    * e - edge of flight line flag
    * d - direction of scan flag.

参照：
https://liblas.org/utilities/txt2las.html


###  *.lasからDSM作成
ではこのlasデータをDSMに変換してみましょう。
こんどはQGISからFUSIONを使います。QGISのプロセッシングツールボックスにあるFUSION→Grid surface createでまず*.dtmファイル形式でDSMを作成します。

なお、ここでいう*.dtm形式は地形標高モデルのDSM/DTMのDTMや宅録のDTMとは違いますので、ご注意を。

![](https://i.imgur.com/UEazdwi.png)|
|:-:|

材料lasデータ、セルサイズ（1とか0.5とか）、XYZの単位（今回はMeter）、出力先を設定します。

![](https://i.imgur.com/pQ2JUBg.png)|
|:-:|

実行を押すと、DSMファイルが生成されます

### DSMをラスタ形式にしてみる

さきほど作成されたDSMファイルをFUSIONツール内のConvension→DTMtoASCIIで変換します(DTMtoTIFというツールもありますがD高さ(m)が256階調に変えられてしまうので今回は使いません)

![](https://i.imgur.com/TKxfLk0.png)|
|:-:|

とりあえず入力ファイル(DSM)と出力ファイルのみを入れましょう
![](https://i.imgur.com/Fi3SMKB.png)


するとQGISのマップキャンバスに変換されたTIFが追加されます
<img src="https://i.imgur.com/lZyXAa2.png" width=80%>

…生成されたGeoTiffには投影法が付与されないので、[地理院地図を追加後](https://qiita.com/Yfuruchin/items/735f786a86380a7c24ac)、マップキャンバスを強引に元データの投影法であるEPSG：2444にしてみると、お、位置合ってます（ラスタメニューから`投影法の設定`でもOKです）
<img src="https://i.imgur.com/hcUI6za.png" width=80%>
出典：地理院タイル

国土地理院写真レイヤを重ねて[QGIS２Threejs](https://qgis2threejs.readthedocs.io/ja/docs/Tutorial.html)で３次元化すると妙にココだけリアリティのあるDSMが表現されています。左の畑の畦もくっきり出ています。
<img src="https://i.imgur.com/vbW6bml.png" width=60%>
出典：地理院タイル

こっちは標準で読める地理院さんのPNG標高タイルから読んできた同じ地域ですが、DEMであることやセルサイズも大きめなのでlasデータと比較するとそれぞれの特徴が一目瞭然ですね。
<img src="https://i.imgur.com/mSv1HXY.png" width=60%>
出典：地理院タイル

###  *.lasからDTM作成

樹冠を算出するために必要なDTMが必要です、上記にもありますが、DTMとは地表面のことで、表面の高さデータであるDSMでは樹木などが邪魔になります。そのためこの算出には[巨神兵で樹木を焼き払って](https://www.google.com/search?q=%E7%81%AB%E3%81%AE%E4%B8%83%E6%97%A5%E9%96%93&sxsrf=ACYBGNR1-aMI7LjFy8kg455624Ofbw3UWQ:1577580884704&source=lnms&tbm=isch&sa=X&ved=2ahUKEwjNyvLh0tnmAhUAyosBHd4XDX4Q_AUoAXoECA8QAw&biw=1680&bih=853)測量し直したり、GISで決められたセル範囲内で中央値などを使ったり、LiDAR生データからファースト/ラストパルスを出したり、機械学習などを用いる方法などいろいろあるようです。

![](https://i.imgur.com/DdYrc9z.png)|
|:-:|
出典：[航空レーザ測量の仕組み（国土地理院）](https://www.gsi.go.jp/kankyochiri/Laser_senmon.html)

地域や建造物・樹種による調整など凝り始めたらLiDAR沼にハマるようなので、今回は比較的シンプルな手法でFUSIONから算出します。
沼にどっぷり浸かりたい方や論文執筆のためにアルゴリズムのリファレンスが必要な方は[FUSIONのドキュメント](http://forsys.cfr.washington.edu/Software/FUSION/FUSION_manual.pdf)や[GRASSのドキュメント](https://grasswiki.osgeo.org/wiki/LIDAR)なども参考になるかとおもいます。

ではQGISのFUSIONから、Point Cloud analysis→Ground filterを選択します

<img src="https://i.imgur.com/t6L0U9T.png" width=50%>


- Input LAS layer…上記の生txtから作成したlasを指定
- Cellsize for intermediate surface（intermediate surface モデル※１における検索半径)
…とりあえず()`５`
- Create .dtm surface(あとの処理で必要になるDTMの*.dtmファイル（名前ややこしい）)…`✓`つける
- Output ground LAS file…フィルタ済み*.dtmファイルの出力先※２。

<img src="https://i.imgur.com/rKoAjr3.png" width=60%>

`実行`を押すと、単に地表面と推定されるデータを抽出したアウトプットと、さらに平滑化をしてくれた`出力ファイル名_ground_surface.dtm`が出力されているはずです。ビューで比較すると、まあまあ樹木などがいいかんじでフィルタされているぽいですね、巨神兵に外注出さなくてもよさそうです。

![](https://i.imgur.com/y2yAELq.gif)
上…DSM　左下…フィルタDTM　右下…フィルタ＋平滑化DTM

余裕があればでよいのですが、DSMを使って[CSMapmaker](https://www.geospatial.jp/ckan/dataset/cs-tool/resource/889f7a44-9212-4f9e-adaa-cb03c3cbd1cf)でLiDAR結果から下のようなCS立体図も作ることが可能ですし、水文解析もやってみるとより詳細な結果が得られるかもしれませんね。
<img src="https://i.imgur.com/TyedyND.png" width=80%>
出典：地理院タイル

※１ [FUSIONのドキュメント](http://forsys.cfr.washington.edu/Software/FUSION/FUSION_manual.pdf)P93によると
Kraus and Pfeifer(1998) やKraus and Mikhail(1972) から導き出したとのこと
※２ もっとがっちり様々なパラメーターを駆使してDTMをこさえたいLiDAR沼な方は、FUSIONですとSurface->Ground surface createなどがあります。




###  DSMからCHMの算出

樹冠を示すCHM(Canopy Heght Model)を算出してみましょう。CHMはFUSIONの`Surface→Canopy model`ツールを使います。(FUSIONはとりあえずDSMのみからCHMぽいものを算出し、最高樹冠高ポイントを出す際にDTMを用いるようです。なのでツール内では `FUSION canopy height model`と呼んでいます)

![](https://i.imgur.com/ITXkQz5.png)

- Input LAS layer(DSMのlasファイル) …最初に作ったDSM
- Cellsize(出力セルサイズ)…なんとなく最終結果の良かった`0.2`mで…（ここは試行錯誤必要かも）
- XY Units(XY方向の単位)…メートルなので`m`
- Z units(Z方向の単位)…メートルなので`m`
- Add an ASCII output　（ASCIIラスタファイル作るか）…確認用に`✔`
- Output surface（出力ファイル）…CHM計算結果（*.dtm）保存先
<img src="https://i.imgur.com/ODzv5cR.png" width=60%>


これを実行するとCHMの*.dtmが算出されます。同時にASCIIラスタも`出力ファイル名.asc`でフォルダに出力されていますので、このファイルをQGISにドラッグア＆ドロップして陰影表示にすると樹冠のモコモコが優先的に抽出されているはずです。
<img src="https://i.imgur.com/OiAtNEa.png" width=60%>
出典：地理院タイル

###  CHMから樹冠ポイントの抽出

では、このCHMからFUSIONの`CanopyMaxima`というツールを用いて樹冠ポイント(Treetop)を算出します。このツールでは樹冠とみられる高さデータを持ったパッチの最高点を樹冠ポイントとして認識します。
もちろん樹形を考慮すると照葉樹などの頂点は必ずしも樹冠ポイント＝樹木そのものとは言えない場合が多いことに注意が必要です（佐々木2012）。
この算出に関する詳細なアルゴリズムは、[FUSIONのドキュメント](http://forsys.cfr.washington.edu/Software/FUSION/FUSION_manual.pdf)P26によるとKraus and Pfeifer(1998) やKraus and Mikhail(1972) から得られています。

ではFUSIONの`Point cloud analysis→Canopy maxima`ツールを選択します。

![](https://i.imgur.com/dqGYvA3.png)

設定画面ではとりあえずこのように埋めます

- Input FUSION canopy height model …上記で作ったCHM
- Input ground .dtm layer…上記で作ったDTMの*.dtm(ややこしい)
- Height threshold(最低の高さの閾値)…今回は5mとしました
- Output…出力先とファイル名。拡張子dtmが付与されるが結果はcsvで出てくる。
<img src="https://i.imgur.com/lczeZ0u.png" width=60%>

その他に関しては詳細な算出アルゴリズムパラメーターのようですが、今回はデフォルトで実行してみます。するとこのようなCSVが指定フォルダに生成されます。
<img src="https://i.imgur.com/JnShI88.png" width=60%>


XY座標列のあるCSVファイルですので、QGISにちょっと慣れてる方ならすぐ地図化できますよね。また`Height`列もあるので、高さごとに色なども変えれそうです。

*XY座標情報の入ったCSVからQGISに表示させるやり方は[このあたりをご参照ください](https://groups.google.com/forum/#!topic/qgisshitumon01/YJrYZ-mraMs)

さすればQGISで可視化をやってみましょう（今回のメイン投影法はEPSG:2444ですのでお間違えなく）
<img src="https://i.imgur.com/aO9tBSh.png" width=80%>
出典：地理院タイル

※max422.9ｍの世界最大級の高木がでているのは、lasデータ周囲のエラーですのでお気になさらずに。

このポイントからボロノイ多角形なども発生させてみると、そう！森林なGISの研究パンフなどに掲載されてるアレが再現できます（この処理する意味あんまりわかってないのですが）
<img src="https://i.imgur.com/wSfqc2i.png" width=60%>

これは林班図SHPデータなどと重ねて、植林などの情報とともにより詳細な検証ができるかもしれない！と思い、熊本県のオープンデータサイトに飛んで「[くまもとの林業」をクリックすると…](https://www.pref.kumamoto.jp/hpkiji/pub/List.aspx?c_id=3&class_set_id=1&class_id=5501)

<img src="https://i.imgur.com/4LBYZan.png" width=60%>

**くまモンと学ぶ熊本の森林・林業・木材に関する副読本**　ぐらいしかありませんでしたorz
これでは、最後の検証ができません。

ガッチリ検証をされたい方は、林班や毎木調査結果などフィールドデータやオープンデータの揃っている場所や地域で遂行されることをおすすめします。

## もっと学びたい・トライしたい方に

最後がなんかアレでしたが、本記事はほんのトライアルですので、パラメータも適当です。手順も間違ってたらごめんなさい。

もしこの雑い記事でQGISとオープンデータを使ってLiDARによる森林解析をもっと極めてみたいと思われたかた向けに、ネットの海から拾い上げたいくつかの知恵を共有したいと思います。

たとえば、[YouTubeでもLiDAR解析だけみると1600本ものレクチャー動画があります。](https://www.google.co.jp/search?q=Lidar+CHM+DSM&sxsrf=ACYBGNQDpc3ykcZrnQeGmOt014H_F63U3g:1577351731967&source=lnms&tbm=vid&sa=X&ved=2ahUKEwiAwK2N_dLmAhVqyYsBHZgsD_8Q_AUoBHoECAsQBg&biw=1922&bih=820&dpr=0.9)、また、Google Scholarで[Lidar CHM DSM Opensource]をあたると[過去五年で250本ほどの論文](https://scholar.google.co.jp/scholar?as_ylo=2015&q=Lidar+CHM+DSM+opensource&hl=ja&as_sdt=0,5)がでています。

これらの全部の動画観て論文読め！というのも何ですので、いくつかピックアップしますと、操作解説動画はGeoRGBアカウントによる[Lesson 10.1. LiDAR LAStools. 3D model. Hillshade. SASplanet.](https://www.youtube.com/watch?v=T4z6gFeDDbc)はQGIS＋Lastoolsのデモが見れて参考になりました。設定で自動翻訳字幕もできます（比較的はっきり話されてるので翻訳もそんな変でないです）。
<img src="https://i.imgur.com/TvaQMld.png" width=60%>

日本語のドキュメントなら

- [RのForestToolsライブラリにより樹頂点を検出しQGIS上で表示する](https://qiita.com/nireiro_py/items/9fcb7d869b5c515f74c6) by @nireiro_py さん
- [国土地理院(2015)航空レーザ測量データを用いた樹高等のデータ作成](https://www.gsi.go.jp/chirijoho/chirijoho40069.html)
- [森林GISフォーラム 東京シンポジウム資料](https://fgis.jp/wordpress/wp-content/uploads/2018/12/20190220_FGIS_PPT.pdf)
- [佐々木剛 航空機LiDARを用いた森林構造の推定 景観生態学 17.2 (2012): 43-55.](https://ci.nii.ac.jp/naid/130003377443/)

あたりが林学の授業もなにも取ってなかったわたしにとって、どれもイメージや処理手順を掴むのにとてもよかったです。ありがとうございます。
森林GISフォーラム資料の府立大田中先生資料の後半ではオープン化されたLiDARでの活用についての提言がされていますし、佐々木剛さんのレビューに関してはLiDARの基礎から近年の手法まで丁寧に解説がされており、現場でちょっと知ってる方にもおすすめします。

## まとめ

長々とお読みくださり誠にありがとうございました。気がついたら2万字こえててまして、いちど人に見てもらったところ、`オカンがギチギチに編み込んだずっしりしたセーターみたいだね`といかようにでも解釈できるありがたいコメントをいただきました。

 さて、今回は

- LiDARのデータはQGISでもいじれる
- 樹冠・樹高もQGISで出せた（かも）
- 検証をすすめたらさまざまな現場でも活用できそう

といった感触を得ました。

上でも書きましたが、植林されているエリアは比較的データのノイズが少ないでしょうから、より高精度の処理や抽出が可能でしょう。

また、今回はLiDARからの処理でしたが、ドローンなどの写真からSfM処理を行って作成されたDSMでもこのような抽出などもできるかもれませんね。

文頭でも述べましたがQGISを使った3Dデータの可視化や処理の事例はほとんどありません。このような情報が国内でもっと増えると、様々な分野でより効率化や高収益化がすすむかもしれませんし、オープンデータとしての3Dデータの価値も上がるでしょう。今後いい検証や事例ができたらぜひシェアしていただけると嬉しいです🌲


## この記事のライセンス

[![クリエイティブ・コモンズ・ライセンス](https://i.creativecommons.org/l/by/4.0/88x31.png)](http://creativecommons.org/licenses/by/4.0/deed.ja)
この記事は[CC BY 4.0（クリエイティブ・コモンズ 表示 4.0 国際 ライセンス）](http://creativecommons.org/licenses/by/4.0/deed.ja)の元で公開します、どうぞどうぞ。