気になるメモメモφ(..)

日常気になる事を取り留めもなくメモメモしてます。

pythonで上高地、奥穂高付近の航空写真を立体写真地図にしてみる。

前回は、pythonのMatplotlibを使い、富士山の斜面勾配図を立体表示してみましたが、

 

memomemokun.hateblo.jp

 

 

plot_surfaceのfacecolorsパラメーターの色指定に、そのまま航空・衛星写真の色情報を渡してあげれば、地理院地図3Dと同じようにpythonでも立体写真地図を作れるんじゃない(?)

 

ということで、早速、国土地理院上高地、奥穂高付近のシームレス航空写真を使い、試してみることに!

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

 

結果は以下で。

 

上高地の北東上空から南東方向へ。奥上高地の横尾山荘、蝶ケ岳、七段ノ滝付近北東上空から飛行機に乗って上高地方向を見渡してる感じで、梓川沿いに上高地、明神池・河童橋方向へ、右手に涸沢カール、奥穂高岳、北穂高。左手に長塀山辺りを俯瞰した感じで立体写真地図を描いてみました。

 

今回のミソは、国土地理院の標高データをもとに作成した立体地図の表面を、load_imgColorsでGETしてきた航空写真の色情報をピクセル単位にRGB数値データとして取り出し、それをplot_surfaceのfacecolorsパラメーターに渡してあげて、立体地形図を航空写真と同じ色で表面を塗りつぶしてやろうというもの。

 

大雑把なズームレベルのタイル1枚の航空写真を立体化してみたのがこちら。

 

 

f:id:memomemokun:20181208181909j:plain

上の立体写真地図を作る元になった標高データと国土地理院さんの航空写真データは以下。

・標高データ
   http://cyberjapandata.gsi.go.jp/xyz/dem/12/3614/1604.txt
・航空写真
   http://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/12/3614/1604.jpg

f:id:memomemokun:20181208210854j:plain

 

また、こちらは、同じ場所の電子国土基本図を遊び半分で立体化してみたもの。

奥上高地から上高地方向への立体校区写真図(2)

 

で、より細かいズームレベルの標高データと航空写真から同じ場所を立体化してみたものが以下。

奥上高地から上高地方向への立体校区写真図(3)

 

それで、こちらは 地理院地図3Dから、今回描画した辺りの航空写真3D画像。倍率や方向、描画範囲が若干異なりますが、地理院地図3Dと同じような立体画像を作成するのはpythonでも難しくないようですね。

f:id:memomemokun:20181208191121j:plain

 

 

ということで!

地図に限らず、他の写真でも、写真データと写真の各ピクセルに対する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()

 

 

前回記事と一緒に、さらっと以下も目を通していただくと、上のソースでやってることの理解は早いと思います。

 

memomemokun.hateblo.jp

 

memomemokun.hateblo.jp