OpenStreetMapの地図画像で立体地図をつくってみました。

Googleマップにしろ、国土地理院の地図にしろ、その画像をネット上で自由に加工して配信したいとなると、その利用には何かと著作権が絡むので、みんなで使えるフリーマップOpenStreetMapを加工し、立体地図をPythonで作成してみました。

OpenStreetMap Japan | 自由な地図をみんなの手に/The Free Wiki World Map

もともとのOpenStreetMapの画像に

『pythonで富士山の斜面勾配図ぽいものを描いてみる』の記事で作成した方法で斜面勾配図を作成して合成すると

こんな感じになります。

 

ソースはこんな感じ

なお、以下のソースで画像のRGB値などを0..1に正規化したりしてるのは、地形を描くためのplot_surface関数のパラメータfacecolorsが受け付ける値が、RGB値を0..1に正規化した値のためなど、matplotlibのデータの扱い方に合わせるためで、それ以上の意味はありません。

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

from scipy import misc
import cv2


標高タイルを読み込み連結して１個の標高タイルとして返す
def load_gis(urlFormat, z, x1, x2, y1, y2):

    for x in range(x1, x2+1):
        for y in range(y1, y2+1):

            #地形データを読み込む
            url = urlFormat.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 y == y1:
                 gis_v = Z
            else:
                 #gis_v = np.append(gis_v,Z,0)
                 gis_v = np.concatenate((gis_v, Z), axis = 0) #縦
                
        # 標高タイルを横に連結
        if x == x1:
            gis = gis_v
        else:
            #gis = np.append(gis,gis_v,1)
            gis = np.concatenate((gis,gis_v), axis = 1) #横
            
    return gis
    


# 地図画像データを読み込み各ピクセルをRGB値に変換した配列を返す
def load_imgColors(urlFormat, z, x1, x2, y1, y2):

    for x in range(x1, x2+1):
        for y in range(y1, y2+1):

            #地図画像データを読み込む
            url = urlFormat.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 = misc.imread(StringIO(response.content))
                 
            #　画像タイルを縦に連結
            if y == y1:
                 im_v = img
            else:
                 #im_v = cv2.vconcat([im_v, img])
                 #im_v = np.append(im_v,img,0)
                 im_v = np.concatenate((im_v, img), axis = 0) #縦
                
        # 画像タイルを横に連結
        if x == x1:
            im = im_v
        else:
            #im = cv2.hconcat([im, im_v])
            #im = np.append(im,im_v,1)
            im = np.concatenate((im,im_v), axis = 1) #横
            
    return im



z=10
x1=906
x2=906
y1=404
y2=404


# 標高タイルのURLフォーマット
urlFormat = 'http://cyberjapandata.gsi.go.jp/xyz/dem/{0}/{1}/{2}.txt'

# 標高タイルを読み込み連結して１個の標高タイルとして返す
Z = load_gis(urlFormat, 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))


# OpenStreetMapのタイル画像のURLフォーマット
urlFormat = 'http://a.tile.openstreetmap.org/{0}/{1}/{2}.png'

# OpenStreetMapのタイル画像を読み込み連結して１枚の画像として返す
imgColors = load_imgColors(urlFormat, z, x1, x2, y1, y2) 

# 地図画像RGBデータは256で割って0..1に正規化しておく
imgColors = imgColors/256.


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

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

# Y軸方向の勾配をRGB値が0..1に正規されたグレイスケール画像に変換する
projectionIMG = cm.binary(Zgradient_norm)

# グレースケールの色彩を2.5で割って弱くしてみる
projectionIMG /= 2.5

#透過情報はカット
projectionIMG = projectionIMG[:,:,:3]

# 地図画像と射影印影図を足して２で割りゃ画像が合成される
imgColors = (imgColors + projectionIMG) / 2


# 合成画像の輝度値の標準偏差を32,平均を80になんとなく変更して0..1に正規化
imgColors = (imgColors-np.mean(imgColors))/np.std(imgColors)*32+80
imgColors = (imgColors-imgColors.min())/(imgColors.max()-imgColors.min())



# 立体地図を作成する
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=-60)


# 横を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=1.0)


ax.set_xlim(0, xlim)
ax.set_ylim(0, ylim)

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

また、 上のソースにこんなコードを追加して

#起伏を示した地図
urlFormat2 = 'https://cyberjapandata.gsi.go.jp/xyz/relief/{0}/{1}/{2}.png'

#起伏を示した地図画像をロード
imgColors2 = load_imgColors(urlFormat2, z, x1, x2, y1, y2) 
imgColors2 = imgColors2/256.

国土理智院の起伏を示した地図を読み込み

色彩を弱めてからOpenStreetMapの画像に合成してあげると

# 2で割って起伏を示した地図の色彩を弱くする
imgColors2 /= 2

imgColors = (imgColors + imgColors2) / 2

こんな感じになります。

なを、最初に掲載した画像はズームレベル11、X座標1812〜1814、Y座標807〜809 のタイル座標でOpenStreetMapと国土地理院の標高データから合成して作成してみた、より広範囲の立体地形図になります。

1860x1200ピクセルの大画像となってますので、画像を新規ウィンドウで開いて拡大してどんな感じか確認してみてください。

以上、OpenStreetMapの地図画像から立体地図を作成する例でしたが、OpenStreetMapの地図画像に斜面勾配図を合成して、上のソースで立体地図を描いているところを単純に

fig = plt.figure(dpi=72, facecolor='w', edgecolor='k')
plt.imshow(imgColors)
plt.savefig('1.jpg', dpi=72)
plt.show()

と変えてあげて2D平面地図画像にするだけでもOpenStreetMapの地図画像が以下の感じでより立体的になります。（実際には以下の画像は、ガンマ補正とか、コントラスト補正とかのロジックも入れ込んで若干カラー補正もして作図したものではありますが。。）

ということで、今回はOpenStreetMapの地図画像と国土地理院の標高データを使って立体地図を作成しましたが、OpenStreetMapの地図画像とこちらも著作権フリーの標高データ、スペースシャトルのSRTMを使って立体地図を作成すれは、画像の片隅にOpenStreetMapのコピーライト表示さえしておけば、誰に憚ることなく気兼ねなく自由にネットで地図画像を配信できますね^^/

memomemokun.hateblo.jp