はじめに
様々な場所で点群データが公開されるようになりましたが、その提供範囲は図郭単位であることが多いと思います。色々いじっていると、もっと小さな範囲で切り出して、後続作業に用いたくなることもあるでしょう。
例えば下のような点群があったとして、競技場の部分がだけ欲しい時の切り出し方法を考えてみます。
任意の範囲の点群を切り出す方法
点群データの出典はオープンナガサキです。
まずは、欲しい範囲のポリゴンを作成します（赤線）。
QGIS上でgeopackegeとして作成しました。
lasデータの操作は、QGISを離れてpdalというものが使えます。
pdalはpython上でも使えるので、miniconda環境にインストールしました。
切り出しは１回の処理で出来ないようなので、pipelineを使っての処理になります。
python上ではこんな感じで動きました↓
import pdal
json = """{
"pipeline":
[
"01_las/point_cloud.las",
{
"type":"filters.overlay",
"dimension":"Classification",
"datasource":"vector.gpkg",
"layer":"polygon",
"column":"id"
},
{
"type":"filters.expression",
"expression":"Classification==2"
},
{
"type":"writers.las",
"filename":"01_las/output.las",
"forward":"all"
}
]
}"""
pipeline = pdal.Pipeline(json)
count = pipeline.execute()
arrays = pipeline.arrays
metadata = pipeline.metadata
log = pipeline.log
入力する点群は複数図郭にまたがっていたので、事前に１つにマージしておきました。
切り取りに使うポリゴンは、適当な値を属性に持たせておく必要があります。
- 入力LAS：01_las/point_cloud.las
- 入力ポリゴン：vector.gpkg｜polygon（属性"id"＝2)
やってることのイメージとしては、
①lasを読み込む
②ポリゴン内の点にClassificationとして標識を付ける
③標識が付いた点のみを抽出する
（リンク先で使われていたfilters.rangeが廃止予定とのことだったので、filters.expressionで代替しています）
④lasを書き出す
というような感じだと思います。"dimension"の"Classification"は、どうやらlasのclassificationと同じものらしいので、分類済みのlasを使う場合は、別のdimensionを使った方が良いなと思いました。
pdalはpython上でpipelineが利用できますが、上記のようにその部分だけjsonで記述し、python上では文字列として扱います。
jsonで記述しなくても、そのまま｜で処理を繋げていく方法もあるようです。
結果は、、、
見事に欲しい点群だけが切り取れました。
注意：変換に伴う座標値の精度低下について
lasの出力ファイル名を指定する記述の以下の部分ですが、
{
"type":"writers.las",
"filename":"01_las/output.las",
"forward":"all"
}
以下にような簡易的な記述でも可能です。
"01_las/output.las"
但し、前者だと入出力で点群の座標値が変わらないのに対して、後者の記述だと切り出し前後で同じ点の座標値が異なることが確認されました。
(左：後者で切り出し、右：元データ)
lasファイルでは、各ポイントの座標値はそのままの値では無く、scale値とoffset値で変換されて格納されています。
格納されている値 = (実際の値 - offset) / scale
両値はlasファイルのヘッダーに記載されていますが、この値が元データと出力データで異なることで、実際の座標値も影響を受け、異なる値になってしまいました。
対策として、ヘッダー情報を入出力で揃える、ということがあります。
入力データのヘッダー情報を維持したい場合は、"forward"で指定が可能です。
ヘッダーのどの値を維持するか細かく指定が可能ですが、"forward":"all"で全て維持してくれます。（scale値やoffset値、lasのversionなど）
なお、"writers.las"のデフォルトでは、scale値とoffset値はそれぞれ0.01と0になっています。
複数の範囲でそれぞれ切り出す方法
続いて、複数箇所で切り出し、それぞれを別のlasとして保存したい場合はどうすれば良いでしょうか。
先ほどの陸上競技場付近のいくつかの建物を囲うようにポリゴンを作成しました。
属性には"id"として1種類の値、"oid"として連番を持たせています。
↓こんな感じのコードでoidの順に点群の切り出しと出力が出来ました。
import geopandas as gpd
import pdal
import json
json_s = """{
"pipeline":
[
"INPUT_LAS",
{
"type":"filters.overlay",
"dimension":"Classification",
"datasource":"polygons.shp",
"query":"QUERY",
"column":"id"
},
{
"type":"filters.expression",
"expression":"Classification==2"
},
{
"type":"writers.las",
"filename":"OUTPUT_LAS",
"forward":"all"
}
]
}"""
json_d = json.loads(json_s)
gdf = gpd.read_file("polygons.shp")
for i in range(len(gdf)):
INPUT_LAS = "01_las/point_cloud.las"
QUERY = "SELECT * FROM polygons WHERE oid={}".format(i+1)
OUTPUT_LAS = "01_las/output_{}.las".format(i+1)
json_d["pipeline"][0] = INPUT_LAS
json_d["pipeline"][1]["query"] = QUERY
json_d["pipeline"][3]["filename"] = OUTPUT_LAS
json_s2 = json.dumps(json_d)
pipeline = pdal.Pipeline(json_s2)
count = pipeline.execute()
arrays = pipeline.arrays
metadata = pipeline.metadata
log = pipeline.log
json.loads()で文字列だったJSONを辞書として扱うことができます。
辞書とするとfor文の中でパラメーターが置き換えやすくなりました。
最後は、json.dumps()でもとに戻します。
1つずつポリゴンを処理に使えるように、filter.overlayにQUERY＝を追加しています。
属性"oid"の連番を使ったのですが、geopackageでは文法が分からずうまくクエリをかけれませんでした。
shpだとできたので、shpを使っています。
結果は、狙い通り６つのlasファイルが出力されました。
【おまけ】QGISでも点群処理できます
実は、最近のバージョンのQGISは、点群操作にも対応しています。
ver3.32以降です。
裏ではpdalが使われているようです。
用いるポリゴンの数が少ない場合は、GUI上の簡単な操作で処理が可能です。
QGISは点群の表示もできるので、単純作業だったらこっちのほうが簡単ですね。