今天看到有人已经在讨论如何获取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代码如下:
1
function
GetQuadtreeURL_(lon, lat)
2
![ExpandedBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
{
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![ExpandedSubBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
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代码如下:
1
function
getSatURL(zoom: integer; X, Y: double):
string
;
2
var
3
wx, wy, cx, cy: double;
4
tid:
string
;
5
i: integer;
6
begin
7
cx :
=
0
;
8
cy :
=
0
;
9
wx :
=
180
;
10
wy :
=
180
;
11
tid :
=
'
t
'
;
12
![](https://www.cnblogs.com/Images/OutliningIndicators/None.gif)
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;
43
end
;
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
![ExpandedBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
{
10
//trtqst
11
string tileid = "t";
12
double halflat = row;
13
double locxmin, locxmax, locymin, locymax, locxmoy, locymoy;
14![](https://www.cnblogs.com/Images/OutliningIndicators/InBlock.gif)
15
locxmin = 0; locxmax = Math.Pow(2, zoom);
16
locymin = 0; locymax = Math.Pow(2, zoom);
17![](https://www.cnblogs.com/Images/OutliningIndicators/InBlock.gif)
18
for (int i = 0; i < zoom; i++)
19![ExpandedSubBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
20
locxmoy = (locxmax + locxmin) / 2;
21
locymoy = (locymax + locymin) / 2;
22
if ((halflat < locymin) ||
23
(halflat > locymax) ||
24
(col < locxmin) ||
25
(col > locxmax))
26![ExpandedSubBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
27
return ("transparent");
28
}
29
if (halflat < locymoy)
30![ExpandedSubBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
31
locymax = locymoy;
32
if (col < locxmoy)
33![ExpandedSubBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{ //q quadrant (top left)
34
tileid += "q";
35
locxmax = locxmoy;
36
}
37
else
38![ExpandedSubBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{//r quadrant (top right)
39
tileid += "r";
40
locxmin = locxmoy;
41
}
42
}
43
else
44![ExpandedSubBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
45
locymin = locymoy;
46
if (col < locxmoy)
47![ExpandedSubBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{ //t quadrant (bottom right)
48
tileid += "t";
49
locxmax = locxmoy;
50
}
51
else
52![ExpandedSubBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{ //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![](https://www.cnblogs.com/Images/OutliningIndicators/InBlock.gif)
15
locxmin = 0
16
locxmax = Math.Pow(2, zoom)
17
locymin = 0
18
locymax = Math.Pow(2, zoom)
19![](https://www.cnblogs.com/Images/OutliningIndicators/InBlock.gif)
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进制
2
function
GoogleQRSTtoDegree(x)
3
![ExpandedBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
{
4
// 23° 27′ 30"
5
var ret = "";
6
if (x < 0)
7![ExpandedSubBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
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
}