QGIS3.32を利用してLASからCOPCとSTAC仕様に基づいたVirtual Point Clouds（仮想点群）を生成する

QGIS3.32（2023/07時点最新リリース）ではクラウドファウンディングにより得た資金でLutraConsulting社が点群データ処理を大幅に改良し、超パワーアップしました！

Native point cloud processing in QGIS

詳しくは上の記事を見てほしいですが、本当にいろんな機能が追加されているので、列挙します。

• LAS ⇆ LAZの変換
• グリッド化、ラスターへのエクスポート
• ベクトルデータ（ポイント）へのエクスポート
• 投影法割り当て
• 仮想点群の構築
• ベクターレイヤーでのクリップ
• COPCの作成
• メタデータの表示
• 点群データのマージ
• 点群データの投影変換
• 点群データの間引き
• 点群データからのポリゴン生成
• 点群データの境界作成
• 点群データから密度ラスターの作成
• 点群データの式によるフィルタリング

恐ろしいですね。

上記機能は、元々PythonなどからOpen3D・PyntcloudやPDALを利用することで自動化することなど出来ましたが、これがGUIでサクッとできるようになったのは本当にすごいと思います。

今回は、この中で個人的に激推しである「仮想点群（VPC）の構築」と「COPCの作成」について触れていきたいと思います。

データをダウンロード

フリーな点群データといえば「VIRTUAL SHIZUOKA」か「オープンナガサキ」が有名ですが、今回は静岡の方を使っていきましょう。

VIRTUAL SHIZUOKA 静岡県 中・西部 点群データから「LPデータ オリジナルデータ」をいくつか取得していきましょう。

静岡市周辺でいくつかデータをダウンロードします。

urlリストはこちら。

https://gic-shizuoka.s3.ap-northeast-1.amazonaws.com/2022/p/LP/LAS/08nd7793.zip
https://gic-shizuoka.s3.ap-northeast-1.amazonaws.com/2022/p/LP/LAS/08nd8703.zip
https://gic-shizuoka.s3.ap-northeast-1.amazonaws.com/2022/p/LP/LAS/08nd8713.zip
https://gic-shizuoka.s3.ap-northeast-1.amazonaws.com/2022/p/LP/LAS/08nd8714.zip
https://gic-shizuoka.s3.ap-northeast-1.amazonaws.com/2022/p/LP/LAS/08nd8704.zip
https://gic-shizuoka.s3.ap-northeast-1.amazonaws.com/2022/p/LP/LAS/08nd7794.zip
https://gic-shizuoka.s3.ap-northeast-1.amazonaws.com/2022/p/LP/LAS/08nd7795.zip
https://gic-shizuoka.s3.ap-northeast-1.amazonaws.com/2022/p/LP/LAS/08nd8705.zip
https://gic-shizuoka.s3.ap-northeast-1.amazonaws.com/2022/p/LP/LAS/08nd8715.zip

解答すると以下のようになりました。

280MB × 9ファイルですので、大体2.5GB程度ですね。
さほど大きくはないです。

投影法の設定

QGIS3.32を開き、静岡県なので、プロジェクトの投影法は事前に「EPSG:6676」に設定した上で、点群データをQGISに読み込ませておきましょう。

これだけで点群データが表示されると思いますが、点群データ自体にも投影法の設定も行ってあげましょう。

適当なレイヤのプロパティを見るとCRSが「不明」になっていると思います。
（日本で公開されているデータはほぼ全てCRSの設定がなされていません。）

レイヤプロパティを開いて、CRSを「EPSG:6676」に設定してあげましょう。

仮想点群の作成

プロセシングツールボックスの検索窓にVPCと入力すると「仮想点群（VPC）を作成」が表示されると思います。

レイヤを全て選択します。

各種パラメータは以下のような意味です。

• 境界ポリゴンを計算: データの正確な境界を表示できるようになる
• 統計量を計算: 属性の値の範囲計算する
• 概観点群をビルド: 元のデータの1000個ごとの点のみを使用しCOPCを生成する

今回はお試しなので、全てのオプションをオフにしておきます。
ただし、仮想点群はファイルとして出力するようにしておきましょう。

一瞬で処理が完了すると思います。

その後、作成された仮想点群にもCRSを設定しましょう。

現在は入力された点群データがCOPCではない場合は、点群の範囲のみが表示されるようですので、元データをCOPCに変換していきましょう。

LASをCOPCに変換

COPCと検索すると作成アルゴリズムが表示されます。

対象9レイヤを全て選択しましょう。

その後、適当なフォルダを選択して「実行」します。

※それなりに時間がかかります。

作成されたCOPCを読み込んでみましょう。

高速に描画されるようになっているはずです。
LASと同様にCRSの設定も行っていきましょう。

その後、作成したVPCをレイヤーに追加し、ズームしていくと、分類された点群データが表示されます。

（スタイル設定などで分類する前の生データなどの表示も行えると思いますが、今回は試していません。）

.vpcを確認してみる

最後に、生成されたデータを確認してみましょう。

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "stac_version": "1.0.0",
      "stac_extensions": [
        "https://stac-extensions.github.io/pointcloud/v1.0.0/schema.json",
        "https://stac-extensions.github.io/projection/v1.1.0/schema.json"
      ],
      "id": "08nd7793.copc",
      "geometry": {
        "coordinates": [
          [
            [
              -10800.0,
              -114000.0,
              9.79
            ],
            [
              -10800.0,
              -113700.0,
              9.79
            ],
            [
              -10400.0,
              -113700.0,
              109.67
            ],
            [
              -10400.0,
              -114000.0,
              109.67
            ],
            [
              -10800.0,
              -114000.0,
              9.79
            ]
          ]
        ],
        "type": "Polygon"
      },
      "bbox": [
        -10800.0,
        -114000.0,
        9.79,
        -10400.0,
        -113700.0,
        109.67
      ],
      "properties": {
        "datetime": "0001-01-01T00:00:00Z",
        "pc:count": 7787332,
        "pc:encoding": "?",
        "pc:schemas": [
          {
            "name": "X",
            "size": 8,
            "type": "floating"
          },
          {
            "name": "Y",
            "size": 8,
            "type": "floating"
          },
          {
            "name": "Z",
            "size": 8,
            "type": "floating"
          },
          {
            "name": "Intensity",
            "size": 2,
            "type": "unsigned"
          },
          {
            "name": "ReturnNumber",
            "size": 1,
            "type": "unsigned"
          },
          {
            "name": "NumberOfReturns",
            "size": 1,
            "type": "unsigned"
          },
          {
            "name": "ScanDirectionFlag",
            "size": 1,
            "type": "unsigned"
          },
          {
            "name": "EdgeOfFlightLine",
            "size": 1,
            "type": "unsigned"
          },
          {
            "name": "Classification",
            "size": 1,
            "type": "unsigned"
          },
          {
            "name": "ScanAngleRank",
            "size": 4,
            "type": "floating"
          },
          {
            "name": "UserData",
            "size": 1,
            "type": "unsigned"
          },
          {
            "name": "PointSourceId",
            "size": 2,
            "type": "unsigned"
          },
          {
            "name": "GpsTime",
            "size": 8,
            "type": "floating"
          },
          {
            "name": "ScanChannel",
            "size": 1,
            "type": "unsigned"
          },
          {
            "name": "ClassFlags",
            "size": 1,
            "type": "unsigned"
          },
          {
            "name": "Red",
            "size": 2,
            "type": "unsigned"
          },
          {
            "name": "Green",
            "size": 2,
            "type": "unsigned"
          },
          {
            "name": "Blue",
            "size": 2,
            "type": "unsigned"
          }
        ],
        "pc:type": "lidar",
        "proj:bbox": [
          -10800.0,
          -114000.0,
          9.79,
          -10400.0,
          -113700.0,
          109.67
        ],
        "proj:geometry": {
          "coordinates": [
            [
              [
                -10800.0,
                -114000.0,
                9.79
              ],
              [
                -10800.0,
                -113700.0,
                9.79
              ],
              [
                -10400.0,
                -113700.0,
                109.67
              ],
              [
                -10400.0,
                -114000.0,
                109.67
              ],
              [
                -10800.0,
                -114000.0,
                9.79
              ]
            ]
          ],
          "type": "Polygon"
        },
        "proj:wkt2": ""
      },
      "links": [],
      "assets": {
        "data": {
          "href": "./08nd7793.copc.laz",
          "roles": [
            "data"
          ]
        }
      }
    },
...
  ]
}

データは上記のようなGeoJSONになっています。

ただし、これはただのGeoJSONではなく、「STAC Item」という仕様に則ったGeoJSONになっています。

STACについては以下の記事などが参考になります。

• AWSのオープンデータとPythonを使って、宇宙から地球を覗いてみよう！
• 静的な STAC カタログを stac-ruby でつくる
• SpatioTemporal Asset Catalogとは

おわりに

今回はQGISでCOPCへの変換と仮想点群の生成を行っていきました！

QGISの3.32では今回紹介した以外にも以下のような点群操作ができるようになっているので、皆さん使ってみてください！