墨卡托投影原理及瓦片公式推导

墨卡托投影

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

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

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

xmax=ymaxx_{max}=y_{max} xmax=ymax

即：

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

墨卡托投影防畸变处理

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

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

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

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

或者由反双曲正弦函数arcsinh⁡x=ln⁡(x+x2+1)arc\sinh x=\ln(x+\sqrt{x^2+1})arcsinhx=ln(x+x2+1)可将y(ϕ)y(\phi)y(ϕ)变形为：
y(ϕ)=Rln⁡1+sin⁡ϕcos⁡ϕ=Rln⁡(tan⁡ϕ+1cos⁡ϕ)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)=Rln⁡(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ϕ)=Rln⁡(tan⁡ϕ+tan⁡2ϕ+1)=R\ln(\tan\phi+\sqrt{\tan^2\phi+1})=Rln(tanϕ+tan2ϕ+1)=Rarcsinh⁡(tan⁡ϕ)\color{red}=Rarc\sinh (\tan\phi)=Rarcsinh(tanϕ)

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

解出：

ϕmax≈±85.05∘\phi_{max}\approx\pm85.05^\circ ϕmax≈±85.05∘

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

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

最大最小经度

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

最大最小纬度

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

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

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

nyn=πR−y(latnmin)2πR=πR−Rarcsinh⁡(tan⁡latnmin)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)
其中瓦片编号为nyn_yny，从而很容易解出：
latnmin=arctan⁡(sinh⁡π(1−2yn))\color{red}lat_{nmin}=arc\tan(\sinh\pi(1-\frac{2y}{n})) latnmin=arctan(sinhπ(1−n2y))
最大经纬度直接将nxn_xnx、nyn_yny替换成nx+1n_x+1nx+1、ny+1n_y+1ny+1即可