pythonで富士山の斜面勾配図を立体表示してみる

先日、pythonのnumpy.gradientで遊んでみながら、富士山の斜面勾配図ぽいものを描いてみましたが、

memomemokun.hateblo.jp

2Dで表示しても立体的に見える斜面勾配図ですので、斜面勾配図を3Dで描画したらより立体的に見えるのではと、今回は、斜面勾配図を3D表示して遊んでみました。

plot_surfaceで3D地形を描き、地形の配色は、前回作成した斜面勾配図と同様、勾配に応じてfacecolorsで3D地形を塗りつぶすというもの。

結果は以下。（配色をもう少し工夫すれば結構使えるかも）

なお、facecolorsのカラーッマップのパラメーターとし、勾配を以下の形で最小値0、最大値1に正規化して渡してます。

Min-Max Normalization式

ソースは以下

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
import matplotlib.colors as colors
from matplotlib import cm

def load_gis(z, x1, x2, y1, y2):

    a=np.array([])
    for x in range(x1, x2+1):
        b=np.array([])
        for y in range(y1, y2+1):
            print  x, y

            #地形データを読み込む
            url = 'http://cyberjapandata.gsi.go.jp/xyz/dem/{0}/{1}/{2}.txt'.format(z,x,y)
            print url
            response = requests.get(url)
            if response.status_code == 404:
                 Z = np.zeros((256, 256))
            else:
                 #標高値がない区画はとりあえず0mに置換する
                 maptxt = string.replace(response.text, u'e', u'0.0')
                 Z = pd.read_csv(StringIO(maptxt), header=None)
                 Z = Z.values
                 #Z = Z.transpose() 
            print Z
            
            if len(b)==0:
                b = Z
            else:
                b = np.append(b,Z,0)
            
    
        if len(a)==0:
            a = b
        else:
            a = np.append(a,b,1)
            
    return a
    



z=12
x1=3624
x2=3628
y1=1615
y2=1619

        
        
Z = load_gis(z, x1, x2, y1, y2)        

xlim = Z.shape[1]
ylim = Z.shape[0]

#X,Y軸のグリッドを生成
X, Y = np.meshgrid(np.linspace(0,xlim-1,xlim), np.linspace(ylim-1,0,ylim))

# 勾配を求める
 (Zy,Zx) = np.gradient(Z)

# Y軸方向の勾配を最小値0、最大値1に正規化
Zgradient_norm = (Zy-Zy.min())/(Zy.max()-Zy.min())

# 立体地図を作成する
fig = plt.figure(figsize=(20, 10), dpi=72, facecolor='w', edgecolor='k')
plt.subplots_adjust(left=0, right=1, bottom=0, top=1)
ax = fig.add_subplot(111, projection='3d')
ax.view_init(elev=65, azim=-70)

# 横を1/0.5=2倍長く設定
ax.set_aspect(0.5, adjustable='box') 

# 配色は勾配に応じて地形を塗りつぶす 影(shade)とパッチの境界線(antialiased)は無効化
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0, shade=False, antialiased=False, facecolors=plt.cm.binary(Zgradient_norm), alpha=0.80)


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

今回の方法を拡張すれば平面航空写真をカラーマップにマッピングしてfacecolorsにパラメーターとして渡せば平面航空写真を立体表示できるかもしれませんね。。