緯度経度と地図タイル XYZ 座標の相互変換 Python スクリプトを書きました。

この変換の記事は N 番煎じですが、「なぜその計算かの説明 + Python 実装 + 実装の正しさの検証」が揃った記事は少ないかもしれないので投稿します。

参考文献

地理院地図｜地理院タイルについて
- このページでタイルの XYZ 座標の規則がわかりますが、メルカトル投影に疎いため最初「北緯及び南緯約85.0511度」の数値の意味がわかりませんでした。[2] を合わせて読んだ後、全世界のメルカトル投影図を正方形になるよう上下をトリミングしたときの上下端だと理解しました。
メルカトルとは何かを計算式から見る
- このページがメルカトル図法の目指すところについてわかりやすかったです。
- 積分計算もしてくださっており、なので私は数式を書くのをさぼります。
メルカトル図法 - Wikipedia
- ここの逆変換式を参照しました。

相互変換 Python スクリプト

個々の計算の説明を直後の print() に書いています。不要なら print() を削除ください。
2025/11/30 スクリプト内の説明は短いので、「計算の説明」を記事の最後に追加しました。

script.py

import math


def lonlat_to_xyz(lon, lat, z):
    print(f'横幅360, 縦幅180, 原点中央の長方形でここは ({lon:10.6f}, {lat:10.6f})')
    x = math.radians(lon)  # ラジアンへ
    y = math.radians(lat)  # ラジアンへ
    print(f'横幅2pi, 縦幅 pi, 原点中央の長方形でここは ({x:10.6f}, {y:10.6f})')
    y = math.log(math.tan(0.25 * math.pi + 0.5 * y))  # [2] の積分結果
    print(f'赤道から離れるほど横にのびる分縦ものばすと ({x:10.6f}, {y:10.6f})')
    x /= 2 * math.pi
    y /= 2 * math.pi
    print(f'2piで割って横幅1に正規化した長方形でここは ({x:10.6f}, {y:10.6f})')
    x = x - (-0.5)
    y =  - (y - 0.5)
    print(f'(-0.5, 0.5)を原点にしてy軸の向きを逆にして ({x:10.6f}, {y:10.6f})')
    n = 2 ** z
    x *= n
    y *= n
    print(f'ズームレベル{z}のグリッド({n:8d})を切ると ({x:10.3f}, {y:10.3f})')
    x = math.floor(x)
    y = math.floor(y)
    path = f'/{z}/{x}/{y}.png'
    print(f'したがってこの地点を含むタイルのXYZパスは  {path}\n')
    return x, y, path


def xyz_to_lonlat(x, y, z):
    n = 2 ** z
    print(f'ズームレベル{z}のグリッド({n:8d})でここは ({x:10.3f}, {y:10.3f})')
    x = (x + 0.5) / n
    y = (y + 0.5) / n
    print(f'ズームレベル 0のグリッド換算でタイル中心は ({x:10.6f}, {y:10.6f})')
    x = x - 0.5
    y =  - (y - 0.5)
    print(f'(0.5, 0.5)を原点にしてy軸の向きを逆にして  ({x:10.6f}, {y:10.6f})')
    x *= 2 * math.pi
    y *= 2 * math.pi
    y = math.asin(math.tanh(y))  # [3]
    print(f'横も縦も 2πにして縦はのばしていた分縮めて  ({x:10.6f}, {y:10.6f})')
    lon = math.degrees(x)
    lat = math.degrees(y)
    print(f'弧度法を度数法に直してタイル中心の経緯度は ({lon:10.6f}, {lat:10.6f})\n')
    return lon, lat


if __name__ == '__main__':
    z = 14
    lat, lon = 43.044706, 144.194578  # 釧路空港
    x, y, _ = lonlat_to_xyz(lon, lat, z)
    xyz_to_lonlat(x, y, z)

    lat, lon = 85.0511, 0.0  # XYZ 方式タイルがぎりぎりカバーできる地点
    x, y, _ = lonlat_to_xyz(lon, lat, z)
    xyz_to_lonlat(x, y, z)

    lat, lon = 85.0512, 0.0  # XYZ 方式タイルがぎりぎりカバーできない地点
    x, y, _ = lonlat_to_xyz(lon, lat, z)  # タイルの y が -1 になり対応できない
    xyz_to_lonlat(x, y, z)  # 変換・逆変換自体はできる

検証

どこでもよいですが釧路空港の座標で検証してみます。釧路空港のウィキペディアからリンクされている GeoHack から緯経度をコピーして変換と逆変換を実行すると以下になります。

lat, lon = 43.044706, 144.194578  # 釧路空港

横幅360, 縦幅180, 原点中央の長方形でここは (144.194578,  43.044706)
横幅2pi, 縦幅 pi, 原点中央の長方形でここは (  2.516670,   0.751272)
赤道から離れるほど横にのびる分縦ものばすと (  2.516670,   0.833908)
2piで割って横幅1に正規化した長方形でここは (  0.400540,   0.132721)
(-0.5, 0.5)を原点にしてy軸の向きを逆にして (  0.900540,   0.367279)
ズームレベル14のグリッド(   16384)を切ると ( 14754.455,   6017.506)
したがってこの地点を含むタイルのXYZパスは  /14/14754/6017.png

ズームレベル14のグリッド(   16384)でここは ( 14754.000,   6017.000)
ズームレベル 0のグリッド換算でタイル中心は (  0.900543,   0.367279)
(0.5, 0.5)を原点にしてy軸の向きを逆にして  (  0.400543,   0.132721)
横も縦も 2πにして縦はのばしていた分縮めて  (  2.516687,   0.751274)
弧度法を度数法に直してタイル中心の経緯度は (144.195557,  43.044805)

特定されたタイル /14/14754/6017.png は OSMFJ タイルサーバだと https://tile.openstreetmap.jp/14/14754/6017.png になります。
この URL をブラウザからみるとちゃんと釧路空港があることがわかります。

おまけとして私が手元で立てたタイルサーバのタイル /14/14754/6017.png も以下に貼っておきます (出典： OpenStreetMap データ)。

逆に、タイル /14/14754/6017.png の中心経緯度が (144.195557, 43.044805) で正しいかは、以下で 43.044805, 144.195557 を検索して上図と照らし合わせてください (スクリプトは経度, 緯度を返しますが検索時は緯度, 経度にしてください; 適宜直してください)。
https://www.openstreetmap.org/search?query=43.044805%2C+144.195557&zoom=14&minlon=144.15564537048343&minlat=43.02567075676123&maxlon=144.23546791076663&maxlat=43.063934102821534#map=14/43.04481/144.19556

備考

XYZ 方式タイルがぎりぎりカバーできる、北緯 85.0511 度の地点を変換・逆変換すると以下になります。タイルの Y 座標が 0 ですがぎりぎり対応しています。OSMFJ タイルサーバでこのタイルをみられます。まあみてもただの海ですが、タイルはあります。

横幅360, 縦幅180, 原点中央の長方形でここは (  0.000000,  85.051100)
横幅2pi, 縦幅 pi, 原点中央の長方形でここは (  0.000000,   1.484422)
赤道から離れるほど横にのびる分縦ものばすと (  0.000000,   3.141587)
2piで割って横幅1に正規化した長方形でここは (  0.000000,   0.499999)
(-0.5, 0.5)を原点にしてy軸の向きを逆にして (  0.500000,   0.000001)
ズームレベル14のグリッド(   16384)を切ると (  8192.000,      0.015)
したがってこの地点を含むタイルのXYZパスは  /14/8192/0.png

ズームレベル14のグリッド(   16384)でここは (  8192.000,      0.000)
ズームレベル 0のグリッド換算でタイル中心は (  0.500031,   0.000031)
(0.5, 0.5)を原点にしてy軸の向きを逆にして  (  0.000031,   0.499969)
横も縦も 2πにして縦はのばしていた分縮めて  (  0.000192,   1.484406)
弧度法を度数法に直してタイル中心の経緯度は (  0.010986,  85.050181)

XYZ 方式タイルがぎりぎりカバーできない、北緯 85.0512 度の地点を変換・逆変換すると以下になります。タイルの Y 座標が -1 になってしまい、XYZ 方式タイルのサポート範囲を超えてしまったことがわかります。ただそれはあくまで XYZ 方式タイルの仕様の話であり、上記スクリプトの変換・逆変換自体は北緯/南緯 90 度に近づくまでは数学的に成立します。

横幅360, 縦幅180, 原点中央の長方形でここは (  0.000000,  85.051200)
横幅2pi, 縦幅 pi, 原点中央の長方形でここは (  0.000000,   1.484423)
赤道から離れるほど横にのびる分縦ものばすと (  0.000000,   3.141607)
2piで割って横幅1に正規化した長方形でここは (  0.000000,   0.500002)
(-0.5, 0.5)を原点にしてy軸の向きを逆にして (  0.500000,  -0.000002)
ズームレベル14のグリッド(   16384)を切ると (  8192.000,     -0.038)
したがってこの地点を含むタイルのXYZパスは  /14/8192/-1.png

ズームレベル14のグリッド(   16384)でここは (  8192.000,     -1.000)
ズームレベル 0のグリッド換算でタイル中心は (  0.500031,  -0.000031)
(0.5, 0.5)を原点にしてy軸の向きを逆にして  (  0.000031,   0.500031)
横も縦も 2πにして縦はのばしていた分縮めて  (  0.000192,   1.484439)
弧度法を度数法に直してタイル中心の経緯度は (  0.010986,  85.052076)

計算の説明

数字のきりのよさのため、北緯35°, 東経135° の地点 (兵庫県西脇市) で上記スクリプトを実行すると以下になります。

横幅360, 縦幅180, 原点中央の長方形でここは (135.000000,  35.000000)  # 【ア】
横幅2pi, 縦幅 pi, 原点中央の長方形でここは (  2.356194,   0.610865)  # 【イ】
赤道から離れるほど横にのびる分縦ものばすと (  2.356194,   0.652837)  # 【ウ】
2piで割って横幅1に正規化した長方形でここは (  0.375000,   0.103902)  # 【エ】
(-0.5, 0.5)を原点にしてy軸の向きを逆にして (  0.875000,   0.396098)  # 【オ】
ズームレベル14のグリッド(   16384)を切ると ( 14336.000,   6489.667)  # 【カ】
したがってこの地点を含むタイルのXYZパスは  /14/14336/6489.png

この【ア】【イ】【ウ】【エ】【オ】【カ】の意味を順に追うと以下になります。

【ア】地球上の地点を地図上に表すために、まずは経緯度をそのまま2次元直交座標系としてみてみると「横幅360, 縦幅180, 原点が中央にある長方形」内の座標といえます。
- この長方形で、兵庫県西脇市は (135, 35) にあります。
  - ちなみに原点 (0, 0) は、グリニッジ天文台から真南に向かって赤道とぶつかった地点 (ガーナ沖の海上) にあります。
【イ】扱いやすさのために経緯度をラジアンに変換して「横幅2π, 縦幅 π, 原点が中央にある長方形」内の座標に変換します。
- この長方形で、西脇市は ((135/180)π, (35/180)π) ≒ (2.356194, 0.610865) にあります。
- この座標系は、地球の球面を赤道から離れるほどに横に引き伸ばします。なぜなら、地球の半径 R とすると、赤道の円周は 2πR で、地球を北緯 35°で輪切りにした切り口の円周は 2πRcos((35/180)π) ですが、このどちらの円周もこの座標系では同じ長さになるためです。なので、日本にある 100m 四方の公園は赤道上にある 100m 四方の公園に比べて 1/cos((35/180)π) ≒ 1.2 倍も横に引き伸ばされます。
- よって、この座標系ではあらゆる緯度の正方形の公園を正方形として表示できません。赤道上の公園に合わせるなら日本の正方形の公園は横長になってしまい、日本の公園に合わせるなら赤道上の公園は縦長になってしまいます。
【ウ】そこで、赤道から離れるほど横に引き伸ばされる分、縦にも引き伸ばすことにします (＝メルカトル図法)。そうすればあらゆる緯度の正方形の公園が正方形になります。
- ただ、赤道から離れるほど引き伸ばしが大きくなるので、アイスランド (北緯65°) などにある公園は縦にも横にも 1/cos((65/180)π) ≒ 2.4 倍くらい引き伸ばされますが、許容することにします。
- ではこのとき西脇市の y 座標はどうなるか考えると、西脇市の真南にある赤道上の点から西脇市まで北上する間のどの地点も緯度 θ に対して 1/cos(θ) 倍に引き伸ばすので、これを積分する必要があります。よって、西脇市の y 座標は ∫_{0}^{(35/180)π} 1/cos(θ) dθ になります。ここで、解析的に ∫_{0}^{φ} 1/cos(θ) dθ = log(tan(π/4 + φ/2)) なので、北緯 φ ラジアンの地点の y 座標は log(tan(π/4 + φ/2)) になります。
- よって、この長方形 (といいつつ縦方向にはどこまでも伸びていますが) で、西脇市は ((135/180)π, log(tan(π/4 + (35/360)π))) ≒ (2.356194, 0.652837) にあります。
【エ】これで地球上の地点を具合よく平面にできたので、扱いやすさのため、縦も横も 2π で割って横幅を 1 に正規化します。
- 西脇市は (135/360, log(tan(π/4 + (35/360)π))/(2π)) ≒ (0.375, 0.103902) にあります。
【オ】ここからは地図タイルとしての運用管理面を考えていきます。まず、現時点で長方形は縦方向にどこまでも伸びています (北極点/南極点に近づくほどいくらでも縦に引き伸ばされるため)。これは管理しづらいので、縦幅を [-0.5, 0.5] の範囲でトリミングします。これで世界が辺の長さ 1 の正方形になります。トリミングによって緯度 85.0511° を超える地点はサポート範囲外になってしまいましたが、それほどの極地の地図はあまり利用しないものとして許容します。さらに、タイル番号を自然数にして左上隅から描画するのに便利なため、原点を正方形の中央から左上に移動して y 軸の向きを逆にします。
- この正方形で、西脇市は (0.375 - (-0.5), - (0.103902 - 0.5)) = (0.875, 0.396098) にあります。
  - ちなみに左上に移動した原点 (0, 0) は経緯度でいうと (180°, 85.0511°) で、ロシアのウランゲリ島から 1500km くらい真北の地点にあります。1500km という南北差は北海道の北端から九州の南端より短いですが、Google マップではウランゲリ島から表示が切れる北端まで日本が何個も入りそうなほど海が広がっていて、極地ほど拡大されていることがわかります。
【カ】次にいよいよ、必要な範囲の地図画像 (のための情報) だけを取り寄せられるように、グリッドを切ります。いま世界が辺の長さ 1 の正方形になっていますが、これを 2×2 のタイルに切ることをズームレベル 1 とし、左上, 右下, 左下, 右下のタイル番号を (0, 0), (1, 0), (0, 1), (1, 1) とします。同様に、世界を 2^{z}×2^{z} のタイルに切ることをズームレベル z とします。
- 赤道の長さ 40000km を 2^14 = 16384 で割ると 2.44km となり、日本の緯度では 1 タイルに収まる範囲が約 2km 四方になるので、ズームレベル 14 だと多くの細かい道も表示しやすいです。正方形の辺の長さを 1 から 2^14 = 16384 にすると兵庫県西脇市は (14336.000, 6489.667) にあります。よって、ズームレベル 14 における西脇市のタイルは (14336, 6489) 及び周辺になり、タイルパスは /14/14336/6489.png になります。