【面】空间索引技术之Google S2

S2算法采用正方体投影的方式将地球展开,然后利用希尔伯特分形曲线将展开后的二维地球进行填充,完成了对三位地球的降维和分形,从而得到空间坐标点与希尔伯特分形曲线的函数关系,即将球面经纬度坐标转换成球面 xyz 坐标,再转换成正方体投影面上的坐标,最后变换成修正后的坐标在坐标系变换,映射到 [0,2^30^-1]区间,最后一步就是把坐标系上的点都映射到希尔伯特曲线上。最终,映射到希尔伯特曲线上的点成为 Cell ID,即是空间坐标点的索引。

S2的最大的优势在于精度高。Geohash有12级,从5000km到3.7cm,中间每一级的变化比较大。有时候可能选择上一级大了一些,选择下一级又小了一些。而S2有30级,从0.7cm^2^到 85,000,000km^2^,中间每一级的变化都比较平缓,接近于4次方的曲线。所以选择精度时不会出现Geohash选择困难的问题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
class S2:
def __init__(self):
self.hilbert_map = {
'a': {(0, 0): (0, 'd'), (0, 1): (1, 'a'), (1, 0): (3, 'b'), (1, 1): (2, 'a')},
'b': {(0, 0): (2, 'b'), (0, 1): (1, 'b'), (1, 0): (3, 'a'), (1, 1): (0, 'c')},
'c': {(0, 0): (2, 'c'), (0, 1): (3, 'd'), (1, 0): (1, 'c'), (1, 1): (0, 'b')},
'd': {(0, 0): (0, 'a'), (0, 1): (3, 'c'), (1, 0): (1, 'd'), (1, 1): (2, 'd')},
}
self.un_hilbert_map = {
'a': { 0: (0, 0,'d'), 1: (0, 1,'a'), 3: (1, 0,'b'), 2: (1, 1,'a')},
'b': { 2: (0, 0,'b'), 1: (0, 1,'b'), 3: (1, 0,'a'), 0: (1, 1,'c')},
'c': { 2: (0, 0,'c'), 3: (0, 1,'d'), 1: (1, 0,'c'), 0: (1, 1,'b')},
'd': { 0: (0, 0,'a'), 3: (0, 1,'c'), 1: (1, 0,'d'), 2: (1, 1,'d')}
}
self.lon_range = [-180.0, 180.0]
self.lat_range = [-90.0, 90.0]

def encode(self, lon, lat, order=16):
'''
S2编码(cartesian_to_hilbert)
'''
current_square = 'a'
position = 0

for i in range(order - 1, -1, -1):
position <<= 2

lng_mid = (self.lon_range[0] + self.lon_range[1])/2
lat_mid = (self.lat_range[0] + self.lat_range[1])/2

if lon >= lng_mid :
quad_x = 1
self.lon_range[0] = lng_mid
else:
quad_x = 0
self.lon_range[1] = lng_mid

if lat >= lat_mid :
quad_y = 1
self.lat_range[0] = lat_mid
else:
quad_y = 0
self.lat_range[1] = lat_mid

quad_position,current_square = self.hilbert_map[current_square][(quad_x, quad_y)]
position |= quad_position

return position

def decode(self, d, order=16):
'''
S2解码(hilbert_to_cartesian)
'''
current_square = 'a'
lon = lat = lon_mid = lat_mid = 0

for i in range(order - 1, -1, -1):
lon_mid = (self.lon_range[0] + self.lon_range[1] ) / 2
lat_mid = (self.lat_range[0] + self.lat_range[1] ) / 2

mask = 3 << (2*i)
quad_position = (d & mask) >> (2*i)

quad_x, quad_y, current_square= self.un_hilbert_map[current_square][quad_position]

if quad_x:
self.lon_range[0] = lon_mid
else:
self.lon_range[1] = lon_mid

if quad_y:
self.lat_range[0] = lat_mid
else:
self.lat_range[1] = lat_mid

lon = self.lon_range[0]
lat = self.lat_range[0]

return lon, lat

if __name__ == '__main__':
s = S2()
d = s.encode(-50.555443, 77.776655,36)
print(d)
lon, lat = s.decode(d, 36)
print(lon, lat)
0%