編譯命令：g++ main.cpp -lgdal

呼叫命令：./a.out  輸出shp名稱 操作選項

註釋：操作選項（1：多邊形A - 多邊形B，2：B - A，3：A和B的交集部分）

#include "ogrsf_frmts.h"
#include <iostream>
using namespace std;
int main(int argc, char* argv[])
{
    const char *pszDriverName = "ESRI Shapefile";
    char *pShpName = argv[1];
    GDALDriver *poDriver;
    GDALAllRegister();
    poDriver = GetGDALDriverManager()->GetDriverByName(pszDriverName );
    if( poDriver == NULL )
    {
        printf( "%s driver not available.\n", pszDriverName );
        exit( 1 );
    }
    GDALDataset *poDS;
    char cShpName[20];
    sprintf(cShpName, "%s.shp", pShpName);
    poDS = poDriver->Create( cShpName, 0, 0, 0, GDT_Unknown, NULL );
    if( poDS == NULL )
    {
        printf( "Creation of output file failed.\n" );
        exit( 1 );
    }
    OGRLayer *poLayer;
    poLayer = poDS->CreateLayer( pShpName, NULL, wkbPolygon, NULL );
    if( poLayer == NULL )
    {
        printf( "Layer creation failed.\n" );
        exit( 1 );
    }
    OGRFieldDefn oField( "Name", OFTString );
    oField.SetWidth(32);
    if( poLayer->CreateField( &oField ) != OGRERR_NONE )  
    {
        printf( "Creating Name field failed.\n" );
        exit( 1 );
    }
    char szName[33] = "hello geos";

    OGRFeature *poFeatureG1, *poFeatureG2, *poFeatureG3;
    poFeatureG3 = OGRFeature::CreateFeature(poLayer->GetLayerDefn());
    poFeatureG3->SetField("Name", szName);

    // Geometry1
    poFeatureG1 = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
    OGRPolygon* pPolygon1 = new OGRPolygon();
    OGRLinearRing* pLinearRing = new OGRLinearRing();
    pLinearRing->addPoint(0, 0);
    pLinearRing->addPoint(2, 0);
    pLinearRing->addPoint(2, 2);
    pLinearRing->addPoint(0, 2);
    pLinearRing->addPoint(0, 0);
    pPolygon1->addRing(pLinearRing);
        
    poFeatureG1->SetGeometry( pPolygon1 ); 
    OGRGeometry* baseGeo = poFeatureG1->GetGeometryRef();

    // Geometry2
    poFeatureG2 = OGRFeature::CreateFeature(poLayer->GetLayerDefn());
    OGRPolygon* pPolygon2 = new OGRPolygon();
    OGRLinearRing* pLinearRing2 = new OGRLinearRing();
    pLinearRing2->addPoint(1, 0);
    pLinearRing2->addPoint(3, 0);
    pLinearRing2->addPoint(3, 3);
    pLinearRing2->addPoint(1, 3);
    pLinearRing2->addPoint(1, 0);
    pPolygon2->addRing(pLinearRing2);
    poFeatureG2->SetGeometry(pPolygon2);
    OGRGeometry* otheGeo = poFeatureG2->GetGeometryRef();

    int choose = atoi(argv[2]);
    OGRGeometry* pRet = NULL;
    switch(choose)
    {
        case 1:
            pRet = baseGeo->Difference(otheGeo);   // baseGeo - otheGeo
            break;
        case 2:
            pRet = otheGeo->Difference(baseGeo);   // otheGeo - baseGeo
            break; 
        case 3:
            pRet = baseGeo->Intersection(otheGeo); // baseGeo 和 otheGeo 相交
            break;
        default:
            break;
    }
    poFeatureG3->SetGeometry(pRet);
    if( poLayer->CreateFeature( poFeatureG3 ) != OGRERR_NONE )
    {
        printf( "Failed to create feature in shapefile.\n" );
        exit( 1 );
    }
    OGRFeature::DestroyFeature( poFeatureG1 );
    OGRFeature::DestroyFeature( poFeatureG2 );
    OGRFeature::DestroyFeature( poFeatureG3 );
    GDALClose( poDS );
}