点群データやテキストファイルのDEMをGeoTiffにする

兵庫県 50cmメッシュ DEM（2021～2022年度） - G空間情報センター

「DEM」は、航空レーザ測量成果のうち、建物、植生等を除去した地表面の数値標高モデルです。データは50cm間隔の格子状で、XYZ座標値を示したテキストデータです。 ※兵庫県が実施した50cmメッシュの航空レーザ測量データを使用して作成.....

VIRTUAL SHIZUOKA 静岡県中・西部点群データ - G空間情報センター

航空レーザ測量（LP）および移動計測車両（MMS）により取得し、統合して活用できる3次元点群データです。各ダウンロードページより、図郭単位で、LAS形式をZIPまたは7z圧縮したファイルのダウンロードができます。データの座標参照系は、日本...

「DEMをダウンロードしたらテキスト形式だった！」とか、「点群をラスタライズしたいけど、わざわざQGISを起動するのもなあ」とか「gdal_translaeteでなぜかエラーが出る！」思うことがありませんか？
私は何度もあります。

久しぶりの投稿はGMTからは脱線しますが、点群データをコマンドラインでGeotiffにする方法です。

準備
テキストデータからlasにする
lasからDEM（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）にする

writers.gdal — 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 を使ってください。

なぜこれだとうまくいくか分からないのですが、今回はここまで。