编译命令：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);// Geometry1poFeatureG1 = 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();// Geometry2poFeatureG2 = 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 - otheGeobreak;case 2:pRet = otheGeo->Difference(baseGeo);   // otheGeo - baseGeobreak; 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 );
}