自ら環境を構築して執筆するブログ

やり方は色々あるのですが、今回は、「重ね合わせた地図を画像ファイルに出力する」ことを目的としました。

今回の手法

手法としては、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 言語でやりたい。