也说Google卫星地图的URL地址的qrts编码算法

      今天看到有人已经在讨论如何获取google卫星图片，见http://www.cnblogs.com/tangf/archive/2006/07/23/457902.html?login=1&CommentID=1507040#Post的一篇博客，里面抄了JavaScript和delphi的代码过来，我在此再重抄一遍免得大家找来找去。

      算法思想如下：Google卫星图片服务器，由不同层次的256x256大小的jpeg图片无缝拼接而成，其编码方式是按照qrst编码方法进行索引：

       zoom=1时，全球只有一个256x256的图片，它的中心经纬度为(0,0)，其URL为“http://kh.google.com/kh?v=3&t=t”(大家可以把改URL复制到浏览器看下)，其范围是地球按等角纵切圆柱投影后，左右为从西径180度到东径180度，上下范围为从南180度到北180度（这里并不是完全按地球南北只有90度进行划分），中点为赤道与中央子午线的交点，其编码为t。
      当zoom=2时进行第二级编码，即从中点开始上下左右从中分成相等的四等份，从左上开始按顺时针方向分别编为左上q,右上r,右下s,左下t，每一块的编码就为：tq,tr,ts,tt。
      依此类推，每增大一级编码，就放大一倍，每一块都从中分为四块进行下一级编码，编码在前组编码的基础上再分别加上q,r,s,t。即一级编码由一个字母组成，二级编码由两个字母组成，三级由三个字母组成，其它级依次类推，不同地区提供下载的图层级数不尽相同，最多可分到21级。
        对于全球全部卫片的管理来讲，这种编码方法是最好的，是典型的金字塔索引编码方法，采用这种编码要得到某一个图块编号时，看似要对编号进行字串生成计算，其实不然，这种编码方法对于编号是可以用数值加或减就能直接计算，然后做一个数值与字符的转换即可。这对于查找、平移、定位等操作非常快速方便。

JavaScript代码如下：

 1function GetQuadtreeURL_(lon, lat)
 2{
 3    var PI = 3.1415926535897;
 4    var digits = 18; // how many digits precision
 5    // now convert to normalized square coordinates
 6    // use standard equations to map into mercator projection
 7    var x = (180.0 + parseFloat(lon)) / 360.0;
 8    var y = -parseFloat(lat) * PI / 180; // convert to radians
 9    
10    y = 0.5 * Math.log((1+Math.sin(y)) / (1 - Math.sin(y)));
11    
12    y *= 1.0/(2 * PI); // scale factor from radians to normalized
13    y += 0.5; // and make y range from 0 - 1
14    
15    var quad = "t"; // google addresses start with t
16    var lookup = "qrts"; // tl tr bl br
17    
18    while (digits–)    {
19        // make sure we only look at fractional part
20        x -= Math.floor(x);
21        y -= Math.floor(y);
22        quad = quad + lookup.substr((x >= 0.5 ? 1 : 0) + (y >= 0.5 ? 2 : 0), 1);
23        // now descend into that square
24        x *= 2;
25        y *= 2;
26    }
27    return quad;
28}

Delphi代码如下：

 1function getSatURL(zoom: integer; X, Y: double): string;
 2var
 3  wx, wy, cx, cy: double;
 4  tid: string;
 5  i: integer;
 6begin
 7  cx := 0;
 8  cy := 0;
 9  wx := 180;
10  wy := 180;
11  tid := 't';
12
13  for i := 1 to zoom-1 do
14  begin
15    if (x >= cx) and (y >= cy) then
16    begin
17      tid := tid + 'r';
18      cx := cx + wx / 2;
19      cy := cy + wy / 2;
20    end
21    else if (x >= cx) and (y < cy) then
22    begin
23      tid := tid + 's';
24      cx := cx + wx / 2;
25      cy := cy - wy / 2;
26    end
27    else if (x < cx) and (y < cy) then
28    begin
29      tid := tid + 't';
30      cx := cx - wx / 2;
31      cy := cy - wy / 2;
32    end
33    else
34    begin
35      tid := tid + 'q';
36      cx := cx - wx / 2;
37      cy := cy + wy / 2;
38    end;
39    wx := wx / 2;
40    wy := wy / 2;
41  end;
42  result := 'http://kh.google.com/kh?v=2&t=' + tid;
43end;

C#代码如下：

 1        /**//// <summary>
 2        /// 得到Google earth地图图片编码的qrst算法
 3         /// </summary>
 4        /// <param name="row">row:lon</param>
 5        /// <param name="col">col:lat</param>
 6        /// <param name="zoom">zoon:地图级别</param>
 7        /// <returns></returns>
 8        public string getSatTileId(int row, int col, int zoom)
 9        {
10            //trtqst
11            string tileid = "t";
12            double halflat = row;
13            double locxmin, locxmax, locymin, locymax, locxmoy, locymoy;
14
15            locxmin = 0; locxmax = Math.Pow(2, zoom);
16            locymin = 0; locymax = Math.Pow(2, zoom);
17
18            for (int i = 0; i < zoom; i++)
19            {
20                locxmoy = (locxmax + locxmin) / 2;
21                locymoy = (locymax + locymin) / 2;
22                if ((halflat < locymin) ||
23                    (halflat > locymax) ||
24                    (col < locxmin) ||
25                    (col > locxmax))
26                {
27                    return ("transparent");
28                }
29                if (halflat < locymoy)
30                {
31                    locymax = locymoy;
32                    if (col < locxmoy)
33                    { //q quadrant (top left)
34                        tileid += "q";
35                        locxmax = locxmoy;
36                    }
37                    else
38                    {//r quadrant (top right)
39                        tileid += "r";
40                        locxmin = locxmoy;
41                    }
42                }
43                else
44                {
45                    locymin = locymoy;
46                    if (col < locxmoy)
47                    { //t quadrant (bottom right)
48                        tileid += "t";
49                        locxmax = locxmoy;
50                    }
51                    else
52                    { //s quadrant (bottom left)
53                        tileid += "s";
54                        locxmin = locxmoy;
55                    }
56                }
57            }
58            return tileid;
59        }

VB.net代码如下：

 1    /**/''' <summary>
 2    ''' google earth地图图片地址的qrst算法
 3    ''' </summary>
 4    ''' <param name="row">row:lon</param>
 5    ''' <param name="col">col:lat</param>
 6    ''' <param name="zoom">zoon:地图级别</param>
 7    ''' <returns></returns>
 8    Public Function getSatTileId()Function getSatTileId(ByVal row As Integer, ByVal col As Integer, ByVal zoom As Integer) As String
 9        'trtqst
10        Dim tileid As String = "t"
11        Dim halflat As Double = row
12        '
13        Dim locxmin As Double, locxmax As Double, locymin As Double, locymax As Double, locxmoy As Double, locymoy As Double
14
15        locxmin = 0
16        locxmax = Math.Pow(2, zoom)
17        locymin = 0
18        locymax = Math.Pow(2, zoom)
19
20        For i As Integer = 0 To zoom - 1
21            locxmoy = (locxmax + locxmin) / 2
22            locymoy = (locymax + locymin) / 2
23            If (halflat < locymin) OrElse (halflat > locymax) OrElse (col < locxmin) OrElse (col > locxmax) Then
24                Return ("transparent")
25            End If
26            If halflat < locymoy Then
27                locymax = locymoy
28                If col < locxmoy Then
29                    'q quadrant (top left)
30                    tileid += "q"
31                    locxmax = locxmoy
32                Else
33                    'r quadrant (top right)
34                    tileid += "r"
35                    locxmin = locxmoy
36                End If
37            Else
38                locymin = locymoy
39                If col < locxmoy Then
40                    't quadrant (bottom right)
41                    tileid += "t"
42                    locxmax = locxmoy
43                Else
44                    's quadrant (bottom left)
45                    tileid += "s"
46                    locxmin = locxmoy
47                End If
48            End If
49        Next
50        Return tileid
51    End Function

注：每块图片的URL格式为“http://kh.google.com/kh?v=3&t=tsrrtsqrssqqttsstrst”

反向递推算法如下（由tsrrtsqrssqqttsstrst编码得到经纬度）：

 1//根据google地图的字符地址得到经纬度的度值，Sexagesimal:60进制
 2function GoogleQRSTtoDegree(x)
 3{
 4    // 23° 27′ 30"
 5    var ret = "";
 6    if (x < 0)
 7    {
 8        ret += '-';
 9        x = -x;
10    }
11    ret += Math.floor(x);
12    ret += '°';
13    
14    x = (x - Math.floor(x)) * 60;
15    ret += Math.floor(x);
16    ret += "'";
17    
18    x = (x - Math.floor(x)) * 60;
19    ret += Math.floor(x);
20    ret += '"';
21    
22    return ret;
23}