Google Map瓦片圖演算法分析

轉自去年看到的一篇文章 http://www.earthplayer.com/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 360longitude=180+longitude;

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

//find the tile coordinatesint 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-*-

importmath

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

defgetMercatorLatitude(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))

returnresult;

defgetTile(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]

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

printgetTile(114.28,30.55,5)

#測試結果http://mt2.google.cn/mt?v=cn1.6&hl=zh-CN&x=3348&y=1682&z=12&s=Galile
推薦閱讀：