1up noob: 9月 2021

国土地理院の基盤地図情報（数値標高モデル）などのDEMデータをその他のデータと組み合わせてmatplotlibでプロットしたい時に使えそうなプログラムが出来たので公開します。
例で使用したDEMデータは、国土地理院からダウンロードしたデータを予め結合したGeoTIFFを使用しています。

また、GeoTIFFの読み込みには、GDALを使用しています。

DEMデータの結合方法やGDALの導入については、色々なサイトで公開されているので、検索してください。本記事では、割愛します。

matplotlibでの標高塗り分けが本記事の中心部分となります。matplotlibのterrainでの塗り分けで満足出来ない場合にご活用いただければ幸いです。

スクリプト例

	import requests
	import numpy as np
	from osgeo import gdal
	from scipy.interpolate import griddata
	import matplotlib.pyplot as plt
	from matplotlib.colors import Normalize
	from matplotlib.colors import ListedColormap

	# 標高モデル
	ds = gdal.Open('./merge.tif', gdal.GA_ReadOnly)
	elevation = ds.GetRasterBand(1).ReadAsArray()

	# 切り出したい領域を設定
	new_x_coord = np.linspace(139.5, 140.5, 910)
	new_y_coord = np.linspace(35, 35.5, 550)
	xx, yy = np.meshgrid(new_x_coord, new_y_coord)
	nrows, ncols = elevation.shape
	x0, dx, dxdy, y0, dydx, dy = ds.GetGeoTransform()
	x1 = x0 + dx * ncols
	y1 = y0 + dy * nrows
	ele_xx, ele_yy = np.meshgrid(np.linspace(x0, x1, ncols), np.linspace(y1, y0, nrows))
	ele_zz = griddata(points=np.stack([ele_xx.reshape(-1,), ele_yy.reshape(-1,)], 1), values=np.flipud(elevation).reshape(-1,), xi=(xx, yy), method='nearest')

	# カラーマップを取得して、設定
	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]
	newcmp = ListedColormap(np.array(clist)[:, 1:]/255)

	fig = plt.figure(figsize=(4, 3), dpi=300)
	plt.rcParams['font.size'] = 5
	ax1 = plt.subplot(111)
	# 標高モデル
	ax1.pcolormesh(xx, yy, np.where(ele_zz > -9999, ele_zz, np.nan), norm=Normalize(vmin=0, vmax=400), cmap=newcmp, shading='auto')
	ax1.contour(xx, yy, ele_zz, levels=[0, 100, 200, 300, 400], colors=['black'], linewidths=0.1)

	# 目盛設定
	ax1.set_xticks(np.linspace(139.5, 140.5, 11))
	xlabels = [str(round(label, 3)) + '°E' for label in np.linspace(139.5, 140.5, 11).tolist()]
	ax1.set_xticklabels(xlabels)
	ax1.set_xlim(139.5, 140.5)
	ax1.set_yticks(np.linspace(35, 35.5, 6))
	ylabels = [str(round(label, 3)) + '°N' for label in np.linspace(35, 35.5, 6).tolist()]
	ax1.set_yticklabels(ylabels)
	ax1.set_ylim(35, 35.5)
	plt.show()

view raw matplotlib_dem.py hosted with ❤ by GitHub

出力例

matplotlibのカラーマップに好みのものがなかったので、ArcGISっぽい?カラーマップを外部サイトから読み込んで使用しています。

カラーマップのファイルは、ローカルから直接読み込んで使用した方が良いと思いますので、上記のコードを利用される際はファイルをダウンロードして、コードを書き換えてください。

	import requests
	import numpy as np
	from matplotlib.colors import ListedColormap
	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]
	newcmp = ListedColormap(np.array(clist)[:, 1:]/255)
	newcmp

view raw create_schwarzwald_cmap.py hosted with ❤ by GitHub

関連記事

1up noob

2021/09/14

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

スクリプト例

出力例