PythonからPDALを呼び出してDEM作成

はじめに

PDALと呼ばれる航空機LiDARで取得した点群向けのソフトウェアがあります。このPDALはコマンドライン上で引数を渡して点群処理を行うのが一般的な使用方法ですが、Pythonから呼び出せます。
今回は、PythonからPDALを呼び出して点群からデジタル標高モデル(DEM)を作成する方法を共有します。

コマンドライン上での処理方法

PDALはjsonに処理対象のファイルや、点群処理の内容を記述するのが一般的なやり方です。
今回は、点群から地面の点群だけを抽出してDEMを作成します。
この場合のJSONは以下のように作成します。

pipeline.json

{
    "pipeline": [
        "./08ME3552.las",
        {
            "type":"filters.expression",
            "expression":"Classification == 2"
        },
        {
            "filename":"./08ME3552.tif",
            "gdaldriver":"GTiff",
            "output_type":"mean",
            "resolution":"1.0",
            "type": "writers.gdal"
        }
    ]
}

このJSONの意味は、Classificationが2(地面)の点のみを抽出し、その後、画像としてエクスポートする処理を記述しています。エクスポートする際の"resolution":"1.0"を変更すると任意の解像度を指定できます。
このJSONを使用して、コマンドライン上でpdal pipeline ./gdal.jsonと実行すると上記の処理が走ります。
点群が一つの場合はよいのですが、複数の点群を対象にする場合にJSONをいちいち編集するのは面倒です。
そこで、Pythonを使ってもう少し楽をしたいと思います。

PythonからPDALを呼び出して実行する方法

Pythonでは、PDALに渡すJSONを文字列として記述します。
以下のように処理対象の点群や解像度を定義してpythonを実行するとコマンドライン上で行った結果と同じことが実現できます。

main.py

import pdal

json = """{
    "pipeline": [
        "./08ME3552.las",
        {
            "type":"filters.expression",
            "expression":"Classification == 2"
        },
        {
            "filename":"./08ME3552.tif",
            "gdaldriver":"GTiff",
            "output_type":"min",
            "resolution":"1.0",
            "type": "writers.gdal"
        }
    ]
}"""

pipeline = pdal.Pipeline(json)
count = pipeline.execute()
arrays = pipeline.arrays
metadata = pipeline.metadata
log = pipeline.log

ディレクトリに格納されている点群を全なめしてDEMを作りたい場合などは、JSONを文字列として記述した点に着目してどうにかします。
辞書型になっていないので、雑な実装ですが、ファイル名を文字列置換したり print(json.replace("08ME3552", "08ME3553")) 解像度を文字列置換 print(json.replace("""1.0""", """0.5""")) したりで対応できます。