みなさん、こんにちは！
PDAL使ってますか！！！
PDALは、点群データを操作するための強力なCLIツールかつ、ライブラリです。
と言うのはこちらの記事でいろいろ書いていますので、ぜひご覧ください。
で、今回はこのPDALをPythonから呼び出して利用してみようと思います。
PDALはC++で書かれたライブラリですが、Pythonからも利用できるようにラッパーが用意されています。
PDALにはパイプラインという概念があり、JSON形式でルールに従って記載することで、点群データの一連の処理をまとめて実行することができます。
このため、CLIからJSONを指定して実行することが多いのですが、Pythonからも同様にJSONを指定して実行することができます。
Pythonから利用できると、例えばWebアプリケーションからPDALを利用したり、Pythonのライブラリ（numpyや今回利用するOpen3Dなど）と組み合わせて利用することが容易になります。
この記事ではPDALでDEMを作成する・地形のglTFを作成する、といった利用方法について解説していきます。
この記事のコードはこちらに載せているので、参考にしてください！
データの準備
今回は、LAS形式のデータを変換します！
どんなデータを利用してもいいのですが、今回はこちらの記事で紹介した東京都の点群データを利用します。
どこでもいいのですが、適当な場所のデータをダウンロードしておいてください。
必要なパッケージのインストール
uvの利用を前提とするので、uvがインストールされていない場合は、uvをインストールしてください。
その後、必要なパッケージをインストールします。
uv add pdal open3d rasterio matplotlib jupyter
こんな感じのpyproject.tomlができていればOKです。
[project]
name = "pointcloud-processing-sample"
version = "0.1.0"
description = "Add your description here"
readme = "README.md"
requires-python = ">=3.12"
dependencies = [
"jupyter>=1.1.1",
"matplotlib>=3.10.1",
"open3d>=0.19.0",
"pdal>=3.4.5",
"rasterio>=1.4.3",
]
DEMの作成
PDALを利用して、DEMを作成してみます。
jupyter labを利用しますが、ただのPythonスクリプトでも問題ありません。
まずは必要なライブラリをインポートします。
import json
import pdal
import rasterio
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource
ファイルのパスを指定しますが、ダウンロードしたパスに合わせて変更してください。
input_path = "./data/input/09LD1895.las"
output_path = "./data/output/09LD1895_dem.tif"
定義するパイプラインはこの程度の記述でOKです。
今回はPythonで書いていますが、同様の内容をJSON形式で書けばCLIから実行することもできます。
pipeline_json = {
"pipeline": [
{
"type": "readers.las",
"filename": input_path
},
{
"type": "filters.smrf",
"window": 50, # 除去したいオブジェクトの大きさ（m）
"slope": 0.2, # 傾斜の閾値
"threshold": 0.45, # 高度の閾値（m）
"cell": 1 # 解像度（m）
},
{
"type": "filters.range",
"limits": "Classification[2:2]"
},
{
"type": "writers.gdal",
"filename": output_path,
"output_type": "mean",
"resolution": 1.0,
"radius": 1.0,
"nodata": -9999
}
]
}
pipeline_definition = json.dumps(pipeline_json)
このパイプラインは、以下の処理を行います。
-
readers.lasでLAS形式のデータを読み込みます -
filters.smrfで地面を抽出します -
filters.rangeで地面の点群を抽出します -
writers.gdalでDEMを出力します
パラメータは、デフォルト値を参考にしながら調整してください。
特に、地表面抽出のsmrfは目的によってシビアかつパラメーターが多いので、色々試してみてください！
ではこのパイプラインを実行してみましょう。
pipeline = pdal.Pipeline(pipeline_definition)
pipeline.execute()
実行が終わったら、出力されたDEMを確認してみましょう。
with rasterio.open(output_path) as src:
dem_data = src.read(1, masked=True) # nodataをマスクして読み込む
transform = src.transform
するとこんな感じのDEMができているはずです。
QGISでも確認してみましょう。
いい感じですね！
地形のglTFの作成
利用用途はさほどないんですが、こんなこともできるよーという紹介として、地形のglTFを作成してみます。
まずは必要なライブラリをインポートします。
import json
import pdal
import open3d as o3d
ファイルのパスを指定しますが、ダウンロードしたパスに合わせて変更してください。
pointcloud_path = "./data/input/09LD1895.las"
ply_path = "./data/output/09LD1895_mesh_terrain.ply"
glb_path = "./data/output/09LD1895_mesh_terrain.glb"
パイプラインは以下のように記述します。
mesh_pipeline = {
"pipeline": [
{
"type": "readers.las",
"filename": pointcloud_path
},
{
"type": "filters.smrf",
"window": 50,
"slope": 0.2,
"threshold": 0.45,
"cell": 1
},
{
"type": "filters.range",
"limits": "Classification[2:2]"
},
{
"type": "filters.poisson",
"depth": 8,
},
{
"type": "writers.ply",
"filename": ply_path,
"faces": True
}
]
}
pipeline_obj = pdal.Pipeline(json.dumps(mesh_pipeline))
このパイプラインは、以下の処理を行います。
-
readers.lasでLAS形式のデータを読み込みます。 -
filters.smrfで地面を抽出します。 -
filters.rangeで地面の点群を抽出します。 -
filters.poissonでメッシュを作成します。 -
writers.plyでply形式のメッシュを出力します。
パラメータは、デフォルト値を参考にしながら調整してください。
以下のように実行します。
pipeline_obj.execute()
これでply形式のメッシュができました。
このファイルをOpen3DをglTF形式に変換してみます。
mesh_terrain = o3d.io.read_triangle_mesh(ply_path)
o3d.io.write_triangle_mesh(glb_path, mesh_terrain)
これでglTF形式のメッシュができました。
glTF Viewerで確認してみましょう。
おわりに
今回はPDALをPythonから利用して、DEMを作成したり、地形のglTFを作成したりしてみました。
PDALはCLIからも利用できるので、Pythonから利用する必要がない場合はCLIから実行しても問題ありません。
PDALは非常に強力なツールなので、ぜひ活用してみてください！