Claude Codeで点群の出来形をLASから計算 — 設計面との離れと点密度を判定

出来形管理の面管理では、設計面と点群の各点との離れを計算し、規格値に対する超過の有無と点密度で合否を見ます。計算そのものは標高の引き算と集計です。ところが入力がLASで数百万点になると、手元の表計算では扱えず、帳票ソフトの結果を別経路で確かめる手段も欲しくなります。ここでは、LASを読んで離れの分布・棄却点・点密度を出すPythonスクリプトをClaude Codeに書かせ、答えが分かっている点群で検算する流れを示します。

面管理が点群に求める3つの数字

国土交通省の「３次元計測技術を用いた出来形管理要領（案）」は令和８年３月版が掲載されています。面管理は、1m間隔以下（1点/㎡以上）の点密度で出来形を計測し、3次元設計データと各ポイントとの離れから良否を面的に判定する手法と定義されています。

要領が出来形管理資料に求める出力は次のとおりです。

• 平均値、最大値、最小値（棄却点を除く）
• データ数（計測密度の下限値と評価面積から計算した値以上であること）
• 評価面積
• 棄却点数（規格値を外れた点の数。データ数の0.3％以内）

標高較差は、各ポイントの標高値と、同じ平面座標にある設計面上の設計標高値との差として算出します。土工では、法肩や法尻などの変化点から水平方向に±50mm以内の計測点を評価から除外できるとも書かれています。このスクリプトでは単純化のため、平らな天端だけを対象にします。

様式の作成例（道路土工の路体盛土工、天端）では、標高較差の平均値の規格値が±50mm、個々の最大・最小の規格値が±150mmです。ここに挙げた数字は様式の作成例に載っている値で、実際の工事では、国土交通省の施工管理基準及び規格値（案）の工種区分に従って規格値を決めます。スクリプトは規格値を定数として外に出しておくのが前提です。

環境と、答えが分かっている点群を作る

LASの読み書きにはlaspyを使います。laspyはLAS・LAZを読み書きし、座標をNumPy配列で扱えるライブラリです。座標のX・Y・Zは、スケールとオフセットを掛けた後の倍精度の値としてlas.xなどで取り出せます。

python3 -m venv .venv
.venv/bin/pip install laspy numpy
.venv/bin/python -c "import laspy; print(laspy.__version__)"

この手順で入った版は2.7.0でした。

実測の点群を渡す前に、答えを自分で決めた合成点群で検算します。25m×40mの天端を縦断勾配1％の設計面とし、標高の誤差を平均−11mm・標準偏差20mmの正規乱数にします。さらに2か所に細工をします。

• 4m×5mの範囲を設計面より140mm掘りすぎた窪みにする
• 3m×3mの範囲は点を抜いて欠測にする

import numpy as np, laspy
rng = np.random.default_rng(7)
W, L = 25.0, 40.0                      # 幅25m x 延長40m = 1000 m2
n = 4200                                # 約4点/m2
x = rng.uniform(0, W, n) + 1000.0
y = rng.uniform(0, L, n) + 2000.0
# 設計面: 縦断勾配 1%、基準 100.000m
z_design = 100.0 + 0.01 * (y - 2000.0)
dz = rng.normal(-0.011, 0.020, n)       # 平均 -11mm、標準偏差 20mm
dz[(x > 1010) & (x < 1014) & (y > 2020) & (y < 2025)] += -0.140  # 掘りすぎの窪み
keep = ~((x > 1003) & (x < 1006) & (y > 2005) & (y < 2008))     # 欠測 3m x 3m
x, y, dz = x[keep], y[keep], dz[keep]
z = 100.0 + 0.01 * (y - 2000.0) + dz
hdr = laspy.LasHeader(point_format=3, version="1.2")
hdr.offsets = [1000.0, 2000.0, 100.0]
hdr.scales = [0.001, 0.001, 0.001]
las = laspy.LasData(hdr)
las.x, las.y, las.z = x, y, z
las.write("survey.las")
print(len(x), "points written")

.venv/bin/python make_cloud.py

出力は4161 points writtenでした。点数は乱数の種で決まるので、種を変えると変わります。

間引き・離れ・棄却点を1本のスクリプトにする

判定のスクリプトは次の3段です。

1. 1mメッシュごとに、離れが中央値に最も近い実計測点を1点だけ残す
2. 設計面との標高較差を点ごとに計算する
3. 規格値と照らして棄却点・平均・最大・最小・点密度を集計する

間引きでは座標を書き換えません。要領は、点群密度の変更方法について、路面切削工と河川浚渫工の記述で「座標値を変更するような処理をとってはならない」とし、鉛直方向の最下点や中央値を抽出するのはよいが平均処理はいけないと例示しています。メッシュ内の平均標高を新しい点にする処理は避け、実在する点を選ぶ形にしています。

import sys
import numpy as np
import laspy
 
las = laspy.read(sys.argv[1])
area = float(sys.argv[2])            # 評価面積 (m2)
MEAN_LIM, EACH_LIM = 0.050, 0.150    # 平均 ±50mm / 個々 ±150mm (単位 m)
 
xyz = np.asarray(las.xyz)
x0, y0 = xyz[:, 0].min(), xyz[:, 1].min()
# 設計面: 縦断勾配 1%、起点標高 100.000m の平面
dz = xyz[:, 2] - (100.0 + 0.01 * (xyz[:, 1] - 2000.0))
 
# 1m メッシュごとに、離れが中央値に最も近い実計測点を1点だけ残す
cell = (np.floor(xyz[:, 0] - x0).astype(int) * 100000
        + np.floor(xyz[:, 1] - y0).astype(int))
order = np.argsort(cell, kind="stable")
edges = np.r_[0, np.flatnonzero(np.diff(cell[order])) + 1, len(cell)]
keep = []
for a, b in zip(edges[:-1], edges[1:]):
    idx = order[a:b]
    keep.append(idx[np.argmin(np.abs(dz[idx] - np.median(dz[idx])))])
d = dz[np.array(keep)]
 
rej = np.abs(d) > EACH_LIM           # 棄却点 = 個々の規格値を外れた点
ok = d[~rej]
need = int(np.ceil(area))            # 1点/m2 × 評価面積
print(f"raw points   : {len(dz)} ({len(dz) / area:.2f} pts/m2)")
print(f"thinned      : {len(d)} (required >= {need})")
print(f"mean         : {ok.mean() * 1000:+.0f} mm (limit ±50)")
print(f"max / min    : {ok.max() * 1000:+.0f} / {ok.min() * 1000:+.0f} mm"
      " (limit ±150)")
print(f"rejected     : {rej.sum()} pts = {rej.mean() * 100:.2f}%"
      f" (limit 0.3% = {int(len(d) * 0.003)} pts)")
print(f"within 80%   : {(np.abs(d) <= 0.8 * EACH_LIM).mean() * 100:.1f}%")
print(f"within 50%   : {(np.abs(d) <= 0.5 * EACH_LIM).mean() * 100:.1f}%")

.venv/bin/python check.py survey.las 1000

第2引数の1000は評価面積（㎡）で、合成点群では25m×40mから決まります。実測の点群では、設計データの管理範囲から求めた面積を渡します。点群の外形から面積を逆算すると、欠測の穴が面積に含まれず、点密度の不足を見逃します。

合成点群での出力は次のとおりです。

raw points   : 4161 (4.16 pts/m2)
thinned      : 974 (required >= 1000)
mean         : -13 mm (limit ±50)
max / min    : +37 / -146 mm (limit ±150)
rejected     : 10 pts = 1.03% (limit 0.3% = 2 pts)
within 80%   : 98.2%
within 50%   : 98.0%

出力から読み取れる不合格の理由

平均は−13mmで規格値の範囲内です。注入した平均−11mmに近く、窪みの分だけ負側へ寄っています。最大・最小も範囲内ですが、これは棄却点を除いた値です。

不合格になる理由は2つあります。1つ目は棄却点で、10点（1.03％）が±150mmを外れ、上限の2点を超えています。窪みの部分は誤差の平均が約−151mmなので、半数近くが規格値の外へ出る計算になります。2つ目は点密度で、間引き後に974点しか残らず、1000㎡×1点/㎡の要件に26点足りません。欠測の穴は9㎡ほどで、残りの空きメッシュは乱数でたまたま点が落ちなかった所です。

間引く前の4161点を数え直すと、±150mmを外れる点は36点で、すべて窪みの範囲に入っていました。窪みの中の75点の平均は約−150mmです。間引きでメッシュごとに1点へ絞ると、この36点が10点に減ります。棄却点の数は点密度に左右されるため、間引きの前後で何点になるかを両方控えておくと、あとで説明がつきます。

ここで見るべきは、生の点数が4.16点/㎡あるのに判定は不合格という点です。密度は全体平均ではなく、評価の単位となるメッシュごとの充足で決まります。平均の標高較差が良くても、窪みと欠測が別の理由で不合格を作ります。

要領の出来形管理図表には、離れの規格値に対する割合を−100％〜＋100％のヒートマップで平面図に示す分布図が含まれ、±50％と±80％の前後が区別できる色で表示することとされています。上のスクリプトのwithin 80%とwithin 50%は、この色分けに対応する割合を数値で出したものです。分布図の描画はここでは扱いません。

Claude Codeに任せる範囲と、検算の仕掛け

点群の座標や現場の名称は、外に出したくない情報であることが多いデータです。Claude Codeには点群そのものを読ませず、スクリプトの作成と実行結果の確認を任せる形が合います。LASは大きなバイナリなので、ファイルを直接開かせない指示をCLAUDE.mdに書いておきます。

# 点群の出来形チェック
 
- LAS/LAZは直接開かない。必ず check.py の出力だけを読む
- 規格値は check.py 冒頭の定数。要領や基準の値を推測で書かない
- 間引きは実在点の選択のみ。平均や補間で座標を作らない
- 変更後は合成点群 survey.las で実行し、出力を貼って報告する

CLAUDE.mdの長さは200行以内が目安とされています。長くなったら、CLAUDE.mdが長すぎるという警告の記事にあるように分割します。

スクリプトの実行は、権限ルールで許可しておくと毎回の確認が減ります。ルールはdeny、ask、allowの順に評価され、最初に一致したものが結果を決めます。次のように、実行を許すコマンドを絞って書きます。

{
  "permissions": {
    "allow": [
      "Bash(.venv/bin/python check.py *)",
      "Bash(.venv/bin/python make_cloud.py)"
    ]
  }
}

検算のやり方は、危険物の倍数を扱った指定数量の倍数計算の記事と同じ型です。答えが決まった入力を作り、出力が決めた答えと合うまで直させます。合成点群では、次の点を目安に見ます。

• 平均が注入した−11mmの近くに出るか
• 窪みの範囲が棄却点として数えられるか
• 欠測の穴の分だけ、点密度が足りなくなるか

数値の食い違いを検算で見つける発想は、有限要素法の入力ファイルを手計算で確かめたCalculiXの片持ち梁の記事にも通じます。

つまずきやすい点

LASの座標は整数で保存されている。laspyのX・Y・Zはスケールとオフセットを掛ける前の整数で、掛けた後の倍精度がx・y・zです。ミリ単位の差を見る計算で整数側を取ると、単位がずれます。ヘッダのスケールが0.001より粗いLASは、そもそもミリ単位の差を表せません。

点群全体を一度に読む前に、点数を見る。laspyのドキュメントには、ヘッダの点数が1億点未満のときだけ全体を読み込む例があります。大きいファイルは、laspy.openでヘッダを先に見て、チャンク読みに切り替えます。

PDALで先に外形を確かめる。PDALをインストールしてあれば、pdal info --summary survey.lasで点数や範囲を、--boundaryで点群の外形を確認できます。点群が設計範囲を覆っているかの目視確認に向きます。評価面積そのものは設計データから決めます。

規格値は工種で変わる。平均値の規格値は、路体盛土工の天端、法面、舗装、河川浚渫工で別々です。上のスクリプトは天端の標高較差だけを見ます。水平較差を含む法面の判定や、法肩・法尻の±50mm除外は組み込んでいません。

提出物の代わりにはならない。要領が求める出来形管理図表には、帳票の様式と分布図が含まれます。このスクリプトは手元の確認用です。提出する資料は、要領の様式に沿って作ります。

次に足すと役立つ拡張

既存のスクリプトに足していける材料は次のとおりです。

• 標高較差を平面にプロットし、ヒートマップのPNGを出す
• 設計面を平面でなく、LandXMLなどの3次元設計データから読む
• 法面を別の規格値で判定する
• 棄却点の座標をCSVに出し、現場で確認する箇所を絞る

どれも、合成点群に対応する細工を足してから実装させると、結果を検算できます。