GEO LBS服务

Geohash

Geohash算法就是将经纬度编码,将二维变一维,给地址位置分区的一种算法。

基本原理:GeoHash是一种地址编码方法。他能够把二维的空间经纬度数据编码成一个字符串

我们知道,经度范围是东经180到西经180,纬度范围是南纬90到北纬90,我们设定西经为负,南纬为负,所以地球上的经度范围就是[-180, 180],纬度范围就是[-90,90]。如果以本初子午线、赤道为界,地球可以分成4个部分。

如果纬度范围[-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
posted on 2021-08-03 15:41  duanxz  阅读(212)  评论(0编辑  收藏  举报