pythonで複数の標高タイルを結合して富士山の3D立体地形図を描いてみる

先日は、国土地理院の標高データを、MySQLに一旦格納してからSQLで指定した範囲の標高データを抽出して富士山周辺の立体地形図を描いてみましたが、

memomemokun.hateblo.jp

実際にメッシュ各点を１レコードとして標高データをデータベースに突っ込んでみたところ、ズームレベルが上がるとちょっとした範囲の区画でも、数千万レーコドのデータにあっという間になってしまったので、実際問題、日本全国のデータを扱いたいとなってくるとこの方法では実用的ではないということで、

ズームレベルと範囲を指定して、指定範囲の複数の国土地理院の標高タイルを直接読み込み1タイルに結合して富士山の3D立体地形図を描いてみることに。

国土地理院の標高タイルデータについては、以下を見ていただくとして

memomemokun.hateblo.jp

ズームレベル12。

x座標3624〜3628、y座標1615〜1619の標高タイルを結合して富士山を俯瞰的に描いたのが以下。
ax.view_init(elev=80, azim=-70)

ズームレベル14 。

x座標14503〜14508、y座標6467〜6471の標高タイルを結合して富士山を描いたのが以下。
ax.view_init(elev=65, azim=-70)

ソースは以下。

load_gisという関数で指定範囲の標高タイルを読み込んで、np.appendでnumpyの配列を縦横結合して、標高データの入った１つのnumpy行列として返す形に。

ここまで、出来れば、指定の緯度経度の範囲に最適のズームレベルと、x座標とy座標の範囲を求める関数でも作れば、緯度経度の範囲を指定すれば、日本中の立体地形図を作成するなんてのも簡単に出来ますね。

ということで。

 
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.ticker as ticker
import matplotlib.colors as colors
from matplotlib.mlab import bivariate_normal

# 指定範囲の国土地理院標高タイルを読込み結合し結果を返す
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):

            #地形データを読み込む
            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:
                 #メッシュデータが無い区画は全ての点が0mのメッシュとしてとりあえず用意する
                 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()             
            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=14
x1=14503
x2=14508
y1=6467
y2=6471

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

    
# 指定範囲の標高タイルを読み込む    
Z = load_gis(z, x1, x2, y1, y2)   

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

#メッシュを間引く時はスライスする
#X=X[::2,::2]
#Y=Y[::2,::2]
#Z=Z[::2,::2]

# 立体地図を作成する
fig = plt.figure(figsize=(24, 12), 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)
#ax.view_init(elev=80, azim=-70)

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


# 標高に応じて塗りつぶす
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0, alpha=.7, cmap='YlOrRd')


# 標高50間隔で等高線を描く
elevation = range(-10000,10000,50)
cont = ax.contour(X, Y, Z, levels=elevation, cmap='binary', alpha=0.9)

# 水深・標高バーをつける
#cb = fig.colorbar(surf, shrink=0.4, aspect=20, extend='both', ticks=bounds, spacing='proportional')
#cb.set_label('elevation')

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

プログラム日記φ(．．)

おもにPython関係のプログラムメモ

pythonで複数の標高タイルを結合して富士山の3D立体地形図を描いてみる