緯度経度をXYZ方式のタイル座標に変換する

2024年4月より、国交省が不動産情報ライブラリという素敵なツールとそのAPIを公開しています。これにより、不動産の価格がわかる取引データや、用途地域や災害危険区域や駅ごとの乗降客数など、不動産に関するデータに簡単にアクセスできるようになりました。

しかし、用途地域など地理情報系のAPIは入力するパラメーターにXYZ方式のタイル座標を指定する必要があり、「この座標、どうやって求めればいいんだ？」と思った方も多いのではないでしょうか（私は思いました）

本記事では、緯度経度の座標をタイル座標に変換する方法を共有いたします。

XYZ方式とは

地図タイル

GoogleマップのようなWeb上の地図は、メルカトル図法の地図の一部を切り取ったタイルの集まりで表現し、必要なタイルだけを都度読み込んで高速化しています。

タイルはズームレベル（z）と横方向の番号（x）、縦方向の番号（y）によって座標が表現されます。例えば国土地理院の標準地図だと

https://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png

の形式でURLが構成されています。

例えば座標がz=6, x=57, y=23だと、 https://cyberjapandata.gsi.go.jp/xyz/std/6/57/23.png となるわけです。

なお、他の地図タイルも同様のURLの構造になっておりまして、OpenStreetMap（オープンデータの地図）だと https://tile.openstreetmap.org/6/57/23.png となります。

ちなみに国土地理院のタイルの確認には国土地理院のこちらのページが便利です。

maps.gsi.go.jp

座標の仕組み

この座標がどうやって決まるかというと、メルカトル図法の地図全体を1つの正方形に収めた状態をzoom level = 0として、縦横半分に切って4つに分割した状態をzoom level = 1とし、同様に4分割するたびにzoom levelが上がっていくものとみなして座標を割り振っています。

このあたりは国土地理院のこのページがわかりやすく参考になります。

maps.gsi.go.jp

座標の変換方法

さて、本題です。

緯度経度からXYZ座標への変換については、実はOpenStreetMapのWiki

wiki.openstreetmap.org

に計算式が紹介されています。

この数式をPythonで再現すると以下のようになります。

from math import radians, cos, tan, log, pi, floor

def latlon_to_tile(lat, lon, zoom):
    n = 2 ** zoom
    lat_rad = radians(lat)
    x = floor(n * ((lon + 180) / 360))
    y = floor(n * (1 - (log(tan(lat_rad) + 1 / cos(lat_rad)) / pi)) / 2)
    return x, y

例

六本木駅の座標でズームレベル15のタイルを見てみます。

z = 15
lat, lon = 35.6638271, 139.7316455 # 六本木駅
x, y = latlon_to_tile(lat, lon, zoom=z)
print(f"https://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png")

このコードの出力は https://cyberjapandata.gsi.go.jp/xyz/std/15/29102/12905.png になります。ちゃんと駅周辺が写っています。