FOSS4G Advent Calendar 2021 12日目の記事です。

追記
本記事はQGIS バージョン3.14を使用して書いております。
本記事で紹介しているSAGAツールはQGIS 3.16以降では使えないようですのでご注意ください。

QGISの地図上にお好きな3Dモデルを生やすことができるんですよ。

いきなりですが、そうなんですよ。知っていましたか？

方法はQgis2threejsプラグインを使います。

QGIS2threejsの使い方

プラグインを実行すると、別ウインドウが立ち上がります。
左パネルのレイヤの一覧から「GSI Elevation Tile」にチェックを入れると立体地形の地図が描画されます。（自身で追加しているDEMのラスタデータにチェックを入れてもOKです）

3Dモデルを配置するポイントデータをレイヤに追加しておき、プラグイン画面でポイントのプロパティを開きます。

設定箇所は以下3点

Object type をModel Fileにする
Z coordinateのModeを「Relative to "GSI Elevation Tile"」（もしくは、自身で用意したDEMラスター）にする
Model fileで配置する3Dモデルのファイルを選択する。（ローカルにあるファイルを選択しても良いし、ネット上のURLを入力しても良い。

以上、設定してOKを押すと

もうサクッと、モデルが配置されるわけなんですね。
これでみなさんも、お好きなモデルを配置してグリグリ動かしてみてください

以上、「QGISで楽しく3D表現してみる」でした。

おわり

ここからは余談

もっと木とかいろんな3Dモデル配置して、実際の風景に近づけたいなぁ

記事の本編は以上になりますが、さすがにこれで終わると味気ないのでもうちょっと遊んでみようと思います。

せっかく3Dモデルを配置できることがわかったので、樹林や建物などいろんな種類のモデルを配置して実際の景色を再現したりできると、楽しそうだと思いませんか？
思いますよね。
うん、思ってきましたね。

例えば本編で紹介した風車のモデルで言うと、障害物がないため風車の足元までしっかり見えていますが、実際には周辺の樹木などによって風車の足元までしっかり見えるということはあまりないかと思います。そのため、このままだと3Dモデルで風車を見せても現実味が薄れてしまいます。

そこで、ここからは樹木などいろいろなモデルを描画させて、さらに現実味の高い3D地図を作成していきたいと思います。

完成イメージとしてこんな感じ

QGISを使って「おうとまちっく」に樹林のポイントを作ってみる。

樹木のモデルを配置するにしても、本編で紹介したようにその地点のポイントデータが必要になります。つまり樹木1本1本の位置にポイントデータを作成する必要があります。

国土地理院などが公開している航空写真は高解像のものも多く、目で見て1本1本の木がわかるものもあります。
航空写真から判読して手作業でポイントを落としていくのもいいのですがそれだと芸がないので、航空写真からある程度自動で判読して、1本1本の木の位置をポイント化ししていきたいと思います。

（作業イメージ）

材料

対象エリアの航空写真（高解像度であればあるほど良い）

樹木のポイントデータを生成するための航空写真を用意します。

今回は国土地理院からタイル配信されているこちらのデータを使います。

令和3年7月1日からの大雨正射画像熱海伊豆山地区（7/6撮影）
URL：https://cyberjapandata.gsi.go.jp/xyz/20210705oame_0706do/{z}/{x}/{y}.png

土地利用データ

木のモデルを生やすのはもちろん樹林地だけにしたいので、土地利用から樹林地かそうでないかの判別を行うために使います。
今回は、JAXAが提供している高解像度土地利用土地被覆図を使用していきます。

(出典：　JAXA 高解像度土地利用土地被覆図)

樹木高さデータ

せっかくなら樹木の高さに合わせてモデルの大きさを変えられたら面白いなと思うので、樹木の地面からの高さ（DHM：Digital Height Model)も算出していきます。

DHMは、建物等の含んだ地物の表層の高さ（DSM）から地盤高の高さ(DTM)を差し引いて求めることができます。
（DHM、DTM、DHMの用語の説明はこちらをご参照ください）

DSMはJAXA 第一宇宙技術部門地球観測研究センター（EORC）より提供されている、ALOS全球数値地表モデルを、
DEMは一般財団法人宇宙システム開発利用推進機構(Japan Space Systems)より提供されているASTER全球3次元地形データ／ASTER GDEMを活用しました。

##調理方法

樹木のポイントデータ作成

タイルを読み込んで、任意の場所の地図を表示する。

xyz tileの読み込み方はこちらをご参照ください。
QGIS3でウェブ地図利用が劇的に変化した

タイルを読み込んだらお好きな場所を表示させます。
広範囲すぎると樹木のポイントデータが大量にできてしまって、いざ3Dモデルを描画しようとした際に時間がかかってしまうため、狭いエリアで試してみましょう。

地図を画像にエクスポートする

好きなエリアを表示したら、画像として保存します。

プロジェクト　＞インポートとエクスポート ＞　　地図を画像にエクスポートを選択します。
エクスポートの設定は特にいじらなくて良いですが、地理参照情報を追加にはチェックを入れておいてください

樹冠の判読

エクスポートした画像の樹冠を判別していきます。
樹冠の判別方法は様々あるかと思いますが、今回は流域解析(watershed)を使って簡易的に行います。

(watershedを使った樹冠判読はこちらなどを参考にしました。)
[Wan Shafrina Wan Mohd Jaafar et al(2018), Improving Individual Tree Crown Delineation and Attributes Estimation of Tropical Forests, Forests 2018,9(12)](Using Airborne LiDAR Datahttps://www.mdpi.com/1999-4907/9/12/759)

WatershedはQGISのプロセッシング（正確にはSAGAツール）にも備わっております。
（プロセッシング > SAGA > Watershed segmentation）

設定は以下の通り

グリッド：先ほどエクスポートした画像
出力：SegmentID
方法：Maxima
結合セグメント：seed to saddle difference　
閾値：5

接合セグメントと閾値は、樹冠と樹冠の間をどの様に設定するかのパラメータです。結合セグメントはseed(樹冠のトップ)と樹間のヘリのピクセル値の差が閾値で設定した値より大きいところを樹冠の境界にする。というイメージです。
閾値が大きいと１つ１つの樹冠が大きく判読され、小さいと細かい樹冠として判読されます。

処理結果にはセグメント、seed point、境界の３つが出力されますが、必要なのはセグメントのみです。
セグメント：樹冠ごとの領域を示したラスターデータ
seed：それぞれの樹冠ポイントデータ
境界：樹間の境界を示したラスターデータ

ポイントデータ(seed)も出力されますが、このポイントの位置は、閾値などの設定が反映されていないので、segmentから改めてポイントデータを作ります。

ポイントデータの生成

それぞれ以下のプロセッシングで、segmentラスターデータをポリゴン化し、そのポリゴンの重心にポイントを作ります。
segmentラスターデータをポリゴン化：ラスタのベクタ化
ポリゴンの重心にポイントを生成：重心

結果、こんな感じで、なんとなく1本1本の木にポイントが落ちて、（なんならソーラーパネルにもポイントができて）いい感じになりました。

樹高を求める

ラスタ計算機でDSM-DTMを計算してDHMを求めます。

処理結果はこんな感じ。これだけみると何が何だが。。

ベクタレイヤにラスタ値を付与ポイントデータの属性に、DHMの値を付与します。

ポイントデータに土地利用の情報を付与する

JAXAの高解像度土地利用土地被覆図を追加します。
ラスタには以下の凡例にあるように土地利用のコードが記録されています。

DHMと同様に「ベクタレイヤにラスタ値を付与」を使って、今度はこのコードの値をポイントデータの属性に付与します。

(ポイントデータの属性。heiht_1がDHMの値、landuse_1が土地利用のコード)

ポイントデータを土地利用ごとに色分けしたらこんな感じ

ポイントデータの属性に3DモデルのURLを付与する。

それぞれのポイントの地点に土地利用に適した3Dモデルを描画するために、属性テーブルを１列追加し、土地利用に合わせたモデルへのリンクを入力していきます。

3Dモデルは自身で用意してもいいですが、今回はこちらに公開されているモデルを使用していきます。
https://market.pmnd.rs/

こちらのサイトで公開されている3Dモデルはダウンロードして使用することもできるし、モデルへのダイレクトリンクも取得できるのでリンクにアクセスして使用することもできます。
（利用の際はそれぞれのモデルのCCライセンスに従ってください。）

例えばこちらのTree Beech

広葉樹林（land useが6か8のポイント）のポイントには、このモデルを生やしたいと思います。
Copy direct linkをクリックして、このモデルのURLをコピーします。

landuseが6と8は広葉樹林なので、3DモデルURLのカラムにコピーしたTree Beechへのリンク入力します。
（属性による選択　→ フィールド計算機）

同様に針葉樹林としてTree spruceのモデルのリンクをコピーして、landuseが7と９（針葉樹林）の地物にURLを入力します。

３D表現してみる

とりあえず広葉樹林と針葉樹林の用意ができたので、3D描画してみたいと思います。
Qgis2threejsを起動し、ポイインとデータのプロパティで以下の項目を設定します。　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　

Object type :　Model File
Z coordinate　Mode : Relative to "GSI Elevation Tile"
Model file ：Model path or URLにし、属性テーブルから３DモデルのURLのカラムを指定する。
Scale：樹高などの数値によってモデルの大きさを変えるときにはこちらをいじります。ひとまずデフォルトの１のままやってみます。