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測量）等の成果品として一般的なグリッドデータ[^1] を想定します。

[^1]: [航空レーザ測量製品仕様書応用スキーマ](https://psgsv2.gsi.go.jp/koukyou/public/seihinsiyou/data/%E2%91%A1%E8%88%AA%E7%A9%BA%E3%83%AC%E3%83%BC%E3%82%B6%E6%B8%AC%E9%87%8F_%E5%BF%9C%E7%94%A8%E3%82%B9%E3%82%AD%E3%83%BC%E3%83%9E.pdf)

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

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

GDAL で事足りるという話かもしれないし、GIS ソフトで変換できるという話かもしれないですが、ここでは python ライブラリの [geopandas](https://geopandas.org/) と [rasterio](https://rasterio.readthedocs.io/en/latest/index.html) を使っていきます。


## 処理の流れ

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

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

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

1. numpy の ndarray に z を格納

1. GeoTIFF 形式で出力


## アファイン変換について

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

[wikipedia](https://ja.wikipedia.org/wiki/%E3%82%A2%E3%83%95%E3%82%A3%E3%83%B3%E5%86%99%E5%83%8F) には「一般に、アフィン変換は線型変換（回転、拡大縮小、剪断（せん断））と平行移動の組み合わせである。」と載っています。

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

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

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


## サンプル

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

[^2]: [兵庫県_全域_標高ラスター／CS立体図（2010年度～2018年度）、兵庫県](https://www.geospatial.jp/ckan/dataset/2010-2018-hyogo-geo-dem)

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

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

![図郭.JPG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1409923/b6a666d5-e529-844e-9179-6f58b23c4bb2.jpeg)

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

ダウンロードした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 です。


## コード

```py
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](https://rasterio.readthedocs.io/en/latest/api/rasterio.transform.html)  モジュールにはアファイン変換行列を生成する関数がいくつあります。）

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

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



## 結果

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

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

`作成した GeoTIFF（DEM_05PF151_1g.tif）の段彩図`
![レイアウト_DEM_05PF151_1g.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1409923/cad4613d-0a78-2f1f-0c64-a9232ce6accb.png)

`ダウンロードした TIFF（alt_05PF151.tif）の段彩図`
![レイアウト_alt_05PF151.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1409923/c80e5ea7-d14b-0e51-0fde-8654f19df8bf.png)

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

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

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


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

`作成した GeoTIFF（DEM_05PF151_1g.tif）のプロパティ`
![DEM_05PF151_1g_info.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1409923/14ec2984-2c25-58c8-af85-771e63fc69d0.png)


`ダウンロードした TIFF（alt_05PF151.tif）のプロパティ`
![alt_05PF151_info.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1409923/3fa4300a-918c-d4c9-e7f1-db9cdf00bb2f.png)

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

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

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

ということで、

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

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

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

## 最後に

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

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

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

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