Python版PDALを使用してLasファイルをGeoPackage形式に変換する - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

本記事はGitHub、nbviewer、Colabでも公開しています。

※Colabの解説記事はこちら

導入

Lasファイルを処理する場合、PDAL - Point Data Abstraction Libraryを使用することが多いです。PDALはC++で書かれたライブラリですが、Python バインディングも提供されています。ここでは、Python版PDALを使用してLasファイルをGeoPackage形式に変換する方法を紹介します。

Pyhon版PDALはオプション一覧をJSON形式のテキストで渡すというPythonらしからぬ書き方で正直かなり使いにくいですが、後処理含めてすべてPythonで完結するのでメリットはあります。（GDALのように誰かがラッパーを書いてくれると嬉しいのですが。）

単に変換する方法に加えて、GeoPandasを用いた後処理の例も示しています。

仮想環境構築

他のパッケージへの依存もあるので、新しい仮想環境を作成することをお勧めします。また、仮想環境はcondaで作成することをお勧めします。

以下のコマンドで必要なパッケージをインストールした仮想環境を作成できます。

conda create -n pdalws -c conda-forge jupyterlab python-pdal gdal entwine matplotlib geopandas

Pythonコード

基本版：PDALのみを使用してLasファイルをGeoPackage形式に変換する

PDALのオプションはpiplineとしてJSON形式で指定します。オプションの詳細は、公式サイトをご覧ください。

import pdal
import json

input_path = "tmp.las"

# PDALのパイプラインをJSON形式で定義
pipeline_json ={
  "pipeline": [
        {"type": "readers.las", "filename": input_path},
    { "type": "filters.reprojection",
      "in_srs": "EPSG:6675",           
      "out_srs": "EPSG:6675" 
     }, # 入力ファイルに設定されている場合は省略可
    {
      "type": "writers.ogr",
      "filename": "output.gpkg",
      "ogrdriver": "GPKG",
      "attr_dims": "Intensity,Classification",
    #   "attr_dims": "all",
      "ogr_options": "SPATIAL_INDEX=YES,FID=FID"
    }
  ]
}

# パイプラインの実行
p = pdal.Pipeline(json.dumps(pipeline_json))
p.execute()

応用版：PDALでLasファイルを読み込みGeoPandas形式に変換する。

基本編を元に、PDALでLasファイルを読み込みGeoPandas形式に変換する方法を示します。

import pdal
import json

import numpy as np
from shapely.geometry import Point
import geopandas as gpd

input_path = "tmp.las"

# PDALのパイプラインをJSON形式で定義
pipeline_json = {
    "pipeline": [
        {"type": "readers.las", "filename": input_path},
    ]
}

# パイプラインの実行
p = pdal.Pipeline(json.dumps(pipeline_json))
p.execute()

# 結果の取得
arrays = p.arrays

# 配列を結合
data = arrays[0] if len(arrays) == 1 else np.concatenate(arrays)

# 必要フィールド取り出し（存在チェック）
def get_col(name, default=None):
    return data[name] if name in data.dtype.names else default

x = get_col("X"); y = get_col("Y"); z = get_col("Z")
cls = get_col("Classification"); intensity = get_col("Intensity")
ret = get_col("ReturnNumber"); nret = get_col("NumberOfReturns")

geoms = [Point(float(xi), float(yi)) for xi, yi in zip(x, y)]
gdf = gpd.GeoDataFrame(
    {
        "Z": z.astype(float),
        **({"Classification": cls} if cls is not None else {}),
        **({"Intensity": intensity} if intensity is not None else {}),
        **({"ReturnNumber": ret} if ret is not None else {}),
        **({"NumberOfReturns": nret} if nret is not None else {}),
    },
    geometry=geoms
)

# 座標参照系の設定
d = gdf.set_crs("EPSG:6675", inplace=True)

ここまででGeoPandasのGeoDataFrameが作成されましたので、さまざまな後処理が可能です。

# GeoPackage形式で保存する場合
gdf.to_file("output2.gpkg", layer="points", driver="GPKG")

GitHub

github.com

参考サイト

jupyterのTeXをはてブロのmarkdownに変換するpythonスクリプト:2024/1/3時点 - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning