点群データ Advent Calendar 2025

クラス分類された3次元点群データからDEMを作成する方法について

Last updated at 2026-01-25Posted at 2025-12-07

この記事は点群データ Advent Calendar 2025 7日目の記事です。

はじめに

本記事では、クラス分類された航空レーザ測量の 3 次元点群データ（オリジナルデータ（las））から、0.5m DEM（数値標高モデル）を作成する一連のワークフローを整理する。
埼玉県が公開している山間部の3次元点群データを対象に、
- PDAL と laspy による前処理
- GMT（triangulate）によるTIN内挿補間
- GDAL による VRT モザイクと GeoTIFF 出力

といった一般的なDEM生成フローを明確に示すことを目的とする。

なお、DEM作成時に参考とした文献は以下のとおりである。

参考文献:
国土地理院「航空レーザ測量による数値標高モデル（ＤＥＭ）作成マニュアル（案）（平成18年4月）」
https://www.gsi.go.jp/common/000258818.pdf

3次元点群データ及びDEMの概要

3次元点群データ（オリジナルデータ）の特徴

本記事で用いるデータは、埼玉県がオープンデータとして提供する航空レーザ測量に基づく3次元点群データである。LAS形式で提供され、
- XYZ座標
- 分類情報（Classification）
- RGB値（Red / Green / Blue）

などの属性を含む。

座標系は JGD2011 / 平面直角座標系 9系（EPSG:6677）。
データは国土基本図図郭地図情報レベル2500（2km×1.5km）単位で提供される。
Classification 値は、測技協による LAS Specification 1.4 日本語版に基づき定義されており、地面（Ground）は Classification=2 である。
本稿ではこの Class 2 の点群のみを抽出して DEM を作成する。

出典: 国土基本図の図郭とは？図郭コードの意味と読み方をわかりやすく解説

出典: 測技協 LAS Specification 1.4 – R15 (ASPRS) の日本語版解説書

DEM とは

DEM（Digital Elevation Model: 数値標高モデル）は、地表面標高を格子状のラスターで表したものであり、地形解析・土木設計・災害評価など幅広い用途がある。今回は0.5m解像度で生成する。

出典: 中学生のページ：DEM　https://www.gsi.go.jp/KIDS/KIDS16.html

PC処理環境

OS：Windows 11 Pro(WSL2/Ubuntu)
CPU：AMD Ryzen 7 5700X 8-Core Processor
メモリ：64GB
SSD：4TB（空き容量 1TB）

※ TINの生成では大量にメモリを消費するため、大規模なデータを処理する場合は高スペックなPC処理環境が必要となる。

使用ツール

使用する点群データ

埼玉県がオープンデータとして公開している山間部の 3 次元点群データ（オリジナルデータ）を使用する。
埼玉県 GIS ポータルサイトからダウンロード可能である。
https://portal-pref-saitama.hub.arcgis.com/

あるいはG空間情報センターからもデータのダウンロードが可能である（こちらのほうがダウンロードしやすい）。
https://www.geospatial.jp/ckan/dataset/mountainous-area-pointcloud-saitama

座標参照系：JGD2011 / 平面直角座標系 9 系（EPSG: 6677）
ライセンス：CC BY 4.0
国土基本図図郭（地図情報レベル2500）単位で提供される。
使用する 3 次元点群データ（オリジナルデータ）の図郭（9つ）：

国土基本図図郭とは？
上の図で示したとおり、国土基本図とは、国土地理院が統一された図式に基づいて作成・整備している大縮尺の基本図のことである。また、国土基本図の図郭（ずかく）とは、国土基本図を一定の基準で区切った個別の地図範囲のことである。図郭は平面直角座標系を使用し、原点からの距離（メートル）を基準にして区分される。
詳細は下記を参照。
https://club.informatix.co.jp/?p=1293

DEM作成範囲は4つの図郭（赤枠の範囲）だが、図郭単位のDEMをマージする際に図郭境界の整合性を確保するため、周辺を含めた9つの図郭（青枠の範囲）の点群データを使用する。
図郭単位のデータは、DEM生成時に図郭境界でギャップが生じやすいため、隣接図郭を含めて処理範囲を拡張する方法を採用する。
具体には、1つの図郭のDEMの作成に対しては、右、下、右下の図郭を含めた計4つの図郭を使用して、DEM作成範囲を右、下にそれぞれ100mオフセットした範囲で切り抜く形で作成する。

点群処理〜DEM生成の詳細手順

DEM 生成の全体フロー

DEM生成の全体フローは以下のとおりである。

点群（LAS）の内容確認
LASのマージ
グラウンド点群の抽出
TIN内挿補間（GMT triangulate）
VRT作成（gdal_translate）
VRTモザイク（gdalbuildvrt）
GeoTIFF出力（gdalwarp）

以下で、具体的な手順を詳述する。

オリジナルデータ（las）の中身を確認

コマンドプロンプトからWSLを起動して、pdal info 09kc131.las -p 0を実行する。
3次元点群データの各点には以下の情報が付与されている。
XYZ座標、分類情報（Classification）、RGB値（Red / Green / Blue）等が確認できる。

$ pdal info 09kc131.las -p 0
{
  "file_size": 1952036277,
  "filename": "09kc131.las",
  "now": "2025-11-30T11:32:43+0900",
  "pdal_version": "2.6.2 (git-version: Release)",
  "points":
  {
    "point":
    {
      "Blue": 81,
      "Classification": 5,
      "EdgeOfFlightLine": 0,
      "GpsTime": 269774.5172,
      "Green": 78,
      "Intensity": 65,
      "KeyPoint": 0,
      "NumberOfReturns": 2,
      "Overlap": 0,
      "PointId": 0,
      "PointSourceId": 9,
      "Red": 74,
      "ReturnNumber": 2,
      "ScanAngleRank": 14,
      "ScanDirectionFlag": 1,
      "Synthetic": 0,
      "UserData": 25,
      "Withheld": 0,
      "X": -68000,
      "Y": -3130.81,
      "Z": 369.83
    }
  },
  "reader": "readers.las"
}

オリジナルデータ（las）のマージ

オリジナルデータ（las）は、図郭単位となっているが、1つの図郭の中でもlasが計測エリア（発注業務単位）ごとに分割されているケースがあるため（例：09kc232.lasはchichibuとhannoに跨る）、lasをマージする必要がある。

las
├─chichibu
│      09kc131.las
│      09kc132.las
│      09kc133.las
│      09kc134.las
│      09kc141.las
│      09kc143.las
│      09kc231.las
│      09kc232.las
│      09kc241.las
│
└─hanno
        09kc232.las
        09kc241.las

lasのマージは、以下のスクリプト（pdal merge）を用いて行う。
コマンドプロンプトからWSLを起動して、以下のスクリプト（1_las-merge.sh）を実行する。

#!/bin/bash

# 実行方法
# chmod +x las-merge.sh
# dos2unix las-merge.sh
# ./las-merge.sh

# lasのマージ
find las -name "*.las" -printf "%f\n" | sort | uniq | while read fname; do
    files=$(find las -name "$fname")
    count=$(echo "$files" | wc -l)

    if [ $count -gt 1 ]; then
        echo "Merging $fname"
        pdal merge $files las_merge/$fname
    else
        echo "Copying $fname"
        cp $files las_merge/
    fi
done

グラウンドデータの作成

グラウンドデータの作成は、以下のスクリプト（laspy）を用いて行う。
コマンドプロンプトからWSLを起動して、以下のスクリプト（2_las2grd.py）を実行する。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
las_merge 配下の *.las から Classification=2 (ground) の点だけ抜き出して、
grd 配下に 〈元ファイル名〉_grd.txt を出力する。

使い方:
    python las2grd.py
"""

from pathlib import Path
import laspy  # pip install laspy

# ===== ここが「埋め込んだ引数」相当 =====
INPUT_DIR  = Path("./las_merge")   # -i ./las_merge
OUTPUT_DIR = Path("./grd")         # -o ./grd
CLASSVAL   = 2                     # ground = 2
PRECISION  = 2                     # x,y,z の小数点以下桁数
# ===================================


def export_grd_from_las(las_path: Path, out_txt: Path):
    """単一の LAS から class=CLASSVAL の点だけ id,x,y,z で書き出す"""
    las = laspy.read(str(las_path))

    cls = las.classification
    mask = (cls == CLASSVAL)

    out_txt.parent.mkdir(parents=True, exist_ok=True)

    if mask.sum() == 0:
        # class=2 が一つもない場合は空ファイルだけ作る
        out_txt.write_text("", encoding="utf-8")
        print(f"EMPTY (no class={CLASSVAL}): {las_path.name}")
        return

    idx = mask.nonzero()[0]
    x = las.x[mask]
    y = las.y[mask]
    z = las.z[mask]

    fmt = f"{{:d}},{{:.{PRECISION}f}},{{:.{PRECISION}f}},{{:.{PRECISION}f}}\n"
    with out_txt.open("w", encoding="utf-8", newline="") as f:
        for i, xi, yi, zi in zip(idx, x, y, z):
            f.write(fmt.format(int(i), xi, yi, zi))

    print(f"DONE: {out_txt}")


def main():
    OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

    las_files = sorted(INPUT_DIR.glob("*.las"))
    if not las_files:
        print(f"No LAS files found in: {INPUT_DIR}")
        return

    print(f"Input dir : {INPUT_DIR}")
    print(f"Output dir: {OUTPUT_DIR}")
    print(f"Files     : {len(las_files)}")

    for las_path in las_files:
        basename = las_path.stem
        out_txt = OUTPUT_DIR / f"{basename}_grd.txt"
        if out_txt.exists():
            print(f"SKIP (exists): {out_txt}")
            continue
        export_grd_from_las(las_path, out_txt)

    print("Finished.")


if __name__ == "__main__":
    main()

グリッドデータの作成

グリッドデータの作成は、以下のGMTコマンド（gmt triangulate）を用いて行う。
GMT Command PromptからWSLを起動して、以下のコマンドを実行する。

gmt triangulate 09kc131_grd.txt 09kc133_grd.txt 09kc132_grd.txt 09kc134_grd.txt -i1,2,3 -Gkc131-trn.grd -I0.5 -R-67900/-65900/-4600/-3100 -E-9999 -V
gmt triangulate 09kc132_grd.txt 09kc134_grd.txt 09kc141_grd.txt 09kc143_grd.txt -i1,2,3 -Gkc132-trn.grd -I0.5 -R-65900/-63900/-4600/-3100 -E-9999 -V
gmt triangulate 09kc133_grd.txt 09kc231_grd.txt 09kc134_grd.txt 09kc232_grd.txt -i1,2,3 -Gkc133-trn.grd -I0.5 -R-67900/-65900/-6100/-4600 -E-9999 -V
gmt triangulate 09kc134_grd.txt 09kc232_grd.txt 09kc143_grd.txt 09kc241_grd.txt -i1,2,3 -Gkc134-trn.grd -I0.5 -R-65900/-63900/-6100/-4600 -E-9999 -V

-i1,2,3：使用する列番号
1列目＝x、2列目＝y、3列目＝z を読み込む指定。

-Gkc131-trn.grd：出力GRDファイル
TIN 内挿補間した DEM を kc131-trn.grd に保存。

-I0.5：解像度（セルサイズ）
出力 DEM の解像度を 0.5 m に設定。

-R-67900/-65900/-4600/-3100：処理範囲

-Rxmin/xmax/ymin/ymax は GMT が処理する空間領域を示す。

xmin：左端（西端）の X 座標
xmax：右端（東端）の X 座標
ymin：下端（南端）の Y 座標
ymax：上端（北端）の Y 座標

本記事の例では、元の図郭（KC131）の平面直角座標の範囲をもとに、
境界付近の継ぎ目を抑えるために上下左右を 100 m ずらした範囲を指定している。

たとえば -R-67900/-65900/-4600/-3100 は、

X 方向：図郭の左端・右端からそれぞれ 100 m 内側に入った範囲
Y 方向：図郭の上下端から 100 m ずつずらした範囲

といったイメージで設定している。（実際の値は図郭の座標に依存する）

VRTファイルの作成

VRTファイルの作成は、以下のGDALコマンド（gdal_translate）を用いて行う。
GMT Command PromptからWSLを起動して、以下のコマンドを実行する。

gdal_translate -of VRT -a_srs EPSG:6677 -a_nodata -9999 kc131-trn.grd kc131.vrt
gdal_translate -of VRT -a_srs EPSG:6677 -a_nodata -9999 kc132-trn.grd kc132.vrt
gdal_translate -of VRT -a_srs EPSG:6677 -a_nodata -9999 kc133-trn.grd kc133.vrt
gdal_translate -of VRT -a_srs EPSG:6677 -a_nodata -9999 kc134-trn.grd kc134.vrt

VRTファイルのマージ

VRTファイルのマージは、以下のGDALコマンド（gdalbuildvrt）を用いて行う。
GMT Command PromptからWSLを起動して、以下のコマンドを実行する。

# vrtファイルのフルパスリストを作成
dir /s /b *.vrt > list.txt

# モザイクVRTを作成
gdalbuildvrt -input_file_list list.txt mosaic.vrt

GeoTIFFの作成

GeoTIFFの作成は、以下のGDALコマンド（gdalwarp）を用いて行う。
GMT Command PromptからWSLを起動して、以下のコマンドを実行する。

gdalwarp ^
  -t_srs EPSG:6677 ^
  -tr 0.5 0.5 -tap ^
  -r bilinear ^
  -srcnodata -9999 -dstnodata -9999 ^
  -co BIGTIFF=YES -co TILED=YES -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 ^
  mosaic.vrt merged_05m.tif

作成した0.5mDEM（GeoTIFF）をQGISで表示して確認

QGISに0.5mDEM（GeoTIFF）をレイヤーとして追加し、レイヤプロパティ＞シンポロジ＞レンダリングタイプ（単バンド擬似カラー）で好みのカラーランプを適用して確認する。
なお、下記では、0.5mDEM（GeoTIFF）をもうひとつレイヤーとして追加し、レイヤプロパティ＞シンポロジ＞レンダリングタイプ（陰影図(hillshade)）を30%の不透明度に設定して重ねている。

また、上記の0.5mDEM（GeoTIFF）とQGISプラグイン「CSMap Plugin」を用いて、CS立体図を作成できる。

課題と今後の展望

課題

本ワークフローは動作するが、実務的に改善すべき点が複数存在する。

① 計算負荷が大きい

triangulate（TIN内挿補間）は点数が多いとメモリ消費が非常に大きく、数千万点規模では処理が不安定になりやすい。
分割処理等が必要だが、手順が煩雑。

② Class=2 の信頼性問題

公開点群の分類精度は発注業務に依存し、地面が誤って分類されているケースもある。
特に急峻地形・樹林境界などではノイズや誤分類によりDEMに「凹凸」が生じることがある。

③ 処理フローがスクリプト＋手作業混在

PDAL, laspy, GMT, GDAL が別々のステップで動作しており、再現性確保が難しい。
一括処理（バッチ化）が不十分。

④ 0.5m DEM の最適性検証が未実施

点群密度に対して0.5m格子が統計的に適切かどうかの考察が未実施。
特に「TIN内挿補間でどの程度の精度が確保できるか」の検証フェーズが必要。

今後の展望

上記の課題を踏まえ、今後取り組むべき改善方向をまとめる。

① PDALフィルタや地形抽出アルゴリズムの活用

PDAL には ground segmentation（SMRF, PMF など）の地面抽出アルゴリズムがある。
オリジナルの Class=2 に依存しない「再分類」フローを構築可能か検証。

より高品質なDEMを生成できる可能性がある。

② 自動DEM生成フレームワーク化

グラウンド抽出
TIN補間
モザイク

を自動で行うバッチ処理を構築する。

③ 内挿補間方法の選択肢拡大

triangulate（TIN）以外にも
- IDW
- 最近隣法
- Kriging
- 平均法

などの内挿補間の手法との比較評価を行い、地形に応じた最適手法を選択できるようにする。

④ 品質評価（RMSE等）の導入

国土地理院の基準（DM・DSM の精度評価）に基づき、統計的な誤差評価を行う。
これにより生成DEMの品質を客観的に説明可能となる。

⑤ 0.5m → 0.25m DEM の比較検証

点群密度から見て
- 0.25mが妥当か
- 0.5mで十分か

の比較検証を行う。

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up