こんにちは
地図のwebメルカトル投影計算式をまとめました（投影は、経度 $\lambda$ ⇔ $x$、緯度 $\varphi$ ⇔ $y$）。web地図タイル座標の計算でも用いられます。下記の間の変換です（真球モデル）1。
- EPSG:4326 WGS 84 geographical
- EPSG:3857 WGS84 Web Mercator (auxiliary sphere)
緯度
緯度 $\varphi$ ⇔ $y$ の変換・逆変換は（$R$ は地球長半径、$\operatorname{gd} x$ は Gudermannian functionです）、
\frac{y}{R} = \operatorname{gd}^{-1} \varphi = \operatorname{arctanh} \left( \sin\varphi \right) = 2 \operatorname{arctanh} \left( \tan \frac{\varphi}{2} \right)\\
\varphi = \operatorname{gd} \frac{y}{R} = \arcsin \left( \tanh \frac{y}{R} \right) = 2 \arctan \left( \tanh \frac{y}{2 R} \right)
\operatorname{arctanh} x = \frac{1}{2} \ln \left| \frac{1+x}{1-x} \right|
なお Gudermannian function は、上記のような Braun投影 （$\tan \frac{\varphi}{2}$）を経る形を含めて、様々な形で表現可能です2。
経度
経度 $\lambda$ ⇔ $x$ については、
\lambda = \frac{x}{R}
スケール係数
EPSG:3857（真球モデル）のスケール係数は（東西EWおよび南北NS方向）、地球が回転楕円体（離心率 $e$）であることを考慮する必要があり、
k_\text{EW} = \sec \varphi \ \sqrt{1-e^2 \sin^2 \varphi} \\
k_\text{NS} = \sec \varphi \ \frac{\left(1-e^2 \sin^2 \varphi\right)^\frac{3}{2}}{1-e^2} \\
\sec \varphi = \frac{1+\tan^2 \frac{\varphi}{2}}{1-\tan^2 \frac{\varphi}{2}} = \cosh \frac{y}{R} \\
\sin \varphi = \frac{2 \tan \frac{\varphi}{2}}{1+\tan^2 \frac{\varphi}{2}} = \tanh \frac{y}{R}
EPSG:3395（回転楕円体モデル）とは異なり、$k_\text{EW} = k_\text{NS}$ は成立しません。ただし、近似計算では、$e \simeq 0$ として計算することも多いです。
k_\text{EW} \simeq k_\text{NS} \simeq \sec \varphi