投影
對于地圖制圖:原面為地球的旋轉橢球面,是三維的;承受面(對瓦片地圖而言為瓦片)為二維平面的。如何在原面與承受面之間建立點、線、面的一一對應關系是地圖制圖的必須過程,這一過程通常稱之為:地圖投影。
目前地圖投影方法通常采用數學分析法,也即:在原面和承受面之間建立點與點的函數關系。地球面上的點用地理坐标經緯度(L,B)來确定,平面上的點一般多用平面直角坐标系(x,y)确定。投影即建立函數關系:
式中
,
是單值而連續的函數,當地球面上的點連續移動時,平面上的對應點也同時連續移動。
投影方法很多,按投影變形性質分:
- 等角投影:投影承受面上任意兩方向線夾角同地球面相應夾角;
- 等面積投影:投影承受面上的有限面積與地球面上相應面積相等;
- 任意投影:即不等角又不等面積的投影,包括等距投影。
按正軸投影經緯網形狀分:
- 方位投影:緯線投影為同心圓,經線投影為同心圓的直徑,兩經線間夾角與相應經差相等;
- 圓柱投影:緯線投影為平行線,經線投影為與緯線垂直而且間隔相等的平行線;
- 圓錐投影:緯線投影為同心圓弧,經線投影為同心圓半徑,兩經線間夾角與相應經差成正比;
- 還有多圓錐投影、僞方位投影、僞圓柱投影、僞圓錐投影,不再一一介紹。
按投影面與地球面位置關系分:
- 正軸投影:正軸方位投影,投影面與地軸垂直;正軸圓柱和正軸圓錐投影,圓柱軸或圓錐軸與地軸重合;
- 橫軸投影:橫軸方位投影,投影平面與地球赤道一直徑垂直;橫軸圓柱或橫軸圓錐投影,圓柱軸或圓錐軸同赤道一直徑重合;
- 斜軸投影:斜軸方位投影,投影平面與地軸和赤道直徑以外的任意大圓直徑垂直;斜軸圓柱或斜軸圓錐投影,圓柱軸或圓錐軸與地軸和赤道直徑以外的任意大圓直徑重合。
投影是個大課題,專業程度要求高,推薦孫達、蒲英霞老師的《地圖投影》。對于瓦片地圖而言,由于要覆寫整個世界範圍且要保證制圖美觀,能選擇的投影方法并不多,是以瓦片地圖的投影部分要簡化很多。目前主流瓦片地圖一般采用Web墨卡托投影、經緯度直投及百度地圖獨有的百度墨卡托投影。
Web墨卡托投影
Web墨卡托投影方法由谷歌提出,并最先應用于谷歌地圖。說起web墨卡托投影,就不得不先說說墨卡托投影。
墨卡托投影由荷蘭地圖學家墨卡托(G.Mercator)于1569年創立。假想一個與地軸方向一緻的圓柱切于地球,按等角條件将經緯網投影到圓柱面上,将圓柱面展為平面後,即為墨卡托投影。
為了更好地了解墨卡托投影,上圖如圖1:
圖1 墨卡托投影示例
對于墨卡托投影,可以想象為地球被圍在一個中空的圓柱裡,其赤道與圓柱相接觸,然後再假設地球中心有一盞燈,把球面上的圖形投影到圓柱體上,再把圓柱體展開,這就是一幅标準緯線為零度(即赤道)的“墨卡托投影”繪制出的世界地圖。
墨卡托投影為等角正軸圓柱投影,由于該投影具有等角性質、經緯網構成簡單和等角航線投影為直線的有限,至今仍廣泛應用于航海和航空領域。它的缺點也非常明顯:面積畸變,而且越往地球南北兩級面積畸變越大。從上圖中可見格陵蘭島與非洲的面積差不多大小,事實上格陵蘭島面積僅為非洲面積的1/14。
但由于墨卡托投影能滿足全世界範圍内投影,等角性質能保證地物形狀不會被改變,且面積畸變大的兩級地區也非人類生産生活集中區域,瓦片地圖的先行者谷歌依然選擇了墨卡托投影作為谷歌地圖的投影方法。
不過為了投影計算簡單,降低二次開發難度,谷歌對墨卡托投影做了簡化:将地球視為球體而非橢球體。為了區分,通常将谷歌簡化後的墨卡托投影稱為:Web墨卡托投影(亦稱:球面墨卡托投影)。
Web墨卡托投影将地球當做球體而非橢球體,在數學意義上并非一種嚴密的投影方法。正因為如此,在Web墨卡托投影剛被推出時并不被業界所認可,EPSG也隻是戲谑性地給了它一個非正式編号900913(挺有趣的編号,跟google是不是很像)。
但采用Web墨卡托投影的谷歌一經推出,後續的瓦片地圖,如:國内的高德地圖、騰訊地圖、超圖、360地圖、搜狗地圖等, 國外的ArcGIS online、MapQuest、BingMaps等紛紛采用該投影方法。Web墨卡托投影已成為瓦片地圖事實上的投影标準,後來EPSG也順勢而為為它賦予了正式編号3857。
Web墨卡托投影面坐标系
Web墨卡托投影面展開後,以赤道作為标準緯線,本初子午線作為中央經線,兩者交點為坐标原點,向東向北為正,向西向南為負建立投影面直角坐标系,機關為米。
X軸取值範圍:赤道半徑r取為6378137.0,則赤道周長為2*π*r = 2*20037508.3427890167,是以X軸的取值範圍:[-20037508.3427890167,20037508.3427890167]。
Y軸取值範圍:由墨卡托投影方法可知,當緯度接近兩極即±90°時,Y值趨向于無窮,無法明确Y軸取值範圍。為了解決這個為題,也為了瓦片切片友善,谷歌工程師們把Y軸取值範圍也限定在[-20037508.3427890167,20037508.3427890167],如此既能明确Y軸取值範圍,也能保證第0級瓦片為正方形,友善各級瓦片切片。
在确定投影面坐标系X/Y軸取值範圍後,反算易知對應地理坐标系經度取值範圍[-180, 180],緯度取值範圍[-85.05112877980659, 85.05112877980659]。
你一定發現了墨卡托投影并未覆寫全球範圍,不過沒關系企鵝和北極熊們并不在乎。
如下,給出Web墨卡托投影公式:
其中r為地球球體半徑,取值為6378137.0m;L為經度,機關為°;L為緯度,機關為°。
經緯度直投
經緯度直投可以說是一種最簡單的投影方法,看看投影公式就知道了。
經緯度直投的特點是:經緯網在投影平面上按等間距依次排開,沒有複雜的數學變換(簡直是碼農的福音)。
但經緯度直投屬于任意投影,既不等角也不等面積。在低緯度地區長度、角度、面積、形狀變化比較小,但越往高緯度,水準長度畸變越嚴重,很小的緯圈都變得和赤道一樣長。受此影響,高緯度地區地物要素變形較為嚴重,十字路口變成了X型路口,矩形變成了上寬下窄的倒梯形。
上兩張來自邁高圖的截圖,如圖2、圖3(截圖地點位于哈爾濱城區某地)。直覺感受下經緯度直投在高緯度地區的畸變。
圖2 墨卡托投影下的十字路口
圖3 經緯度直投下的十字路口
兩張截圖均來自同一地點,可見墨卡托投影下的标準十字路口,在經緯直投下變成了X型路口。饒是如此,仍然有瓦片地圖采用經緯度直投作為制圖投影方法,如:天地圖,谷歌地球(二維)。
百度墨卡托投影
百度地圖采用自家獨有投影方法,其核心還是墨卡托投影,具有等角特性。
但它的投影計算方法與Web墨卡托投影公式不同 ,具體可參見百度地圖Js API V1.3。這就導緻了百度地圖的投影與Web墨卡托投影有一定的誤差(投影平面坐标系原點與Web墨卡托投影平面原點不重合,而且每個地區不重合偏移都不同),是以百度地圖的投影計算隻能通過抽取、轉譯Js API V1.3的相關代碼來實作。可參考以下Qt代碼。
/*h file*/
#include <QPoint>
class BaiduMercatorProj
{
public:
BaiduMercatorProj();
~BaiduMercatorProj();
QPoint convertLL2MC(const QPoint& gcp);
QPoint convertMC2LL(const QPoint& pcp);
private:
QPoint convertor(const QPoint& cM, const double cN[10]);
double getRange(double val, double minVal, double maxVal);
double getLoop(double val, double minVal, double maxVal);
private:
const double MCBAND[6] = { 12890594.86, 8362377.87, 5591021, 3481989.83,1678043.12, 0 };
const double LLBAND[6] = { 75, 60, 45, 30, 15, 0 };
const double MC2LL[6][10] = {
{1.410526172116255e-8, 0.00000898305509648872,
-1.9939833816331, 200.9824383106796,
-187.2403703815547, 91.6087516669843,
-23.38765649603339, 2.57121317296198,
-0.03801003308653, 17337981.2
},
{-7.435856389565537e-9, 0.000008983055097726239,
-0.78625201886289, 96.32687599759846,
-1.85204757529826, -59.36935905485877,
47.40033549296737, -16.50741931063887,
2.28786674699375, 10260144.86
},
{-3.030883460898826e-8, 0.00000898305509983578,
0.30071316287616, 59.74293618442277,
7.357984074871, -25.38371002664745,
13.45380521110908, -3.29883767235584,
0.32710905363475, 6856817.37
},
{-1.981981304930552e-8, 0.000008983055099779535,
0.03278182852591, 40.31678527705744,
0.65659298677277, -4.44255534477492,
0.85341911805263, 0.12923347998204,
-0.04625736007561, 4482777.06
},
{3.09191371068437e-9, 0.000008983055096812155,
0.00006995724062, 23.10934304144901,
-0.00023663490511, -0.6321817810242,
-0.00663494467273, 0.03430082397953,
-0.00466043876332, 2555164.4
},
{2.890871144776878e-9, 0.000008983055095805407,
-3.068298e-8, 7.47137025468032,
-0.00000353937994, -0.02145144861037,
-0.00001234426596, 0.00010322952773,
-0.00000323890364, 826088.5}
};
const double LL2MC[6][10] = {
{-0.0015702102444, 111320.7020616939,
1704480524535203, -10338987376042340,
26112667856603880, -35149669176653700,
26595700718403920, -10725012454188240,
1800819912950474, 82.5
},
{0.0008277824516172526, 111320.7020463578,
647795574.6671607, -4082003173.641316,
10774905663.51142, -15171875531.51559,
12053065338.62167, -5124939663.577472,
913311935.9512032, 67.5
},
{0.00337398766765, 111320.7020202162,
4481351.045890365, -23393751.19931662,
79682215.47186455, -115964993.2797253,
97236711.15602145, -43661946.33752821,
8477230.501135234, 52.5
},
{0.00220636496208, 111320.7020209128,
51751.86112841131, 3796837.749470245,
992013.7397791013, -1221952.21711287,
1340652.697009075, -620943.6990984312,
144416.9293806241, 37.5
},
{-0.0003441963504368392, 111320.7020576856,
278.2353980772752, 2485758.690035394,
6070.750963243378, 54821.18345352118,
9540.606633304236, -2710.55326746645,
1405.483844121726, 22.5
},
{-0.0003218135878613132, 111320.7020701615,
0.00369383431289, 823725.6402795718,
0.46104986909093, 2351.343141331292,
1.58060784298199, 8.77738589078284,
0.37238884252424, 7.45
}
};
};
/*cpp file*/
#include "BaiduMercatorProj.h"
BaiduMercatorProj::BaiduMercatorProj()
{
}
BaiduMercatorProj::~BaiduMercatorProj()
{
}
QPoint BaiduMercatorProj::convertLL2MC(const QPoint& pt)
{
double cN[10];
QPoint cL(getLoop(pt.x(), -180, 180), getRange(pt.y(), -74, 74));
const staitc int len = 6;//sizeof LLBAND / sizeof(double);
int cNok = false;
for (int i = 0; i < len; i++)
{
if (cL.y() >= LLBAND[i])
{
memcpy_s(cN, 80, LL2MC[i], 80);
cNok = true;
break;
}
}
if (!cNok) {
for (int i = len - 1; i >= 0; i--)
{
if (cL.y() <= -LLBAND[i]) {
memcpy_s(cN, 80, LL2MC[i], 80);
break;
}
}
}
return convertor(pt, cN);
}
QPoint BaiduMercatorProj::convertMC2LL(const QPoint& pt)
{
double cO[10];
QPoint cM(fabs(pt.x()), fabs(pt.y()));
const static int len = 6;// sizeof MCBAND / sizeof(double);
for (int i = 0; i < len; i++)
{
if (cM.y() >= MCBAND[i])
{
memcpy_s(cO, 80, MC2LL[i], 80);
break;
}
}
return convertor(pt, cO);
}
Point BaiduMercatorProj::convertor(const Point& cM, const double cN[10])
{
double T = cN[0] + cN[1] * fabs(cM.x());
double cL = fabs(cM.y()) / cN[9];
double cO = cN[2] + cN[3] * cL + cN[4] * cL * cL + cN[5] * cL *
cL * cL + cN[6] * cL * cL * cL * cL + cN[7] * cL *
cL * cL * cL * cL + cN[8] * cL * cL * cL * cL *
cL * cL;
T *= (cM.x() < 0 ? -1 : 1);
cO *= (cM.y() < 0 ? -1 : 1);
return Point(T, cO);
}
double BaiduMercatorProj::getLoop(double val, double minVal, double maxVal)
{
double loopVal = 0.0;
while (val > maxVal) {
loopVal -= maxVal - minVal;
}
while (val < minVal) {
loopVal += maxVal - minVal;
}
return loopVal;
}
double BaiduMercatorProj::getRange(double val, double minVal, double maxVal)
{
double rangeVal = 0.0;
rangeVal = max(val, minVal);
rangeVal = min(val, maxVal);
return rangeVal;
}
主流瓦片地圖的投影彙總
主流瓦片地圖的地理坐标系可分以下幾類:
投影 | ||
1 | Web墨卡托投影 | 谷歌地圖、高德地圖、騰訊地圖、搜狗地圖、必應地圖、OpenStreetMap、ArcGIS online、MapQuest、超圖、360地圖 |
2 | 經緯度直投 | 谷歌地球、天地圖(經緯直投) |
3 | 百度墨卡托 | 百度地圖 |