第3步，按照Base32进行编码

Base32编码表的其中一种如下，是用0-9、b-z（去掉a, i, l, o）这32个字母进行编码。具体操作是先将上一步得到的合并后二进制转换为10进制数据，然后对应生成Base32码。需要注意的是，将5个二进制位转换成一个base32码。上例最终得到的值为

  wx4g0ec1

Geohash比直接用经纬度的高效很多，而且使用者可以发布地址编码，既能表明自己位于北海公园附近，又不至于暴露自己的精确坐标，有助于隐私保护。

• GeoHash用一个字符串表示经度和纬度两个坐标。在数据库中可以实现在一列上应用索引（某些情况下无法在两列上同时应用索引）
• GeoHash表示的并不是一个点，而是一个矩形区域
• GeoHash编码的前缀可以表示更大的区域。例如wx4g0ec1，它的前缀wx4g0e表示包含编码wx4g0ec1在内的更大范围。 这个特性可以用于附近地点搜索

编码越长，表示的范围越小，位置也越精确。因此我们就可以通过比较GeoHash匹配的位数来判断两个点之间的大概距离。

 

问题

geohash算法有两个问题。首先是边缘问题。

MongoDb geohash

1、关于MongoDB

在众多NoSQL数据库，MongoDB是一个优秀的产品。其官方介绍以下： 
MongoDB (from "humongous") is a scalable, high-performance, open source, document-oriented database.html

看起来，十分诱人！值得说明的是，MongoDB的document是以BSON（Binary JSON）格式存储的，彻底支持Schema Free。这对地理空间数据是十分友好的。由于有著名的GeoJSON可供使用。另外OGR库也支持将Geometry类型导出为JSON格式。mongodb

本文将尝试使用OGR库把Shapefile导入到MongoDB存储，而后创建空间索引，进行空间查询。数据库

著名的Foursquare使用了MongoDB数据库。json

2、开发环境

MongoDB+Python+Pymongo+GDAL for Python数组

关于MongoDB和Python安装，本文不作介绍。关于GDAL for Python的安装，可见个人另外一篇博文：http://blog.3sdn.net/311.html服务器

在继续本文以前，请先启动你的MongoDB服务器。本文默认采用以下服务器参数： 
Server:localhost 
Post:27017 
Database Name:gisdbapp

3、将shapefile导入到MongoDB

这里我直接提供代码，代码中已经有比较详尽的注释了。代码基本源于“引文1”，只是作了些改动，将MongoDB的Geometry的存储格式由wkt改为json。你可直接复制并运行下面的代码，固然须要修改一下Shapefile路径和MongoDB服务器相关参数。性能

import os 
import sys 
import json 
from pymongo import json_util 
from pymongo.connection import Connection 
from progressbar import ProgressBar 
from osgeo import ogr测试

def shp2mongodb(shape_path, mongodb_server, mongodb_port, mongodb_db, mongodb_collection, append, query_filter): 
        """Convert a shapefile to a mongodb collection""" 
        print ‘Converting a shapefile to a mongodb collection ‘ 
        driver = ogr.GetDriverByName(‘ESRI Shapefile’) 
        print ‘Opening the shapefile %s…’ % shape_path 
        ds = driver.Open(shape_path, 0) 
        if ds is None: 
                print ‘Can not open’, ds 
                sys.exit(1) 
        lyr = ds.GetLayer() 
        totfeats = lyr.GetFeatureCount() 
        lyr.SetAttributeFilter(query_filter) 
        print ‘Starting to load %s of %s features in shapefile %s to MongoDB…’ % (lyr.GetFeatureCount(), totfeats, lyr.GetName()) 
        print ‘Opening MongoDB connection to server %s:%i…’ % (mongodb_server, mongodb_port) 
        connection = Connection(mongodb_server, mongodb_port) 
        print ‘Getting database %s’ % mongodb_db 
        db = connection[mongodb_db] 
        print ‘Getting the collection %s’ % mongodb_collection 
        collection = db[mongodb_collection] 
        if append == False: 
                print ‘Removing features from the collection…’ 
                collection.remove({}) 
        print ‘Starting loading features…’ 
        # define the progressbar 
        pbar = ProgressBar(maxval=lyr.GetFeatureCount()).start() 
        k=0 
        # iterate the features and access its attributes (including geometry) to store them in MongoDb 
        feat = lyr.GetNextFeature() 
        while feat: 
                mongofeat = {} 
                geom = feat.GetGeometryRef() 
                mongogeom = geom.ExportToJson() 
                # store the geometry data with json format 
                mongofeat[‘geom‘] = json.loads(mongogeom,object_hook=json_util.object_hook)
                # iterate the feature’s  fields to get its values and store them in MongoDb 
                feat_defn = lyr.GetLayerDefn() 
                for i in range(feat_defn.GetFieldCount()): 
                        value = feat.GetField(i) 
                        if isinstance(value, str): 
                                value = unicode(value, "gb2312") 
                        field = feat.GetFieldDefnRef(i) 
                        fieldname = field.GetName() 
                        mongofeat[fieldname] = value 
                # insert the feature in the collection 
                collection.insert(mongofeat) 
                feat.Destroy() 
                feat = lyr.GetNextFeature() 
                k = k + 1 
                pbar.update(k) 
        pbar.finish() 
        print ‘%s features loaded in MongoDb from shapefile.’ % lyr.GetFeatureCount() 
        
        
input_shape = ‘/home/evan/data/map/res4_4m/XianCh_point.shp’ 
mongodb_server = ‘localhost’ 
mongodb_port = 27017 
mongodb_db = ‘gisdb’ 
mongodb_collection = ‘xqpoint’ 
filter = ”spa

print ‘Importing data to mongodb…’ 
shp2mongodb(input_shape, mongodb_server, mongodb_port, mongodb_db, mongodb_collection, False, filter)

4、MongoDB中空间数据的存储格式

在MongoDB的Shell中执行： 
>db.xqpoint.findOne() 
结果以下：

{ 
    "_id" : ObjectId("4dc82e7f7de36a5ceb000000"), 
    "PERIMETER" : 0, 
    "NAME" : "漠河县", 
    "PYNAME" : "Mohe Xian", 
    "AREA" : 0, 
    "ADCODE93" : 232723, 
    "CNTYPT_ID" : 31, 
    "CNTYPT_" : 1, 
    "geom" : { 
        "type" : "Point", 
        "coordinates" : [ 
            122.53233, 
            52.968872 
        ] 
    }, 
    "ID" : 1031, 
    "PN" : 1, 
    "CLASS" : "AI" 
} 

这即是一个document，使用JSON格式，一目了然。其中的"geom"即为Geometry类型的数据，即地理空间数据，也是采用JSON格式存储，这样后续的空间索引与空间查询将十分方便。

MongoDB原生地支持了空间索引与空间查询，这一点比PostgreSQL方便，再也不须要使用PostGIS进行空间扩展了。至于性能，我还没测试，在此不敢妄加评论。

5、在MongoDB中创建空间索引

>db.xqpoint.ensureIndex({‘geom.coordinates’:’2d’})

是否是十分简单？其它参数及用法请自行查看MongoDB手册。

6、在MongoDB中进行空间查询

>db.xqpoint.find({"geom.coordinates":[122.53233,52.968872]})

便可查询到上述“莫河县”这个点。固然，像这种精确查询，实际应用并很少。实际应用的空间查询大多为范围查询。MongoDB支持邻域查询（$near），和范围查询（$within）。

1. 邻域查询($near)

>db.xqpoint.find({"geom.coordinates":{$near:[122,52]}}) 
上述查询语句查询点[122,52]附近的点，MongoDB默认返回附近的100个点，并按距离排序。你也能够用limit()指定返回的结果数量， 如：>db.xqpoint.find({"geom.coordinates":{$near:[122,52]}}).limit(5)

另外，你也能够指定一个最大距离，只查询这个距离内的点。 
>db.xqpoint.find({"geom.coordinates":{$near:[122,52],$maxDistance:5}}).limit(5)

MongoDB的find()方法可很方便的进行查询，同时MongoDB也提供了geoNear命令，用于邻域查询。 
>db.runCommand({geoNear:"xqpoint",near:[122,56],num:2}) 
上述语句用于查询[122,56]点附近的点，并只返回2个点。结果以下：

{ 
    "ns" : "gisdb.xqpoint", 
    "near" : "1110011000111101111100010000011000111101111100010000", 
    "results" : [ 
        { 
            "dis" : 3.077515616588727, 
            "obj" : { 
                "_id" : ObjectId("4dc82e7f7de36a5ceb000000"), 
                "PERIMETER" : 0, 
                "NAME" : "漠河县", 
                "PYNAME" : "Mohe Xian", 
                "AREA" : 0, 
                "ADCODE93" : 232723, 
                "CNTYPT_ID" : 31, 
                "CNTYPT_" : 1, 
                "geom" : { 
                    "type" : "Point", 
                    "coordinates" : [ 
                        122.53233, 
                        52.968872 
                    ] 
                }, 
                "ID" : 1031, 
                "PN" : 1, 
                "CLASS" : "AI" 
            } 
        }, 
        { 
            "dis" : 4.551319677334594, 
            "obj" : { 
                "_id" : ObjectId("4dc82e7f7de36a5ceb000001"), 
                "PERIMETER" : 0, 
                "NAME" : "塔河县", 
                "PYNAME" : "Tahe Xian", 
                "AREA" : 0, 
                "ADCODE93" : 232722, 
                "CNTYPT_ID" : 66, 
                "CNTYPT_" : 2, 
                "geom" : { 
                    "type" : "Point", 
                    "coordinates" : [ 
                        124.7058, 
                        52.340332 
                    ] 
                }, 
                "ID" : 1059, 
                "PN" : 1, 
                "CLASS" : "AI" 
            } 
        } 
    ], 
    "stats" : { 
        "time" : 0, 
        "btreelocs" : 85, 
        "nscanned" : 85, 
        "objectsLoaded" : 4, 
        "avgDistance" : 3.814417646961661, 
        "maxDistance" : 4.551319677334594 
    }, 
    "ok" : 1 
}

固然，咱们也可附加条件查询条件，如查询[122,56]附近的且"PYNAME"为"Tahe Xian"的点： 
>db.runCommand({geoNear:"xqpoint",near:[122,56],num:2,query:{"PYNAME":"Tahe Xian"}) 
返回结果以下：

{ 
    "ns" : "gisdb.xqpoint", 
    "near" : "1110011000111101111100010000011000111101111100010000", 
    "results" : [ 
        { 
            "dis" : 4.551319677334594, 
            "obj" : { 
                "_id" : ObjectId("4dc82e7f7de36a5ceb000001"), 
                "PERIMETER" : 0, 
                "NAME" : "塔河县", 
                "PYNAME" : "Tahe Xian", 
                "AREA" : 0, 
                "ADCODE93" : 232722, 
                "CNTYPT_ID" : 66, 
                "CNTYPT_" : 2, 
                "geom" : { 
                    "type" : "Point", 
                    "coordinates" : [ 
                        124.7058, 
                        52.340332 
                    ] 
                }, 
                "ID" : 1059, 
                "PN" : 1, 
                "CLASS" : "AI" 
            } 
        } 
    ], 
    "stats" : { 
        "time" : 45, 
        "btreelocs" : 2095, 
        "nscanned" : 2096, 
        "objectsLoaded" : 2096, 
        "avgDistance" : 4.551319677334594, 
        "maxDistance" : 4.551319677334594 
    }, 
    "ok" : 1 
}

2. 范围查询($within)

MongoDB的$within操做符支持的形状有$box（矩形）,$center（圆形）,$polygon（多边形，包括凹多边形和凸多边形）。全部的范围查询，默认是包含边界的。

查询一个矩形范围，须要指定矩形的左下角和右上角两个坐标点，以下： 
> box = [[80,40],[100,50]] 
> db.xqpoint.find({"geom.coordinates":{$within:{$box:box}}})

查询一个圆形范围，须要指定圆心坐标和半径，以下： 
> center = [80,44] 
> radius =5 
> db.xqpoint.find({"geom.coordinates":{$within:{$center:[center,radius]}}})

查询一个多边形范围，须要指定多边形的各个顶点，能够经过一个顶点数组或一系列点对象指定。其中，最后一个点是默认与第一个点链接的。以下： 
> polygon1 = [[75,35],[80,35],[80,45],[60,40]] 
> db.xqpoint.find({"geom.coordinates":{$within:{$polygon:polygon1}}}) 
或者 
> polygon2 = ｛a:{75,35},b:{80,35},c:{80,45},d:{60,40}} 
> db.xqpoint.find({"geom.coordinates":{$within:{$polygon:polygon2}}})

注意：MongoDB 1.9及以上版本才支持多边形范围查询。

P.S. MongoDB还支持复合索引，球面模型（可简单理解为投影吧），多位置文档（Multi-location Documents，即一个文档中包括多个Geometry），可参见“引文2”或MongoDB手册。