geohash算法

来源：互联网发布：cdn网宿和阿里云cdn 编辑：程序博客网时间：2024/06/07 03:03

今天在点饿了吗外卖的时候，跳出了个geohash不能为空，导致下不了单，于是百度了一下！原来是没有定位，选择配送地址！

1.算法背景

Geohash的初衷是如何用尽量短的URL来标志地图上的某个位置，而地图上的位置一般是用经纬度来表示，问题就转化为如何把经纬度转化为一个尽量短的URL。

Geohash的算法描述请参考:http://en.wikipedia.org/wiki/Geohash ,本文的主要目的是更加细致地解释该算法的原理及实用场景。

2.算法

算法的主要思想是对某一数字通过二分法进行无限逼近，比如纬度的区间是[-90,90]，假如给定一个纬度值：46.5，可以通过下面算法对46.5进行无限逼近:

（1）把区间[-90,90]进行二分为[-90,0),[0,90]，称为左右区间，可以确定46.5属于右区间[0,90]

（2）递归上述过程46.5总是属于某个区间，无论第几次迭代46.5总是属于某个区间[a,b]。随着每次迭代区间[a,b]总在缩小，根据极限可知[a,b]会收敛到46.5，用δ来描述就是，任意给定一个 ε，总存在一个N使得： δ=|x-a/2^N |< ε,x为任意给定的纬度

（3）上述分析过程保证了算法收敛性的同时，再记录一下收敛的过程：如果给定的纬度x（46.5）属于左区间，则记录0，如果属于右区间则记录1，这样随着算法的进行会产生一个序列1011100，序列的长度跟给定的收敛次数N相关。

反过来，如果我们知道了序列1011100，我们就可以分别能确定纬度x（46.5）属于哪个更小的迭代区间，也就是说该算法是可逆的

（4）算法的精度：显然的是，不可能让计算机执行无穷计算，加入执行N此计算，则x属于的区间长度为(b-a)/2^N+1 ,以纬度计算为例，则为180/2^N ，误差近似计算为:err= 180/2^N+1 =90/2^N ,如果N=20，则误差最大为：0.00009。但无论如何这样表明Geohash是一种近似算法。

3.编码

在对纬度产生了序列1011100后，在对经度做相同的算法也会产生一个序列，比如0011101。根据偶数位放经度，奇数位放纬度（0被视为偶数），把2个串合二为一，产生一个新串：01001111110010，对该串进行Base32编码，则可获得一个ASIIC码的字符串，关于Base32编码，请参考：http://en.wikipedia.org/wiki/Base32

4.解码

解码的过程相对比较简单

（1）对拿到的字符串进行Base32解码

（2）根据奇偶位取出纬度、经度

（3）根据序列反向得到每个区间，并取中间值（0为左区间，1为右区间）

5.应用

该算法目前主要用在地图的地址搜索，有了该算法可以为数据库中的地址建立索引，极大提高地图数据检索的速度。

仔细观察，该算法还有另为一个特点，对相近的x,y，会得到相同前缀的序列，原因是相近的x,y，在递归的绝大多数时间会处在同一个区间，故而，逼近的轨迹是一致的，这又可以解决地图中“离我最近的搜索“的问题，同时，对进行hash的模糊检索也有一定的启发作用。

6.代码

下面的实现代码很精美，引用自：http://bloggermap.org/rss/readblog/70907，格式不太好，大家自己整理一下即可

[cpp] view plaincopy
#define BASE32 "0123456789bcdefghjkmnpqrstuvwxyz"  
static void encode_geohash(double latitude, double longitude, int precision, char *geohash) {  
int is_even=1, i=0;  
double lat[2], lon[2], mid;  
char bits[] = {16,8,4,2,1};  
int bit=0, ch=0;  
lat[0] = -90.0; lat[1] = 90.0;  
lon[0] = -180.0; lon[1] = 180.0;  
while (i < precision) {  
if (is_even) {  
mid = (lon[0] + lon[1]) / 2;  
if (longitude > mid) {  
ch |= bits[bit];  
lon[0] = mid;  
} else  
lon[1] = mid;  
} else {  
mid = (lat[0] + lat[1]) / 2;  
if (latitude > mid) {  
ch |= bits[bit];  
lat[0] = mid;  
} else  
lat[1] = mid;  
}  
is_even = !is_even;  
if (bit < 4)  
bit++;  
else {  
geohash[i++] = BASE32[ch];  
bit = 0;  
ch = 0;  
}  
}  
geohash[i] = 0;  
}  

------------------------------------------------------------------------------------------------------------------------------------------

我一直在奇怪geohash实现周边查询的原理是什么，毕竟地理数据可是二维的坐标，而geohash的结果只是一个简单的字符串，要说通过简单字符串的比较就能找到周边的点，无论如何我也不能相信，因此我研究了一下geohash的算法，果然发现我以前的做法其实是做不到精确的周边查询的，不过我也不得不承认geohash确实是一个很好的索引模式，下面简单的介绍一下。
    具体的geohash算法，应该去参考Geohash - Wiki，我仅仅做一下简单的原理介绍。
    geohash将一个二维的坐标转化为一个简单的可排序可比较的字符串的方法是这样的：先把经纬度范围(纬度范围-90到90，经度范围-180到180)当作一个纯平面的矩形，这样就得到一个二维的平面坐标系。怎么将这个二维坐标转化为一个一维的坐标呢，请看下面的示意图：
    如图所示，平面二维坐标(经纬度)和1维的坐标(geohash之后的字符串)转化方式如下(我绘图的本领不高，将就看看吧)：

    也就是说geohash先将平面坐标平分为4块，按照以下的顺序进行编号，之后再对每一个子块以同样的方式进行划分，得到另一个编号，一直按照这个规则划分下去，最终编号越来越多，而方格越来越小，直到能表达该坐标需要的精度为止。而在划分之中的得到的每个编号对应一维坐标之中的一段，这些编号最后可以编码成一个字符串，具体的编码过程请看前面提到的wiki的网址，我这里只是表述二维坐标向一维坐标的转化过程。
    知道了这个hash的原理之后，我们来进行下面的讨论：
    1.两个离的越近，geohash的结果相同的位数越多，对么？这一点是有些用户对geohash的误解，虽然geo确实尽可能的将位置相近的点hash到了一起，可是这并不是严格意义上的(实际上也并不可能，因为毕竟多一维坐标)，例如在上图大方格4的左下部分的点和大方格1的右下部分的点离的很近，可是它们的geohash值一定是相差的相当远，因为头一次的分块就相差太大了，很多时候我们对geohash的值进行简单的排序比较，结果貌似真的能够找出相近的点，并且似乎还是按照距离的远近排列的，可是实际上会有一些点被漏掉了。
    2.既然不能做到将相近的点hash值也相近，那么geohash的意义何在呢？我觉得geohash还是相当有用的一个算法，毕竟这个算法通过无穷的细分，能确保将每一个小块的geohash值确保在一定的范围之内，这样就为灵活的周边查找和范围查找提供了可能。
    我看到过一些其他实现周边查找的方法(这里不讨论直接通过组合sql语句进行查询的模式，因为这种模式实际上是很耗性能的)，基本上都是建立一个单一的块索引，例如将经纬度以0.1为单位将地图划分为若干个单元格，建立每个坐标所在单元格编号的索引，这样做虽然直观，可能只能支持一种特定的查询方式，如果需要更改为实现以0.01为单位的查询，就需要重新建立单元格索引，这显然会带来一些问题。
    如果使用geohash则不用担心这个问题，hash的结果可以支持多个级别范围的查询，每个级别之间查询的单元范围大小是2倍的关系，这实际已经可以实现灵活的查询了，例如要实现指定范围的查找，我们可以指定一个深度(深度越大，每个单元块就越小)，然后将此地理范围所覆盖的所有单元块的hash范围全部找出来，然后逐个在这些范围内查找。最后，假如有必要的话，再对结果进行筛选，去除不在此范围内的点。
    下面我贴出我通过geohash实现一个指定范围查找的例子，我的网站也是在这个例子的基础上实现了周边查找，如果不想研究具体的实现，可以考虑直接使用此代码，如果需要详细的明白实现的过程建议先去通过前面的wiki链接了解geohash的原理。

class GetNearStations(BaseService):
2

def getIndexByTab(self,(t,n),bl):
3

#根据单元格的二维序号计算1维序号
4

num=0
5

#先考虑经度多分一层的情况
6

num+=n%2
7

n=(n-num)/2
8

#再逐次计算出num
9

for i in range(0,bl):
10

num+=(int(t%2))<<(2*i+1)
11

t=(t-(t%2))/2
12

num+=(int(n%2))<<(2*i+2)
13

n=(n-(n%2))/2
14

return num
15

def getHashByIndex(self,index,bl):
16

#根据单元格的一维序号计算单元格的起始哈希值
17

if(index>=pow(2,bl*2+1)):
18

return "zzzzzzzzzzzzzz";
19

BASE_32 = "0123456789bcdefghjkmnpqrstuvwxyz"
20

index=index<<(64-bl*2)
21

return "".join([BASE_32[(long(index) >> ((13-i-1)*5)) & 31] for i in range(0,13)])
22

def searchBounds(self,bounds,bl):
23

#根据bounds(minLat,minLon,maxLat,maxLon)查找范围内的点
24

#此函数应确保在此范围内的点都被返回，但不确保返回的一定严格在此范围内
25

spans=[]
26

#单元格大小
27

precision=180*pow(2,-bl)
28

for t in range(int(floor((bounds[0]+90)/precision)),int(ceil((bounds[2]+90)/precision))):
29

if t<0 or t>=pow(2,bl): continue
30

for n in range(int(floor((bounds[1]+180)/precision)),int(ceil((bounds[3]+180)/precision))):
31

if n<0 or n>=pow(2,bl+1): continue
32

#调试时显示方格范围
33

#print "%.5f,%.5f,%.5f,%.5f" % (t*precision-90,n*precision-180,(t+1)*precision-90,(n+1)*precision-180)
34

index=self.getIndexByTab((t,n),bl)
35

spans.append([index,index+1])
36

#合并相邻的内容
37

spans.sort(lambda x,y: cmp(x[0],y[0]))
38

for i in range(len(spans)-1,0,-1):
39

if spans[i-1][1]==spans[i][0]:
40

spans[i-1][1]=spans[i][1]
41

spans.pop(i)
42

return [(self.getHashByIndex(span[0],bl),self.getHashByIndex(span[1],bl)) for span in spans]
43

def get(self):
45

#此方法不完整，仅仅显示如何实现周边查找
46

#参数7表示深度，越大搜索单元范围越小，搜索的结果越精确，需要搜索的次数越多
47

spans=self.searchBounds((lat-0.5*span,lon-0.5*span,lat+0.5*span,lon+0.5*span),7)
48

stations=[]
49

for span in spans:
50

stations.extend(db.GqlQuery("select * from Train_stations where geohash >= '%s' and geohash<= '%s' order by geohash desc"%span).fetch(100))
51

self.outputJson(stations)

0 0