简明版墨卡托投影坐标系 (原理到实现)

[toc]

前言

FYI 本博客初稿完成于 2017 年,内容更新于 个人网站 - 简明版墨卡托投影坐标系,请移步阅读最新内容。

本文讲述了墨卡托投影坐标系的基本原理和实现,但是因为地球非标准椭圆,经纬度和米坐标的转换复杂,本文提供的算法存在较大误差,仅适用于初步验证。

墨卡托投影坐标系

墨卡托投影 (Mercator Projection) 是一种 “等角正切圆柱投影”,荷兰地图学家墨卡托 (Mercator) 在 1569 年拟定:假设地球被围在一个中空的圆柱里,其赤道与圆柱相接触,然后再假想地球中心有一盏灯,把球面上的图形投影到圆柱体上,再把圆柱体展开,这就是一幅标准纬线为零度(即赤道)的 “墨卡托投影” 绘制出的世界地图。
墨卡托投影在今天对于航海事业起着极为重要的作用,目前世界各国绘制海洋地图时仍广泛使用墨卡托投影,国际水路局 (IHB) 规定:“除特殊情况外,各国都要用墨卡托投影绘制海图”。国际水路局发行的《大洋水深总图》是把全世界分成 24 幅编辑的,在南北纬 72 度之间就是使用墨卡托投影绘成的。

墨卡托投影性质

由于墨卡托投影的经纬线离开赤道逐渐以相同倍数伸长,所以又称为渐长投影,由于它是具有等角性质的圆筒投影,所以也叫做等角圆筒投影。注意:这种投影不适合高纬地区,通常纬度 60 度以上区域,不用此投影

墨卡托投影有一个特别的特性:所有罗盘等角线,或称斜航线(就是与所经过的所有经线形成相同角度的航线,也称恒向航线)在墨卡托投影下都是直线。这使得在航海领域这个投影非常重要。
注意:经纬线的伸长与纬线的正割成比例变化,随纬度增高极具拉伸,到极点成为无穷大;面积的扩大更为明显,在 60 度的地方面积要扩大四倍。如下图所示,地理上等半径圆在高纬度面积明显扩大。

墨卡托投影是按等角条件修改透视圆筒投影而得到的投影,等角(也称为保形) 是指当地图上任何一点的各方向具有相同的比例,称为局部保形,透视圆筒投影如图 1 所示。从墨卡托投影图上可以看出,经线间隔的经度如果相等,则经线是等距平行的直线, 纬线也是平行的直线,而且经纬线是相互垂直的。墨卡托投影对透视圆筒投影改造点:要使圆筒投影称为等角的性质,必须使由赤道向两极经线逐渐伸长的倍数与经线上各点相应的纬度扩大的倍数相同。

透视圆筒投影

墨卡托投影方程式

墨卡托投影以整个世界范围,赤道作为标准纬线,本初子午线作为中央经线,两者交点为坐标原点,向东向北为正,向西向南为负。南北极在地图的正下、上方,而东西方向处于地图的正右、左。由于墨卡托投影在两极附近是趋于无限值,因此它并没完整展现了整个世界,地图上最高纬度是 85.05 度(通过纬度取值范围 ys 反解计算可得到纬度值为 85.05112877980659)。为了简化计算,我们采用球形映射,而不是椭球体形状。

公式推导具体见文献 墨卡托投影与大圆投影的构成及其在定航线计算航程与航向方面的应用_程光举


利用等角条件 m=n 来讨论具体公式,具体分为三步:

  1. 根据 m=n 得到地球表面投影到平面上的微积线段的关系式。
  2. 把地球视为球体:
    设地球表面 A 点经纬坐标为(λ,Φ),对应的投影坐标为(x,y), 基准纬线设置为赤道,则 R 为地球半径;

墨卡托投影方程式为:

  1. 把地球视为旋转椭球体

墨卡托投影正反解公式:

公式推导具体见文献 墨卡托投影与大圆投影的构成及其在定航线计算航程与航向方面的应用_程光举

程序实现

最新版见 Geography Coordinate Transform Lite/Mercator

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
// Lite版,仅适用于初步验证
//把地球视为球体实现经纬度和墨卡托投影的函数
typedef struct Point
{
double x;
double y;
}WayPoint;

//经纬度转墨卡托
WayPoint lonLat2Mercator(WayPoint lonLat)
{
WayPoint mercator;
double x = lonLat.x * 20037508.34 / 180;
double y = log(tan((90 + lonLat.y) * Pi / 360))/(Pi / 180);
y = y * 20037508.34 / 180;
mercator.x = x;
mercator.y = y;
return mercator;
}

//墨卡托转经纬度
WayPoint Mercator2lonLat(WayPoint mercator)
{
WayPoint lonLat;
double x = mercator.x / 20037508.34 * 180;
double y = mercator.y / 20037508.34 * 180;
y = 180 / Pi * (2 * atan(exp(y * Pi / 180)) - Pi / 2);
lonLat.x = x;
lonLat.y = y;
return lonLat;
}

Reference