【Python】DEMデータをvtkに変換して、ParaViewで表示する

ParaViewは、2次元・3次元データの可視化のためのオープンソースソフトウェアです。可視化機能が強力でフリーで使えるため人気があります。
Pythonの可視化ライブラリには、matplotlibやplotlyなどがありますが、基本的には2次元までデータの可視化向きで、3次元データの可視化にはParaViewを使った方が動作も軽いです。

本記事では、以前の記事と同様に、国土地理院の基盤地図情報（数値標高モデル）のDEMデータから変換したGeoTIFFを用いました。DEMデータの結合方法やGDALの導入については、本記事でも割愛します。
また今回は、3次元表示に当たり軸の単位が揃っていた方が都合が良いため、GeoTIFFへの変換時にUTM座標に変換したデータを用いました。

実行環境：MacBook Air（M2、2022、メモリ：16GB）

ParaViewでは、GeoTIFF形式のDEMデータも読み込めるようですが、下記の例では、DEMデータの指定範囲での切り出しや解像度変更、vtk(Visualization Toolkit)形式への変換を試しました。

スクリプト例：

	import numpy as np
	import pyvista as pv
	from osgeo import gdal
	from pyproj import Transformer
	from scipy.interpolate import griddata
	# UTM座標に変換済の標高モデルを読み込む
	ds = gdal.Open('./merge_utm.tif', gdal.GA_ReadOnly)
	zz = ds.GetRasterBand(1).ReadAsArray()
	nrows, ncols = zz.shape
	x0, dx, dxdy, y0, dydx, dy = ds.GetGeoTransform()
	x1 = x0 + dx * ncols
	y1 = y0 + dy * nrows
	xx, yy = np.meshgrid(np.linspace(x0, x1, ncols), np.linspace(y1, y0, nrows))
	# 緯度経度(JGD2011)=>UTM座標系(JGD2011)
	epsg6668_to_epsg6690 = Transformer.from_crs('epsg:6668', 'epsg:6690')
	# 緯度経度で切り出し範囲を指定する
	lim_x0, lim_y0 = epsg6668_to_epsg6690.transform(34.685, 133)
	lim_x1, lim_y1 = epsg6668_to_epsg6690.transform(35.96, 135)
	# 切り出した範囲の解像度を指定する
	new_ncols = 1400
	new_nrows = 1800
	# 指定した設定でメッシュを構築する
	new_xx, new_yy = np.meshgrid(np.linspace(lim_x0, lim_x1, new_ncols), np.linspace(lim_y1, lim_y0, new_nrows))
	new_zz = griddata(points=np.stack([xx.reshape(-1,), yy.reshape(-1,)], 1), values=np.flipud(zz).reshape(-1,), xi=(new_xx, new_yy), method='nearest')
	# geotiffへの変換の際に海域の標高を-9999としているため、NaNに置換する
	new_zz[new_zz < -999] = np.nan
	# vtk形式に格納する
	points = np.stack((new_xx.ravel(order='F'), new_yy.ravel(order='F'), new_zz.ravel(order='F')), axis=1)
	grid = pv.StructuredGrid()
	grid.points = points
	grid.dimensions = 1, new_nrows, new_ncols
	grid['Elevation'] = new_zz.ravel(order='F')
	# NaNを含むため、Binary形式で保存する
	grid.save('./test_DEM.vtk', binary=True)

view raw dem_vtk_convert.py hosted with ❤ by GitHub

出力例：Scaleの設定からZ軸方向を10倍に強調

参考：

ParaViewはVTK ASCIIではなくBinaryでないとNaNを読み込めない - Qiita

NaNを含むデータをParaViewで表示するにはVTKライブラリを使ってVTK ASCII形式でファイルに出力すると「nan」と表示されます。このファイルをParaViewで開くと

また、以下のようにカラーマップを作成すると好きなカラーマップを追加できます。

スクリプト例：前回の記事と同様にwiki-schwarzwald-contを設定

	import requests
	pg = requests.get('http://soliton.vm.bytemark.co.uk/pub/cpt-city/wkp/schwarzwald/wiki-schwarzwald-cont.pg')
	clist = [list(map(int, txt.split())) for txt in pg.text.splitlines()][::-1]
	with open('./wiki-schwarzwald-cont.xml', mode='w') as f:
	f.write('<ColorMaps>\n')
	f.write('<ColorMap name="wiki-schwarzwald-cont" space="RGB">\n')
	for c in clist:
	point = '<Point x="'+str(c[0]/1500)+'" o="1" r="'+str(c[1]/255)+'" g="'+str(c[2]/255)+'" b="'+str(c[3]/255)+'"/>\n'
	f.write(point)
	f.write('</ColorMap>\n')
	f.write('</ColorMaps>')

view raw paraview_colormap.py hosted with ❤ by GitHub

出力例：

関連記事

【Python】matplotlibでDEMデータを表示する