初めに

地図タイルダウンロードツール作成にあたり、緯度経度をタイル座標に変換する必要があり、計算式等を調べたが理解できる解説がなかったので、自分で式を作りました。もっと簡易な方法があるかもしれません

本解説では、緯度経度を一度ＸＹ座標へ変換していますがこれは個人的に理解しやすかったからです

変換方法

ざっくりと次の手順で計算している。

北緯約８５度をＹ座標０としたときの南緯約８５度のＹ座標を算出する
西経１８０度をＸ座標０としたときの東経１８０度のＸ座標を算出する
ズームレベル「Ｚ」におけるＸＹ座標それぞれのタイルの最大数を算出する
１、２をそれぞれ３で割り、タイル１枚あたりの間隔を算出する
変換したい座標のうち緯度のＹ座標を算出する
変換したい座標のうち経度のＸ座標を算出する
５、６をそれぞれ４で割る
７の小数点以下を切り捨てた値がタイルのＸ（経度）Ｙ（緯度）座標となる

各手順を説明する前に

地図タイルとは主にＷｅｂメルカトルで用いられる地図を表すためのパーツのようなもの
詳しい説明は国土地理院を参照

前提

・測地系はＷＧＳ８４
　赤道半径(m)：$a = 6,378,137$
・ 地球は真球として扱う
・今回は緯度35.685227、経度139.753347(皇居)のタイル座標を求めてみる。

手順１、２：ＸＹ座標の最大値を算出

緯度南端

$S = -85.05112877980659$

Ｙ座標南端

$Ys = a * Atanh(Sin(S * π / 180))$
$Ys = -20,037,508.34278924$

赤道を基準とした値から北緯約８５度を基準とした値に変換する
$Ys=40,075,016.68557848$

経度東端

$E = 180$

Ｘ座標東端

$Xe = a * E * π / 180$
$Xe = 20,037,508.34278924$

本初子午線を基準とした値から西経180度を基準とした値に変換する
$Xe=40,075,016.68557848$

手順３：タイル最大数

タイルはズームレベルが１つ上がる毎に２倍されていく
従ってズームレベルｎにおけるＸＹ座標それぞれのタイル最大数は$2^n$

手順４：タイル間隔

ズームレベルｎにおけるタイル間隔は以下の通り
$40,075,016.68557848/2^n$

手順５、６：ＸＹ座標の算出

手順１、２の最大値を求める緯度経度を今回タイル座標を求める緯度経度へ置き換える

求める緯度のＹ座標

$N = 35.685227$
$Yn = a * Atanh(Sin(N * π / 180))$
$Yn = 4,257,395.08697871$

赤道を基準とした値から北緯約８５度を基準とした値に変換する
$Yn=20,037,508.34278924 - 4,257,395.08697871$
$Yn=$$15,780,113.25581053$

求める経度のＸ座標

$W = 139.753347$
$Xw = a * W * π / 180$
$Xw = 15,557,271.42469566$

本初子午線を基準とした値から西経180度を基準とした値に変換する
$Xw=20,037,508.34278924 + 15,557,271.42469566$
$Xw=$$35,594,779.7674849$

手順７、８

手順５、６の値を手順４の値で割って、小数点以下を切り捨てた値がタイル座標となる
緯度35.685227、経度139.753347(皇居)のズームレベル14のタイル座標は次のとおり

ズームレベル：$n = 14$
タイル１枚の間隔：$40,075,016.68557848/2^n=2,445.98490512$
Ｙ座標：$15,780,113.25581053$
Ｘ座標：$35,594,779.7674849$

タイル座標Ｙ

$15,780,113.25581053 / 2,445.98490512 = 6,451.43525733$
小数点以下を切り捨てて$6,451$

タイル座標Ｘ

$35,594,779.7674849 / 2,445.98490512 = 14,552.33010350$
小数点以下を切り捨てて$14,552$

最後に

タイル座標は次の式で求められる。

赤道半径：$a=6,378,137$
ズームレベル：$n$
緯度：$Lat$
経度：$Lon$
Ｙ座標の最大値：$Ys=40,075,016.68557848$
Ｙ座標の中央値：$Yq=20,037,508.34278924$
Ｘ座標の最大値：$Xe=40,075,016.68557848$
Ｘ座標の中央値：$Xq=20,037,508.34278924$

$タイル座標Ｙ=(Yq - a * Atanh(Sin(Lat * π / 180)))/(Ys / 2^n)$
$タイル座標Ｘ=(Xq + a * Lon * π / 180)/(Xe / 2^n)$