matplotlib でグラフに地図タイルを載せる

**実行環境**
Windows10
Python 3.9.5
matpotlib 3.7.1
numpy 1.24.2
opencv 4.7.0
pyproj 3.5.0


## はじめに
[matplotlib](https://matplotlib.org/) の Axes に地図タイルを載せる方法について書きます。

matplotlib はあくまでもデータ可視化のためのライブラリで、地図を描くライブラリではないと人は言うかもしれません。

[folium](https://python-visualization.github.io/folium/) でいいではないかと。

それでも私は matplotlib に地図タイルを載せたいのです。

インタラクティブな地図を作りたいのではないのです。


## サンプルユースケース
[東京スカイツリー](https://www.tokyo-skytree.jp/) を中心とした東西南北 0.01 度の範囲に、地理院タイルの [標準地図](https://maps.gsi.go.jp/development/ichiran.html#std2) を載せた地図（グラフ）を作ります。

地図タイルのズームレベルは15、地図の座標系はブラウザでの表示と同じく Web メルカトルとします。

東京スカイツリーの座標は、地理院地図で「東京スカイツリー」を検索したときの以下の値を使います。

緯度 35.709529
経度 139.810713

なお地理院地図で東京スカイツリー付近をズームレベル 15 で表示すると、[こんな感じ](https://maps.gsi.go.jp/#15/35.709529/139.810713) です。


## 座標変換に関する関数

登場する座標系は、以下の３つです。

+ タイル座標（いわゆるxyz）
+ 緯度経度（WGS84）
+ Webメルカトル座標

地理院タイルの測地系は「世界測地系（JGD2011）」です [^1] が、WGS84 と実質同じとします。

[^1]: [地図投影法](https://maps.gsi.go.jp/development/siyou.html) 参照

これらの座標系間を変換する関数を定義します。

といっても、**deg2num()** と **num2deg()** は [wiki](https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames) からの引用で、緯度経度とタイル座標間の変換をする関数です。

**lonlat2m()** は、緯度経度をメートル単位の Web メルカトル座標に変換する関数です。

[pyproj](https://pyproj4.github.io/pyproj/stable/) というライブラリを使っています。

「4326」は WGS84 の EPSG コード、「3857」は Webメルカトル座標系の EPSG コードです。

```python
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](https://numpy.org/) の np.ndarray の配列として扱いますが、[OpenCV](https://opencv.org/) の仕様上、RGBの順ではなくBGRの順になります。

**concat_tile()** はタイル画像を取得し結合する関数で、[こちら](https://sorabatake.jp/7325/) [^2] を参考にしています。

[^2]: [get_combined_image()](https://sorabatake.jp/7325/) 参照

[urllib.request](https://docs.python.org/ja/3/library/urllib.request.html) を使っているのは、[requests](https://requests.readthedocs.io/en/latest/) だと ssl エラーがでてしまったためです。

**clip_tile()** は、結合した画像を切り抜く関数です。

内部で、先に定義した関数を使っています。

いろいろな座標値をごちゃごちゃと計算していますが、一言でいうと np.ndarray のスライス処理になります。

```python
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 の背景に設定し、画像出力します。

<details><summary>こちらはインポート部と先に載せた関数部です。</summary>

```python
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
```
</details>

ここから続きです。

```python
# 東京スカイツリーの座標
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')
```

出力画像がこちらです。

![map_skytree.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1409923/5b6cde39-a9b0-856b-ae17-a778cda863b8.png)

スカイツリー付近の地図を切り抜くことができました。

長方形になっているのは、投影の事情 [^3] になります。

[^3]: たとえば [メルカトル図法とは](https://club.informatix.co.jp/?p=1183) 参照

軸は、桁が大きいですが、Webメルカトル座標系でメートル単位です。


## GISで確認
GIS ソフトで同じことをして、それぞれの結果を比べてみます。

ポイントとポリゴンを WGS84 で座標値を与えて作成し、Web メルカトルのマップに載せました。

下の画像は、対象範囲より少し広い範囲を [QGIS](https://qgis.org/ja/site/) で画像出力したものです。

![レイアウト_skytree.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1409923/c100f209-fed6-c199-932a-f25185d1414f.png)

東京スカイツリー（シアンの点）が川の中ですが、matplotlib での出力と同じ地点です。

太い黒枠は切り抜いた範囲で、これも matplotlib での出力と一致しています。

これなら、matplotlib で地図タイルを精度よくグラフに載せられたと言ってよさそうです。


## おわりに
タイル画像をつなげて、切って、貼るという、原始的な方法ではありますが、matplotlib でグラフに地図タイルを載せる方法について書きました。

[geopandas](https://geopandas.org/en/stable/) というライブラリでは GIS ファイルを読み込んで matplotlib で出力できるので、今回の方法を利用すると、より地図らしいグラフができるかもしれません。