geohash

在上次写了文章《在GAE之中使用基于地理位置的查询》，之后，我一直在奇怪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的原理。

 1class 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                
44    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)