Geohash
Geohash算法就是将经纬度编码，将二维变一维，给地址位置分区的一种算法。
基本原理：GeoHash是一种地址编码方法。他能够把二维的空间经纬度数据编码成一个字符串
如果纬度范围[-90°, 0°)用二进制0代表，（0°, 90°]用二进制1代表，经度范围[-180°, 0°)用二进制0代表，（0°, 180°]用二进制1代表，那么地球可以分成如下4个部分
如果在小块范围内递归对半划分呢？
可以看到，划分的区域更多了，也更精确了。geohash算法就是基于这种思想，划分的次数更多，区域更多，区域面积更小了。通过将经纬度编码，给地理位置分区
Geohash算法
Geohash算法一共有三步。
首先将经纬度变成二进制。
比如这样一个点（39.923201, 116.390705）
纬度的范围是（-90，90），其中间值为0。对于纬度39.923201，在区间（0，90）中，因此得到一个1；（0，90）区间的中间值为45度，纬度39.923201小于45，因此得到一个0，依次计算下去，即可得到纬度的二进制表示，如下表：
10111000110001111001
同理可以得到经度116.390705的二进制表示为：
11010010110001000100
第2步，就是将经纬度合并。
经度占偶数位，纬度占奇数位，注意，0也是偶数位。
11100 11101 00100 01111 00000 01101 01011 00001第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算法有两个问题。首先是边缘问题。
如图，如果车在红点位置，区域内还有一个黄点。相邻区域内的绿点明显离红点更近。但因为黄点的编码和红点一样，最终找到的将是黄点。这就有问题了。
要解决这个问题，很简单，只要再查找周边8个区域内的点，看哪个离自己更近即可。
另外就是曲线突变问题。
本文第2张图片比较好地解释了这个问题。其中0111和1000两个编码非常相近，但它们的实际距离确很远。所以编码相近的两个单位，并不一定真实距离很近，这需要实际计算两个点的距离才行。
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手册。
参考：https://www.jianshu.com/p/2fd0cf12e5ba
参考：https://www.shangmayuan.com/a/ddb0b7d19cec4ab4aecc18f5.html