Python Basemapの背景地図をOpenStreetMapで描画してみました。
Python Basemapの背景地図をOpenStreetMapにしてみた。
Pythonには地図作図関係のライブラリーに
basemapとcartopyの代表的な２つのライブラリーが存在しますが、
- Welcome to the Matplotlib Basemap Toolkit documentation — Basemap Matplotlib Toolkit 1.1.0 documentation
- Introduction — cartopy 0.17.0 documentation
今回は、Pythonの地図作成ライブラリーのなかでも歴史が長いbasemapを使用し、basemapの背景地図をOpenStreetMapを使用し描画していました。
描画の投影法は、Googleマップをはじめ、Web上のオンライン地図サービスで標準的に使われているEPSGコード3857。測地系WGS84の球面メルカトル図法の地図にOpenStreetMapの地図タイル画像を流し込んでみることに。
以下、basemap上の緯度経度とOpenStreetMapの地図画像を読み込んで来て描画した地図画像の間にズレが生じてないかの確認に、basemapのdrawcoastlines関数で海岸線も描いてますが、地図画像の輪郭とdrawcoastlines関数で描いた海岸線の輪郭を比べてみると、実用上、特に問題ないようで
パブリックドメインの地図で画像を作成できると、はてなブログを始めTwitterたFacebook、Instagramなどなど、著作権フリーで大腕を振って地図画像を色々加工したりして貼り付けることが出来るのでとても便利です。
なお、今回も。OpenStreetMapの地図画像に国土地理院の日本全国の標高メッシュデータをもとに作成した斜面勾配図を合成して、山岳部などは若干立体的な陰影を付けるようにしてみましたが、国土地理院の標高メッシュデータは著作権的にいろいろ制限があるので、著作権フリーで地形の陰影付き地図を作りたいときは、ETOPO1の地形データなどを使うのがいいかもしれませんね。
また、以下では、地図に必要なタイル画像をOpenStreetMapから毎回ネットでダウンロードしてくる形にしてますが、常用で使うような場合は、読み込んできたタイル画像は、ファイルなどに保存しておくなどして、新規のタイル区画だけオンラインで読み込むようにして、以前利用しているタイルはファイルから読み込むなどの工夫をしてOpenStreetMapのサーバーに不必要に負荷をかけないようにするのがベストかもしれませんね。
ズームレベル5
X軸のタイル座標27〜29、Y軸のタイル座標11〜13で日本全体を描いてみたのが以下。
ズームレベル7
X軸のタイル座標112〜114、Y軸のタイル座標49〜51で本州中央部を描いてみたのが以下。
ズームレベル9
X軸のタイル座標453〜456、Y軸のタイル座標200〜203で関東地方を描いてみたのが以下。
ズームレベル11
X軸のタイル座標1765〜1767、Y軸のタイル座標841〜843で屋久島を描いてみたのが以下。
ズームレベル13
X軸のタイル座標7067〜7070、Y軸のタイル座標3336〜3339で桜島を描いてみたのが以下。
ズームレベル15
X軸のタイル座標29103〜29105、Y軸のタイル座標12901〜12903で皇居付近を描いてみたのが以下。
上の地図を作図したコードは以下になります
以下ソースで使用しているピクセル座標から緯度経度への変換は以下参照。
import requests from StringIO import StringIO import string import pandas as pd import numpy as np import matplotlib.pyplot as plt import matplotlib.colors as colors from matplotlib import cm from scipy import misc from mpl_toolkits.basemap import Basemap from osgeo import gdal from math import pi from math import tanh from math import sin from math import asin from math import exp from numpy import arctanh # ピクセル座標を緯度経度に変換する def fromPointToLatLng(pixelLat,pixelLon, z): L = 85.05112878 lon = 180 * ((pixelLon / 2.0**(z+7) ) - 1) lat = 180/pi * (asin(tanh(-pi/2**(z+7)*pixelLat + arctanh(sin(pi/180*L))))) return lat, lon # 地形データを読み込む 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 = cv2.vconcat([gis_v, Z]) #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 = cv2.hconcat([gis, gis_v]) #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=5 x1=27 x2=29 y1=11 y2=13 # ズームレベルに応じた経緯度線のステップ coordinateLineStep = { 5 : 10.0, 6 : 10.0, 7 : 5.0, 8 : 2.5, 9 : 1.0, 10 : 0.5, 11 : 0.25, 12 : 0.2, 13 : 0.1, 14 : 0.05, 15 : 0.02 } # 標高タイルのURLフォーマット urlFormat = 'http://cyberjapandata.gsi.go.jp/xyz/dem/{0}/{1}/{2}.txt' # 標高タイルを読み込み連結して１枚の標高タイルとして返す Z = load_gis(urlFormat, z, x1, x2, y1, y2) # 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軸方向の勾配をグレイスケール化する projectionIMG = cm.binary(Zgradient_norm) #透過情報はカット projectionIMG = projectionIMG[:,:,:3] # 地図画像と射影印影図を足して２で割りゃ画像が合成される imgColors = (imgColors*1.2 + projectionIMG) / 2 # 合成画像の輝度値の標準偏差を32,平均を80になんとなく変更 imgColors = (imgColors-np.mean(imgColors))/np.std(imgColors)*32+80 imgColors = (imgColors-imgColors.min())/(imgColors.max()-imgColors.min()) # 地図の描画範囲の緯度経度を求める xlim = Z.shape[1] ylim = Z.shape[0] minLat, minLon = fromPointToLatLng(y1 * 256 + ylim, x1 * 256, z) maxLat, maxLon = fromPointToLatLng(y1 * 256, x1 * 256 + xlim, z) # 地図を作成する fig = plt.figure(figsize=(12,9.8)) map = Basemap(epsg=3857, lat_ts=38, resolution='h', llcrnrlon=minLon, llcrnrlat=minLat, urcrnrlon=maxLon, urcrnrlat=maxLat) # 海岸線を描画 map.drawcoastlines( linewidth=1.0, color='k' ) # 経緯度線を描画 map.drawparallels(np.arange(10.0, 120.0, coordinateLineStep[z]), labels = [1,0,0,0], fontsize=18) map.drawmeridians(np.arange(50.0, 180.0, coordinateLineStep[z]), labels = [0,0,0,1], fontsize=18) # 背景地図画像を描画 mx0, my0 = map(minLon, minLat) mx1, my1 = map(maxLon, maxLat) extent = (mx0, mx1, my0, my1) plt.imshow(imgColors, vmin = 0, vmax = 255, extent = extent) #X,Y軸のグリッドを生成 #x = [fromPointToLatLng(y1* 256, x1* 256 + i, z)[1] for i in range(0, xlim)] #y = [fromPointToLatLng(y1* 256 + i, 0, z)[0] for i in range(ylim, 0, -1)] #x, y = map(x, y) #x = np.linspace(0, map.urcrnrx, Z.shape[1]) #y = np.linspace(0, map.urcrnry, Z.shape[0]) #X, Y = np.meshgrid(x, y) #標高250m間隔で等高線を描く #elevation = range(0,4000,250) #cont = map.contour(X, Y, Z, levels=elevation,linewidth=1, linestyles = 'solid', cmap='PiYG') #cont = map.contour(X, Y, Z, levels=elevation,linewidth=1, colors = '#c0d0d0') #cont = map.contour(X, Y, Z, levels=elevation,linewidth=1, colors = 'k', linestyles = 'solid') #cont.clabel(fmt='%1.1fm', fontsize=6) plt.savefig('1.jpg', dpi=72) plt.show()
