matplotlibで国土地理院標高タイルから等高線を描いてみる

国土地理院の地形データからpythonのグラフ描画ライブラリmatplotlibで等高線を描いてみました。

国土地理院の地形データは、国土地理院標高タイルとして、標高データ（単位はｍ）が入った 256 × 256 サイズの数値標高モデル（DEM：Digital Elevation Model) データが公開されますが、

地理院地図｜標高タイルの詳細仕様

• 1行にカンマ区切りで256個の標高値を表す数値データが格納されている。
• 1枚の標高タイルは256行からなる。 
• 数値データは対応する地図タイルのピクセル座標における標高値を表す。座標系はタイル座標と同一である。つまり、東方向がX、南方向がYになっている。 
• 標高値が存在しない画素には「e」の文字が格納されている。 
• 標高データは小数点第二位までデータとして入っている（単位はm）。

その国土地理院標高タイルのデータを使い

http://maps.gsi.go.jp/development/ichiran.html#dem

富士山周辺の地形の等高線をmatplotlibで描いてみることに。

やり方は簡単で、

１．国土地理院のタイル座標確認ページを開き

http://maps.gsi.go.jp/development/tileCoordCheck.html

２．富士山周辺の「レイヤ番号/X/Y」を確認

以下は、ズームレベル12で富士山周辺を表示した場合ですが

上を見ると、富士山山頂を含む標高タイルの「レイヤ番号/X/Y」は「12/3626/1617」なので、以下のアドレスにアクセスすると

http://cyberjapandata.gsi.go.jp/xyz/dem/12/3626/1617.txt

カンマ区切りの標高データの入ったテキストデータが取得できます。

３．標高データの入ったアドレスがわかれば、あとは簡単。

matplotlibで等高線を描いてみる

国土地理院の航空写真と比べてみると、富士山西斜面の大沢崩れとか、南斜面の宝永山の噴火口あたりの等高線が中に入り込んでいるのがわかります。

ソースはこんな感じ

import requests
from io import StringIO
import string
import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
from matplotlib import cm


#地形データを読み込む
url = 'http://cyberjapandata.gsi.go.jp/xyz/dem/12/3626/1617.txt'
response = requests.get(url)
#標高値がない区画はとりあえず0mに置換する
maptxt = string.replace(response.text, u'e', u'-0.0')
Z = pd.read_csv(StringIO(maptxt), header=None)


#X,Y軸のグリッドを生成
X, Y = np.meshgrid(np.linspace(0,255,256), np.linspace(255,0,256))

#標高100m間隔で等高線を描く
fig = plt.figure(figsize=(14, 10), dpi=80, facecolor='w', edgecolor='k')
elevation = range(0,4100,100)
cont = plt.contour(X, Y, Z, levels=elevation, cmap='spring_r')
cont.clabel(fmt='%1.1f', fontsize=14)

#ラベルをつける
cb = plt.colorbar(cont, shrink=0.5, aspect=10)

#plt.savefig('map.jpg', dpi=72)
plt.show()

標高により背景を塗りつぶして見たいときは、こんな感じ

ソース

import requests
from io import StringIO
import string
import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
from matplotlib import cm


#地形データを読み込む
url = 'http://cyberjapandata.gsi.go.jp/xyz/dem/12/3626/1617.txt'
response = requests.get(url)
#標高値がない区画はとりあえず0mに置換する
maptxt = string.replace(response.text, u'e', u'-0.0')
Z = pd.read_csv(StringIO(maptxt), header=None)


#X,Y軸のグリッドを生成
X, Y = np.meshgrid(np.linspace(0,255,256), np.linspace(255,0,256))

fig = plt.figure(figsize=(14, 10), dpi=80, facecolor='w', edgecolor='k')

#標高に応じて塗りつぶす
pcolor = plt.pcolor(X, Y, Z, cmap='hot_r')

#標高100m間隔で等高線を描く
elevation = range(0,4100,100)
cont = plt.contour(X, Y, Z, levels=elevation, cmap='binary_r')
cont.clabel(fmt='%1.1f', fontsize=14)

#ラベルをつける
cb = plt.colorbar(pcolor, shrink=0.5, aspect=10)

#plt.savefig('map.jpg', dpi=72)
plt.show()