実行環境
Windows10
Python 3.9.5
matpotlib 3.7.1
numpy 1.24.2
opencv 4.7.0
pyproj 3.5.0
はじめに
matplotlib の Axes に地図タイルを載せる方法について書きます。
matplotlib はあくまでもデータ可視化のためのライブラリで、地図を描くライブラリではないと人は言うかもしれません。
folium でいいではないかと。
それでも私は matplotlib に地図タイルを載せたいのです。
インタラクティブな地図を作りたいのではないのです。
サンプルユースケース
東京スカイツリー を中心とした東西南北 0.01 度の範囲に、地理院タイルの 標準地図 を載せた地図（グラフ）を作ります。
地図タイルのズームレベルは15、地図の座標系はブラウザでの表示と同じく Web メルカトルとします。
東京スカイツリーの座標は、地理院地図で「東京スカイツリー」を検索したときの以下の値を使います。
緯度 35.709529
経度 139.810713
なお地理院地図で東京スカイツリー付近をズームレベル 15 で表示すると、こんな感じ です。
座標変換に関する関数
登場する座標系は、以下の３つです。
- タイル座標（いわゆるxyz）
- 緯度経度（WGS84）
- Webメルカトル座標
地理院タイルの測地系は「世界測地系（JGD2011）」です 1 が、WGS84 と実質同じとします。
これらの座標系間を変換する関数を定義します。
といっても、deg2num() と num2deg() は wiki からの引用で、緯度経度とタイル座標間の変換をする関数です。
lonlat2m() は、緯度経度をメートル単位の Web メルカトル座標に変換する関数です。
pyproj というライブラリを使っています。
「4326」は WGS84 の EPSG コード、「3857」は Webメルカトル座標系の EPSG コードです。
import math
import pyproj
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.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
return (xtile, ytile)
def num2deg(xtile, ytile, zoom):
n = 2.0 ** zoom
lon_deg = xtile / n * 360.0 - 180.0
lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
lat_deg = math.degrees(lat_rad)
return (lat_deg, lon_deg)
def lonlat2m(lon: float, lat: float) -> tuple[float, float]:
""" 座標変換: WGS84[度] -> Webメルカトル[m]
"""
transformer = pyproj.Transformer.from_crs(4326, 3857, always_xy=True)
x, y = transformer.transform(lon, lat)
return x, y
タイル画像の処理に関する関数
複数枚のタイル画像を取得・結合し、境界範囲に合わせて切り抜きます。
画像データは numpy の np.ndarray の配列として扱いますが、OpenCV の仕様上、RGBの順ではなくBGRの順になります。
concat_tile() はタイル画像を取得し結合する関数で、こちら 2 を参考にしています。
urllib.request を使っているのは、requests だと ssl エラーがでてしまったためです。
clip_tile() は、結合した画像を切り抜く関数です。
内部で、先に定義した関数を使っています。
いろいろな座標値をごちゃごちゃと計算していますが、一言でいうと np.ndarray のスライス処理になります。
import time
import urllib.request
import cv2
import numpy as np
def concat_tile(
urlformat: str,
zoom: int,
left_tile: int,
right_tile: int,
top_tile: int,
bottom_tile: int
) -> np.ndarray:
""" タイル画像を取得し結合
"""
# requestsライブラリだとsslエラーになってしまったので、
# urllib.requestを使っている
rows = []
for ytile in range(top_tile, bottom_tile + 1):
row = []
for xtile in range(left_tile, right_tile + 1):
url = urlformat.format(z=zoom, x=xtile, y=ytile)
# タイル画像1枚分 (256, 256, 3) の多次元行列を取得
img = cv2.imdecode(
buf=np.array(
bytearray(urllib.request.urlopen(url).read()),
dtype=np.uint8
),
flags=-1
)
time.sleep(1)
row.append(img)
# 横方向に結合
rows.append(np.hstack(row))
# 縦方向に結合
arr = np.vstack(rows)
return arr
def clip_tile(
urlformat: str,
zoom: int,
lat_n: float,
lat_s: float,
lon_w: float,
lon_e: float
) -> tuple[np.ndarray, float, float, float, float]:
""" 地図タイルを指定の範囲で切り抜く
Parameters
----
urlformat : str
地図タイルのURL
zoom : int
ズームレベル
lat_n : float
北端緯度 [度]
lat_s : float
南端緯度 [度]
lon_w : float
西端経度 [度]
lon_e : float
東端経度 [度]
Returns
----
arr : np.ndarray
BGRの多次元配列 (RGBではない)
north_m : float
北端y座標 [m]
south_m : float
南端y座標 [m]
west_m : float
西端x座標 [m]
east_m : float
東端x座標 [m]
"""
# --境界範囲の情報--
# 各方位のxy[m]座標を取得
west_m, north_m = lonlat2m(lon_w, lat_n)
east_m, south_m = lonlat2m(lon_e, lat_s)
# --境界範囲にかかるタイルの情報--
# タイル座標を取得
left_tile, top_tile = deg2num(lat_deg=lat_n, lon_deg=lon_w, zoom=zoom)
right_tile, bottom_tile = deg2num(lat_deg=lat_s, lon_deg=lon_e, zoom=zoom)
# タイル境界の緯度経度を取得
top_deg, left_deg = num2deg(left_tile, top_tile, zoom)
bottom_deg, right_deg = num2deg(right_tile + 1, bottom_tile + 1, zoom)
# タイル境界のxy[m]座標を取得
left_m, top_m = lonlat2m(left_deg, top_deg)
right_m, bottom_m = lonlat2m(right_deg, bottom_deg)
# 幅と高さを[m]で取得
width_m = right_m - left_m
height_m = top_m - bottom_m
# 画素数を取得
px_h = (right_tile - left_tile + 1) * 256
px_v = (bottom_tile - top_tile + 1) * 256
# 境界範囲にかかるタイルを結合
arr_concat = concat_tile(
urlformat=urlformat,
zoom=zoom,
left_tile=left_tile,
right_tile=right_tile,
top_tile=top_tile,
bottom_tile=bottom_tile
)
# 切り抜き境界[pixel] 誤差許容
px_top = round((top_m - north_m) / height_m * px_v)
px_bottom = round((top_m - south_m) / height_m * px_v)
px_left = round((west_m - left_m) / width_m * px_h)
px_right = round((east_m - left_m) / width_m * px_h)
# 切り抜き
arr_trim = arr_concat[px_top: px_bottom, px_left: px_right]
return arr_trim, north_m, south_m, west_m, east_m
matplotlib で画像出力
切り抜いた画像を matplotlib の axes の背景に設定し、画像出力します。
こちらはインポート部と先に載せた関数部です。
import math
import time
import urllib.request
import cv2
import matplotlib.pyplot as plt
import numpy as np
import pyproj
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.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
return (xtile, ytile)
def num2deg(xtile, ytile, zoom):
n = 2.0 ** zoom
lon_deg = xtile / n * 360.0 - 180.0
lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
lat_deg = math.degrees(lat_rad)
return (lat_deg, lon_deg)
def lonlat2m(lon: float, lat: float) -> tuple[float, float]:
""" 座標変換: WGS84[度] -> Webメルカトル[m]
"""
transformer = pyproj.Transformer.from_crs(4326, 3857, always_xy=True)
x, y = transformer.transform(lon, lat)
return x, y
def concat_tile(
urlformat: str,
zoom: int,
left_tile: int,
right_tile: int,
top_tile: int,
bottom_tile: int
) -> np.ndarray:
""" タイル画像を取得し結合
"""
# requestsライブラリだとsslエラーになってしまったので、
# urllib.requestを使っている
rows = []
for ytile in range(top_tile, bottom_tile + 1):
row = []
for xtile in range(left_tile, right_tile + 1):
url = urlformat.format(z=zoom, x=xtile, y=ytile)
# タイル画像1枚分 (256, 256, 3) の多次元行列を取得
img = cv2.imdecode(
buf=np.array(
bytearray(urllib.request.urlopen(url).read()),
dtype=np.uint8
),
flags=-1
)
time.sleep(1)
row.append(img)
# 横方向に結合
rows.append(np.hstack(row))
# 縦方向に結合
arr = np.vstack(rows)
return arr
def clip_tile(
urlformat: str,
zoom: int,
lat_n: float,
lat_s: float,
lon_w: float,
lon_e: float
) -> tuple[np.ndarray, float, float, float, float]:
""" 地図タイルを指定の範囲で切り抜く
Parameters
----
urlformat : str
地図タイルのURL
zoom : int
ズームレベル
lat_n : float
北端緯度 [度]
lat_s : float
南端緯度 [度]
lon_w : float
西端経度 [度]
lon_e : float
東端経度 [度]
Returns
----
arr : np.ndarray
BGRの多次元配列 (RGBではない)
north_m : float
北端y座標 [m]
south_m : float
南端y座標 [m]
west_m : float
西端x座標 [m]
east_m : float
東端x座標 [m]
"""
# --境界範囲の情報--
# 各方位のxy[m]座標を取得
west_m, north_m = lonlat2m(lon_w, lat_n)
east_m, south_m = lonlat2m(lon_e, lat_s)
# --境界範囲にかかるタイルの情報--
# タイル座標を取得
left_tile, top_tile = deg2num(lat_deg=lat_n, lon_deg=lon_w, zoom=zoom)
right_tile, bottom_tile = deg2num(lat_deg=lat_s, lon_deg=lon_e, zoom=zoom)
# タイル境界の緯度経度を取得
top_deg, left_deg = num2deg(left_tile, top_tile, zoom)
bottom_deg, right_deg = num2deg(right_tile + 1, bottom_tile + 1, zoom)
# タイル境界のxy[m]座標を取得
left_m, top_m = lonlat2m(left_deg, top_deg)
right_m, bottom_m = lonlat2m(right_deg, bottom_deg)
# 幅と高さを[m]で取得
width_m = right_m - left_m
height_m = top_m - bottom_m
# 画素数を取得
px_h = (right_tile - left_tile + 1) * 256
px_v = (bottom_tile - top_tile + 1) * 256
# 境界範囲にかかるタイルを結合
arr_concat = concat_tile(
urlformat=urlformat,
zoom=zoom,
left_tile=left_tile,
right_tile=right_tile,
top_tile=top_tile,
bottom_tile=bottom_tile
)
# 切り抜き境界[pixel] 誤差許容
px_top = round((top_m - north_m) / height_m * px_v)
px_bottom = round((top_m - south_m) / height_m * px_v)
px_left = round((west_m - left_m) / width_m * px_h)
px_right = round((east_m - left_m) / width_m * px_h)
# 切り抜き
arr_trim = arr_concat[px_top: px_bottom, px_left: px_right]
return arr_trim, north_m, south_m, west_m, east_m
ここから続きです。
# 東京スカイツリーの座標
lat_skytree = 35.709529
lon_skytree = 139.810713
clip_deg = 0.01 # 切り抜く範囲
file_png = 'map_skytree.png' # 出力ファイル名
zoom = 15 # ズームレベル
lat_n = lat_skytree + clip_deg / 2 # 北端緯度
lat_s = lat_skytree - clip_deg / 2 # 南端緯度
lon_w = lon_skytree - clip_deg / 2 # 西端経度
lon_e = lon_skytree + clip_deg / 2 # 東端経度
# 地理院タイル 標準地図の URL
urlformat = 'https://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png'
# 境界範囲で地図タイルを切り抜き
arr, north_m, south_m, west_m, east_m = clip_tile(
urlformat=urlformat,
zoom=zoom,
lat_n=lat_n,
lat_s=lat_s,
lon_w=lon_w,
lon_e=lon_e
)
# ここから matplotlib 処理
fig = plt.figure()
ax = fig.add_subplot()
# 色の順番をBGRからRGBに変更
# 範囲をメルカトル座標系の境界座標に指定
ax.imshow(
cv2.cvtColor(arr, cv2.COLOR_BGR2RGB),
extent=[west_m, east_m, south_m, north_m]
)
# 東京スカイツリーをプロット
x_skytree, y_skytree = lonlat2m(lon_skytree, lat_skytree)
ax.scatter(x_skytree, y_skytree, color='cyan', edgecolors='black')
# 調整
ax.grid(True, color='black')
ax.tick_params(axis='x', labelrotation=90)
# 出力
fig.savefig(file_png, bbox_inches='tight')
出力画像がこちらです。
スカイツリー付近の地図を切り抜くことができました。
長方形になっているのは、投影の事情 3 になります。
軸は、桁が大きいですが、Webメルカトル座標系でメートル単位です。
GISで確認
GIS ソフトで同じことをして、それぞれの結果を比べてみます。
ポイントとポリゴンを WGS84 で座標値を与えて作成し、Web メルカトルのマップに載せました。
下の画像は、対象範囲より少し広い範囲を QGIS で画像出力したものです。
東京スカイツリー（シアンの点）が川の中ですが、matplotlib での出力と同じ地点です。
太い黒枠は切り抜いた範囲で、これも matplotlib での出力と一致しています。
これなら、matplotlib で地図タイルを精度よくグラフに載せられたと言ってよさそうです。
おわりに
タイル画像をつなげて、切って、貼るという、原始的な方法ではありますが、matplotlib でグラフに地図タイルを載せる方法について書きました。
geopandas というライブラリでは GIS ファイルを読み込んで matplotlib で出力できるので、今回の方法を利用すると、より地図らしいグラフができるかもしれません。