使用GDAL/GEOS求面特徵的並集

存在這樣一個示例的矢量文件，包含了兩個重疊的面特徵：

一個很常見的需求是求取這個矢量中全部面元素的並集，經過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打開顯示以下：