PDALで点群データを処理

1年を振り返って

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

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

点群データ

静岡県の点群データ公開は話題になりましたね。
その他にも土木学会インフラデータチャレンジにおいて、国土交通省　那賀川河川事務所・福井河川国道事務から九頭竜川、那賀川のLASデータが公開されていますし、国土交通データプラットフォームからICT土工の工事成果品として取得した3次元データのうち233工事が提供されています。
簡単にデータが手に入るようになってきました。

PDAL

みなさんどうやって点群データ扱ってるの？

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

まずは点群データを用意

G空間情報センターから那賀川ALBデータ(LAS)をダウンロードしましょう。
ダウンロードページに表示されるメッシュ枠をクリックすると、ダウンロードURLが表示されます。

ダウンロードしたzipを解凍すると、数面のlasデータが入っています。

表示して確認するならCloudCompareが便利です。
04GG5490に含まれる3ファイルを表示してみたところです。

PDALのインストール

ここは省略。PDAL開発チームにしたがってCondaを使って行いましょう。

点群の情報をみる

さてインストールが終わって、手元に点群データもありますので、まずは点群データの情報をみてみます。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で処理を記述しておき実行させることができます。

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を使用して書くこともできます。

merge.json

[
    "*las",
    "merge.las"
]

pdal pipeline merge.json

CloudCompareでマージしたlasを表示してみます。

おわり

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