墨卡托投影

墨卡托投影将地球球面投影到一个圆柱体柱面上，将地球看作一个正球体时，以 O O O为地球球心，从球心向外辐射射线，与地球外接圆柱面交与 P ′ P' P′。
设纬度为 ϕ \phi ϕ，经度为 λ \lambda λ，其中：

ϕ ∈ ( − π 2 , π 2 ) \phi\in(-\frac{\pi}{2},\frac{\pi}{2}) ϕ∈(−2π​,2π​) λ ∈ ( − π , π ) \lambda\in(-\pi,\pi) λ∈(−π,π)

对于经纬度为( ϕ \phi ϕ， λ \lambda λ)的坐标点，投影到圆柱面的坐标为( R tan ⁡ ϕ R\tan\phi Rtanϕ， R λ R\lambda Rλ)，坐标轴分别与柱面高以及柱面弧平行，设经度轴为 x x x轴，纬度轴为 y y y轴。
投影后， x x x最大取值为赤道长，即 x m a x = 2 π R x_{max}=2\pi R xmax​=2πR，而 y m a x = 2 R tan ⁡ ϕ m a x y_{max}=2R\tan\phi_{max} ymax​=2Rtanϕmax​，其中 R R R为球半径。
由于投影将圆柱面投影成正方形，长宽相等，则：

x m a x = y m a x x_{max}=y_{max} xmax​=ymax​

即：

2 π R = 2 R tan ⁡ ϕ m a x ( 1 ) 2\pi R=2R\tan\phi_{max} (1) 2πR=2Rtanϕmax​(1)

墨卡托投影防畸变处理

墨卡托投影将地球映射成了一个经纬度均匀分布的坐标系，由于将球体直接投影到柱面时经度是均匀的而纬度是不均匀的，所以地图需要对纬度进行防畸变处理。下面讨论纬度 ϕ \phi ϕ与畸变处理后地图纵坐标 y y y的函数关系。

P P P点在地球上的经度为 λ \lambda λ，纬度为 ϕ \phi ϕ，设 Q Q Q点为地球上与 P P P点无限接近的一点，则 Q Q Q点经纬度分别为 λ + Δ λ \lambda+\Delta\lambda λ+Δλ、 ϕ + Δ ϕ \phi+\Delta\phi ϕ+Δϕ。
P P P、 M M M所在纬线圈的半径
r = R cos ⁡ ϕ r=R\cos\phi r=Rcosϕ
其中 R R R为赤道圆半径。

由于 P P P、 M M M在同一纬线圈上，所以 P M PM PM间的弧长为它们纬线圈半径与两点经度差之积，即：
P M = r ∗ Δ λ = R Δ λ cos ⁡ ϕ PM=r*\Delta\lambda=R\Delta\lambda\cos\phi PM=r∗Δλ=RΔλcosϕ Q M = R ∗ Δ ϕ QM=R*\Delta\phi QM=R∗Δϕ

P ′ P' P′、 Q ′ Q' Q′为地球上点 P P P、 Q Q Q由墨卡托投影变换得到的点，此坐标系下经度、纬度线的分布都是均匀的， x x x、 y y y为地图坐标系，不再代表经纬度。所以 P ′ M ′ = Δ x P'M'=\Delta x P′M′=Δx Q ′ M ′ = Δ y Q'M'=\Delta y Q′M′=Δy
因为投影后将每一个纬线圈都放大到与赤道圆相同大小，所以投影后 Δ x \Delta x Δx长度就是 P ′ P' P′、 M ′ M' M′经度夹角与赤道圆半径的乘积，即：
Δ x = R Δ λ \Delta x=R\Delta\lambda Δx=RΔλ
由于投影后需要保证航线角相同，所以：
P M Q M = P ′ M ′ Q ′ M ′ \frac{PM}{QM}=\frac{P'M'}{Q'M'} QMPM​=Q′M′P′M′​即 R Δ λ cos ⁡ ϕ R Δ ϕ = Δ x Δ y \frac{R\Delta\lambda\cos\phi}{R\Delta\phi}=\frac{\Delta x}{\Delta y} RΔϕRΔλcosϕ​=ΔyΔx​ ⇒ Δ y Δ ϕ = Δ x Δ λ cos ⁡ ϕ \Rightarrow\frac{\Delta y}{\Delta\phi}=\frac{\Delta x}{\Delta\lambda\cos \phi} ⇒ΔϕΔy​=ΔλcosϕΔx​ ⇒ Δ y Δ ϕ = R Δ λ Δ λ cos ⁡ ϕ = R cos ⁡ ϕ \Rightarrow\frac{\Delta y}{\Delta\phi}=\frac{R\Delta\lambda}{\Delta\lambda\cos\phi}=\frac{R}{\cos\phi} ⇒ΔϕΔy​=ΔλcosϕRΔλ​=cosϕR​
根据导数定义，上式即： y ′ ( ϕ ) = R cos ⁡ ϕ y^\prime(\phi)=\frac{R}{\cos\phi} y′(ϕ)=cosϕR​
下面积分求 y ( ϕ ) y(\phi) y(ϕ)就可以了:
y ( ϕ ) = ∫ y ′ ( ϕ ) d ϕ = R ∫ 1 cos ⁡ ϕ d ϕ y(\phi)=\int y^\prime(\phi)d\phi=R\int\frac{1}{\cos\phi}d\phi y(ϕ)=∫y′(ϕ)dϕ=R∫cosϕ1​dϕ = R ∫ 1 cos ⁡ 2 ϕ d sin ⁡ ϕ = R ∫ 1 1 − sin ⁡ 2 ϕ d sin ⁡ ϕ =R\int\frac{1}{\cos^2\phi}d\sin\phi=R\int\frac{1}{1-\sin^2\phi}d\sin\phi =R∫cos2ϕ1​dsinϕ=R∫1−sin2ϕ1​dsinϕ
设 t = sin ⁡ ϕ t=\sin\phi t=sinϕ，
y ( ϕ ) = R ∫ 1 1 − t 2 d t = R 2 ∫ ( 1 1 + t + 1 1 − t ) d t y(\phi)=R\int\frac{1}{1-t^2}dt=\frac{R}{2}\int(\frac{1}{1+t}+\frac{1}{1-t})dt y(ϕ)=R∫1−t21​dt=2R​∫(1+t1​+1−t1​)dt
易得：
y ( ϕ ) = R 2 ln ⁡ ∣ 1 + t 1 − t ∣ = R 2 ln ⁡ ∣ 1 + sin ⁡ ϕ 1 − sin ⁡ ϕ ∣ y(\phi)=\frac{R}{2}\ln|\frac{1+t}{1-t}|=\frac{R}{2}\ln|\frac{1+\sin\phi}{1-\sin\phi}| y(ϕ)=2R​ln∣1−t1+t​∣=2R​ln∣1−sinϕ1+sinϕ​∣ = R 2 ln ⁡ ∣ ( 1 + sin ⁡ ϕ ) 2 1 − sin ⁡ 2 ϕ ∣ = R ln ⁡ 1 + sin ⁡ ϕ cos ⁡ ϕ =\frac{R}{2}\ln|\frac{(1+\sin \phi)^2}{1-\sin^2\phi}|=R\ln\frac{1+\sin \phi}{\cos \phi} =2R​ln∣1−sin2ϕ(1+sinϕ)2​∣=Rlncosϕ1+sinϕ​
由反双曲正切函数 a r c tanh ⁡ x = 1 2 ln ⁡ 1 + x 1 − x arc\tanh x=\frac{1}{2}\ln\frac{1+x}{1-x} arctanhx=21​ln1−x1+x​可知， y ( ϕ ) y(\phi) y(ϕ)还可以表示成：
y ( ϕ ) = R a r c tanh ⁡ ( sin ⁡ ϕ ) \color{red}y(\phi)=Rarc\tanh(\sin\phi) y(ϕ)=Rarctanh(sinϕ)

或者由反双曲正弦函数 a r c sinh ⁡ x = ln ⁡ ( x + x 2 + 1 ) arc\sinh x=\ln(x+\sqrt{x^2+1}) arcsinhx=ln(x+x2+1 ​)可将 y ( ϕ ) y(\phi) y(ϕ)变形为：
y ( ϕ ) = R ln ⁡ 1 + sin ⁡ ϕ cos ⁡ ϕ = R ln ⁡ ( tan ⁡ ϕ + 1 cos ⁡ ϕ ) y(\phi)=R\ln\frac{1+\sin \phi}{\cos \phi}=R\ln(\tan\phi+\frac{1}{\cos\phi}) y(ϕ)=Rlncosϕ1+sinϕ​=Rln(tanϕ+cosϕ1​) = R ln ⁡ ( tan ⁡ ϕ + sin ⁡ 2 ϕ + cos ⁡ 2 ϕ cos ⁡ 2 ϕ ) =R\ln(\tan\phi+\frac{\sqrt{\sin^2\phi+\cos^2\phi}}{\sqrt{\cos^2\phi}}) =Rln(tanϕ+cos2ϕ ​sin2ϕ+cos2ϕ ​​) = R ln ⁡ ( tan ⁡ ϕ + tan ⁡ 2 ϕ + 1 ) =R\ln(\tan\phi+\sqrt{\tan^2\phi+1}) =Rln(tanϕ+tan2ϕ+1 ​) = R a r c sinh ⁡ ( tan ⁡ ϕ ) \color{red}=Rarc\sinh (\tan\phi) =Rarcsinh(tanϕ)

经过畸变后，(1)式可变为：
2 π R = 2 ∗ R a r c tanh ⁡ ( sin ⁡ ϕ m a x ) 2\pi R=2*Rarc\tanh(\sin\phi_{max}) 2πR=2∗Rarctanh(sinϕmax​)

解出：

ϕ m a x ≈ ± 85.0 5 ∘ \phi_{max}\approx\pm85.05^\circ ϕmax​≈±85.05∘

则墨卡托投影出来的地图纬度只能表示到 ± 85.0 5 ∘ \pm85.05^\circ ±85.05∘，此时经度取值范围为 ( − π , π ) (-\pi,\pi) (−π,π)。

地图瓦片最大最小经纬度计算

最大最小经度

设瓦片最小经度为 l o n n x m i n lon_{n_{xmin}} lonnxmin​​，最大经度为 l o n n x m a x lon_{n_{xmax}} lonnxmax​​，每行或列瓦片个数为 n n n瓦片经度编号为 n x n_x nx​，由于经度的分布是均匀的，所以（按照弧度计算）：
l o n n x m i n − ( − π ) 2 π = n x n \frac{lon_{n_{xmin}}-(-\pi)}{2\pi}=\frac{n_x}{n} 2πlonnxmin​​−(−π)​=nnx​​
很容易解出
l o n n x m i n = n x ∗ 2 π n − π \color{red}lon_{n_{xmin}}=\frac{n_x*2\pi}{n}-\pi lonnxmin​​=nnx​∗2π​−π

最大最小纬度

找出两组等比例关系放在等式左右两边：

1. 地球球面上北极到瓦片最小纬度的弧长，与地球半球（半个经线圈）弧长，即图中红色弧长与半圆总弧长之比
2. 经过墨卡托投影及放畸变处理后红色弧长投影与地图纵向总长度之比，即图中红色线段与线段总长度只比

列出等式，注意由于地图是正方形的，所以纵向总长度为 2 π R 2\pi R 2πR：

n y n = π R − y ( l a t n m i n ) 2 π R = π R − R a r c sinh ⁡ ( tan ⁡ l a t n m i n ) 2 π R \frac{n_y}{n}=\frac{\pi R-y(lat_{nmin})}{2\pi R}=\frac{\pi R-Rarc\sinh (\tan lat_{nmin})}{2\pi R} nny​​=2πRπR−y(latnmin​)​=2πRπR−Rarcsinh(tanlatnmin​)​
其中瓦片编号为 n y n_y ny​，从而很容易解出：
l a t n m i n = a r c tan ⁡ ( sinh ⁡ π ( 1 − 2 y n ) ) \color{red}lat_{nmin}=arc\tan(\sinh\pi(1-\frac{2y}{n})) latnmin​=arctan(sinhπ(1−n2y​))
最大经纬度直接将 n x n_x nx​、 n y n_y ny​替换成 n x + 1 n_x+1 nx​+1、 n y + 1 n_y+1 ny​+1即可