0. はじめに

ある土地を計測した三次元点群データの利用方法として、地表面の高さを表現したDigital Surface Model (DSM)を作成するか、標高を表現したDegital Elevation Model (DEM)を作成するなどがあります。DSMは点群の高さを用いることで作成できますが、DEMは地面 (Ground)とそれ以外の点群に分ける必要があります。これをGround Filterといいます。

Ground Filterにはいくつかの方法がありますが、その方法をまとめたソースはあまりないです。このサイトでは、Ground Filterができる私の知る限りの各ツールの手法をまとめております。
　
※あくまでも、手法のまとめとコードの解説が目的であるため、アルゴリズムの内容については参考になるサイトを記載するのみにします。

使用するデータセット

本サイトでは、三次元点群データデータとして静岡県のG空間情報センターからランダムに得た土地のデータを使用します。

実行環境

R: 4.40
CloudCompare: 2.13.1
QGIS: 3.36.3
python: 3.11
VS Code: 1.93.1
Anaconda: 2024.06
pdal: 3.2.0

1. Rを使った方法

〇メリット
・出力結果に自動で地面ラベルが付与される。RでDEMを作成するためには、地面ラベルが必要であるため、手間が解消される
・扱いやすい
・自作のアルゴリズムを導入できる

〇デメリット
・処理時間が長い。1千万点を超えると、半日以上かかる

初めに、必要な準備をします。以下のコードでデータを読み込みます。

library(lidR) #使用するライブラリ
setwd(Directory_Path) #読み込みたいlasファイルがあるディレクトリのパス
las = readLAS(File_Name.las) #読み込みたいlasファイルの名前

1.1 Progressive Morphological Filter (PMF)

pmf = classify_ground(las, algorithm = pmf(ws = 5, th = 3))
pmf_ground = filter_ground(pmf) #地面点群の抽出
plot(pmf_ground) #可視化

1.2 Cloth Simulation Function (CSF)
CSFはRにインストールされていない可能性があるため、インストールのコマンドを記載しています。

install.packages("RCSF") #インストールコマンド。一度すれば、以降は不要
csf = classify_ground(las, algorithm = csf())
csf_ground = filter_ground(csf)

1.3 Multiscale Curvature Classification (MCC)
MCCはRにインストールされていない可能性があるため、インストールのコマンドを記載しています。

install.packages("RMCC") #インストールコマンド。一度すれば、以降は不要
mcc = classify_ground(las, mcc(1.5,0.3))
mcc_ground = filter_ground(mcc)

保存方法

writeLAS(pmf_ground, File_Name.las) #保存したいファイル名

引用
lidR 3 Ground Classification

2. CloudCompareを使用した方法

CloudCompareで使用できるアルゴリズムは、CSFです。Rでも実装されています。
〇メリット
・クリックだけで使用できる
・地形によってパラメータを変更できる。選択が簡単。

〇デメリット
・精度があまりよくない可能性がある (森林地帯の点群データで試すと、満足な出力にならないことが多い)

使い方は、ウィンドウ上の[Plugins]→[CSF]で使用できます。選択すると以下のウィンドウが出てきます。

計測した地形を選択することで、自動で必要なパラメータが決定します。

3. QGISを使用した方法

〇メリット
・パラメータを設定しやすい
・処理が速い
・出力結果を見て、修正しやすい

〇デメリット
・良いパラメータを見つけるまでに時間がかかる

QGISは、プラグインのFUSIONをインストールすることで点群データを処理できるようになります。
FUSIONの実装と使用については、以下のYfuruchin様のページで詳細に記されているため、参照してください。
https://qiita.com/Yfuruchin/items/7b9ea3e823824d2e4e86#fusion%E3%81%AE%E3%82%BB%E3%83%83%E3%83%88%E3%82%A2%E3%83%83%E3%83%97

実装した後、Processing Toolsの中にFUSIONが追加されます。
[FUSION]→[Point Cloud Analysis]→[Ground Filter]を選択することでウィンドウを開けます。
パラメータを変更しながら、良い出力になるパラメータを探索します。
[Cellsize for intermediate surface]により、Ground Filterのパラメータを変更できます。

4. pythonを使用した方法

〇メリット
・コードの組み方次第で完全自動化しやすい
・気合を入れれば、自分でもアルゴリズムを作れる

〇デメリット
・少し遅い

pythonで使用できるのは、私が確認した中ではRやCloudCompareで使用できるCSFでした。
しかし、pythonはコードが組みやすいので効率化を図るにはもってこいです。
下記のいずれかでCSFのライブラリをインストールできます。

py -m pip install cloth-simulation-filter
pip install cloth-simulation-filter

加えて、lasファイルを扱うためのライブラリも必要となります。

py -m pip install laspy

次のコードで実行できますが、引用をコピペしているだけです。詳しく見たい場合は、引用にあるURLを閲覧してください。

# coding: utf-8
import laspy
import CSF
import numpy as np

inFile = laspy.read(r"in.las") # read a las file
points = inFile.points
xyz = np.vstack((inFile.x, inFile.y, inFile.z)).transpose() # extract x, y, z and put into a list

csf = CSF.CSF()

# prameter settings
csf.params.bSloopSmooth = False
csf.params.cloth_resolution = 0.5
# more details about parameter: http://ramm.bnu.edu.cn/projects/CSF/download/

csf.setPointCloud(xyz)
ground = CSF.VecInt()  # a list to indicate the index of ground points after calculation
non_ground = CSF.VecInt() # a list to indicate the index of non-ground points after calculation
csf.do_filtering(ground, non_ground) # do actual filtering.

outFile = laspy.LasData(inFile.header)
outFile.points = points[np.array(ground)] # extract ground points, and save it to a las file.
out_file.write(r"out.las")

ただし、もしかしたら出力するためのコードが動かない可能性はあります。
そうなった場合は、以下のコードを参考にしてください。

new_las=laspy.create(point_format=inFile.header.point_format)
new_las.header.scales=inFile.header.scales
new_las.header.vlrs=inFile.header.vlrs
new_las.header.offsets=inFile.header.offsets
new_las.header.add_crs(inFile.header.parse_crs())
new_las.points=points[np.array(ground)]
new_las.classification=np.full(len(new_las.X),2) # 地面点群にラベルづけ。慣例的に2を与える。
new_las.write(r"out.las")

引用
https://github.com/jianboqi/CSF

5. pdalを使用した方法

〇メリット
・精度が高い

〇デメリット
・ユーザーが少なく、情報が少ない

pdalは点群データを処理するためのライブラリで、pythonで使用可能です (gdalの点群データ用のイメージ)。python上でも動かせますが、anaconda (仮想環境)での構築が簡単です。

pdalのインストールは、以下のコマンドをanaconda promptに入力すれば完了します。
※anaconda promptはWindowsボタンを押して、anacondaと打てば出てきます。

conda create -n new_env python=3.9 #新たな仮想環境の構築 (new_envは仮想環境名)
conda activate new_env #仮想環境に移る
conda install -c conda-forge python-pdal #仮想環境にpdalをインストール

実行するときにインストールの許可を求めてきますが、Enterを押しておけば問題ないです。
pdalの実行方法には二種類あります。

jsonファイルを作成する方法
pythonのコード内に組み込む方法

今回は1を試します。

1. jsonファイルを作成する方法
こちらの方法ではjsonファイルに実行コマンドを記載し、anaconda promptからファイルを実行するという手順です。
まず、以下がjsonファイルの中身です。
#以降は解説です。実際に使用する際は削除してください。

{
    "pipeline":[
        "C:/Users/User/Data/PointCloud.las", #実行対象のlasファイルのパス
    {
    "type":"filters.assign", #GroundFilterの対象とするラベルを選択
    "assignment":"Classification[:]=0" #今回は0のラベルがつけられている点群を対象とする
    },
    {
    "type":"filters.elm" #ELMというアルゴリズムにより、ノイズ点群を見つけます。7のラベルが与えられます
    },
    {
    "type":"filters.outlier" #極端に離れている点をノイズとします。7のラベルが与えられます
    },
    {
    "type":"filters.smrf", #SMRFというアルゴリズムにより、地面点群を抽出します。
    "returns":"last", #last returnの点を対象に
    "ignore":"Classification[7:7]", #7がつけられた点群を無視して、
    "slope":0.2, #以下のパラメータに従って分類されます (公式サイトのデフォルト値を使用)
    "window":16,
    "threshold":0.45,
    "scalar":1.2
    },
    {
    "type":"filters.range", #抽出した点群に2のラベルを与える
    "limits":"Classification[2:2]"
    },
    "C:/Users/User/Data/Ground.las" # 出力するファイルのパス
    ]
}

上記のファイルをanaconda promptにより実行する。
仮に、jsonファイル (ground_filter.jsonというファイル名)がDownloadのファイルにあるとした場合

conda activate new_env
cd Downloads #ファイルが存在するディレクトリに移動
pdal pipeline ground_filter.json #実行コマンド

により、ground filterが実行されます。

引用
https://pdal.io/en/latest/tutorial/ground-filters.html

おまけ

ここまで紹介した方法により抽出した地面点群からDEMを作成することで、より使いやすい情報にできます。
Rによる作成方法とpdalによる作成方法をご紹介いたします。

DEMの作成 with R

library(lidR) #lasファイル用のライブラリ
library(raster) #ラスタ用のライブラリ
setwd("C:/Users/User/directory_of_your_file")
las=readLAS("Ground.las") #2のラベルを持っている点群でなければならない
#resはDEMの解像度を指定
#algorithmは、DEM作成に使用するアルゴリズムを指定
#tin, knnidw, krigingなどがある
dem = rasterize_terrain(las, res = 1, algorithm = tin()) #DEMの作成
writeRaster(dem, "DEM.tif", overwrite=TRUE) #overwrite=Trueで同名ファイルを上書きできる

引用
https://r-lidar.github.io/lidRbook/dtm.html#:~:text=Generating%20a%20Digital%20Terrain%20Model%20(DTM)%20is%20usually%20the%20second

DEMの作成 with pdal
pdalを用いた場合、単純に点群をラスタに変換するため、後処理が必須になります。
そういった面では、RのほうがきれいなDEMを作成できるかもしれません。
pdal公式から後処理が簡単に紹介されておりますので、参考にしてください。
※紹介された方法以外では、pythonでrasterioやQGISを用いることで効率的に後処理ができます。

{
    "pipeline": [
		"C:/Users/User/directory_of_your_file/Ground.las",
        {
            "filename":"C:/Users/User/directory_of_your_file/dem.tif",
            "gdaldriver":"GTiff",
            "output_type":"all",
            "resolution":"0.5",
            "type": "writers.gdal"
        }
    ]
}

引用
https://pdal.io/en/2.4.3/workshop/exercises/analysis/dtm/dtm.html#:~:text=gdaldem,%20which%20powers%20the%20QGIS%20DEM%20tools,%20is%20a%20very

地面ラベルの付与

Rとpdal以外の手法で地面点群を取得した場合、地面ラベルが付与されていないため、Rを使用したDEMの作成ができません。
そこで、以下にpythonで簡単に地面ラベルを付与できるプログラムを記載します。
使用するlaspyはデフォルトでインストールされていないため、以下でインストールします。

py -m pip install laspy

import laspy
import os
import numpy as np
path="Directory_Path" #lasファイルが保存されているディレクトリのパス
files=[f for f in os.listdir(path) if f.endswith(".las")] #ディレクトリ内のlasファイルのみを選択
for f in files:
    las=laspy.read(path+"/"+f) #laspyによるファイルの読み込み
    las.classification=np.full(len(las.X), 2) #lasファイルに地面ラベル(2)を与える
    las.write(path+"/"+f) #読み込んだファイルを上書き保存

まとめ

本ページでは、無料で知っている限りの方法のみに焦点を当てました。
有料のもの (例えば、DF LAT)などもございます。学会で紹介いただいた際には、かなり精度が高いと報告を頂きました。
使用する地形条件や点群密度によって、適している方法も異なっております。
色々試して、精度が良さそうなものまたは、加工しやすいものを選択する手助けになるといいと思います。
もし、ご指摘やご質問、新たな手法がございましたらコメントを頂けたら幸いです。