开源软件在地图数据处理中的应用
作为开源软件的受益者,在享受开源带来的技术便利同时,我们也积极拥抱开源,同时也回馈开源。城市交通指数(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来说明:
图5 3度分带图
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();//添加几何对象不需要按照顺序,只要道路首尾坐标点重合即可lineMerger.add(geom_1);lineMerger.add(geom_2);lineMerger.add(geom_3);Collection mergedLineStrings = lineMerger.getMergedLineStrings();System.out.println(mergedLineStrings.toString());打印结果:[LINESTRING (116.96832 36.64882, 116.96849 36.64882, 116.96862 36.64882, 116.96878 36.64881)]
- GeoTools
功能和GDAL类似,空间拓扑算法使用JTS来实现。
- GeoMesa
GeoMesa是一个开源的进行时空数据处理的工具包,可以支持大数据场景下的地理信息分析和分布式计算。能较好的兼容大数据处理框架,如HBase、Spark,公司大数据架构部提供了完整的GeoMesa解决方案,可参考内网相关信息
作者:张深圳
现在注册滴滴云,有机会可得30元无门槛滴滴出行券
新购云服务1月5折 3月4.5折 6月低至4折
滴滴云使者招募,推荐最高返佣50%