1. 程式人生 > 其它 >真正的墨卡託投影轉換詳解

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

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

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]

這跳的也跳快了,為啥緯度取值範圍也要半周長,反計算過程又在哪?


墨卡託投影的模型:地球近似為正球體,球體中心有個發亮的光點,光點將球面上的每個點都投影到正切赤道的圓柱內表面上,將圓柱體內表面展開就是一張世界地圖。

形象的理解:有個特殊的乒乓球(地球),這個乒乓球球體內中心點會發光,乒乓球卡在塑料水管內部,把塑料水管水管按有乒乓球陰影部分鋸開碾平後水管內表面就是世界地圖。

![](data:image/svg+xml;utf8,)

在WebGIS中,經常有已知WGS84(EPSG:4326)經緯度[lng, lat],在Canvas2D上要顯示該點就要換算成墨卡託投影(EPSG:3857),其實只要知道地球變成平面地圖時的轉換比率,那就能通過經緯度求出墨卡託的值。

Openlayers以及Mapbox/sphericalmercator都有轉換函式,但是這個計算公式怎麼來的呢?

假設有個Q點離P點無限近,地球半徑為 ,如下圖

![](data:image/svg+xml;utf8,)

P的經緯度為 ( 和 都為弧度),相鄰點Q的經緯度為 ,那麼:

平行赤道的 緯線所在圓的半徑: ,弧度 ,如下示意圖

![](data:image/svg+xml;utf8,)


與緯線平行轉換比率:
與經線平行轉換比率:

由於經線轉化後未發生形變,那麼根據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,如有侵權,請聯絡刪除。