兵庫県 50cmメッシュ DEM（2021～2022年度） - G空間情報センター
「DEM」は、航空レーザ測量成果のうち、建物、植生等を除去した地表面の数値標高モデルです。 データは50cm間隔の格子状で、XYZ座標値を示したテキストデータです。 ※兵庫県が実施した50cmメッシュの航空レーザ測量データを使用して作成.....
www.geospatial.jp
VIRTUAL SHIZUOKA 静岡県 中・西部 点群データ - G空間情報センター
航空レーザ測量（LP）および移動計測車両（MMS）により取得し、統合して活用できる3次元点群データです。 各ダウンロードページより、図郭単位で、LAS形式をZIPまたは7z圧縮したファイルのダウンロードができます。データの座標参照系は、日本...
www.geospatial.jp
「DEMをダウンロードしたらテキスト形式だった！」とか、「点群をラスタライズしたいけど、わざわざQGISを起動するのもなあ」とか「gdal_translaeteでなぜかエラーが出る！」思うことがありませんか？
私は何度もあります。
久しぶりの投稿はGMTからは脱線しますが、点群データをコマンドラインでGeotiffにする方法です。
準備
以下の内容を実行するには、GDAL, PDAL, LAStools無料版が必要です。
私だったらcondaで入れちゃいます。
#!/bin/bash
conda install gdal python-pdal lastools -c conda-forgeテキストデータからlasにする
実はこれが鍵です。いったんlasに変換させます。
これは、LAStoolsの txt2las を使えばOKです。
たとえば、兵庫県の50cmDEMはテキストデータなので、これをlasにするならば、
カレントディレクトリにテキストデータがあることを想定すれば、
#!/bin/bash
find . -type f -name "*.txt" | sed 's/\.txt//g' | xargs -P5 -I{} txt2las {}.txt -o {}.lasで終わりです。実行したディレクトリに*.lasが大量にできているはずです。
lasからDEM（Geotiff）にする
Page not found — Point Data Abstraction Library (PDAL)
pdal.io
lasからDEM（Geotiff）には、PDALを用います。
詳細なパラメータは公式ドキュメントを確認してください。
ここでは、カレントディレクトリに配置してある、*.las を一括でGeoTiffに変換します。
resolution は適当に1mにしてあります。50cmDEMだから0.5mにしてもいいのかもしれません。
import json
from osgeo import gdal
import pdal
import os
import glob
# Las ファイルのパス
input_files = glob.glob("./*.las")
for input_file in input_files:
output_file = os.path.splitext(input_file)[0] + ".tif"
pipeline_json = """
{
"pipeline":[
"%s",
{
"type":"writers.gdal",
"output_type":"mean",
"filename":"%s",
"resolution": 1
}
]
}
""" % (input_file, output_file)
pipeline = pdal.Pipeline(json.dumps(json.loads(pipeline_json)))
pipeline.execute()こうすれば、カレントディレクトリにGeoTiffが作成されているはずです。
マージしたい人は gdal_merge.py を使ってください。
なぜこれだとうまくいくか分からないのですが、今回はここまで。
コメント