存在這樣一個示例的矢量檔案，包含了兩個重疊的面特征：
一個很常見的需求是求取這個矢量中所有面元素的并集，通過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打開顯示如下：
轉載請註明出處，本文鏈接：https://www.uj5u.com/qiye/3959.html
標籤：GIS
上一篇：使用GDAL/OGR讀寫矢量檔案