PythonとPDALを使って点群（LAS）データをDEMに変換したり、Open3DでglTFに変換してみる！

![1745409134711.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/b7672eb4-9122-4125-8956-2261979327ab.png)


みなさん、こんにちは！
PDAL使ってますか！！！

PDALは、点群データを操作するための強力なCLIツールかつ、ライブラリです。

と言うのはこちらの記事でいろいろ書いていますので、ぜひご覧ください。

https://qiita.com/nokonoko_1203/items/fb842d163cc9c3e56a4b

で、今回はこのPDALをPythonから呼び出して利用してみようと思います。

PDALはC++で書かれたライブラリですが、Pythonからも利用できるようにラッパーが用意されています。

PDALにはパイプラインという概念があり、JSON形式でルールに従って記載することで、点群データの一連の処理をまとめて実行することができます。
このため、CLIからJSONを指定して実行することが多いのですが、Pythonからも同様にJSONを指定して実行することができます。

Pythonから利用できると、例えばWebアプリケーションからPDALを利用したり、Pythonのライブラリ（numpyや今回利用するOpen3Dなど）と組み合わせて利用することが容易になります。

この記事ではPDALでDEMを作成する・地形のglTFを作成する、といった利用方法について解説していきます。

この記事のコードはこちらに載せているので、参考にしてください！

https://github.com/nokonoko1203/pointcloud-processing-sample

## データの準備

今回は、LAS形式のデータを変換します！

どんなデータを利用してもいいのですが、今回は[こちら](https://qiita.com/nokonoko_1203/items/847cb7fcd6fea4b36133)の記事で紹介した東京都の点群データを利用します。

https://www.geospatial.jp/ckan/dataset/tokyopc-23ku-2024

どこでもいいのですが、適当な場所のデータをダウンロードしておいてください。

## 必要なパッケージのインストール

uvの利用を前提とするので、uvがインストールされていない場合は、uvをインストールしてください。

その後、必要なパッケージをインストールします。

```bash
uv add pdal open3d rasterio matplotlib jupyter
```

こんな感じの`pyproject.toml`ができていればOKです。

```toml
[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スクリプトでも問題ありません。

まずは必要なライブラリをインポートします。

```Python
import json
import pdal
import rasterio
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource
```

ファイルのパスを指定しますが、ダウンロードしたパスに合わせて変更してください。

```Python
input_path = "./data/input/09LD1895.las"
output_path = "./data/output/09LD1895_dem.tif"
```

定義するパイプラインはこの程度の記述でOKです。
今回はPythonで書いていますが、同様の内容をJSON形式で書けばCLIから実行することもできます。

```Python
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)
```

このパイプラインは、以下の処理を行います。

1. `readers.las`でLAS形式のデータを読み込みます
2. `filters.smrf`で地面を抽出します
3. `filters.range`で地面の点群を抽出します
4. `writers.gdal`でDEMを出力します

パラメータは、デフォルト値を参考にしながら調整してください。
- [readers.las](https://pdal.io/en/stable/stages/readers.las.html)
- [filters.smrf](https://pdal.io/en/stable/stages/filters.smrf.html)
- [filters.range](https://pdal.io/en/stable/stages/filters.range.html)
- [writers.gdal](https://pdal.io/en/stable/stages/writers.gdal.html)

特に、地表面抽出のsmrfは目的によってシビアかつパラメーターが多いので、色々試してみてください！

ではこのパイプラインを実行してみましょう。

```Python
pipeline = pdal.Pipeline(pipeline_definition)
pipeline.execute()
```

実行が終わったら、出力されたDEMを確認してみましょう。

```Python
with rasterio.open(output_path) as src:
    dem_data = src.read(1, masked=True) # nodataをマスクして読み込む
    transform = src.transform
```

するとこんな感じのDEMができているはずです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/bb99021d-e825-4223-9990-dca63dbe6dbe.png)

QGISでも確認してみましょう。

![1745415753004.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/7e2ec35c-1f7b-4d25-aa25-83b51501c583.png)

いい感じですね！

## 地形のglTFの作成

利用用途はさほどないんですが、こんなこともできるよーという紹介として、地形のglTFを作成してみます。

まずは必要なライブラリをインポートします。

```Python
import json
import pdal
import open3d as o3d
```

ファイルのパスを指定しますが、ダウンロードしたパスに合わせて変更してください。

```Python
pointcloud_path = "./data/input/09LD1895.las"
ply_path = "./data/output/09LD1895_mesh_terrain.ply"
glb_path = "./data/output/09LD1895_mesh_terrain.glb"
```

パイプラインは以下のように記述します。

```Python
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))
```

このパイプラインは、以下の処理を行います。
1. `readers.las`でLAS形式のデータを読み込みます。
2. `filters.smrf`で地面を抽出します。
3. `filters.range`で地面の点群を抽出します。
4. `filters.poisson`でメッシュを作成します。
5. `writers.ply`で`ply`形式のメッシュを出力します。

パラメータは、デフォルト値を参考にしながら調整してください。

- [filters.poisson](https://pdal.io/en/stable/stages/filters.poisson.html)
- [writers.ply](https://pdal.io/en/stable/stages/writers.ply.html)

以下のように実行します。

```Python
pipeline_obj.execute()
```

これで`ply`形式のメッシュができました。
このファイルをOpen3DをglTF形式に変換してみます。

```Python
mesh_terrain = o3d.io.read_triangle_mesh(ply_path)
o3d.io.write_triangle_mesh(glb_path, mesh_terrain)
```

これでglTF形式のメッシュができました。
[glTF Viewer](https://gltf-viewer.donmccurdy.com/)で確認してみましょう。

![1745409134711.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/e4338e5d-bb8c-46e0-bc05-af5993e5dc8a.png)

## おわりに

今回はPDALをPythonから利用して、DEMを作成したり、地形のglTFを作成したりしてみました。
PDALはCLIからも利用できるので、Pythonから利用する必要がない場合はCLIから実行しても問題ありません。
PDALは非常に強力なツールなので、ぜひ活用してみてください！