坐标参考

最近在用webgl手撕一个三维GIS引擎，引擎内部使用4326的参考系，切片使用的3857的数据源，出现的效果就是经度正常，但是纬度方向被压缩了，如下图所示。

这是一个可以预见的情况，但是也引出了一些长期困扰自己的问题，趁这次机会深入理清一下。

一、Mercator投影

mercator投影是一种”等角正切圆柱投影”，以地球为例，假设地球被围在一中空的圆柱里，其基准纬线与圆柱相切（赤道）接触，然后再假想地球中心有一盏灯，把球面上的图形投影到圆柱体上，再把圆柱体展开，这就是一幅选定基准纬线上的“墨卡托投影”绘制出的地图，如下图：

当然，这是一个对于任何一个giser都知道的概念，在这里重点强调的一点是：墨卡托投影的对象不仅仅是正球，对椭球也一样通用。

更准确的来说，椭球才是标准对象。我们的地球是一种更接近与梨形的球体，如下

每个地区高低不平，所以各个地区都会根据本区域特点，选取一个最适合本区域的椭球体作为参考了对象，往往一个地区的参考椭球体并不使用于另外一个区域。实际上，对于局部小区域，也可用会用圆锥切面投影，或是高斯-克吕格投影。

Mercator 投影坐标系统，全球范围尺度上其基准面可以是 WGS 1984 ，

WGS 1984定义如下：

1
2
3
4
5
6
7
8
9
GCS_WGS_1984
WKID: 4326 Authority: EPSG
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_WGS_1984
Spheroid: WGS_1984
Semimajor Axis: 6378137.0
Semiminor Axis: 6356752.314245179
Inverse Flattening: 298.257223563

墨卡托投影的“圆柱”特性，保证了南北（纬线）和东西（经线）都是平行直线，并且相互垂直。而且经线间隔是相同的，纬线间隔从标准纬线（此处是赤道，也可能是其他纬线）向两级逐渐增大。

既然是投影了，肯定会涉及到投影的换算，它的投影计算公式应该是椭球的投影公司如下：

  • x、y是投影展开成平面后以赤道本初子午线交点为原点的平面坐标系的坐标
  • a 是椭球体长半轴，b是短半轴
  • L是经度（弧度制），B是纬度（弧度制）

二、WebMercator投影

Web Mercator 坐标系使用的投影方法不是严格意义的墨卡托投影，而是一个被 EPSG（European Petroleum Survey Group）称为伪墨卡托的投影方法，这个伪墨卡托投影方法的大名是 Popular Visualization Pseudo Mercator，PVPM。

因为这个坐标系统是 Google Map 最先使用的，或者更确切地说，是Google 最先发明的。在投影过程中，将表示地球的参考椭球体近似的作为正球体处理（正球体半径 R = 椭球体半长轴 a）。这也是为什么在 ArcGIS 中我们经常看到这个坐标系叫 WGS 1984 Web Mercator (Auxiliary Sphere)。Auxiliary Sphere 就是在告知你，这个坐标在投影过程中，将椭球体近似为正球体做投影变换，虽然基准面是WGS 1984 椭球面。

后来，Web Mercator 在 Web 地图领域被广泛使用，这个坐标系就名声大噪。尽管这个坐标系由于精度问题一度不被GIS专业人士接受，但最终 EPSG 还是给了 WKID:3857

三、切片和切片原点

我们知道通过墨卡托投影后，x轴是经度投影的结果，经度区间是[-π，π]；y轴时纬度投影的结果，而纬度区间为[-π/2,π/2]。而我们的3857的切片服务，通常切片原点是[-2.0037508342787E7,2.0037508342787E7]。

这就不好理解了，x轴是以赤道周长来算的，x轴的范围是[-2.0037508342787E7,2.0037508342787E7]可以计算出来，可是为啥y轴的的范围也是[-2.0037508342787E7,2.0037508342787E7]呢？另一个与此相光的问题是为什么在切片的每个级别上分辨率是一个值呢？常规理解应该是x轴一个分辨率，y轴一个分辨率的。并且由于切片是正方形，x轴的分辨率应该是y轴分辨率的2倍才对呀。

其实虽然纬度区间[-π/2,π/2]只是经度区间[-π，π]的二分之一，但是在投影后同一纬度上的经度属于等间距投影，但是纬度却是随着纬度的增加投影距离不断增大，直到纬度为90度时，投影点趋于无穷大。纬度上的距离区间应该远远大于经度区间才对，如下图

此时有两点可以明确：

  • 同一纬度上，经度投影分布是均匀的，线性的
  • 纬度的投影，随着纬度的增加，投影距离是不断增大，是非线性的。

那么如何让纬度纬度线的区间分布也是均匀的呢？这场计算纬度线的投影距离公式是y=f(φ)=R*tan(φ)；有一个聪明的家伙模拟了一个近似函数

1
2
3
		-ln(tan(π/4 + abs(φ/2)))     φ<0时
y=g(φ)= +ln(tan(π/4 + abs(φ/2)))     φ=0时
		+ln(tan(π/4 + abs(φ/2)))     φ>0时

这种投影算法使得赤道附近的纬线较密，极地附近的纬线较稀。极点被投影到无穷远，所以这种投影不适合在高纬度地区使用。Google Maps的选取的范围为 -π<y<π ，这样近似的有 -85°<Φ<85°。这也就得到了y轴上的取值区间和x轴上的取值区间一致，一张方形的世界地图就出来了。

这时就注意了，在做x轴的切片计算时，可以采用二分的方式直接计算，但是对于y轴切片的计算则不能使用二分了，因为y轴不是线性切片的，必须严格按照该界别切片提供的分辨率进行计算。这么做就将问题分离开了，可视化人员只需要按照既定的规则进行切片展示；制图人员可以自定义经纬度两个方向切片的计算方法（只需最终输出满足既定规则的切片即可）；

四、像素坐标

提到分辨率，就需要解释像素坐标了。

分辨率表示一个像素代表的实际长度，可以同过分辨率和像素坐标来计算该像素点所在的实际位置。

将地球表面通过墨卡托投影到一个方形平面时，依据展示内容的精细度，这个方形平面可大可小。通常全球范围内就不需要展示特别精细的内容，只需要轮廓就可以，而在一个工业园区级别就需要将具体的厂房道路展示出来。这样最终精细度的提高，我们这张世界地图的尺寸就越来越大。

不过好在，我们通常在精细化程度较高时，往往只需要看到很小的一块区域，为了切分的方便，一般切分方式都是以2的次幂来切分。这方面的内容此处不再赘述。只强调一点：既定的级别都有自己的分辨率，无论经度还是纬度都可以依据这个分辨率来准确的计算某一像素点的实际位置，而不必再考虑该图是如何将经度和纬度统一起来制图的。

地图学中一个重要的概念，就是比例尺，即地图上的一厘米代表着实际上的多少厘米；到了web地图中我们把比例尺转换成另一个概念，分辨率（Resolution），即图上一像素代表实际多少米。

假设地图的坐标单位是米，整张地图的dpi为96，当前地图在赤道处的比例尺为1:125000000（即图上1米等于实地125000000米），1英寸=2.54厘米； 1英寸=96像素。那么计算可得地图赤道上1像素代表实地距离是 125000000*0.0254/96 = 33072.9166666667米。为什么要强调某一条纬度线（上述例子为赤道）的比例尺？因为根据墨卡托投影的特性，同一张地图中不同纬度线的比例尺是变化的，越靠近两极，图上1米相当于实地的距离越小。

参考：

https://www.jianshu.com/p/778fc3e9f889

https://blog.csdn.net/mr_jianrong/article/details/72625811

https://blog.csdn.net/kikitamoon/article/details/46124935

https://www.maptiler.com/google-maps-coordinates-tile-bounds-projection/

https://blog.csdn.net/qq_35732147/article/details/83856513

https://baike.baidu.com/item/%E5%A2%A8%E5%8D%A1%E6%89%98%E6%8A%95%E5%BD%B1/5477927?fr=kg_qa)

https://www.maptiler.com/google-maps-coordinates-tile-bounds-projection/