pythonで上高地、奥穂高付近の航空写真を立体写真地図にしてみる。
前回は、pythonのMatplotlibを使い、富士山の斜面勾配図を立体表示してみましたが、
plot_surfaceのfacecolorsパラメーターの色指定に、そのまま航空・衛星写真の色情報を渡してあげれば、地理院地図3Dと同じようにpythonでも立体写真地図を作れるんじゃない（？）
ということで、早速、国土地理院の上高地、奥穂高付近のシームレス航空写真を使い、試してみることに！
http://maps.gsi.go.jp/development/ichiran.html#seamlessphoto
結果は以下で。
上高地の北東上空から南東方向へ。奥上高地の横尾山荘、蝶ケ岳、七段ノ滝付近北東上空から飛行機に乗って上高地方向を見渡してる感じで、梓川沿いに上高地、明神池・河童橋方向へ、右手に涸沢カール、奥穂高岳、北穂高。左手に長塀山辺りを俯瞰した感じで立体写真地図を描いてみました。
今回のミソは、国土地理院の標高データをもとに作成した立体地図の表面を、load_imgColorsでGETしてきた航空写真の色情報をピクセル単位にRGB数値データとして取り出し、それをplot_surfaceのfacecolorsパラメーターに渡してあげて、立体地形図を航空写真と同じ色で表面を塗りつぶしてやろうというもの。
大雑把なズームレベルのタイル１枚の航空写真を立体化してみたのがこちら。
上の立体写真地図を作る元になった標高データと国土地理院さんの航空写真データは以下。
・標高データ
http://cyberjapandata.gsi.go.jp/xyz/dem/12/3614/1604.txt
・航空写真
http://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/12/3614/1604.jpg
また、こちらは、同じ場所の電子国土基本図を遊び半分で立体化してみたもの。
で、より細かいズームレベルの標高データと航空写真から同じ場所を立体化してみたものが以下。
それで、こちらは 地理院地図3Dから、今回描画した辺りの航空写真3D画像。倍率や方向、描画範囲が若干異なりますが、地理院地図3Dと同じような立体画像を作成するのはpythonでも難しくないようですね。
ということで！
地図に限らず、他の写真でも、写真データと写真の各ピクセルに対する3D空間位置データがセットであれば、写真を立体化することは可能でしょうね。
今回の作図に使ったソースは以下。
import requests from StringIO import StringIO import string import pandas as pd import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import Image # 地形データを読み込む 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: 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 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 # 地図画像データを読み込み各ピクセルをRGB値に変換した配列を返す def load_imgColors(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 url = 'http://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{0}/{1}/{2}.jpg'.format(z,x,y) #電子国土基本地図画像のURL #url = 'http://cyberjapandata.gsi.go.jp/xyz/std/{0}/{1}/{2}.png'.format(z,x,y) print url #地図画像データを読み込む response = requests.get(url) if response.status_code == 404: #地図画像データが無い区画は白塗りにする colors=np.ones((256, 256, 3), dtype=object) else: #画像READ img = Image.open(StringIO(response.content)) imgArray = np.asarray(img) #地図画像データの各ピクセルをRGBの値に変換する colors=np.empty((256,256), dtype=object) maxcol , maxrow = img.size for i in range(maxrow): for j in range(maxcol): # RGBの値は256で割り、0から1の範囲に変換しておく RGB = [imgArray[i,j][0]/256., imgArray[i,j][1]/256., imgArray[i,j][2]/256.] colors[i,j] = RGB if len(b)==0: b = colors else: b = np.append(b,colors,0) if len(a)==0: a = b else: a = np.append(a,b,1) return a #z=12 #x1=3614 #x2=3614 #y1=1604 #y2=1604 z=14 x1=14456 x2=14459 y1=6416 y2=6419 # 地形データを読み込む 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)) # 地図画像データを読み込み各ピクセルをRGB値に変換する imgColors = load_imgColors(z, x1, x2, y1, y2) # 立体地図を作成する fig = plt.figure(figsize=(15, 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=72, azim=30) # 横を1/0.5=2倍長く設定 ax.set_aspect(0.5, adjustable='box') # 立体地図の表面はシームレス航空写真の配色で塗り潰しながら、地形データから立体地図を描画する surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0, shade=False, antialiased=False, facecolors=imgColors, alpha=0.80) ax.set_xlim(0, xlim) ax.set_ylim(0, ylim) plt.savefig('1.jpg', dpi=72) #plt.show()
前回記事と一緒に、さらっと以下も目を通していただくと、上のソースでやってることの理解は早いと思います。
