瓦片地图原理

经纬度

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

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

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

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

投影

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

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

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

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

瓦片

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

对每张图片，我们将其切分为256x256的图片，称为瓦片（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从下向上。

关于中国的经纬度

国内的经纬度有三套系统：

WGS84：为一种大地坐标系，也是目前广泛使用的GPS全球卫星定位系统使用的坐标系。
GCJ02：又称火星坐标系，是由中国国家测绘局制订的地理信息系统的坐标系统。由WGS84坐标系经加密后的坐标系。
BD09：为百度坐标系，在GCJ02坐标系基础上再次加密。其中bd09ll表示百度经纬度坐标，bd09mc表示百度墨卡托米制坐标。

使用OpenStreetMap的坐标为WGS84；使用高德地图、腾讯地图的坐标为GCJ02；使用百度地图的坐标为BD09；谷歌地图和Bing地图的中国部分采用了高德地图的数据，所以坐标为GCJ02。

WGS84的坐标转化为GCJ02的坐标是单向的，即WGS84的坐标能够准确地变换为GCJ02坐标；但GCJ02坐标转换为WGS84时会存在精度损失。

GCJ02的坐标和BD09的坐标转换是双向的，转换规则可以参考下面的python代码：

import math

x_pi = 3.14159265358979324 * 3000.0 / 180.0

def amapcoor2bmapcoor(amap_lon, amap_lat):
    x = amap_lon
    y = amap_lat
    z = math.sqrt(x * x + y * y) + 0.00002 * math.sin(y * x_pi)
    theta = math.atan2(y, x) + 0.000003 * math.cos(x * x_pi)
    bmap_lon = z * math.cos(theta) + 0.0065
    bmap_lat = z * math.sin(theta) + 0.006
    return (bmap_lon, bmap_lat)

def bmapcoor2amapcoor(bmap_lon, bmap_lat):
    x = bmap_lon - 0.0065
    y = bmap_lat - 0.006;
    z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * x_pi);
    theta = math.atan2(y, x) - 0.000003 * math.cos(x * x_pi);
    amap_lon = z * math.cos(theta);
    amap_lat = z * math.sin(theta);
    return (amap_lon, amap_lat)