はじめに

みなさん、こんにちは。研究開発室の大道です。
前回のブログではQGISで航空写真と数値樹冠高モデルを重ねて表示しましたが、今回はいよいよ数値樹冠高モデルから樹頂点を抽出していきます。
全3パートのうち今回はPart2になります。

• Part1: 航空写真と数値樹冠高モデルを重ねて表示してみよう
• Part2: プロセシングを使って数値樹冠高モデルから樹頂点を抽出してみよう
• Part3: モデルデザイナーを使って樹頂点抽出処理を自動化してみよう

これ以降、便宜上、数値樹冠高モデルを樹冠高モデルと呼称します。

想定読者

• QGISを少し触ったことのある方
• 地理空間情報のデータ活用に興味のある方

実行環境

• Windows10
• QGIS 3.34.6

プロセシングを使って数値樹冠高モデルから樹頂点を抽出してみよう

1. 樹頂点の抽出方法を考える

基本的な考え方は画像処理で樹冠高モデルのピークを抽出することです。
国土地理院の樹木疎密度データの作り方では、3×3の最大値フィルタを用意して、最大値を計算した結果とDCM¹との差が0になる点をピークとしています。

上図の赤枠部分が樹頂点の抽出フローですので、今回はこの手続きを参考にしてQGISで処理していきます。

2. プロセシングを用いた樹頂点抽出処理

2.1 プロセシングとは？

では早速、QGISで樹頂点を抽出していきますが、その前にQGISのプロセシングについて少し説明します。
プロセシングは地理空間情報を処理するアルゴリズムがまとめられているフレームワークです。
GDALやGRASS GISなどのサードパーティのアルゴリズムも呼び出すことが可能です。
基本的にはプロセシングツールボックスから必要なアルゴリズムを選択して処理を行います。

他にも、複数のアルゴリズムをグラフィカルにつなげて処理を行うモデルデザイナー（後日Part3で紹介予定）やバッチ処理なども可能です。

2.2 プロセシングを使う

先程のプロセシングを実際に使ってみましょう。
QGISを起動して、前回と同じ状態からスタートします。

2.2.1 樹冠高モデルのピークを抽出する

まずは、樹冠高モデルの範囲が大きすぎるので一部分だけ切り抜きましょう。
［プロセシングツールボックス］から［GDAL］⇒［ラスタ抽出］⇒［範囲を指定して切り抜く］をダブルクリックします。
（プロセシングツールボックスがない場合は、属性ツールバーの歯車マークをオンにすると表示されます。）

• 入力レイヤ：樹冠高モデル
• 切り抜く範囲：15238608.5725,15241034.3522,4187900.5975,4190076.1989 [EPSG:3857]
• このnodata値を出力バンドに割り当てる［オプション］：0

切り抜く範囲については任意に決めてもよいです。
実行をクリックすると、名古屋城周辺の範囲に切り抜かれた結果が出力されます。
また、高さが0mの所はnodataになります。

次は樹冠高モデルのピークを抽出するために、最大値フィルタを用意して樹冠高モデルをフィルタリングしましょう。
［プロセシングツールボックス］から［GRASS］⇒［Raster］⇒［r.neighbors］をダブルクリックします。

• Input raster layer：切り抜いた樹冠高モデル
• Neighborhood operation［オプション］：maximum
• Neighborhood size(must be odd)［オプション］：7

実行をクリックして、最大値フィルタの計算結果を出力します。
続けて、最大値フィルタの計算結果と樹冠高モデルの差が0となる領域を抽出します。
［プロセシングツールボックス］から［GDAL］⇒［ラスタその他］⇒［ラスタ計算機］をダブルクリックします。

• 入力レイヤ A：最大値フィルタの計算結果
• Aのラスタバンド数：バンド1（Palette）
• 入力レイヤ B［オプション］：切り抜いた樹冠高モデル
• Bのラスタバンド数［オプション］：バンド1（Gray）
• gdalnumeric構文で計算。：A==B
• 出力のnodata値を設定［オプション］：0

gdalnumeric構文の「A==B」は、最大値フィルタの計算結果と切り抜いた樹冠高モデルのピクセルの値が一致しているときは1、そうでなければ0を返します。
それと出力結果が0の場合はnodataに設定するようにします。
実行をクリックすると、樹冠高モデルのピークが黒い領域で出力されます。

2.2.2 ピーク領域から樹頂点を抽出する

樹冠高モデルのピーク領域が取り出せたのでここから樹頂点を抽出しましょう。
［プロセシングツールボックス］から［GDAL］⇒［ラスタ変換］⇒［ラスタをベクタ化（polygonize）］をダブルクリックします。

• 入力レイヤ：樹冠高モデルのピーク領域

実行をクリックして、ラスタをベクタに変換します。
続けて、ベクタ化されたピーク領域から重心を計算して樹頂点を抽出します。
［プロセシングツールボックス］から［ベクタジオメトリ］⇒［重心］をダブルクリックします。

• 入力レイヤ：ベクタ化されたピーク領域

実行をクリックすると、ピーク領域の重心、つまり樹頂点が表示されます。

2.2.3 5m以上の樹頂点を抽出する

2.2.2で取り出した樹頂点は高さ情報を考慮していないので、かなり粗い結果です。
そのため、より適切なデータを得るために5m以上の樹頂点を抽出します。
［プロセシングツールボックス］から［GDAL］⇒［ラスタ解析］⇒［ベクタレイヤにラスタ値を付加］をダブルクリックします。

• 入力レイヤ：樹頂点
• ラスタレイヤ：切り抜いた樹冠高モデル
• ラスタ値を収納するカラム名の接頭辞［オプション］：height_

実行をクリックして、樹頂点に樹冠高モデルの高さ情報を付加します。
続けて、この高さ情報から5m以上の樹頂点を抽出します。
［プロセシングツールボックス］から［ベクタ選択］⇒［属性による抽出］をダブルクリックします。

• 入力レイヤ：高さ情報が付加された樹頂点
• 選択基準になる属性：height_1
• 演算子：≥
• 値［オプション］：5

実行をクリックすると、5m以上の樹頂点のみが出力されます。
以上が樹冠高モデルから樹頂点を抽出する処理の手続きになります。

3. 課題

国土地理院の手続きを参考にしながら樹頂点を抽出しましたが、この方法だけだと以下の課題があります。

1. 樹頂点が車道上などの不適切な場所にある
2. 樹木が密である領域のピークが適切でない

まず1つ目の課題ですが、抽出した樹頂点の位置にはズレがあることです。
この樹冠高モデルは航空写真と比べて位置がシフトしている所があり、樹頂点をそのまま使うと車道上に樹木が配置されてしまいます。

2つ目の課題は主に樹木が密集している場所では、樹木1本1本のピークを取り出せていないことです。
そもそも今回使用している樹冠高モデルはあくまで深層学習の推定結果なので、高さの分布が真とは限りません。
ですので、画像処理による抽出方法だと限界があります。

これら2つの課題に対するアプローチですが、手前味噌になりますがCEDEC+KYUSHU2024の講演資料で位置補正の方法やサンプリングの仕方について触れていますので、よろしければご参照ください。

おわりに

今回はQGISのプロセシングを用いて樹冠高モデルから樹頂点を抽出することができました。
様々なプロセシングを利用して簡単に地理空間情報を処理出来ることがご理解頂ければ幸いです。
ただ同時に、複数のプロセシングを1つ1つ操作するのはとても大変だったかと思います。
ですので次回は、これらの手続きをグラフィカルにつなげて自動化するモデルデザイナーをご紹介する予定です。
お楽しみに！

1. Digital Canopy Modelの略称。樹冠高モデルのことである。 ↩︎