用GDAL实现GIS矢量数据读写与空间分析

前的介绍相信读者已经对于GIS的矢量数据有了一个全面的认识尽管可能不是那么深入。但是没有关系我们可以在实际的开发过程中加深对其的认识。4.4.1 第三方开源库OGR/GDAL在第3章的时候我们小试牛刀通过第三方开源库PROJ/GDAL很简单的就实现了地理空间参考系统的相互转换。在这里同样如此我们可以通过GDAL的OGR组件轻松实现GIS矢量数据的基础开发。OGR曾经是一个独立的矢量IO库表示OpenGIS Simple Features Reference Implementation的意思但其实OGR不完全符合OGC的Simple Feature标准规范因此未被批准作为该规范的参考实现当然也非常接近规范了。从GDAL2.0开始GDAL和OGR组件被集成在一起。OGR/GDAL提供了非常强大的矢量读写能力支持市面上绝大多数矢量数据格式。在第3章我们介绍GDAL的时候就提到过GDAL是一个GIS数据抽象库。所谓数据抽象是指无论是哪一种具体的矢量数据格式GDAL都会将其抽象成数据集对象Dataset从而可以支持读取、写出和处理操作。常用的矢量数据格式如下DXF/DWG应用最广泛的几何图形数据格式在CAD领域内用的非常多但是缺点是缺少地理信息。ESRI ShapefileGIS中最常见的矢量数据格式几乎所有的商业和开源GIS软件都支持。除了几何信息之外还包含空间信息和属性信息。GeoJSON一种通过JSON来表述要素的矢量数据格式因而很容易被JavaScript解析和处理适用于Web端的轻量化应用。KMLGoogle提出的一种基于XML标准来描述地理空间信息的数据格式包括点、线、面、多边形和模型等并且已经被OGC认定为开放地理信息编码标准。在GDAL中对数据格式的支持被封装成驱动dirver包括上述格式在内的数据驱动要么内嵌在GDAL中要么通过另外的第三方库来支持。当然这些我们可以暂时不用关心直接使用本书Github主页代码仓库中的GDAL即可。4.4.2 矢量数据的读取、处理和写入了解了第三方库OGR/GDAL接下来我们就通过其实现关于GIS矢量数据的基础开发。一个很容易理解的常识是包含GIS矢量数据在内的任何数据的操作都离不开三个过程读取、处理和写入。因此这里我们结合前面空间坐标参考转换的实践实现一个简单的实例读取一个面的矢量然后对其进行空间参考的转换最后重新写出一个矢量。具体的实现代码如下#include ogrsf_frmts.h #include Eigen/Eigen #include iostream #include vector using namespace std; using namespace Eigen; using Point Vector3d; using Line vectorPoint; using Polygon vectorLine; vectorPolygon polygonData; OGRSpatialReference srcFileSpatialReference; OGRSpatialReference dstFileSpatialReference; void OgrRing2Line(OGRLinearRing *ogrLinearRing, Line line) { for (int i 0; i ogrLinearRing-getNumPoints(); i) { line.emplace_back(ogrLinearRing-getX(i), ogrLinearRing-getY(i), ogrLinearRing-getZ(i)); } } void OgrPolygon2Polygon(OGRPolygon *ogrPolygon, Polygon polygon) { //外环 Line line; OgrRing2Line(ogrPolygon-getExteriorRing(), line); polygon.push_back(line); //内环 for (int ri 0; ri ogrPolygon-getNumInteriorRings(); ri) { Line line; OgrRing2Line(ogrPolygon-getInteriorRing(ri), line); polygon.push_back(line); } } bool ReadShp() { string srcFile getenv(GISBasic); srcFile srcFile /../Data/Vector/multipolygons.shp; GDALDataset *poDS (GDALDataset *)GDALOpenEx(srcFile.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL); if (!poDS) { printf(无法读取该文件请检查数据是否存在问题); return false; } if (poDS-GetLayerCount() 1) { printf(该文件的层数小于1请检查数据是否存在问题); return false; } //原始数据空间参考 char *pszWKT nullptr; poDS-GetLayer(0)-GetSpatialRef()-exportToWkt(pszWKT); srcFileSpatialReference.importFromWkt(pszWKT); CPLFree(pszWKT); pszWKT nullptr; for (int li 0; li poDS-GetLayerCount(); li) { OGRLayer *poLayer poDS-GetLayer(li); //读取层 poLayer-ResetReading(); //输出字段名 OGRFeatureDefn *poFDefn poLayer-GetLayerDefn(); int n poFDefn-GetFieldCount(); for (int iField 0; iField n; iField) { OGRFieldDefn *OGRFieldDefn poFDefn-GetFieldDefn(iField); cout OGRFieldDefn-GetNameRef() endl; } //遍历特征 OGRFeature *poFeature nullptr; while ((poFeature poLayer-GetNextFeature()) ! nullptr) { OGRGeometry *geometry poFeature-GetGeometryRef(); OGRwkbGeometryType geometryType geometry-getGeometryType(); switch (geometryType) { case wkbPolygon: case wkbPolygonM: case wkbPolygonZM: { OGRPolygon *ogrPolygon dynamic_castOGRPolygon *(geometry); if (!ogrPolygon) { continue; } Polygon polygon; OgrPolygon2Polygon(ogrPolygon, polygon); polygonData.push_back(polygon); break; } case wkbMultiPolygon: case wkbMultiPolygonM: case wkbMultiPolygonZM: { OGRMultiPolygon *ogrMultiPolygon dynamic_castOGRMultiPolygon *(geometry); if (!ogrMultiPolygon) { continue; } for (int gi 0; gi ogrMultiPolygon-getNumGeometries(); gi) { OGRPolygon *ogrPolygon dynamic_castOGRPolygon *(ogrMultiPolygon-getGeometryRef(gi)); if (!ogrPolygon) { continue; } Polygon polygon; OgrPolygon2Polygon(ogrPolygon, polygon); polygonData.push_back(polygon); } break; } default: { printf(未处理的特征类型\n); break; } } //输出每个字段的值 for (int iField 0; iField n; iField) { cout poFeature-GetFieldAsString(iField) ; } cout endl; OGRFeature::DestroyFeature(poFeature); } } GDALClose(poDS); poDS nullptr; return true; } void Convert() { dstFileSpatialReference.importFromEPSG(3857); dstFileSpatialReference.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); srcFileSpatialReference.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); OGRCoordinateTransformation *src2DstTransformation OGRCreateCoordinateTransformation(srcFileSpatialReference, dstFileSpatialReference); if (!src2DstTransformation) { return; } for (auto polygon : polygonData) { for (auto line : polygon) { for (auto point : line) { src2DstTransformation-Transform(1, point.data(), point.data() 1, point.data() 2); } } } } bool CreateField(OGRLayer *poLayer) { // 字符串 OGRFieldDefn oField1(Type, OFTString); oField1.SetWidth(8); if (poLayer-CreateField(oField1) ! OGRERR_NONE) { printf(Creating Name field failed.\n); return false; } // 浮点数 OGRFieldDefn oField2(Area, OFTReal); oField2.SetPrecision(3); if (poLayer-CreateField(oField2) ! OGRERR_NONE) { printf(Creating Name field failed.\n); return false; } // 整型 OGRFieldDefn oField3(VertexCount, OFTInteger); if (poLayer-CreateField(oField3) ! OGRERR_NONE) { printf(Creating Name field failed.\n); return false; } return true; } bool WriteGeoJson() { string dstFile getenv(GISBasic); dstFile dstFile /../Data/Out.geoJson; //创建 GDALDriver *driver GetGDALDriverManager()-GetDriverByName(GeoJSON); if (!driver) { printf(Get Driver GeoJSON Error\n); return false; } GDALDataset *dataset driver-Create(dstFile.c_str(), 0, 0, 0, GDT_Unknown, NULL); OGRLayer *poLayer dataset-CreateLayer( FirstLayer, dstFileSpatialReference, wkbPolygon, NULL); if (!CreateField(poLayer)) { return false; } //创建特征 for (const auto polygon : polygonData) { OGRFeature ogrFeature(poLayer-GetLayerDefn()); OGRPolygon ogrPolygon; int vertexNum 0; for (const auto line : polygon) { OGRLinearRing ogrRing; for (const auto point : line) { ogrRing.addPoint(point.x(), point.y(), point.z()); vertexNum; } ogrPolygon.addRing(ogrRing); } ogrFeature.SetGeometry(ogrPolygon); ogrFeature.SetField(Type, Polygon); ogrFeature.SetField(Area, ogrPolygon.get_Area()); ogrFeature.SetField(VertexCount, vertexNum); if (poLayer-CreateFeature(ogrFeature) ! OGRERR_NONE) { printf(Failed to create feature.\n); return false; } } //释放 GDALClose(dataset); dataset nullptr; return true; } int main() { GDALAllRegister(); CPLSetConfigOption(GDAL_FILENAME_IS_UTF8, NO); //支持中文路径 CPLSetConfigOption(SHAPE_ENCODING, ); //解决中文乱码问题 if (!ReadShp()) { return 1; } Convert(); WriteGeoJson(); }在这里读取的矢量是一个Shapefile文件包含的都是多边形的要素空间参考为WGS84地理坐标系EPSG4326显示如下图4.14所示经过空间参考转换后空间参考为Web墨卡托投影坐标系EPSG3857且生成的是一个GeoJson文件。如下图4.15所示可以看到在空间坐标参考转换之后形状似乎“变窄”了一点这是地图投影变换的特性决定的参考第2章介绍的知识。这段代码的关键在于如下几点读取矢量时矢量数据被读取成数据集对象GDALDatasetGDALDataset管理图层对象OGRLayerOGRLayer管理特征类对象OGRFeatureOGRFeature则包含几何对象类OGRGeometry。OGRGeometry是一个抽象父类有意义的是其表达具体的几何对象的子类比如OGRPolygon代表多边形几何对象。OGRPolygon又是由1个外环和多个内环环线OGRLinearRing组成。矢量就是这样一层层抽象组合表达出来的。在进行数据处理也就是空间参考转换之前我们把读取的数据放到我们自己的数据容器中using Point Vector3d; using Line vectorPoint; using Polygon vectorLine; vectorPolygon polygonData;真实的GIS数据处理开发要求是千变万化的操作的往往不是读取组件这里指OGR定义的数据结构对象。我们可以将其读取在成自定义的数据结构中以方便我们进行更加通用的操作或者进行并行优化。在这里我们遍历了自定义的数据结构对象中的顶点逐点进行空间坐标参考的转换。进行其他的数据处理也是如此根据自己的需要操作自己的数据结构容器。写出矢量时将处理结束的自定义数据结构对象恢复成数据集对象GDALDataset读取矢量的逆操作。另外需要关注的是GDAL的数据驱动GDALDriver其是通过名称来创建对应格式的矢量文件例如这里传入到名称是GeoJson:GDALDriver *driver GetGDALDriverManager()-GetDriverByName(GeoJSON); if (!driver) { printf(Get Driver GeoJSON Error\n); return false; }矢量数据的属性表数据与数据库中的表结构比较类似每个字段具有自己的定义属性例如数据类型、长度、精度等。并且每一个特征都对应了一行记录。因而属性表的操作还是挺复杂的往往与业务操作相关联。这里仅仅只是输出显示了读取的属性表信息写入了自己创建的字段与值。读取和写入的属性表数据如下图4.16、图4.17所示4.4.3 空间拓扑关系的判断对矢量要素进行空间拓扑关系的判断也是GIS开发中常用的功能。在4.3.3节中介绍了OGC的Simple Feature标准规范中拓扑关系OGR/GDAL对其做了具体的实现。任何继承了OGRGeometry的几何对象类如OGRPoint、OGRPolygon都可以调用标准中定义的接口二元谓词判断该几何对象与另外一个几何对象的拓扑关系。在如下实例中判断了空间中某点与某多边形是否存在包含Contains关系#include ogrsf_frmts.h #include iostream using namespace std; int main() { OGRLinearRing linearRing; linearRing.addPoint(268.28, 784.75); linearRing.addPoint(153.98, 600.60); linearRing.addPoint(274.63, 336.02); linearRing.addPoint(623.88, 401.64); linearRing.addPoint(676.80, 634.47); linearRing.addPoint(530.75, 822.85); linearRing.closeRings(); OGRPolygon polygon; polygon.addRing(linearRing); cout 点A是否在多边形内; OGRPoint pointA(407.98, 579.43); cout polygon.Contains(pointA) endl; cout 点B是否在多边形内; OGRPoint pointB(678.92, 482.07); cout polygon.Contains(pointB) endl; }运行的结果如下所示点A是否在多边形内1 点B是否在多边形内0在4.3.3节中介绍过包含(Contains)与内含(Within)是一对相反的关系。在这里我们当然也可以调用OGRPoint的内含(Within)接口会返回与上述同样的结果。不止如此其他空间拓扑关系的接口也是这样调用读者可以根据自己的需求尝试一二这里只是抛砖引玉。