大地經緯度座標系與Web墨卡託座標系的轉換
1. 概述
我在《大地經緯度座標與地心地固座標的的轉換》這篇文章中已經論述了大地座標系/地理座標系的概念,簡單來說就是由經度、緯度以及高程(BLH)確定的座標系,它是一種曲面座標。
然而,在實際使用過程中我們用的最多的還是平面座標,並且單位最好與常用的長度單位(米)一致。所以就產生了從曲面到平面的轉換,這個過程也叫做投影,轉換的結果也就是投影平面座標系。我在《GDAL座標轉換》這篇文章中詳細論述了我們國內常用的三種投影平面座標系:橫軸墨卡託投影,高斯-克呂格投影和UTM投影。本質上來說,高斯-克呂格投影和UTM投影其實都是橫軸墨卡託投影,橫軸墨卡託投影也是用的最為廣泛的地圖投影方式。
但是在GIS,尤其是WebGIS領域中,橫軸墨卡託投影的使用遠沒有Web墨卡託投影方式用的多。最重要的原因是Web墨卡託投影的轉換演算法比橫軸墨卡託投影要簡單很多,符合Web的輕量化的特點。
2. 實現
Web墨卡託投影是橫軸墨卡託投影的特化版,要完全搞清楚Web墨卡託投影就必須得先搞清楚橫軸墨卡託投影,不過橫軸墨卡託投影實在太複雜了,但是我們可以定性地去理解。它的計算過程大概可以這樣理解:
在X方向上,為了保證投影到平面後經線和緯線仍然垂直,那麼每條緯線都會按照赤道周長展開,也就是\(2*PI*r = 2*20037508.3427892\)。由於原點位於平面中心,那麼可以算得X軸的取值範圍:[-20037508.3427892,20037508.3427892]。經度與投影后X長是簡單的線性關係。
在Y方向上,則需要藉助於墨卡託投影公式。為了保證投影的結果是正方形,那麼就把Y軸的取值範圍也取值成[-20037508.3427892,20037508.3427892]之間。這樣做沒什麼道理,純粹是為了希望投影的結果是正方形,便於切片。最後,通過墨卡託投影公式進行反算,得到的經緯度範圍就是[-85.05112877980659,85.05112877980659]。也就是這種投影方式,大於這個範圍是失效的。
參考Cesium的具體實現如下:
#include <iostream> //#include <eigen3/Eigen/Eigen> //#include <osgEarth/GeoData> using namespace std; const double epsilon = 0.000000000000001; const double pi = 3.14159265358979323846; const double d2r = pi / 180; const double r2d = 180 / pi; const double a = 6378137.0; //橢球長半軸 const double f_inverse = 298.257223563; //扁率倒數 const double b = a - a / f_inverse; //const double b = 6356752.314245; //橢球短半軸 const double e = sqrt(a * a - b * b) / a; //墨卡託範圍[-PI, PI]->大地緯度範圍[-PI/2, PI/2] static double mercatorAngleToGeodeticLatitude(double mercatorAngle) { return pi / 2.0 - (2.0 * atan(exp(-mercatorAngle))); //return 2.0 * atan(exp(mercatorAngle)) - pi / 2.0; } //Web墨卡託投影所支援的最大緯度(北和南) static double maximumLatitude = mercatorAngleToGeodeticLatitude(pi); //大地緯度範圍[-PI/2, PI/2]->墨卡託範圍[-PI, PI] static double geodeticLatitudeToMercatorAngle(double latitude) { // Clamp the latitude coordinate to the valid Mercator bounds. if (latitude > maximumLatitude) { latitude = maximumLatitude; } else if (latitude < -maximumLatitude) { latitude = -maximumLatitude; } double sinLatitude = sin(latitude); return 0.5 * log((1.0 + sinLatitude) / (1.0 - sinLatitude)); } void Blh2Wmc(double &x, double &y, double &z) { x = x * d2r * a; y = geodeticLatitudeToMercatorAngle(y * d2r) * a; } void Wmc2Blh(double &x, double &y, double &z) { //var oneOverEarthSemimajorAxis = this._oneOverSemimajorAxis; x = x / a * r2d; y = mercatorAngleToGeodeticLatitude(y / a) * r2d; } int main() { double x = 113.6; double y = 38.8; double z = 100; printf("%.10lf\n", maximumLatitude * r2d); printf("原大地經緯度座標:%.10lf\t%.10lf\t%.10lf\n", x, y, z); Blh2Wmc(x, y, z); printf("Web墨卡託座標:%.10lf\t%.10lf\t%.10lf\n", x, y, z); Wmc2Blh(x, y, z); printf("轉回大地經緯度座標:%.10lf\t%.10lf\t%.10lf\n", x, y, z); }
最終執行的結果:
通過GlobalMapper中的座標轉換工具對照的結果如下:
兩者結果基本一致。