真正的墨卡托投影转换详解

所谓天下文章一大抄，看你会抄不会抄。网上对于墨卡托投影的解释比较多，但是我发现很多都是是漏洞百出，无脑抄。例如：

X轴：由于赤道半径为6378137米，则赤道周长为2*PI*r = 20037508.3427892，因此X轴的取值范围：[-20037508.3427892,20037508.3427892]。

2 * Math.PI * 6378137 = 40075016.0019724

X轴取值范围确实是[-20037508.3427892,20037508.3427892]，但是请不要闭着眼睛说这是周长（实际是半周长）

地理经度的取值范围是[-180,180]，纬度不可能到达90°，通过纬度取值范围为[20037508.3427892,20037508.3427892] ，反计算可得到纬度值为85.05112877980659。因此纬度取值范围是[-85.05112877980659，85.05112877980659]

这跳的也跳快了，为啥纬度取值范围也要半周长，反计算过程又在哪？

---

墨卡托投影的模型：地球近似为正球体，球体中心有个发亮的光点，光点将球面上的每个点都投影到正切赤道的圆柱内表面上，将圆柱体内表面展开就是一张世界地图。

形象的理解：有个特殊的乒乓球（地球），这个乒乓球球体内中心点会发光，乒乓球卡在塑料水管内部，把塑料水管水管按有乒乓球阴影部分锯开碾平后水管内表面就是世界地图。


在WebGIS中，经常有已知WGS84(EPSG:4326)经纬度[lng, lat]，在Canvas2D上要显示该点就要换算成墨卡托投影(EPSG:3857)，其实只要知道地球变成平面地图时的转换比率，那就能通过经纬度求出墨卡托的值。Openlayers以及Mapbox/sphericalmercator都有转换函数，但是这个计算公式怎么来的呢？

假设有个Q点离P点无限近，地球半径为 ，如下图


P的经纬度为 （ 和 都为弧度），相邻点Q的经纬度为 ，那么:

平行赤道的 纬线所在圆的半径： ，弧度 ，如下示意图


与纬线平行转换比率：
与经线平行转换比率：

由于经线转化后未发生形变，那么根据globe和projection对比图（ ）

得出 ，

再利用求导概念， ，那么：

与纬线平行转换比率：
与经线平行转换比率：

由于墨卡托投影为了保证投影后达到等角航线（航向角和经度所成的夹角保持不变），那么可以认定 。

等等，各位看到这里肯定会问，墨卡托设计的为什么需要等角航线？为什么是要和经度夹角而不是维度夹角？其实我写的时候也想了下，我们不妨大胆猜测下：

墨卡托设计出来的投影是为了航海用的，当时航海导航全靠指南针，指南针指向南北两极，与经线平行，那么要根据地图找到航行路线上的陆地的话，一定要根据地图上的所指的角度进行航线行驶（想象下航海电影中主角拿出地图，地图上压着个指南针，然后用尺子画出两点的直线，再量下这条直线和经线的夹角，最后指挥舵手往哪里哪里打几度方向。电影果然没欺骗我们）。

回归正题，由于 ，那么 可得出

那么根据

得 的原函数为 ， 的反函数 称为高德曼函数 , （这里由于为了和截图对应需要把 变为了 ）

在Web中，需要对地图进行缩放，为了更好的计算缩放比例，设定横宽比1。之前由于地图的x的取值范围为 , 那么y的取值范围也为 ，于是可以求得Web下的墨卡托投影的最大纬度值为：

而EPSG:4326转EPSG:3857就很好计算了

// 弧度版
const x = R * lng 
const y = R * Math.log(Math.tan((Math.PI*0.25) + (0.5 * lat)))

// 角度版
const radians = Math.PI / 180
const x = R * lng * radians
const y = R * Math.log(Math.tan((Math.PI*0.25) + (0.5 * lat * radians)))

---

参考资料：

本文转自 https://zhuanlan.zhihu.com/p/326955505，如有侵权，请联系删除。