Google Map瓦片图算法分析

转自去年看到的一篇文章 http://www./bbs/viewthread.php?tid=321 目前已经不能访问了，另外Tile服务器的地址也换了。更新一下，留作备份。

Google Map 瓦片图算法分析
Map tile编码

　　Google map使用两种算法对tile的位置进行编码。

　　对于Google ditu，tile的url地址类似于：http://mt2.google.cn/mt?v=cn1.6&hl=zh-CN&x=6678&y=3557&z=13&s=Galile  使用x和y来设置tile坐标和放大因子。放大因子从17（完全缩小）到0（最大比例）。当放大因子为17时，整个地球在一个tile中显示，此时x=0 ，y=0；放大因子为16时，地球被分为2x2部分，这时0<=x<=1 且0<=y<=1。每放大一次，每个tile被分为4个部分。因此，当放大因子为z时，显示的水平和垂直tile个数为2^(17-z)。

理论上的算法会像下面这样的：

//correct the latitude to go from 0 (north) to 180 (south),
// instead of 90(north) to -90(south)
latitude=90-latitude;

//correct the longitude to go from 0 to 360
longitude=180+longitude;

//find tile size from zoom level
double latTileSize=180/(pow(2,(17-zoom)));
double longTileSize=360/(pow(2,(17-zoom)));

//find the tile coordinates
int tilex=(int)(longitude/longTileSize);
int tiley=(int)(latitude/latTileSize);

上面是假设投影为平面。

实际上googlemap,51ditu在显示时使用了墨卡托投影，因此上述的算法需要进行修改。
在墨卡托投影中，两条纬线间的距离不一定相等，所以纬度对应的titley会依据它的垂直位置。

以下代码为从tile的纬度位置计算的tile的垂直个数

Collapse
/**<summary>Get the vertical tile number from a latitude using mercator ptrojection formula</summary>*/

/**<summary>此函数取得相应纬度对应的瓦片纵轴序号，参数为纬度</summary>*/

        private int getMercatorLatitude(double lati)
        {
            double maxlat = Math.PI;

            double lat = lati;
          
            if (lat > 90) lat = lat - 180;
            if (lat < -90) lat = lat + 180;

            // conversion degre=>radians
            // 转换度数到弧度
            double phi = Math.PI * lat / 180;

            double res;
              //网上其他帖子这个地方有问题，应该为加号
            //double temp = Math.Tan(Math.PI / 4 + phi / 2);
            //res = Math.Log(temp);
            //下面这一句是上面的合并
            res = 0.5 * Math.Log((1 + Math.Sin(phi)) / (1 - Math.Sin(phi)));
            double maxTileY = Math.Pow(2, zoom);
            int result = (int)(((1 - res / maxlat) / 2) * (maxTileY));

            return (result);
        }

 

　　覆盖地带：
　　理论上，纬度范围应该是-90度到90度，但事实上，由于墨卡托投影使得两级无穷大，覆盖的地带小于-90 到 90。

 

最后加上 Python 验证代码

 

# -*- coding: UTF-8 -*-


import math

tileServer = "http://mt2.google.cn/mt?v=cn1.6&hl=zh-CN"

def getMercatorLatitude(lati, zoom):
    
    maxlat = math.pi
    lat = lati
   
    if (lat > 90):
        lat = lat - 180
    if (lat < -90):
        lat = lat + 180

    # 转换度数到弧度
    phi = math.pi * lat / 180.0;

    #temp = math.tan(math.pi / 4.0 + phi / 2.0);
    #res = math.log(temp);
    #下面这一句是上面的合并
    res = 0.5 * math.log((1 + math.sin(phi)) / (1 - math.sin(phi)))
    maxTileY = math.pow(2, 17-zoom)
    result = (int)(((1 - res / maxlat) / 2) * (maxTileY))

    return result;


def getTile(longitude, latitude, zoom):

    longitude=180+longitude
    
    longTileSize=360.0/(pow(2,(17-zoom)))

    tilex =  longitude/longTileSize

    tiley = getMercatorLatitude(latitude, zoom)

    tilex = int(math.floor(tilex))

    tiley = int(math.floor(tiley))

    #return [tilex, tiley]

    return tileServer + "&x=" + str(tilex) + "&y=" + str(tiley) + "&z=" + str(17-zoom) + "&s=Galile"

print getTile(114.28, 30.55, 5)

#测试结果 http://mt2.google.cn/mt?v=cn1.6&hl=zh-CN&x=3348&y=1682&z=12&s=Galile