matplotlibで国土地理院標高タイルから3D地形図を描いてみる

前回、国土地理院の地形データから等高線を描いた時の要領で

memomemokun.hateblo.jp

今回は、国土地理院標高タイルから3D地形図を描いてみました。

場所は、奥穂高岳、涸沢岳、北穂高岳、前穂高岳、西穂高岳などの峰々からなる穂高連峰から、涸沢カール、屏風岩ぐらいにかけての地形図。

だいたい、奥穂高岳、前穂高岳方向から涸沢カール、屏風岩方向を見下ろす感じの地形図ですかね。

涸沢カール - Wikipedia

北穂高岳側から見た涸沢カール（左下：涸沢ヒュッテ、正面：前穂高岳）

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

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

２．穂高連峰辺りの「レイヤ番号/X/Y」を確認

f:id:memomemokun:20181102091913j:plain

穂高連峰が周囲を囲む涸沢カール辺りの「レイヤ番号/X/Y」は「13/7228/3208」なので、以下のアドレスでDEMデータが取得できるので、

http://cyberjapandata.gsi.go.jp/xyz/dem/13/7228/3208.txt

地理院地図3Dから、今回描画した辺りの航空写真3D

f:id:memomemokun:20181102101937j:plain

あとは、matplotlibを使い3Dで地形図を描いてみると！

◆ ワイヤーフレーム（plot_wireframe）

f:id:memomemokun:20181102093940j:plain

ソース

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

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


#地形データを読み込む
url = 'http://cyberjapandata.gsi.go.jp/xyz/dem/13/7228/3208.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=(13, 10), dpi=80, facecolor='w', edgecolor='k')

# プロット中の軸の取得。gca はGet Current Axes の略。
ax = fig.gca(projection='3d')
# 横を1/0.8=1.25倍長く設定
ax.set_aspect(0.8, adjustable='box') 
# 上高地の遙か上空ぐらいから前穂高越しに地形を見下ろす感じに視点を設定
ax.view_init(70, -67)

#wireframeを描く
ax.plot_wireframe(X, Y, Z, rstride=1, cstride=2, linewidth=1)

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

等高線（contour）

25m間隔で等高線を引いてみる。

f:id:memomemokun:20181102125252j:plain

ソース

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

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


#地形データを読み込む
url = 'http://cyberjapandata.gsi.go.jp/xyz/dem/13/7228/3208.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')

# プロット中の軸の取得。gca はGet Current Axes の略。
ax = fig.gca(projection='3d')
# 横を1/0.8=1.25倍長く設定
ax.set_aspect(0.8, adjustable='box') 
# 上高地の遙か上空ぐらいから前穂高越しに地形を見下ろす感じに視点を設定
ax.view_init(70, -67)

#標高25m間隔で等高線を描く
elevation = range(1500,3500,25)
cont = plt.contour(X, Y, Z, levels=elevation, cmap='hot_r')

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

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

表面を塗りつぶす（plot_surface）

f:id:memomemokun:20181102131900j:plain

ソース

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

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


#地形データを読み込む
url = 'http://cyberjapandata.gsi.go.jp/xyz/dem/13/7228/3208.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')

# プロット中の軸の取得。gca はGet Current Axes の略。
ax = fig.gca(projection='3d')
# 横を1/0.8=1.25倍長く設定
ax.set_aspect(0.8, adjustable='box') 
# 上高地の遙か上空ぐらいから前穂高越しに地形を見下ろす感じに視点を設定
ax.view_init(70, -67)

#標高に応じて塗りつぶす
surf = ax.plot_surface(X, Y, Z.values, cstride=1, rstride=1, cmap='hot_r', 
                      linewidth=0, antialiased=False, shade=False)

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

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