开源软件在地图数据处理中的应用

2 年前 · 来自专栏研发提效

滴滴云

作为开源软件的受益者，在享受开源带来的技术便利同时，我们也积极拥抱开源，同时也回馈开源。城市交通指数（TTI）作为公司第16个开源项目，通过盖亚计划对外开放了脱敏数据，下载人员分布于127个高校或科研机构，覆盖了70%的双一流高校。

在地图数据处理中，能经常看到开源软件的身影，常见的有以下几项：

GDAL (Geospatial Data Abstraction Library)

GDAL是一个在X/MIT许可协议下的空间栅格数据转换库，OGR是GDAL项目的一个分支，其提供对矢量数据的支持。栅格数据和矢量数据是地图数据中两种较为常见的数据格式，通俗的理解，栅格数据用像素来表达，矢量数据用坐标点来表达。常见的栅格数据有遥感影像、扫描地图等。常见的矢量数据有各种点、线、面数据，如POI、路网、水系或湖泊。GDAL可以很方便的对栅格或矢量数据进行读写操作。

GDAL读取遥感影像示例代码：

GDALDataset* pDataSet = (GDALDataset*)GDALOpen("/Users/didi/Desktop/test.img",GA_ReadOnly);
//仿射变换6参数
double geoTransform[6] = {0};
pDataSet->GetGeoTransform(geoTransform);
//影像宽
int nWidth = pDataSet->GetRasterXSize();
//影像高
int nHeight = pDataSet->GetRasterYSize();
//像素值矩阵
unsigned char* pPixelValue = (unsigned char *)malloc(sizeof(unsigned char) * nWidth * nHeight);
memset(pPixelValue,0,nWidth * nHeight);
CPLErr err = pDataSet->RasterIO(GF_Read,0,0,nWidth,nHeight,pPixelValue,nWidth,nHeight,GDT_Byte,1,NULL,0,0,0);
free(pPixelValue);
GDALClose(pDataSet);

GDAL中有一些很实用的图像处理算法，如GDALSimpleSURF类可以实现特征点检测以及匹配，GDALSieveFilter可以删除4连通或8连通的像素多边形（具有相同像素值的连通区域），如砂眼噪声。如下图所示

图1 处理前

图2 处理后

OGR读取矢量数据示例代码：

GDALDataset* pDS = (GDALDataset*) GDALOpenEx("/Users/didi/Desktop/test.shp", GDAL_OF_VECTOR | GDAL_OF_READONLY, NULL, NULL, NULL );
OGRLayer* pLayer = pDS->GetLayer(0);
OGRFeature *pFeature = NULL;
pLayer->ResetReading();
while((pFeature = pLayer->GetNextFeature()) != NULL ){
    //获取几何信息
    OGRGeometry* pGeom = pFeature->GetGeometryRef();
    //获取字段Length的值
    double dfLength = pFeature->GetFieldAsDouble("Length");
    OGRFeature::DestroyFeature(pFeature);
//关闭数据集
GDALClose(pDS);

GEOS(Geometry Engine Open Source)
GEOS是非常经典的空间拓扑关系运算库，一些算法性能不亚于商业GIS软件。GDAL读取到矢量数据的几何信息后，在进行空间运算或拓扑判断时，就需要GEOS的支持了。一般情况下，在编译GDAL源码时直接将GEOS编译进去，也可以把GDAL中的几何类转换成GEOS的几何类，这样就可以直接使用GEOS库进行数据处理。使用GEOS可以很方便的计算两个几何对象的差集、交集、对称差、缓冲区。有时候我们需要将多个相邻的多边形合并成一个多边形，常规用法是使用union方法，将其合并，当待合并的多边形个数较多时，效率就会非常的低，这里我们可以使用计算缓冲区的方法进行处理，效率会提升很多。

图3 待合并多边形

图4 合并结果图

示例代码：

//蓝色多边形
char* szWKT_1 = "POLYGON ((113.885 22.6815, 113.9425 22.6585, 113.91 22.7, 113.885 22.6815))";
//橙色多边形
char* szWKT_2 = "POLYGON ((113.91 22.7, 113.9425 22.6585, 113.9675 22.689, 113.91 22.7))";
OGRGeometry* pGeom_1 = NULL;
OGRGeometry* pGeom_2 = NULL;
OGRGeometryFactory::createFromWkt(&szWKT_1, NULL, &pGeom_1);
OGRGeometryFactory::createFromWkt(&szWKT_2, NULL, &pGeom_2);
OGRMultiPolygon* pMultiPolygon = (OGRMultiPolygon*)OGRGeometryFactory::createGeometry(wkbMultiPolygon);
pMultiPolygon->addGeometryDirectly(pGeom_1);
pMultiPolygon->addGeometryDirectly(pGeom_2);
//用Buffer替代Union，缓冲距离设置为0
//pUnion为紫色多边形
OGRGeometry* pUnion = pMultiPolygon->Buffer(0);

PROJ.4

Proj.4是开源GIS中最著名的地图投影库，许多GIS开源软件的投影都直接使用Proj.4的库文件。主要功能有投影坐标与地理坐标的转换，坐标系的转换，基准变换等。WGS84坐标系是我们最常用的坐标系，GPS轨迹点就是地理坐标系。有些时候我们需要基于投影坐标系做一些数学运算，即高斯正算。如长度计算。这里取以3°分带，中央经线为117°的区域内的一条link来说明：

siwei_id:9000025617988

地理坐标：LINESTRING (116.9710037 36.6434771,116.97049 36.64324,116.97014 36.64304,116.96921 36.64244,116.96787 36.64177)

投影坐标：LINESTRING (497406.947036178 4057018.36084719,497361.000282045 4056992.06317016,497329.693821681 4056969.87826856,497246.504857284 4056903.32080414,497126.64614413 4056829.00824338)，在投影坐标下，我们可以直接使用两点之间距离公式来计算这条link的长度，length=339米。

libSpatialindex

libspatialindex是一种高效的C++空间索引库。支持复杂查询，如范围查询、点位置查询、最近邻查询、K邻近查询以及参数化查询。创建内存空间索引示例代码：

IStorageManager* diskfile = StorageManager::createNewMemoryStorageManager();
StorageManager::IBuffer* file = StorageManager::createNewRandomEvictionsBuffer(*diskfile, 10, false);
double fillFactor = 0.7;
uint32_t indexCapacity = 10;
uint32_t leafCapacity = 10;
uint32_t dimension = 2;
RTree::RTreeVariant variant = RTree::RV_RSTAR;
id_type indexIdentifier;
ISpatialIndex* tree = RTree::createNewRTree(*file, fillFactor, indexCapacity, leafCapacity, dimension, variant, indexIdentifier);
x1 = enve.MinX;
y1 = enve.MinY;
x2 = enve.MaxX;
y2 = enve.MaxY;
plow[0] = x1;
plow[1] = y1;
phigh[0] = x2;
phigh[1] = y2;
//r为几何对象的外接矩形，id为唯一的标识符，为长整型
Region r = Region(plow, phigh, 2);
tree->insertData(0, 0, r, id);
delete tree;
delete file;
delete diskfile;

SpatiaLite

SQLite是一款很小巧的关系数据库，经常被集成在各种移动应用程序中，为了扩充SQLite的空间数据存储功能，基于SQLite的内核，增加了空间SQL功能。 SpatiaLite使用RTree作为空间索引，对一张表创建空间索引后，可以进行高效的空间查询处理。创建空间索引示例代码：select CreateSpatialIndex('table_name', 'column_name')

PCL（Point Cloud Library）

激光点云是近些年使用较多一种地图矢量数据，可以高效获取目标的三维坐标。由于其高昂的配套设备价格，只有少数图商能玩的起。对于3D点云处理来说，PCL完全是一个模块化的C++模板库，基于第三方库实现点云相关的获取、滤波、分割、配准、检索、特征提取、识别、追踪、曲面重建、可视化等。

PostGIS

PostGIS只是PostgreSQL的一个插件，但是它将PostgreSQL变成了一个强大的空间数据库。PostGIS具有优异的空间查询性能，如果说Spatialite、libspatialindex等空间数据存储工具是一艘军舰，那么PostGIS就是一艘航空母舰。其中pgRouting增加了路由功能，实现了导航路径规划中的经典算法，如Dijkstra算法、A*算法、旅行商算法。将地图数据整理成符合其规格的格式后，可以很快打造出一个简易的路径规划计算服务。在空间查询中，PostGIS可以实现即查即所得的效果。图6中，如果要获取橙色多边形内的路网数据，普通空间查询结果为绿色查询结果集合，如果要获取多边形内的数据，需要在查询结果集中做二次空间拓扑判断。在PostGIS中，由于其特有的GiST空间索引，一次空间查询即可，大大提高查询效率。

图6 空间查询效果图

GeoServer

GeoServer 是 OpenGIS Web 服务器规范的 J2EE 实现，利用 GeoServer 可以方便的发布地图数据，允许用户对特征数据进行更新、删除、插入操作。其数据源支持PostGIS、MySQL、GeoMesa等。用户编辑SLD文件可以自定义地图绘制风格，制作各种专题地图。下图为使用不同的SLD文件，发布WMS服务示意图。

图7 轨迹图

图8 轨迹热力图

OpenLayer、Mapbox

OpenLayer和MapBox应用场景较为类似，为Web GIS 客户端开发提供的JavaScript 类库包，用于实现标准格式发布的地图数据访问。基于OpenLayer和MapBox可以让前端地图显示的更加漂亮。

QGIS

QGIS是一个免费的、开源的、跨平台(Linux/Windows/Mac/Android)的地理信息系统桌面版软件，基于QT C++开发。它提供强大的空间数据显示、编辑和分析功能。QGIS和很多开源项目一样，使用CMake进行编译，由于QGIS集成了大部分的第三方软件，基于源码编译出可执行的QGIS，这是一个考验耐心的活，如同自己盖一座大楼，每一块砖头都需要自己造。编译出Debug版本的QGIS后，就能调试其源代码，并且可以基于SDK进行二次开发。QGIS的设计架构，很精妙，很值得学习。

Mapnik

Mapnik是一款开源的地图渲染引擎，它能够为PostGIS，Shapefile，Geojson，SQLite等在内的多种数据源提供空间数据计算与可视化服务，包括png瓦片，矢量瓦片，同时它支持自定义渲染样式配置，具有很高的灵活性，提供了C++，python，node接口。Open Street Map主地图层就是用Mapnik渲染得到的。

GRASS GIS

在地理信息系统行业中，如果说ArcGIS（操作系统行当中的windows）是闭源软件中乔峰，GRASS GIS就是开源软件中的慕容复。GRASS可处理矢量、影像数据，并且进行时空分析、空间建模、空间分析、地图可视化等操作。GRASS中的某些影像处理算法，无论从性能还是效果，完全可以媲美专业商业软件。虽然GRASS提供了一些可视化操作界面，但是大部分功能需用命令方式进行交互，对于普通用户来说，使用时有一定困难，比较适合科研人员或高校教师。

JTS

Java版本的GEOS，请参考GEOS部分。在路网操作中，会遇到将首尾相连的多条道路合并成一条道路的情况，使用JTS中的LineMerger类，可以很好的完成这个操作，示例代码：

WKTReader reader = new WKTReader();
Geometry geom_1 = reader.read("LINESTRING (116.96832000000000562 36.64882000000000062, 116.96849000000000274 36.64882000000000062)");
Geometry geom_2 = reader.read("LINESTRING (116.96849000000000274 36.64882000000000062, 116.96862000000000137 36.64882000000000062)");
Geometry geom_3 = reader.read("LINESTRING (116.96862000000000137 36.64882000000000062, 116.96877999999999531 36.64880999999999744)");
LineMerger lineMerger = new LineMerger();
//添加几何对象不需要按照顺序，只要道路首尾坐标点重合即可