LP測量グリッドデータをラスタ化

対象読者
LP測量成果をGISソフトで扱う方
python ユーザ

実行環境
Windows10
Python 3.9.5
geopandas 0.9.0
pandas 1.3.1
rasterio 1.2.7

はじめに

xyz 座標の記載されたテキストファイルから、z 座標を格納した GeoTIFF データを作成します。

入力ファイルは、航空レーザ測量（LP測量）等の成果品として一般的なグリッドデータ¹ を想定します。

x、y、z 座標がひたすら並んだテキストデータです。
（座標値以外の項目を含む場合もあります）

平面的で広域な GIS 処理であれば、las 形式の LiDAR や高密度の TIN サーフェス等のような大容量データよりも、標高値をもたせたラスタデータの方がはるかに軽量で扱いやすいです。

GDAL で事足りるという話かもしれないし、GIS ソフトで変換できるという話かもしれないですが、ここでは python ライブラリの geopandas と rasterio を使っていきます。

処理の流れ

1. テキストファイルを pandas の DataFrame で読込

2. ポイント xy を geometry として geopandas の GeoDataFrame に変換

3. アファイン変換行列を設定

4. numpy の ndarray に z を格納

5. GeoTIFF 形式で出力

アファイン変換について

知っている人には説明不要ですが、線形代数という数学分野をなるべく避けて生きてきた身としては、頭の痛い話です。

wikipedia には「一般に、アフィン変換は線型変換（回転、拡大縮小、剪断（せん断））と平行移動の組み合わせである。」と載っています。

公共測量であれば xy 座標は平面直角座標系になりますが、z 座標を2次元の numpy ndarray に格納して扱うため、その行列の並びと平面直角座標系の対応関係が必要になります。

ここで扱っているアファイン変換は、ピクセル座標から平面直角座標系への変換になります。

GIS ソフトでラスタをジオリファレンスすることがありますが、その処理がラスタのアファイン変換に該当します。
（アファイン変換以外の変換もあります）

サンプル

兵庫県_全域_標高ラスター ² の DEM データを使います。

テキスト形式と TIFF 形式の両方をダウンロードできるので、テキストから変換した GeoTIFF とダウンロードした TIFF を最後に比較したいと思います。

NoData（海域等）のある図郭でテストしたいので、下記の図郭を対象にします。

淡路島北端の明石海峡大橋を含む図郭です。

ダウンロードしたzipファイルを解凍して、以下の３ファイルをフォルダ 05pf15 にいれます。

DEM_05PF151_1g.txt
DEM_05PF153_1g.txt
DEM_05PF154_1g.txt

DEM_05PF151_1g.txt の先頭5行は、こんな感じです。

60000.500 -154499.500 0.00
60001.500 -154499.500 0.00
60002.500 -154499.500 0.00
60003.500 -154499.500 0.00
60004.500 -154499.500 0.10

スペース区切りで xyz 座標が並んでいます。

標高0.00m が並んでいますが、おそらく沿岸部でしょう。

単位m、TP標高、解像度1.0mです。

座標系は平面直角座標系第5系（JGD2000）で、EPSG コードは 2447 です。

コード

import glob
import os.path

import geopandas as gpd
import pandas as pd
import rasterio
import rasterio.transform
import rasterio.features

def get_shape(xmin, xmax, ymin, ymax, res):
	"""
	入力：データ範囲と解像度
	出力：行列サイズ(ny, nx)
	"""
	ny = int(((ymax-ymin)+res) / res)
	nx = int(((xmax-xmin)+res) / res)
	shape = (ny, nx)
	return shape

def get_trans(xmin, ymax, res):
	"""
	入力：北西端のセルの中心座標と解像度
	出力：アファイン変換行列
	"""
	west = xmin - res/2
	north = ymax + res/2
	trans = rasterio.transform.from_origin(
		west, north, res, res
	)
	return trans

# フォルダ
fold = '05pf15'
# 区切り文字
sep = ' '
# 解像度m
res = 1.0
# EPSGコード
crs = 'EPSG:2447'
# NoData値
fill = -9999

files = glob.glob(os.path.join(fold, '*.txt'))

for fp_input in files:
	filename = os.path.splitext(os.path.basename(fp_input))[0]
	fp_output = os.path.join(fold, filename + '.tif')
	# 読込
	df = pd.read_csv(
		fp_input,
		sep=sep,
		names=['x', 'y', 'z'],
		usecols=[0, 1, 2]
		)
	# point 生成
	df['geometry'] = gpd.points_from_xy(df['x'], df['y'])
	# GeoDataFrame に変換
	gdf = gpd.GeoDataFrame(df)
	# 対象範囲を取得
	desc = gdf.describe()
	xmax = desc.loc['max', 'x']
	xmin = desc.loc['min', 'x']
	ymax = desc.loc['max', 'y']
	ymin = desc.loc['min', 'y']
	# 行列サイズ
	shape = get_shape(xmin, xmax, ymin, ymax, res)
	# アフィン変換行列
	trans = get_trans(xmin, ymax, res)
	# np.ndarrayとして標高を格納
	arr = rasterio.features.rasterize(
		[(g, z) for g, z in zip(gdf['geometry'], gdf['z'])],
		out_shape=shape,
		fill=fill,
		transform=trans,
		dtype='float32'
	)
	# 書込
	with rasterio.open(
		fp_output, 'w',
		driver='GTiff',
		height=arr.shape[0],
		width=arr.shape[1],
		count=1,
		dtype='float32',
		crs=crs,
		transform=trans,
		nodata=fill
	) as f:
		f.write(arr, 1)

get_trans 関数内で rasterio.transform.from_origin を使ってアファイン変換行列を生成しています。
（from_origin 以外にも transform モジュールにはアファイン変換行列を生成する関数がいくつあります。）

続いて、rasterio.features.rasterize で numpy ndarray に z 座標を格納しています。

設定している NoData 値は -9999 としていますが、この値はソースに値がない場合に numpy ndarray に格納される値で、出力時には -9999 が NoData に置き換えられます。

結果

作成した GeoTIFF ファイルとダウンロードした Tiff ファイルを、QGIS で開いて比較します。

図郭 05PF151 （図郭の北西範囲）を見てみます。

作成した GeoTIFF（DEM_05PF151_1g.tif）の段彩図

ダウンロードした TIFF（alt_05PF151.tif）の段彩図

段彩図を見てわかる通り、ダウンロードした TIFF データには海域にもデータがあり、セル値は 0.00 となっていました。

作成した GeoTIFF データはコードで指定したとおり、値のない場所は NoData となっています。

陸地の色の付き方は全く同じです。

次にラスタ情報を見てみます。

作成した GeoTIFF（DEM_05PF151_1g.tif）のプロパティ

ダウンロードした TIFF（alt_05PF151.tif）のプロパティ

ラスタ情報からは、作成した GeoTIFF データには座標参照系が定義されていることが分かります。

データのある範囲のみで作成しているため、領域や画像サイズも異なります。

データタイプやピクセルの大きさは一致しています。

ということで、

1. ラスタに座標系を定義したい

2. 欠損セルを NoData 設定したい

といった場合は、今回紹介した方法が有用です。

最後に

rasterio を使ってテキスト形式のグリッドデータを GeoTIFF 形式に変換しました。

グリッドデータの仕様は測量会社やデータ提供元によって異なりますが、xyz の並び順や区切り文字が変わっても、DataFrame のコンストラクタの引数を変えるだけで対応できます。

また、図郭コードから逆算して図郭範囲を割り出したりする手間もありません。

なお本コードは、グリッドではないデータ（オリジナルやグラウンド）には対応していないのでご注意ください。

1. 航空レーザ測量製品仕様書応用スキーマ ↩

2. 兵庫県_全域_標高ラスター／CS立体図（2010年度～2018年度）、兵庫県 ↩