OGR示例:写shp，求面与面的交和差操作

来源：互联网发布：淘宝页头怎么设置编辑：程序博客网时间：2024/04/28 15:06
编译命令：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 );}
1 0