本文节选自作者新书《GIS基础原理与技术实践》第4章，系统讲解矢量数据开发基础，涵盖 OGR/GDAL 使用、坐标转换、拓扑判断等实战内容。
4.4 矢量数据开发基础
通过之前的介绍，相信读者已经对于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 Shapefile：GIS中最常见的矢量数据格式，几乎所有的商业和开源GIS软件都支持。除了几何信息之外，还包含空间信息和属性信息。
- GeoJSON：一种通过JSON来表述要素的矢量数据格式，因而很容易被JavaScript解析和处理，适用于Web端的轻量化应用。
- KML：Google提出的一种基于XML标准来描述地理空间信息的数据格式（包括点、线、面、多边形和模型等），并且已经被OGC认定为开放地理信息编码标准。
在GDAL中，对数据格式的支持被封装成驱动（dirver），包括上述格式在内的数据驱动要么内嵌在GDAL中，要么通过另外的第三方库来支持。当然这些我们可以暂时不用关心，直接使用本书Github主页代码仓库中的GDAL即可。
4.4.2 矢量数据的读取、处理和写入
了解了第三方库OGR/GDAL，接下来我们就通过其实现关于GIS矢量数据的基础开发。一个很容易理解的常识是，包含GIS矢量数据在内的任何数据的操作都离不开三个过程：读取、处理和写入。因此，这里我们结合前面空间坐标参考转换的实践，实现一个简单的实例：读取一个面的矢量，然后对其进行空间参考的转换，最后重新写出一个矢量。具体的实现代码如下：- #include #include #include #include using namespace std;using namespace Eigen;using Point = Vector3d;using Line = vector;using Polygon = vector;vector 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 GetNameRef() GetNextFeature()) != nullptr) { OGRGeometry *geometry = poFeature->GetGeometryRef(); OGRwkbGeometryType geometryType = geometry->getGeometryType(); switch (geometryType) { case wkbPolygon: case wkbPolygonM: case wkbPolygonZM: { OGRPolygon *ogrPolygon = dynamic_cast(geometry); if (!ogrPolygon) { continue; } Polygon polygon; OgrPolygon2Polygon(ogrPolygon, polygon); polygonData.push_back(polygon); break; } case wkbMultiPolygon: case wkbMultiPolygonM: case wkbMultiPolygonZM: { OGRMultiPolygon *ogrMultiPolygon = dynamic_cast(geometry); if (!ogrMultiPolygon) { continue; } for (int gi = 0; gi < ogrMultiPolygon->getNumGeometries(); gi++) { OGRPolygon *ogrPolygon = dynamic_cast(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 GetFieldAsString(iField) 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();}
经过空间参考转换后，空间参考为Web墨卡托投影坐标系（EPSG：3857），且生成的是一个GeoJson文件。如下图4.15所示，可以看到在空间坐标参考转换之后形状似乎“变窄”了一点，这是地图投影变换的特性决定的（参考第2章介绍的知识）。
这段代码的关键在于如下几点：
- 读取矢量时，矢量数据被读取成数据集对象GDALDataset，GDALDataset管理图层对象OGRLayer，OGRLayer管理特征类对象OGRFeature，OGRFeature则包含几何对象类OGRGeometry。OGRGeometry是一个抽象父类，有意义的是其表达具体的几何对象的子类，比如OGRPolygon代表多边形几何对象。OGRPolygon又是由1个外环和多个内环（环线：OGRLinearRing）组成。矢量就是这样一层层抽象组合表达出来的。
- 在进行数据处理也就是空间参考转换之前，我们把读取的数据放到我们自己的数据容器中：
- using Point = Vector3d;using Line = vector;using Polygon = vector;vector polygonData;
- 写出矢量时，将处理结束的自定义数据结构对象恢复成数据集对象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）关系：
[code]#include #include 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 |
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
千斤顶
照妖镜
