地理信息系统之瓦片坐标系

最近抓取了部分百度地图的数据，中间的数据使用的是瓦片经纬度。由于先前对这方面知识没有接触过。今天抽时间整理下，供后续使用。

经纬度与坐标系

地球是一个椭球，Datum是一组用于描述这个椭球的数据集合。最常用的一个Datum是WGS84（World Geodetic System 1984），它的主要参数有：

• 坐标系的原点是地球质心（center of mass）；
• 子午线（meridian），即零度经线，位于格林威治子午线Royal Observatory所在纬度往东5米所对应的的经线圈；
• 椭球截面长轴为a=6378137米；
• 椭圆截面短轴为b=6356752.3142米，可选参数；
• 扁平比例（flattening）f=(a−b)/a=1/298.257223563；
• geoid，海平面，用于定义高度，本文从略。

通过以上参数设定，我们才能对地球上的任意一个位置用经度、纬度、高度三个变量进行描述。所以当我们获取一组经纬度信息时，首先要弄明白这组信息对应的Datum。

WGS84 Datum的信息可以用下图进行概括：

WGS84的坐标转化为GCJ02的坐标是单向的，即WGS84的坐标能够准确地变换为GCJ02坐标；但GCJ02坐标转换为WGS84时会存在精度损失。GCJ02的坐标和BD09的坐标转换是双向的。详细信息见先前的文章地图经纬度及坐标系统转换的那点事。

墨卡托投影

地图是显示在平面上的，因此需要将球面坐标转换为平面坐标，这个转换过程称为投影。最常见的投影是墨卡托（Mercator）投影，它具有等角性质，即球体上的两点之间的角度方位与平面上的两点之间的角度方位保持不变，因此特别适合用于导航。墨卡托于1569年提出的一种地球投影方法，该方法是圆柱投影的一种，又名”等角正轴圆柱投影”。

主要特点：

• 没有角度变形，由每一点向各方向的长度比相等，它的经纬线都是平行直线，且相交成直角，经线间隔相等，纬线间隔从基准纬线处向两极逐渐增大。
• 长度和面积变形明显，但基准纬线处无变形，从基准纬线处向两极变形逐渐增大，但因为它具有各个方向均等扩大的特性，保持了方向和相互位置关系的正确。

在地图上保持方向和角度的正确是墨卡托投影的优点，墨卡托投影地图常用作航海图和航空图，如果循着墨卡托投影图上两点间的直线航行，方向不变可以一直到达目的地，因此它对船舰在航行中定位、确定航向都具有有利条件，给航海者带来很大方便。故广泛用于编制航海图和航空图等

值得注意的是：墨卡托投影并不是一种坐标系，而是为了在二维平面上展示三维地球而进行的一种空间映射。所以在GIS地图和互联网地图中，虽然用户看到的地图经过了墨卡托投影，但依然使用经纬度坐标来表示地球上点的位置。

Web墨卡托投影（又称球体墨卡托投影）是墨卡托投影的变种，它接收的输入是Datum为WGS84的经纬度，但在投影时不再把地球当做椭球而当做半径为6378137米的标准球体，以简化计算。其计算公式推导请参考下图：

参数：

• 球半径：以WGS84椭球的长半轴半径为球半径，即球半径取6378137米。(部门用的Flex地图端测距是6378000米，我认为是不太准确的)
• 赤道周长:2PIr = 2*20037508.3427892
• 投影坐标系：以赤道为标准纬线，本初子午线为中央经线，两者交点为坐标原点，向东向北为正，向西向南为负。
• x：[-20037508.3427892,20037508.3427892]
• y：[-20037508.3427892,20037508.3427892]
• lon：[-180,180]
• lat：[-85.05112877980659，05112877980659]

Web墨卡托投影有两个相关的投影标准，经常搞混：

• EPSG4326：Web墨卡托投影后的平面地图，但仍然使用WGS84的经度、纬度表示坐标；
• EPSG3857：Web墨卡托投影后的平面地图，坐标单位为米。

瓦片

经过Web墨卡托投影后，地图就变为平面的一张地图。考虑到有时候我们需要看宏观的地图信息（如世界地图里每个国家的国界），有时候又要看很微观的地图信息（如导航时道路的路况信息）。为此，我们对这张地图进行等级切分。在最高级（zoom=0），需要的信息最少，只需保留最重要的宏观信息，因此用一张256×256像素的图片表示即可；在下一级（zoom=1），信息量变多，用一张512×512像素的图片表示；以此类推，级别越低的像素越高，下一级的像素是当前级的4倍。这样从最高层级往下到最低层级就形成了一个金字塔坐标体系。

对每张图片，我们将其切分为256×256的图片，称为瓦片（Tile）。这样，在最高级（zoom=0）时，只有一个瓦片；在下一级（zoom=1）时有4个瓦片；在下一级（zoom=2）时有16个瓦片，以此类推。上述过程可以用下图进行总结：

瓦片编码

瓦片生成后，就是一堆图片。怎么对这堆图片进行编号，是目前主流互联网地图商分歧最大的地方。总结起来分为四个流派：

• 谷歌XYZ：Z表示缩放层级，Z=zoom；XY的原点在左上角，X从左向右，Y从上向下。
• TMS：开源产品的标准，Z的定义与谷歌相同；XY的原点在左下角，X从左向右，Y从下向上。
• QuadTree：微软Bing地图使用的编码规范，Z的定义与谷歌相同，同一层级的瓦片不用XY两个维度表示，而只用一个整数表示，该整数服从四叉树编码规则
• 百度XYZ：Z从1开始，在最高级就把地图分为四块瓦片；XY的原点在经度为0纬度位0的位置，X从左向右，Y从下向上。

下图显示了前三个流派在zoom=1层级上的瓦片编号结果：

各大地图服务商都采用了Web Mercator进行投影，瓦片坐标系的不同主要是投影截取的地球范围不同、瓦片坐标起点不同。下表总结了中国主要地图商的瓦片编号流派：

高德地图瓦片坐标

坐标系定义

高德地图瓦片坐标与Google Map、Open Street Map相同。高德地图的墨卡托投影截取了约85.05ºS, 约85.05ºN纬度之间部分的地球，使得投影后的平面地图水平方向和垂直方向长度相等。将墨卡托投影地图的左上角作为瓦片坐标系起点，往左方向为X轴，X轴与北纬85.05º重合且方向向左；往下方向为Y轴，Y轴与东经180º（亦为西经180º）重合且方向向下。瓦片坐标最小等级为0级，此时平面地图是一个像素为256*256的瓦片。在某一瓦片层级Level下，瓦片坐标的X轴和Y轴各有$2^{Level}$个瓦片编号，瓦片地图上的瓦片总数为$2^{Level}\times2^{Level}$。

高德地图Level=2的瓦片坐标编号情况：

如上图所示，此时X方向和Y方向各有4个瓦片编号，总瓦片数为16。中国大概位于高德瓦片坐标的（3，1）中。

坐标转换图解

从高德地图坐标转换图解中可以看出，高德地图的坐标转换具有以下特点：

• 所有坐标转换都在某一瓦片等级下进行，不同瓦片等级下的转换结果不同。
• 经纬度坐标可以直接转换为瓦片坐标和瓦片像素坐标。
• 瓦片像素坐标需要结合其瓦片坐标才能得到该像素坐标的经纬度坐标。

坐标转换公式

经纬度坐标（lng, lat）转瓦片坐标（tileX， tileY）:

$$tileX=\lfloor\frac{lng + 180}{360}\times{2^{Level}}\rfloor$$

$$tileY=\lfloor{(\frac{1}{2}-\frac{\ln(\tan(lat\times\pi/180)+\sec(lat\times\pi/180))}{2\timesπ}})\times{2^{Level}}\rfloor$$

经纬度坐标（lng, lat）转像素坐标（pixelX, pixelY）

$$pixelX=\lfloor\frac{lng + 180}{360}\times{2^{Level}}\times256\%256\rfloor$$

$$pixelY=\lfloor{(1-\frac{\ln(\tan(lat\times\pi/180)+\sec(lat\times\pi/180))}{2\timesπ}})\times{2^{Level}}\times256\%256\rfloor$$

瓦片（tileX， tileY）的像素坐标（pixelX, pixelY）转经纬度坐标（lng, lat）

$$lng=\frac{tileX+\frac{pixelX}{256}}{2^{Level}}\times360-180$$

$$lat=\arctan({\sinh({\pi-2\times\pi\times\frac{tileY+\frac{pixelY}{256}}{2^{Level}}})})\times\frac{180}{\pi}$$

Python代码实现

import math


def num2deg(xtile, ytile, zoom):
    """
    Tile numbers to lon./lat.
    """
    n = 2.0 ** zoom
    lon_deg = xtile / n * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    lat_deg = math.degrees(lat_rad)
    return (lat_deg, lon_deg)


def num2deg(xtile, ytile, zoom):
    """
    Tile numbers to lon./lat.
    """
    n = 2.0 ** zoom
    lon_deg = xtile / n * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    lat_deg = math.degrees(lat_rad)
    return (lat_deg, lon_deg)

参考链接：https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames

百度地图瓦片坐标

坐标系定义

百度地图的瓦片坐标系定义与高德地图并不相同，其墨卡托投影的参数也不同。百度地图瓦片坐标以墨卡托投影地图中赤道与0º经线相交位置为原点，沿着赤道往左方向为X轴，沿着0º经线向上方向为Y轴。百度瓦片坐标定义了另一种二维坐标系，称为百度平面坐标系。百度平面坐标系的坐标原点与百度瓦片坐标原点相同，以瓦片等级18级为基准，规定18级时百度平面坐标的一个单位等于屏幕上的一个像素。平面坐标与地图所展示的级别没有关系，也就是说在1级和18级下，同一个经纬度坐标的百度平面坐标都是一致的。

百度地图Level=2的瓦片坐标编号情况：

此时X方向和Y方向各有4个瓦片编号，但是外围的某些瓦片只有部分区域有地图或完全没有地图。没有地图的区域也可以认为其瓦片是无效的，即百度地图中X方向或Y方向的有效瓦片不一定达到$2^{Level}$个。
中国大概位于百度瓦片坐标的(0, 0)中。

坐标转换图解

从百度地图坐标转换图解中可以看出，百度地图的坐标转换具有以下特点：

• 百度经纬度坐标与百度平面坐标可以直接相互转换，并且与瓦片地图等级无关。
• 经纬度坐标需要先转换为平面坐标，然后才能在某一瓦片等级下转换为瓦片坐标和瓦片像素坐标。
• 瓦片像素坐标需要结合其瓦片坐标才能得到该像素坐标的平面坐标，然后再转换为经纬度坐标。

坐标转换公式

平面坐标（pointX, pointY）转瓦片坐标（tileX， tileY）

$$tileX=\lfloor\frac{pointX\times2^{Level-18}}{256}\rfloor$$

$$tileY=\lfloor\frac{pointY\times2^{Level-18}}{256}\rfloor$$

平面坐标（pointX, pointY）转像素坐标（pixelX, pixelY）

$$pixelX=\lfloor{pointX\times2^{Level-18}-\lfloor\frac{pointX\times2^{Level-18}}{256}\rfloor\times256}\rfloor$$

$$pixelY=\lfloor{pointY\times2^{Level-18}-{\lfloor\frac{pointY\times2^{Level-18}}{256}\rfloor\times256}}\rfloor$$

瓦片（tileX， tileY）的像素坐标（pixelX, pixelY）转平面坐标（pointX, pointY）

$$pointX=\frac{tileX\times256+pixelX}{2^{Level-18}}$$

$$pointY=\frac{tileY\times256+pixelY}{2^{Level-18}}$$

代码实现

百度经纬度坐标与百度平面坐标的相互转换，并没有公开的公式，需要通过百度地图的API实现。

/* 经纬度坐标(lng， lat)转平面坐标（pointX, pointY）*/
lnglatToPoint(longitude, latitude)
{
    let projection = new BMap.MercatorProjection();
    let lnglat = new BMap.Point(longitude, latitude);
    let point = projection.lngLatToPoint(lnglat);
    return {
        pointX: point.x, 
        pointY: point.y
    };
}

/* 平面坐标（pointX, pointY）转经纬度坐标(lng， lat) */
pointToLnglat(pointX, pointY)
{
    let projection = new BMap.MercatorProjection();
    let point = new BMap.Pixel(pointX, pointY);
    let lnglat = projection.pointToLngLat(point);
    return {
        lng: lnglat.lng,
        lat: lnglat.lat
    };
}

更多实现：https://github.com/CntChen/tile-lnglat-transform

参考链接：

• 瓦片地图原理
• https://www.maptiler.com/google-maps-coordinates-tile-bounds-projection/
• https://en.wikipedia.org/wiki/Mercator_projection
• Bing Maps Tile System
• 国内主要地图瓦片坐标系定义及计算原理
• https://github.com/wandergis/coordtransform
• https://en.wikipedia.org/wiki/Web_Mercator_projection
• https://lbs.amap.com/api/javascript-api/reference/map/

坐标拾取：

• 百度：http://api.map.baidu.com/lbsapi/getpoint/index.html
• 腾讯：https://lbs.qq.com/tool/getpoint/index.html
• 高德：https://lbs.amap.com/console/show/picker