MongoDB地理空间数据存储及检索
之前写过MySQL空间索引简单使用，测试也是可用的，当时没有测试效率问题，因为存储的矢量数据都只是一个四点的多边形而已。这次使用mongoDB来做一个行政区划检索的功能，记录一下使用的过程。
参考资料：
1、存入地理数据
MongoDB存储的数据是bson结构，所以只要你的数据符合这个结构都是可以存储的，但是要支持空间索引，就必须按照它的规定来。
早期版本的(2.6之前)仅仅支持简单的点数据的索引，也就是filed:[x,y]这样的结构，这个适用范围太有限了。现在的版本支持GeoJSON形式的数据类型，且支持OCG的空间数据查询模型，使用上非常方便。
GeoJSON数据存入
MongoDB要求把GeoJSON格式的数据以子文档的形式存入，但实际上并不是存入一个完整GeoJSON对象，只需要其中的type和coordinates两个字段就可以了。
下面以存入一个含有地理空间数据的文档为例，把所有支持的GeoJSON对象类型都做个示例。
这里假设存储一个县的信息，数据都以json格式表示。
{
"xian":"潜山县",
"sheng":"安徽"
}
1、Ponit 点数据
现在假设为这个文档添加上中心点位置，那么这个文档就变成了如下样子：
{
"xian":"潜山县",
"sheng":"安徽"
"center":{ "type":"Point","coordinates":[116.45,30.72]}
}
2、LineString 线数据(多段线)
现在加上一个到省会合肥的路径连线，那么文档就变成了如下样子：
{
"xian":"潜山县",
"sheng":"安徽"
"center":{ "type":"Point","coordinates":[116.45,30.72]},
"toShengHui":{
"type":"LineString","
coordinates":[[116.55944824218749,30.58827267102698],
[116.87667846679689,30.791396195188927],[116.96594238281249,31.038815104128687],
[117.18292236328124,31.264465555752835],[117.22412109375,31.819230730326613]]}
}
3、 Polygon 多边形数据
多边形是当前地理信息领域应用的比较多的数据类型。
多边形描述的是一个面对象，其由两部分组成，一个外壳shell和0或多个內洞holes。
外壳和內洞的表示都是一个闭合的线环(LinearRing)，表示一个闭合曲面。外壳包括的区域表示在多边形内的区域，內洞表示的区域则是从外壳表示的区域中排除的区域，如下图所示的样子，蓝色的是外壳，绿色部分的是內洞。
因为行政边界的涉及到的点太多，所以这里就只添加一个外包框作为示例：
因为GeoJSON中使用bbox字段(四个double值的数组)来描述范围外包框，所以这里不能使用bbox来存储一个多边形，否则将报错error inserting documents: location object expected, location array not in correct format。但奇怪的是，这个bbox用于搜索的时候却是无效的。
这里就不使用內洞了，有內洞的情况就是coordinates数组中后面加上线环即可。
{
"xian":"潜山县",
"sheng":"安徽"
"center":{ "type":"Point","coordinates":[116.45,30.72]},
"toShengHui":{
"type":"LineString","
coordinates":[[116.55944824218749,30.58827267102698],
[116.87667846679689,30.791396195188927],[116.96594238281249,31.038815104128687],
[117.18292236328124,31.264465555752835],[117.22412109375,31.819230730326613]]},
"box":{
"type":"Polygon",
"coordinates":[[[116.40701293945311,30.454001045389525],
[116.77505493164062,30.454001045389525],[116.77505493164062,30.76248901825541],
[116.40701293945311,30.76248901825541],[116.40701293945311,30.454001045389525]]]}
}
4、MultiPoint多点、MultiLineString多线、MultiPolygon多多边形
对于多点、多线和多多边形，与单个的区别，也就是将coordinates成员改为一个数组形式，存入多个单个形式的坐标数据。
5、GeometryCollection 几何集合
几何集合就是多个几何对象的集合，就是一个数组里面放多个几何对象。
6、全国区县行政区划入库示例
1、首先下载全国的性质边界矢量数据，这个可以从https://www.gadm.org/download_country_v3.html下载。因为中国的矢量数据中没有台湾和香港澳门的数据，可以下在后合成一个。这份数据还有一些其他的小问题，这里就不提了。这也是我能找到的免费数据中较好的一份。
2、下载的数据可以使用GDAL或QGIS工具将其转换为geojson格式文档，也可以转换，直接写程序来读取。我把转换后的程序再经过了一次简化，因为所有的边界线都是MutilPolygon，而大多数边界是仅仅一个Polygon的，所以我把能转换的都转换成了Polygon。
3、因为MongoDB中存储字段bbox为GeoJSON中的数组形式在查询的时候会有问题，所以我把它改为了多边形。但bbox字段又不支持多边形，所以改为以box字段来存储。因为县级行政区划的边界都比较复杂，点比较多，在查询的时候会比较慢，所以使用$and来先查询box在查询geometry会比较快。
处理好的数据可以在这里下载 链接: https://pan.baidu.com/s/1f2c9FEQhkfDzC1dZHn6XmA 密码:5a2u
4、处理完成之后把所有的json对象按行写入了到了(县级边界.json.txt)文件中，因为每行都是一个json对象，可以使用mongoimport将其导入数据库中。但是直接导入会有一个问题，就是数据量太大，无法一次写入。所以我先拆分成了几个小文件。
split -l 512 县级边界.json.txt -d -a 1 中国县级行政边界_
执行上面命令后生成了5个小文件，然后逐个导入到mongodb即可。
mongoimport --host 192.168.0.28 --db us --collection xzbj --file 中国县级行政边界_0
2018-10-12T14:14:45.165+0800 connected to: 192.168.0.28
2018-10-12T14:14:47.386+0800 imported 512 documents
mongoimport --host 192.168.0.28 --db us --collection xzbj --file 中国县级行政边界_1
mongoimport --host 192.168.0.28 --db us --collection xzbj --file 中国县级行政边界_2
mongoimport --host 192.168.0.28 --db us --collection xzbj --file 中国县级行政边界_3
mongoimport --host 192.168.0.28 --db us --collection xzbj --file 中国县级行政边界_4
bbox只能以数组的形式存在，如图所示，我在处理的时候已经改用box多边形表示了。
2、创建地理索引
MongoDB的空间索引有三种，2dsphere、2d、geoHaystack。
对于某些地理空间查询操作，必须有相应的索引才行。
2.1、2dsphere索引
官网文档：https://docs.mongodb.com/manual/core/2dsphere/
2dsphere索引支持查询球面几何实体对象。2dsphere是MongoDB地理空间索引支持所有查询：用于查询、交点和接近。地理空间查询的更多信息，参见地理空间查询。
2dsphere索引支持点数据（包括传统点数据方式和GeoJSON的Point方式）。2dsphere索引由于是球面索引，所以仅仅支持经纬度数据，坐标系为WGS84。
版本2以后的MongoDB支持附加GeoJSON对象包括MultiPoint、MultiLineString、MultiPolygon和 GeometryCollection，具体的参考官方文档。
创建一个2dsphere索引的语句如下：
db.collection.createIndex( { <location field> : "2dsphere" } )
这里有一个问题，就是创建的时候，有MultiPolygon等不支持的几何类型的时候会出现错误"errmsg" : "Can't extract geo keys...，但看官网文档，这个是可以附加的类型，但目前没有找到相关的文档。
2.2、2d索引
官网文档：https://docs.mongodb.com/manual/core/2d/
2d索引对存储为二维平面上的点的数据使用。 2d索引适用于MongoDB 2.2及更早版本中使用的旧坐标对。
应对仅在下列情况下使用2d索引
您的数据库具有MongoDB 2.2或更早版本的遗留遗留坐标对，以及
您不打算将任何位置数据存储为GeoJSON对象。
创建一个2d索引的语句如下：
db.collection.createIndex( { <location field> : "2d" } )
要在具有非简单排序规则的集合上创建2d索引，必须在创建索引时显式指定{collation：{locale：“simple”}}。
2.3、geoHaystacks索引
官网文档：https://docs.mongodb.com/manual/core/geohaystack/
geoHaystack索引是一种特殊索引，优化小面积内的返回结果。geoHaystack索引可提高在平面进行几何(geometry)查询的性能。
对于使用球面的几何查询，2dsphere索引是一个比geoHaystack索引更好的选择。2dsphere索引允许字段重新排序。
geoHaystack索引要求第一个字段为location字段。此外，geoHaystack索引仅可通过命令使用，因此始终一次返回所有结果。
3、检索地理数据
检索这个也是看官方的文档比较快，这里只是一个简单的介绍。
官方文档：Geospatial Query Operators
3.1地理空间模型
MongoDB的地理空间查询可以实现球面或平面几何实体对象的查询。
2dsphere索引仅仅支持球面查询(即把点坐标数据当做球面经纬度处理)。2d索引支持平面查询(即将点坐标数据当做平面直角坐标系点坐标处理)和一些球面查询，虽然2d索引支持一些球面查询，但是对这些球面查询使用2d索引可能会导致错误，这样的数据尽量优先使用2dsphere索引。
下面列出每个地理空间操作使用的地理空间查询运算符，支持的查询和相关说明：
| 操作符 | 参数类型 | 索引 | 支持查询 | 说明 |
|---|---|---|---|---|
| $near 邻近查询 |
GeoJSON质心点在这个line和下一个line, | 2dsphere | 球面 | 另请参阅$ nearSphere运算符，该运算符在与GeoJSON和2dsphere索引一起使用时提供相同的功能 |
| $near | legacy coordinates | 2d | 平面 | |
| $nearSphere | GeoJSON点 | 2dsphere | 球面 | 提供与使用GeoJSON点和2dsphere索引的$near操作相同的功能 对于球面查询，可能最好使用$nearSphere，它明确指定名称中的球形查询而不是$near运算符 |
| $nearSphere | legacy coordinates | 2d | 球面 | 使用GeoJSON的点来代替 |
| $geoWithin 内部查询 |
GeoJSON几何对象 |
球面 | 查询完全在参数指定地理空间内的文档 | |
| $geoWithin | 2d | 平面 | 只能查询box框住的点 | |
| $geoWithin | 2d | 平面 | 同上 | |
| $geoWithin | 2d | 平面 | 同上，点加半径(弧度值)$geoWithin: { $center: [ [ <x>, <y> ] , <radius> ] } |
|
| $geoWithin | 2d/2dsphere | 球面 | 支持查询框住的GeoJSON对象 | |
| $geoIntersects 相交查询 |
{ $geometry: … } 多边形或多多边形 |
球面 | 查询与参数给定几何对象有相交关系的文档 |
还有$geoNear这个操作符，这里就不摘录了，直接去官网看好了。
因为这些操作符都比较简单，这里只单独介绍一下最有用的$geoIntersects操作符。
$geoIntersects操作符使用
选择地理空间数据与指定GeoJSON对象相交的文档; 即数据和指定对象的交集是非空的。
$geoIntersects运算符使用$geometry运算符来指定一个GeoJSON对象作为参数，使用默认坐标系（CRS）指定GeoJSON多边形或多边形的使用语法如下：
{
空间数据字段名: {
$geoIntersects: {
$geometry: {
type: "<GeoJSON对象类型>" ,
coordinates: [ <coordinates> ]
}
}
}
}
对于$geoIntersects查询,当指定的GeoJSON的几何对象大于半个球面时，使用默认的坐标系(CRS)会导致互补的几何对象在查询结果中。
3.0版中的新功能：要指定单环的GeoJSON多边形使用自定义MongoDB CRS，使用以下语法在$geometry表达式中指定自定义MongoDB CRS：
{
<location field>: {
$geoIntersects: {
$geometry: {
type: "Polygon" ,
coordinates: [ <coordinates> ],
crs: {
type: "name",
properties: { name: "urn:x-mongodb:crs:strictwinding:EPSG:4326" }
}
}
}
}
}
自定义MongoDB CRS使用逆时针顺序包覆，并允许$geoIntersects支持具有单环GeoJSON多边形的查询，该多边形的面积大于或等于单个半球。如果指定的多边形小于单个半球，则带有MongoDB CRS的$ geoIntersects的行为与默认的CRS相同。“Big” Polygons参考
如果指定纬度和经度坐标，请先列出经度然后列出纬度：
有效经度值介于-180和180之间（包括两者）。
有效纬度值介于-90和90之间（包括两者）。
3.2、查询示例（使用全国县级行政边界数据）
上面的数据导入之后，写几个查询的例子来测试一下。
3.2.1、使用$geoIntersects查询相交的区域
$geoIntersects用于查询与给定参数有相交区域的记录，只能查询geojson形式表示的位置字段，有没有索引都可以查，速度与几何对象复杂程度有关。
查询代码如下
db.getCollection('xzbj').find({
"geometry": {
"$geoIntersects": {
"$geometry": {
"type": "Polygon",
"coordinates": [[[116.24633789062499, 40.168380093142446], [116.17492675781251, 40.15998434802335], [116.1199951171875, 40.057052221322], [116.09527587890624, 40.002371935876475], [116.1474609375, 39.890772566959534], [116.10626220703124, 39.70929962338767], [116.3177490234375, 39.69662085337441], [116.57592773437499, 39.7642140375156], [116.6912841796875, 39.86969567045658], [116.69677734375, 39.99605985169435], [116.62261962890624, 40.094882122321145], [116.6143798828125, 40.13899044275822], [116.43310546875, 40.15788524950653], [116.28753662109375, 40.18097176388719], [116.24633789062499, 40.168380093142446]]]
}
}
}
})
这些点是沿着北京六环线画的一个多边形，查询的速度比较慢，耗时达到7.12秒。
下面添加外包框过滤，在进行相交比较，加快查询速度。这里要注意我是要的box字段内存的也是一个Polygon而不是GeoJSON中的bbox。
db.getCollection('xzbj').find({
"$and": [{
"box": {
"$geoIntersects": {
"$geometry": {
"type": "Polygon",
"coordinates": [[[116.0870361328125, 39.69873414348139], [116.71600341796874, 39.69873414348139], [116.71600341796874, 40.17257757632168], [116.0870361328125, 40.17257757632168], [116.0870361328125, 39.69873414348139]]]
}
}
}
},
{
"geometry": {
"$geoIntersects": {
"$geometry": {
"type": "Polygon",
"coordinates": [[[116.24633789062499, 40.168380093142446], [116.17492675781251, 40.15998434802335], [116.1199951171875, 40.057052221322], [116.09527587890624, 40.002371935876475], [116.1474609375, 39.890772566959534], [116.10626220703124, 39.70929962338767], [116.3177490234375, 39.69662085337441], [116.57592773437499, 39.7642140375156], [116.6912841796875, 39.86969567045658], [116.69677734375, 39.99605985169435], [116.62261962890624, 40.094882122321145], [116.6143798828125, 40.13899044275822], [116.43310546875, 40.15788524950653], [116.28753662109375, 40.18097176388719], [116.24633789062499, 40.168380093142446]]]
}
}
}
}]
})
这样查询的速度就大大提升了，仅耗时41毫秒就完成了检索。经过测试，$and数组中的元素顺序对检索速度没有影响，不知道是不是与字段名的排序有关系。
3.2.3 使用$geoWithin查询
使用上和$geoIntersects差不多，查询速度上也差不多。
db.getCollection('xzbj').find({
"geometry": {
"$geoWithin": {
"$geometry": {
"type": "Polygon",
"coordinates": [[[93.69140625, 28.76765910569123], [113.37890625, 28.76765910569123], [113.37890625, 39.30029918615029], [93.69140625, 39.30029918615029], [93.69140625, 28.76765910569123]]]
}
}
}
})
可以指定参数几何对象的坐标系。
db.getCollection('xzbj').find({
"box": {
"$geoWithin": {
"$geometry": {
"type": "Polygon",
"coordinates": [[[93.69140625, 28.76765910569123], [113.37890625, 28.76765910569123], [113.37890625, 39.30029918615029], [93.69140625, 39.30029918615029], [93.69140625, 28.76765910569123]]],
"crs": {
"type": "name",
"properties": {
"name": "urn:x-mongodb:crs:strictwinding:EPSG:4326"
}
}
}
}
}
})
上面使用Polygon作为查询参数，没有索引也可以查。如果使用$center(中心点加半径)和$box查询，则仅有存在2d索引的情况才能查询。
使用$centerSphere作为查询参数的时候，有没有索引都可以查，而且速度很快。可以查询GeoJSON表示的字段，没有类型限制。
$centerShpere的参数也是一个点(先经度后纬度)加一个弧度，示例如下:
db.getCollection('xzbj').find({
"box": {
"$geoWithin": { "$centerSphere": [ [116.3623809,39.9013085] ,0.008 ] }
}
})
上的很快就搜索到了，因为我对box字段建了2dSphere索引。把搜索的字段换为geometry就比较慢了，因为没有建索引。
3.2.2 其他的
其他的这里就不记录了，需要的时候查询即可。
下面是我对QGIS转出的geojson进行提取的代码，放在这里做个存档。因为这是简单的使用一次，没有处理各种错误异常等。
#include <iostream>
#include <rapidjson/filereadstream.h>
#include <rapidjson/rapidjson.h>
#include <rapidjson/document.h>
#include <rapidjson/stringbuffer.h>
#include <rapidjson/writer.h>
#include <Poco/Net/HTTPClientSession.h>
#include <Poco/Net/HTTPRequest.h>
#include <Poco/Net/HTTPResponse.h>
#include <Poco/URI.h>
#include <Poco/StreamCopier.h>
int main()
{
FILE* fp = fopen("./xianjie.geojson","rb");
char buffer[65536];
rapidjson::FileReadStream is(fp,buffer,sizeof(buffer));
rapidjson::Document doc;
doc.ParseStream(is);
if(doc.HasParseError()){
std::cout<<"parse error:"<< doc.GetParseError()<<std::endl;
return 0;
}
rapidjson::Value::ConstMemberIterator iter_features = doc.FindMember("features");
if(iter_features == doc.MemberEnd()){
std::cerr<<"没有找到 iter_features"<<std::endl;
return 0;
}
const auto& features = iter_features->value.GetArray();
// std::cout<<"features size = "<< features.Size()<<std::endl;
for(size_t i =0;i<features.Size();++i){
auto& obj = features[i];
rapidjson::Value::ConstMemberIterator iter_properties = obj.FindMember("properties");
if(iter_properties == obj.MemberEnd()){
std::cout<<"features["<<i<<"] format error"<<std::endl;
continue;
}
std::string sheng = iter_properties->value["NL_NAME_1"].GetString();
size_t pos = sheng.rfind('|');
if(pos != std::string::npos){
sheng = sheng.substr(pos+1);
}
std::string di = iter_properties->value["NL_NAME_2"].GetString();
pos = di.find('|');
if(pos != std::string::npos){
di = di.substr(0,pos);
}
std::string xian = di;
if(iter_properties->value.HasMember("NL_NAME_3") && iter_properties->value["NL_NAME_3"].IsString()){
xian = iter_properties->value["NL_NAME_3"].GetString();
pos = xian.find('|');
if(pos != std::string::npos){
xian = xian.substr(0,pos);
}
}
std::string fullname = sheng + "-" + di + "-" + x2;
rapidjson::Value::ConstMemberIterator iter_geometry = obj.FindMember("geometry");
if(iter_geometry == obj.MemberEnd()){
std::cout<<"Failed "<<fullname<<std::endl;
continue;
}
std::string geometry;
{
rapidjson::StringBuffer buffer;
rapidjson::Writer<rapidjson::StringBuffer> w(buffer);
if(iter_geometry->value["coordinates"].GetArray().Size() == 1){
iter_geometry->value["coordinates"].GetArray()[0].Accept(w);
geometry.append("{\"type\":\"Polygon\",\"coordinates\":");
geometry.append(buffer.GetString(),buffer.GetSize());
geometry.append("}");
}
else{
obj.Accept(w);
geometry.assign(buffer.GetString(),buffer.GetSize());
}
}
double bbox[4];
{
bbox[0] = obj["bbox"].GetArray()[0].GetDouble();
bbox[1] = obj["bbox"].GetArray()[1].GetDouble();
bbox[2] = obj["bbox"].GetArray()[2].GetDouble();
bbox[3] = obj["bbox"].GetArray()[3].GetDouble();
}
std::cout<<"{\"name\":\""<<fullname<<"\",\"type\":\"Feature\","
<<"\"properties\":{\"sheng\":\""<<sheng<<"\",\"di\":\""<<di<<"\",\"xian\":\""<<xian
<<"\"},\"box\":{\"type\":\"Polygon\",\"coordinates\":[["
<<"["<<bbox[0]<<","<<bbox[1]<<"],"
<<"["<<bbox[2]<<","<<bbox[1]<<"],"
<<"["<<bbox[2]<<","<<bbox[3]<<"],"
<<"["<<bbox[0]<<","<<bbox[3]<<"],"
<<"["<<bbox[0]<<","<<bbox[1]<<"]]]}"
/*<<"\"},\"bbox\":["<<bbox[0]<<","<<bbox[1]<<","<<bbox[2]<<","<<bbox[3]<<"]"*/
<<",\"geometry\":"<<geometry<<"}\n";
#ifdef CouchDB_Insert
Poco::Net::HTTPClientSession session("192.168.0.28",5984);
Poco::Net::HTTPRequest request(Poco::Net::HTTPRequest::HTTP_PUT ,"/xzbj/" + fullname);
request.setContentType("application/json");
request.setContentLength((int)geometry.size());
request.set("Authorization","Basic 用户名密码base64");
std::ostream& ss = session.sendRequest(request);
ss.write(geometry.data(),geometry.size());
Poco::Net::HTTPResponse response;
std::istream& rs = session.receiveResponse(response);
std::cout<<"\n\n"<<request.getURI()<<"\t\t"<<response.getStatus()<<std::endl;
Poco::StreamCopier copier;
copier.copyStream(rs,std::cout);
#endif
}
return 0;
}
