使用GDAL/GEOS求面特征的并集
- 2019 年 10 月 10 日
- 筆記
存在这样一个示例的矢量文件，包含了两个重叠的面特征：
一个很常见的需求是求取这个矢量中所有面元素的并集，通过GDAL/GEOS很容易实现这个功能，具体代码如下：
#include <iostream> #include <gdal/ogrsf_frmts.h> using namespace std; bool WritePolygon(string filePath, OGRPolygon *pOgrMerged) { //创建 GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("ESRI Shapefile"); if (!driver) { printf("Get Driver ESRI Shapefile Error！n"); return false; } GDALDataset* dataset = driver->Create(filePath.c_str(), 0, 0, 0, GDT_Unknown, NULL); OGRLayer* poLayer = dataset->CreateLayer("houseType", NULL, wkbPolygon, NULL); //创建特征 { OGRFeature *poFeature = new OGRFeature(poLayer->GetLayerDefn()); poFeature->SetGeometry(pOgrMerged); if (poLayer->CreateFeature(poFeature) != OGRERR_NONE) { printf("Failed to create feature in shapefile.n"); return false; } } //释放 GDALClose(dataset); dataset = nullptr; //GDALDestroyDriverManager(); return true; } int main() { GDALAllRegister(); CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); //支持中文路径 CPLSetConfigOption("SHAPE_ENCODING", ""); //解决中文乱码问题 string filePath = "D:/Work/OSGWork/shpTest/test/src.shp"; GDALDataset *poDS = (GDALDataset*)GDALOpenEx(filePath.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL); if (!poDS) { printf("无法读取该文件，试检查格式是否正确！"); return 1; } if (poDS->GetLayerCount()<1) { printf("该文件的层数小于1，试检查格式是否正确！"); return 1; } OGRLayer *poLayer = poDS->GetLayer(0); //读取层 poLayer->ResetReading(); std::unique_ptr<OGRPolygon> pOgrMerged(new OGRPolygon()); OGRFeature *poFeature; while ((poFeature = poLayer->GetNextFeature()) != NULL) { // OGRGeometry *pGeo = poFeature->GetGeometryRef(); OGRwkbGeometryType pGeoType = pGeo->getGeometryType(); if (pGeoType == wkbPolygon) { OGRPolygon *pPolygon = (OGRPolygon*)pGeo; if (!pPolygon) { continue; } OGRPolygon* pTemp = static_cast<OGRPolygon*>(pOgrMerged->Union(pPolygon)); if (pTemp) { pOgrMerged.reset(pTemp); } } OGRFeature::DestroyFeature(poFeature); } GDALClose(poDS); poDS = nullptr; if (pOgrMerged && pOgrMerged->IsValid() && pOgrMerged->getExteriorRing()) { string path = "D:/Work/OSGWork/shpTest/test/dst.shp"; WritePolygon(path, pOgrMerged.get()); } return 0; }
在这段代码中，遍历了示例矢量文件中的每个面元素，求取了所有面元素的并集，得到最终一个面元素，并将这个面元素保存成新的矢量文件。这个矢量文件用ArcMap打开显示如下：
