PDALで点群データを処理

# 1年を振り返って
今年一年を振り返って、コロナ禍の事務所で1人寂しく何をしていたのかと思い返してみると、ブラウザでは常にMapbox GL JSやDeckGLを使って3D地図を表示していました。

3Dデータを処理する機会が多くなったのですが、その扱いについては情報が共有されてはいないなと感じています。

# 点群データ
静岡県の点群データ公開は話題になりましたね。
その他にも[土木学会インフラデータチャレンジ](https://urbandata-challenge.jp/news/idc2020)において、国土交通省　那賀川河川事務所・福井河川国道事務から[九頭竜川、那賀川のLASデータ](https://www.geospatial.jp/ckan/dataset/doboku-alb)が公開されていますし、[国土交通データプラットフォーム](https://www.mlit-data.jp/platform/)から[ICT土工の工事成果品として取得した3次元データのうち233工事](https://urbandata-challenge.jp/idc/idc-pointcloud)が提供されています。
簡単にデータが手に入るようになってきました。

# PDAL
## みなさんどうやって点群データ扱ってるの？
そんな簡単に手に入るようになった点群データですが、ちょっとしたフォーマット変換や投影変換、間引きしたいなどを行いたい際は、みなさん何で処理してるんでしょう？？
私なんかは「[PDAL](https://pdal.io/)でしょ」とかすぐ思うんですが、なぜかPDALの情報を日本語では見かけないんです...
ということで、PDALの使い方を紹介しておきます。

## まずは点群データを用意
G空間情報センターから[那賀川ALBデータ(LAS)](https://www.geospatial.jp/ckan/dataset/doboku-alb)をダウンロードしましょう。
ダウンロードページに表示されるメッシュ枠をクリックすると、ダウンロードURLが表示されます。
![スクリーンショット 2020-12-07 22.29.01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58748/106ddb52-ea17-c433-8904-03eb19234170.png)
ダウンロードしたzipを解凍すると、数面のlasデータが入っています。

表示して確認するなら[CloudCompare](https://www.danielgm.net/cc/)が便利です。
04GG5490に含まれる3ファイルを表示してみたところです。
![スクリーンショット 2020-12-07 22.32.12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58748/e3439a57-56fb-efcb-0f68-b5e5e4dd3dbe.png)

## PDALのインストール
ここは省略。[PDAL開発チームにしたがってCondaを使って行いましょう。](https://pdal.io/quickstart.html#)

## 点群の情報をみる
さてインストールが終わって、手元に点群データもありますので、まずは点群データの情報をみてみます。`pdal info`コマンドを使用します。

```
pdal info 04GG5490-AB.las
```

とファイル名を指定して実行すると、bounding boxの情報や各カラムの統計情報が表示されます。


```
{
  "file_size": 5652369,
  "filename": "04GG5490-AB.las",
  "now": "2020-12-07T22:44:34+0900",
  "pdal_version": "2.1.0 (git-version: Release)",
  "reader": "readers.las",
  "stats":
  {
    "bbox":
    {
      "native":
      {
        "bbox":
        {
          "maxx": 96399.999,
          "maxy": 102183.51,
          "maxz": 129.79,
          "minx": 96322.921,
          "miny": 102150,
          "minz": 66.881
        },
        "boundary": { "type": "Polygon", "coordinates": [ [ [ 96322.921000000002095, 102150.0, 66.881 ], [ 96322.921000000002095, 102183.51, 66.881 ], [ 96399.998999999996158, 102183.51, 129.79 ], [ 96399.998999999996158, 102150.0, 129.79 ], [ 96322.921000000002095, 102150.0, 66.881 ] ] ] }
      }
    },
    "statistic":
    [
      {
        "average": 96373.39082,
        "count": 148729,
        "maximum": 96399.999,
        "minimum": 96322.921,
        "name": "X",
        "position": 0,
        "stddev": 18.13699802,
        "variance": 328.9506973
      },
....
```

```
pdal info 04GG5490-AB.las --metadata
```
でメタデータだけを表示させることができますし、

```
pdal info 04GG5490-AB.las -p 0
```
のように0個目のポイントの情報だけを表示させる、といったこともできます。

## フォーマット変換してみる

`-drivers`オプションで表示される`readers` `writers`から、対応しているドライバーを確認できます。

```
pdal --drivers
```

例えば、lasからcsvへ変換する場合は`writes.text`を使用します。
下記のように、`入力ファイル名` `出力ファイル名`に続けて`writers.text.format`へ`csv`を指定します。

```
pdal translate 04GG5490-AB.las test.csv --writers.text.format=csv
```

同じく`writes.text.format`へ`geojson`を指定することでgeojsonへの変換することもできます。

```
pdal translate 04GG5490-AB.las test.json --writers.text.format=geojson
```

## フィルタを使ってみる

フォーマット変換に合わせて、投影変換も行ってみましょう。
オプションとして`filters.reprojection`を利用します。

```
pdal translate 04GG5490-AB.las test.csv --writers.text.format=csv --writers.text.precision=9  -f filters.reprojection --filters.reprojection.in_srs="EPSG:2448" --filters.reprojection.out_srs="EPSG:4326"
```

さらにフィルタを重ねることもできます。間引きも一緒に行ってみます。

```
pdal translate 04GG5490-AB.las test.csv --writers.text.format=csv --writers.text.precision=9  -f filters.reprojection --filters.reprojection.in_srs="EPSG:2448" --filters.reprojection.out_srs="EPSG:4326" -f filters.sample --filters.sample.radius=0.00001
```

## pipeline
複数の処理を重ねて行っていく場合は、pipelineを使用するのが便利です。
jsonで処理を記述しておき実行させることができます。

```json:pipeline.json
[
    "04GG5490-AB.las",
    {
        "type":"filters.reprojection",
        "in_srs":"EPSG:2448",
        "out_srs":"EPSG:4326"
    },
    {
        "type":"filters.sample",
        "radius":0.00001
    },
    {
        "type":"writers.text",
        "format":"csv",
        "precision":9,
        "filename":"test.csv"
    }
]
```

コマンド`pipeline`に処理を記述したjsonを渡してあげます。

```
pdal pipeline pipeline.json
```

## merge

折角なので他のコマンドも試しましょう。
複数lasを1つのまとめてみます。
`merge`コマンドへ`複数ファイル名`を与えて、最後を`出力ファイル名`とします。

```
pdal merge 04GG5490-AB.las 04GG5490-BA.las 04GG5490-BB.las merge.las
```

pipelineを使用して書くこともできます。

```json:merge.json
[
    "*las",
    "merge.las"
]
```

```
pdal pipeline merge.json
```

CloudCompareでマージしたlasを表示してみます。
![スクリーンショット 2020-12-08 19.15.17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58748/3351dd64-a81d-e0dd-141c-263f1853162e.png)

# おわり
PDAL活用の情報が共有されてくるとうれしいです。