真正的墨卡托投影轉換詳解

所謂天下文章一大抄，看你會抄不會抄。網上對於墨卡托投影的解釋比較多，但是我發現很多都是是漏洞百出，無腦抄。例如：

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，如有侵權，請聯系刪除。