任意の2地点の緯度経度情報から周辺の地図画像を出力するスクリプト(numpy、matplotlib、pandas)
やり方は色々あるのですが、今回は、「重ね合わせた地図を画像ファイルに出力する」ことを目的としました。
今回の手法
手法としては、2地点の緯度・経度の情報から、その間を埋めるweb上の地図タイル(png画像)をダウンロードして結合し、pythonでグラフを描くときによく用いられる matplotlib を使って(緯度経度を座標軸にして)プロットすることを考えました。
それで書いたのか以下の python スクリプトです
import math
import urllib
from io import BytesIO
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from PIL import Image
# 緯度経度からダウンロードする地図タイルを番号を算出
def deg2num(lat_deg, lon_deg, zoom):
lat_rad = math.radians(lat_deg)
n = 2.0 ** zoom
xtile = int((lon_deg + 180.0) / 360.0 * n)
ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
return (xtile, ytile)
# 地図タイルの左上の角の緯度経度を算出
def num2deg(xtile, ytile, zoom):
n = 2.0 ** zoom
loncal = xtile / n * 360.0 - 180.0
lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
latcal = math.degrees(lat_rad)
return (latcal, loncal)
def getImageCluster(lat_deg, lon_deg, delta_lat, delta_long, zoom):
# ダウンロードするタイルの決定
smurl = r"https://cyberjapandata.gsi.go.jp/xyz/std/{0}/{1}/{2}.png"
xmin, ymax =deg2num(lat_deg, lon_deg, zoom)
xmax, ymin =deg2num(lat_deg + delta_lat, lon_deg + delta_long, zoom)
# 使用するタイルの緯度経度の範囲
latmax_cal, lonmin_cal = num2deg(xmin, ymin, zoom)
latmin_cal0, lonmax_cal0 = num2deg(xmax+1, ymax+1, zoom)
dlat = latmax_cal - latmin_cal0
dlon = lonmax_cal0 - lonmin_cal
latmin_cal = latmax_cal - dlat
lonmax_cal = lonmin_cal + dlon
# 地図タイルのダウンロード、結合
Cluster = Image.new('RGB',((xmax-xmin+1)*256,(ymax-ymin+1)*256) )
for xtile in range(xmin, xmax+1):
for ytile in range(ymin, ymax+1):
try:
imgurl=smurl.format(zoom, xtile, ytile)
print("Opening: " + imgurl)
imgstr = BytesIO(urllib.request.urlopen(imgurl).read())
tile = Image.open(imgstr)
Cluster.paste(tile, box=((xtile-xmin)*256 , (ytile-ymin)*256))
except:
print("Couldn't download image")
tile = None
return (Cluster, latmin_cal, latmax_cal, lonmin_cal, lonmax_cal)
if __name__ == '__main__':
lat1, lon1 = 32.814595,130.871619 # 地点1
lat2, lon2 = 32.793385,130.830152 # 地点2
zoom = 14 # ズームレベル
latmin_deg = min(lat1, lat2)
lonmin_deg = min(lon1, lon2)
delta_lat = abs(lat1 - lat2)
delta_long = abs(lon1 - lon2)
a, latmin_cal, latmax_cal, lonmin_cal, lonmax_cal = getImageCluster(latmin_deg, lonmin_deg, delta_lat, delta_long, zoom)
fig = plt.figure()
fig.patch.set_facecolor('white')
plt.imshow(np.asarray(a),extent=[lonmin_cal, lonmax_cal, latmin_cal, latmax_cal])
plt.savefig('map.png')
ちなみに、このスクリプトは、地点1、地点2 の緯度経度とズームレベルを指定すれば、地図が出力されます。
この地図の上に、プロットしたい地点の緯度経度をcsvファイル(ここでは、location.csv)にまとめておいて、pandas で読み込ませ、matplotlib の scatter でプロットします。
そのために、main の中に、以下の行を追加します。
csvdata= pd.read_csv("location.csv")
x1 = []
y1 = []
y1 = csvdata['緯度']
x1 = csvdata['経度']
plt.scatter(x1, y1, c='red',marker='o')
今回は、国土地理院のものを使用しましたが、タイル地図はopenstreetmapを含め、利用できるものが多くあります。
今後の課題
- 軸の表記等を整理して、もう少しきれいな図にしたい。
- できれば、Julia 言語でやりたい。