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

ERROR: The request could not be satisfied

ERROR: The request could not be satisfied

「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）にする

https://pdal.io/en/2.6.0/stages/writers.gdal.html

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

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