点群からオルソ画像とDSMを作る

このページでは「Pythonプログラムを用いて点群処理を行い、オルソ画像やDSMを作る方法」について解説します。

現在工事中：スタート（2024.06.05）

写真測量ではSfMにより点群やDSM、オルソ画像を生成できますが、研究活動をする上では生成されたデータの生成過程を把握しておく必要があります。またLiDARや衛星画像など、他の計測プラットフォームより得られたデータと重ね合わせを行うためには、自作プログラムにより自分で生成できる方が好ましいです。そこでこのページではPythonのプログラムを用いてDEMとオルソ画像を生成する過程を解説します。

[ 注意書き ]
・参考 Program の環境（python 3.8）
・以下のプログラムをコピーする時は、全角やスペースが含まれていないか注意してください。

1-1. 点群データ

2. Python プログラム

アルゴリズム

（１）モジュールのインポート

（２）I/O の設定

（３）ランダムな点群をボクセルに格納して均質化

（４）水平方向に値のないセルを内挿

（５）DSM画像とオルソ画像を保存

（１）モジュールのインポート

import cv2

import pandas as pd

import numpy as np

from tqdm import tqdm

from osgeo import gdal

from osgeo import osr

from osgeo import gdal_array

import matplotlib.pyplot as plt

# 処理時間計測

import time

# メモリの解放

import gc

#Jupyterでインライン表示するための宣言

%matplotlib inline

（２） I/O の設定

# 対象観測シーンID

SceneID = "Saoka-20220720-1422-X5S"

# パスの設定

Dir1 = SceneID+ "/"

Dir2 = Dir1 + "IntermediateData/"

# I/O

IFile1 = Dir1 + SceneID + "-RandomPointCloud.txt" # 入力ファイル1 - ランダムな点群

OFile1 = Dir2 + SceneID + "-DSMImage.tif" # 出力ファイル1 - DSM画像

OFile2 = Dir2 + SceneID + "-OltImage.tif" # 出力ファイル2 - オルソ画像

# ====================================================================

# 地上解像度(単位:m / ボクセルのセルサイズにも利用)

GSD = 0.20

# ボクセルの範囲（例：日本平面直角座標4系）

xmin,xmax = 19940 , 20440

ymin,ymax = 71440 , 71940

zmin,zmax = 100 , 260

# トリミング範囲

xmin2,xmax2 = 20050 , 20350

ymin2,ymax2 = 71550 , 71850

# XYX軸ごとのセル数を算出

xsize, ysize, zsize = int(round((xmax-xmin)/GSD)),int(round((ymax-ymin)/GSD)),int(round((zmax-zmin)/GSD))

# ランダムな点群データを整形

PointCloud_pd = pd.read_csv(IFile1, sep=' ', skiprows=0, names=("x","y","z","r","g","b","nx","ny","nz"))

PointCloud_pd = PointCloud_pd.drop(["nx","ny","nz"],axis=1) # 法線は使わないのでドロップ

PointCloud_pd["xid"] = round((PointCloud_pd["x"]-xmin)/GSD,4).astype(int) # ボクセル Xid

PointCloud_pd["yid"] = round((ymax-PointCloud_pd["y"])/GSD,4).astype(int) # ボクセル yid

PointCloud_pd["zid"] = round((PointCloud_pd["z"]-zmin)/GSD,4).astype(int) # ボクセル zid

PointCloud_np = PointCloud_pd.values # NumPy に変換

# メモリを解放（点群データは大きいので消しておく）

del PointCloud_pd

gc.collect()

（３）ランダムな点群をボクセルに格納して均質化

3-1. 均質化された点群を格納する箱を作成

# 各XYZ軸に入力する値のリスト(セルの中央座標を使用する)

shift2C = GSD/2 # セルの中央座標へシフトするための値

digitsN = len(str(GSD).split('.')[1]) + 1 # 桁数

list_xv = np.round(np.arange(xmin+shift2C, xmax, GSD), digitsN)

list_yv = np.round(np.arange(ymax-shift2C, ymin,-GSD), digitsN)

list_zv = np.round(np.arange(zmin+shift2C, zmax, GSD), digitsN)

# 均質化された点群が入るボクセル

xv, yv, zv = np.zeros((zsize, ysize, xsize)), np.zeros((zsize, ysize, xsize)), np.zeros((zsize, ysize, xsize))

rv, gv, bv = np.zeros((zsize, ysize, xsize)), np.zeros((zsize, ysize, xsize)), np.zeros((zsize, ysize, xsize))

nv = np.zeros((zsize, ysize, xsize))

# 表層の点群だけが入るボクセル

xzm, yzm, zzm = np.zeros((ysize, xsize)), np.zeros((ysize, xsize)), np.zeros((ysize, xsize))

rzm, gzm, bzm = np.zeros((ysize, xsize)), np.zeros((ysize, xsize)), np.zeros((ysize, xsize))

nzm = np.zeros((ysize, xsize))

print("中央座標へシフトする距離:",shift2C," 座標の小数部桁数",digitsN)

print("セル数と入力値の数が合っているか？","x",xsize,":",len(list_xv)," y",ysize,":",len(list_yv)," z",zsize,":",len(list_zv))

3-2. ランダムな点群の均質化

# 関数：ボクセルによるランダム点群の均質化と表層のセルの抽出

def voxel(PointCloud_np):

for point in tqdm(PointCloud_np):

x, y, z = float(point[0]),float(point[1]),float(point[2])

R, G, B = int(point[3]),int(point[4]),int(point[5])

Xid, Yid, Zid = int(point[6]),int(point[7]),int(point[8])

# （１）均質化された点群：解析エリアに入っている点に絞り込む

if 0 <= Xid < xsize and 0 <= Yid < ysize and 0 <= Zid < zsize:

xv[Zid][Yid][Xid], yv[Zid][Yid][Xid], zv[Zid][Yid][Xid] = list_xv[Xid], list_yv[Yid], list_zv[Zid]

rv[Zid][Yid][Xid] += R

gv[Zid][Yid][Xid] += G

bv[Zid][Yid][Xid] += B

nv[Zid][Yid][Xid] += 1

nzm[Yid][Xid] += 1

# （２）DSM点群：標高が最大値の時に格納する

if zzm[Yid][Xid] < z:

xzm[Yid][Xid], yzm[Yid][Xid], zzm[Yid][Xid] = list_xv[Xid], list_yv[Yid], z

rzm[Yid][Xid], gzm[Yid][Xid], bzm[Yid][Xid] = R, G, B

return xzm,yzm,zzm,rzm,gzm,bzm,nzm

# プログラム実行(ボクセルによるランダム点群の均質化)

start = time.time()

xzm,yzm,zzm,rzm,gzm,bzm,nzm = voxel(PointCloud_np)

print ("プログラム実行時間:{0}".format(time.time()-start) + "[sec]")

3-3. DSMとオルソ画像（点群が無い場合は穴ができている）

# Matplotlib で確認する

plt.figure(figsize=(3, 3), dpi=150)

plt.imshow(zzm,cmap="gray")

plt.show()

# Matplotlib で確認する

img_rgb = cv2.merge((rzm.astype(np.uint8), gzm.astype(np.uint8), bzm.astype(np.uint8)))

plt.figure(figsize=(3, 3), dpi=150)

plt.imshow(img_rgb)

plt.show()

（４）水平方向に値のないセルを内挿

4-1. 内挿処理

# 値の無いセルに周辺から穴埋めした値を入れる箱

X2, Y2, Z2, N2 = np.zeros((ysize,xsize)), np.zeros((ysize,xsize)), np.zeros((ysize,xsize)), np.zeros((ysize,xsize))

R2, G2, B2 = np.zeros((ysize,xsize),np.uint32), np.zeros((ysize,xsize),np.uint32), np.zeros((ysize,xsize),np.uint32)

# 近隣から値を検索する範囲（セル数）

MaxArea = 15

def Interpolation():

for j in tqdm(range(0,ysize)):

for i in range(0,xsize):

if nzm[j][i] > 0: # 点が一つでもあったところはそのまま代入

X2[j][i], Y2[j][i], Z2[j][i] = xzm[j][i], yzm[j][i], zzm[j][i]

R2[j][i], G2[j][i], B2[j][i] = rzm[j][i], gzm[j][i], bzm[j][i]

N2[j][i] = nzm[j][i]

else: # 点が無いところは近傍の平均値を入れる

k = 1

Nsum = 0

while Nsum == 0: # 0以上になるまで周辺から計算する

jmin,jmax = j-k,j+k+1

imin,imax = i-k,i+k+1

jmin = 0 if jmin < 0 else jmin

jmax = ysize if jmax > ysize else jmax

imin = 0 if imin < 0 else imin

imax = xsize if imax > xsize else imax

Zsum, Zn = np.sum(zzm[jmin:jmax,imin:imax]), np.count_nonzero(zzm[jmin:jmax,imin:imax] > 0)

Nsum = np.sum(nzm[jmin:jmax,imin:imax])

if Zsum and Nsum > 0 :

Rsum, Rn = np.sum(rzm[jmin:jmax,imin:imax]), np.count_nonzero(rzm[jmin:jmax,imin:imax])

Gsum, Gn = np.sum(gzm[jmin:jmax,imin:imax]), np.count_nonzero(gzm[jmin:jmax,imin:imax])

Bsum, Bn = np.sum(bzm[jmin:jmax,imin:imax]), np.count_nonzero(bzm[jmin:jmax,imin:imax])

X2[j][i] = round(xmin + i*GSD + shift2C, 4)

Y2[j][i] = round(ymax - j*GSD - shift2C, 4)

Z2[j][i] = round(Zsum / Zn,5)

R2[j][i] = 0 if Rsum == 0 else round(Rsum / Rn,0)

G2[j][i] = 0 if Gsum == 0 else round(Gsum / Gn,0)

B2[j][i] = 0 if Bsum == 0 else round(Bsum / Bn,0)

N2[j][i] = 1

break

elif k == MaxArea:

break

else:

k += 1

return X2,Y2,Z2,R2,G2,B2,N2

X2,Y2,Z2,R2,G2,B2,N2 = Interpolation()

4-2. 内挿した後のDSMとオルソ画像

# Matplotlib で確認する

plt.figure(figsize=(3, 3), dpi=150)

plt.imshow(Z2,cmap="gray")

plt.show()

# Matplotlib で確認する

img_rgb = cv2.merge((R2.astype(np.uint8), G2.astype(np.uint8), B2.astype(np.uint8)))

plt.figure(figsize=(3, 3), dpi=150)

plt.imshow(img_rgb)

plt.show()

（５）DSM画像とオルソ画像を保存(必要な範囲へトリミングする)

5-1. トリミングしてGeoTifとして画像を保存

# トリミング範囲の設定

xmin2num,xmax2num = int(round((xmin2-xmin)/GSD, 0)),int(round((xmax2-xmin)/GSD, 0))

ymin2num,ymax2num = int(round((ymax-ymax2)/GSD, 0)),int(round((ymax-ymin2)/GSD, 0))

xsize2 ,ysize2 = int(round((xmax2-xmin2)/GSD)) ,int(round((ymax2-ymin2)/GSD))

# Nodata の地点を入れるための前処理

SetNodata = -9999

Z3 = np.where(N2==0,SetNodata, Z2)

R3 = np.where(N2==0,SetNodata, R2)

G3 = np.where(N2==0,SetNodata, G2)

B3 = np.where(N2==0,SetNodata, B2)

# DSM画像の保存

driver = gdal.GetDriverByName('GTiff')

dataset = driver.Create(OFile1,xsize,ysize,1,gdal.GDT_Float32) # 画素数と型の設定

dataset.SetGeoTransform((xmin,GSD,0,ymax,0,-GSD)) # 位置情報の設定

srs = osr.SpatialReference() # 空間参照情報

srs.ImportFromEPSG(2446) # 日本測地第４系に座標系を指定

dataset.SetProjection(srs.ExportToWkt()) # 空間情報を結合

dataset.GetRasterBand(1).SetNoDataValue(SetNodata)

dataset.GetRasterBand(1).WriteArray(Z3[ymin2num:ymax2num,xmin2num:xmax2num])

dataset.FlushCache()

dataset = None

# オルソ画像の保存

driver = gdal.GetDriverByName('GTiff')

dataset = driver.Create(OFile2,xsize2,ysize2,3,gdal.GDT_Byte) # 画素数と型の設定

dataset.SetGeoTransform((xmin2,GSD,0,ymax2,0,-GSD)) # 位置情報の設

srs = osr.SpatialReference() # 空間参照情報

srs.ImportFromEPSG(2446) # 日本測地第４系に座標系を指定

dataset.SetProjection(srs.ExportToWkt()) # 空間情報を結合

dataset.GetRasterBand(1).SetNoDataValue(SetNodata)

dataset.GetRasterBand(1).WriteArray(R3[ymin2num:ymax2num,xmin2num:xmax2num])

dataset.GetRasterBand(2).WriteArray(G3[ymin2num:ymax2num,xmin2num:xmax2num])

dataset.GetRasterBand(3).WriteArray(B3[ymin2num:ymax2num,xmin2num:xmax2num])

dataset.FlushCache()

5-2. トリミング結果を確認する

# Matplotlib で確認する

plt.figure(figsize=(3, 3), dpi=150)

plt.imshow(Z3[ymin2num:ymax2num,xmin2num:xmax2num],cmap="gray")

plt.show()

# Matplotlib で確認する

img_rgb = cv2.merge((R3.astype(np.uint8), G3.astype(np.uint8), B3.astype(np.uint8)))

plt.figure(figsize=(3, 3), dpi=150)

plt.imshow(img_rgb[ymin2num:ymax2num,xmin2num:xmax2num])

plt.show()

点群からオルソ画像とDSMを作る

目次

はじめに

1-1. 点群データ

2. Python プログラム

アルゴリズム

（１）モジュールのインポート

（２） I/O の設定

（３）ランダムな点群をボクセルに格納して均質化

3-1. 均質化された点群を格納する箱を作成

3-2. ランダムな点群の均質化

3-3. DSMとオルソ画像（点群が無い場合は穴ができている）

（４）水平方向に値のないセルを内挿

4-1. 内挿処理

4-2. 内挿した後のDSMとオルソ画像

（５）DSM画像とオルソ画像を保存(必要な範囲へトリミングする)

5-1. トリミングしてGeoTifとして画像を保存

5-2. トリミング結果を確認する

点群からオルソ画像とDSMを作る

目次

はじめに

1-1. 点群データ

2. Python プログラム

アルゴリズム

（１） モジュールのインポート

（２） I/O の設定

（３）ランダムな点群をボクセルに格納して均質化

3-1. 均質化された点群を格納する箱を作成

3-2. ランダムな点群の均質化

3-3. DSMとオルソ画像（点群が無い場合は穴ができている）

（４）水平方向に値のないセルを内挿

4-1. 内挿処理

4-2. 内挿した後のDSMとオルソ画像

（５）DSM画像とオルソ画像を保存(必要な範囲へトリミングする)

5-1. トリミングしてGeoTifとして画像を保存

5-2. トリミング結果を確認する

（１）モジュールのインポート