球面墨卡托(Spherical Mercator)

地理信息描述空间位置相关的信息，在空间位置的表达中，需要基于空间参照系来保证数据精度以及不同数据源之间的相互叠加/空间分析操作。自Google Maps与2005年发布以来，电子地图服务与普通民众的日常生活关联日益密切，对于从事GIS行业的人员来说，基于现有地图服务制作相关应用能够达到事半功倍的效果。在本篇文章中主要介绍电子地图数据相关的概念以及构成弟子地图服务基础的Web Mercator投影坐标系的相关内容。

相关概念：

1. Spacial Reference，空间参考
2. EPSG，European Petroleum Survey Group (EPSG)，该组织负责维护并发布坐标参照系统的数据集参数，以及坐标转换描述
3. GCS，地理坐标系，椭球体定义，如WGS84坐标系，世界地理坐标系，EPSG 4326，GPS使用该坐标系，以及国内的54，80坐标系，使用经纬度坐标
4. PCS，投影坐标系，基于GCS椭球体定义的投影方式的定义
5. Mercator投影，由荷兰人Mercator定义的正轴等角切圆柱投影。Google Maps, Yahoo Maps, Microsoft Maps等网络电子地图服务使用基于该投影坐标系的Web Mercator坐标系， EPSG:900913 或 EPSG:3785，使用像素坐标（单位：像素）或者空间直角坐标（单位：米）
6. Groud Resolution，地面精度为地图上的像素表达的地面距离，单位为meters per pixcel
7. Map Scale，地图比例尺为地图距离与地面距离的比值
8. Tile，瓦片表示网络地图中的最小组成部分
9. TMS，Tile Map Service，瓦片地图服务，使用层行列坐标或者Quadkey

坐标转换

 在使用公开的电子地图服务构建应用并叠加数据时，需要用到以下的坐标转换场景：

         LatLon      <->       Meters      <->     Pixels    <->       Tile     

     WGS84 coordinates   Spherical Mercator  Pixels in pyramid  Tiles in pyramid
         lat/lon            XY in metres     XY pixels Z zoom      XYZ from TMS 
        EPSG:4326           EPSG:900913                                         
         .----.              ---------               --                TMS      
        /           <->     |       |     <->     /----/    <->      Google    
              /             |       |           /--------/          QuadTree   
         -----               ---------         /------------/                   
       KML, public         WebMapService         Web Clients      TileMapService

EPSG:900913坐标范围

[-20037508, -20037508, 20037508, 20037508]

let r = 6378137
// 常量20037508根据地球周长（单位米）计算
// 坐标原点位于范围中心
offset = 2 * pi * r  / 2.0

EPSG:900913金字塔缩放级别地面精度(pixels/meter)

// 整个区域在金字塔顶层 (zoom=0)为一张256x256 pixels的 瓦片,往下逐层细分
baseMetersPerPixel = offset * 2 / 256 = 156543

LatLon 与 Meter X/Y互相转换

def YToLat(y)
{
    return pi / 2.0 - 2.0 * atan(exp(-y/r))
} 
def LatToY(lat) 
{ 
    sinLat = sin(lat);
    return 0.5*log((1.0 + sinLat)/(1.0 - sinLat))*r; 
    // return -log(tan(pi / 4.0 - lat / 2.0)) * r; 
}

def LonToX(lon)
{
    return lon * offset / pi
}

def XToLon(x)
{
    return x  * pi / offset
}

其中，LatToY有两种写法，转换原理为：

let t = tan(pi / 4.0 - lat / 2.0)
= tan((pi/2.0–lat) / 2.0)
// 由半角公式
= sin(pi/2.0–lat) / (1+cos(pi/2.0–lat))
// 由sin（π/2-α）= cosα, cos（π/2-α）= sinα
= cos(lat)/(1+sin(lat))
= sqrt((1-sin(lat))/(1+sin(lat)))

Log(t) 
= log(sqrt((1-sin(lat))/(1+sin(lat))))
= -0.5*log((1+sin(lat))/ (1-sin(lat)))

Meter X/Y 与 Pixel x/y 互相转换

此处示范的是Bing Map中的转换方法，设定为：(origin [0,0] in top-left)，如下图所示（Level 3）

def MetersPerPixel2(zoom)
{
    return baseMetersPerPixel / (1 << zoom)
}

def PixelYToY(pixelY,zoom)
{
    return offset - pixelY * MetersPerPixel2(zoom)
}

def YToPixelY(y,zoom)
{
    return (offset - y) / MetersPerPixel2(zoom)
}

def PixelXToX(pixelX, zoom)
{
    return pixelX * MetersPerPixel2(zoom) - offset
}

def XtoPixelX(x,zoom)
{
    return (x + offset) / MetersPerPixel2(zoom)
}

Pixel x/y 与 Tile l/r/c或者quadkey互相转换

此处示范的是Bing Map中的瓦片编号方法，如下图所示（Level 3）

 

def TileToPixel(tile)
{
    east = (tile.X + 1) * 256
    north = tile.Y * 256
    south = (tile.Y+1) * 256
    west = tile.X * 256
}

def PixelToTile(pixelX, pixelY)
{
    tile.X = floor(pixelX / 256)
    tile.Y = floor(pixelY / 256)
}

通过tile的层行列号，转换到Bingmap下的quadkey的方法详见MSDN.

参考：

Bing Maps Tile System

Mercator那些事儿

Global Map Tiles Classes