C++ OGRGeometryCollection::addGeometry方法代码示例

int OGRDB2GeometryValidator::ValidatePolygon(OGRPolygon* poGeom)
{
    int i,j;
    OGRLinearRing* poRing = poGeom->getExteriorRing();
    OGRGeometry* poInteriorRing;

    if (poRing == NULL)
        return FALSE;

    OGRGeometryCollection* poGeometries = NULL;

    if (!ValidateLinearRing(poRing))
    {
        if (poGeom->getNumInteriorRings() > 0)
        {
            poGeometries = new OGRGeometryCollection();
            poGeometries->addGeometryDirectly(poValidGeometry);
        }
    }

    for (i = 0; i < poGeom->getNumInteriorRings(); i++)
    {
        poInteriorRing = poGeom->getInteriorRing(i);
        if (!ValidateLinearRing((OGRLinearRing*)poInteriorRing))
        {
            if (!poGeometries)
            {
                poGeometries = new OGRGeometryCollection();
                poGeometries->addGeometry(poRing);
                for (j = 0; j < i; j++)
                    poGeometries->addGeometry(poGeom->getInteriorRing(j));
            }

            poGeometries->addGeometry(poValidGeometry);
            continue;
        }

        if (poGeometries)
            poGeometries->addGeometry(poInteriorRing);
    }

    if (poGeometries)
    {
        if (poValidGeometry)
            delete poValidGeometry;

        poValidGeometry = poGeometries;
    }

    return (poValidGeometry == NULL);
}

OGRGeometry *OGRGeometryCollection::clone() const

{
    OGRGeometryCollection       *poNewGC;

    poNewGC = new OGRGeometryCollection;
    poNewGC->assignSpatialReference( getSpatialReference() );

    for( int i = 0; i < nGeomCount; i++ )
    {
        poNewGC->addGeometry( papoGeoms[i] );
    }

    return poNewGC;
}

// [[Rcpp::export]]
Rcpp::List CPL_gdal_linestring_sample(Rcpp::List sfc, Rcpp::List distLst) {
	if (sfc.size() != distLst.size())
		throw std::invalid_argument("sfc and dist should have equal length");
	std::vector<OGRGeometry *> g = ogr_from_sfc(sfc, NULL);
	std::vector<OGRGeometry *> out(g.size());
	for (size_t i = 0; i < g.size(); i++) {
		OGRGeometryCollection *gc = new OGRGeometryCollection;
		Rcpp::NumericVector dists = distLst[i];
		for (size_t j = 0; j < dists.size(); j++) {
			OGRPoint *poPoint  = new OGRPoint;
			((OGRLineString *) g[i])->Value(dists[j], poPoint);
			gc->addGeometry(poPoint);
		}
		out[i] = OGRGeometryFactory::forceToMultiPoint(gc);
	}
	Rcpp::List ret = sfc_from_ogr(out, true);
	ret.attr("crs") = sfc.attr("crs");
	return ret;
}

OGRGeometry *OGRGeometryCollection::clone() const

{
    OGRGeometryCollection *poNewGC =
        OGRGeometryFactory::createGeometry(getGeometryType())->
            toGeometryCollection();
    poNewGC->assignSpatialReference( getSpatialReference() );
    poNewGC->flags = flags;

    for( auto&& poSubGeom: *this )
    {
        if( poNewGC->addGeometry( poSubGeom ) != OGRERR_NONE )
        {
            delete poNewGC;
            return nullptr;
        }
    }

    return poNewGC;
}

OGRGeometry *OGRGeometryCollection::clone() const

{
    OGRGeometryCollection       *poNewGC;

    poNewGC = (OGRGeometryCollection*)
            OGRGeometryFactory::createGeometry(getGeometryType());
    if( poNewGC == NULL )
        return NULL;
    poNewGC->assignSpatialReference( getSpatialReference() );
    poNewGC->flags = flags;

    for( int i = 0; i < nGeomCount; i++ )
    {
        if( poNewGC->addGeometry( papoGeoms[i] ) != OGRERR_NONE )
        {
            delete poNewGC;
            return NULL;
        }
    }

    return poNewGC;
}

OGRGeometry* GML_BuildOGRGeometryFromList(const CPLXMLNode* const * papsGeometry,
                                          int bTryToMakeMultipolygons,
                                          int bInvertAxisOrderIfLatLong,
                                          const char* pszDefaultSRSName,
                                          int bConsiderEPSGAsURN,
                                          int bGetSecondaryGeometryOption,
                                          void* hCacheSRS,
                                          int bFaceHoleNegative)
{
    OGRGeometry* poGeom = NULL;
    int i;
    OGRGeometryCollection* poCollection = NULL;
    for(i=0;papsGeometry[i] != NULL;i++)
    {
        OGRGeometry* poSubGeom = GML2OGRGeometry_XMLNode( papsGeometry[i],
                                                          bGetSecondaryGeometryOption,
                                                          0, 0, FALSE, TRUE,
                                                          bFaceHoleNegative );
        if (poSubGeom)
        {
            if (poGeom == NULL)
                poGeom = poSubGeom;
            else
            {
                if (poCollection == NULL)
                {
                    if (bTryToMakeMultipolygons &&
                        wkbFlatten(poGeom->getGeometryType()) == wkbPolygon &&
                        wkbFlatten(poSubGeom->getGeometryType()) == wkbPolygon)
                    {
                        OGRGeometryCollection* poGeomColl = new OGRMultiPolygon();
                        poGeomColl->addGeometryDirectly(poGeom);
                        poGeomColl->addGeometryDirectly(poSubGeom);
                        poGeom = poGeomColl;
                    }
                    else if (bTryToMakeMultipolygons &&
                                wkbFlatten(poGeom->getGeometryType()) == wkbMultiPolygon &&
                                wkbFlatten(poSubGeom->getGeometryType()) == wkbPolygon)
                    {
                        OGRGeometryCollection* poGeomColl = (OGRGeometryCollection* )poGeom;
                        poGeomColl->addGeometryDirectly(poSubGeom);
                    }
                    else if (bTryToMakeMultipolygons &&
                                wkbFlatten(poGeom->getGeometryType()) == wkbMultiPolygon &&
                                wkbFlatten(poSubGeom->getGeometryType()) == wkbMultiPolygon)
                    {
                        OGRGeometryCollection* poGeomColl = (OGRGeometryCollection* )poGeom;
                        OGRGeometryCollection* poGeomColl2 = (OGRGeometryCollection* )poSubGeom;
                        int nCount = poGeomColl2->getNumGeometries();
                        int i;
                        for(i=0;i<nCount;i++)
                        {
                            poGeomColl->addGeometry(poGeomColl2->getGeometryRef(i));
                        }
                        delete poSubGeom;
                    }
                    else if (bTryToMakeMultipolygons &&
                                wkbFlatten(poGeom->getGeometryType()) == wkbMultiPolygon)
                    {
                        delete poGeom;
                        delete poSubGeom;
                        return GML_BuildOGRGeometryFromList(papsGeometry, FALSE,
                                                            bInvertAxisOrderIfLatLong,
                                                            pszDefaultSRSName,
                                                            bConsiderEPSGAsURN,
                                                            bGetSecondaryGeometryOption,
                                                            hCacheSRS);
                    }
                    else
                    {
                        poCollection = new OGRGeometryCollection();
                        poCollection->addGeometryDirectly(poGeom);
                        poGeom = poCollection;
                    }
                }
                if (poCollection != NULL)
                {
                    poCollection->addGeometryDirectly(poSubGeom);
                }
            }
        }
    }
    
    if( poGeom == NULL )
        return NULL;

    std::string osWork;
    const char* pszSRSName = GML_ExtractSrsNameFromGeometry(papsGeometry, osWork,
                                                            bConsiderEPSGAsURN);
    const char* pszNameLookup = pszSRSName;
    if( pszNameLookup == NULL )
        pszNameLookup = pszDefaultSRSName;

    if (pszNameLookup != NULL)
    {
        SRSCache* poSRSCache = (SRSCache*)hCacheSRS;
        SRSDesc& oSRSDesc = poSRSCache->Get(pszNameLookup);
        poGeom->assignSpatialReference(oSRSDesc.poSRS);
        if (oSRSDesc.bAxisInvert && bInvertAxisOrderIfLatLong)
            poGeom->swapXY();
//.........这里部分代码省略.........

void OGRILI1Layer::PolygonizeAreaLayer()
{
    if (poAreaLineLayer == 0) return;

    //add all lines from poAreaLineLayer to collection
    OGRGeometryCollection *gc = new OGRGeometryCollection();
    poAreaLineLayer->ResetReading();
    while (OGRFeature *feature = poAreaLineLayer->GetNextFeatureRef())
        gc->addGeometry(feature->GetGeometryRef());

    //polygonize lines
    CPLDebug( "OGR_ILI", "Polygonizing layer %s with %d multilines", poAreaLineLayer->GetLayerDefn()->GetName(), gc->getNumGeometries());
    OGRMultiPolygon* polys = Polygonize( gc , false);
    CPLDebug( "OGR_ILI", "Resulting polygons: %d", polys->getNumGeometries());
    if (polys->getNumGeometries() != poAreaReferenceLayer->GetFeatureCount())
    {
        CPLDebug( "OGR_ILI", "Feature count of layer %s: %d", poAreaReferenceLayer->GetLayerDefn()->GetName(), GetFeatureCount());
        CPLDebug( "OGR_ILI", "Polygonizing again with crossing line fix");
        delete polys;
        polys = Polygonize( gc, true ); //try again with crossing line fix
    }
    delete gc;

    //associate polygon feature with data row according to centroid
#if defined(HAVE_GEOS)
    int i;
    OGRPolygon emptyPoly;
    GEOSGeom *ahInGeoms = NULL;

    CPLDebug( "OGR_ILI", "Associating layer %s with area polygons", GetLayerDefn()->GetName());
    ahInGeoms = (GEOSGeom *) CPLCalloc(sizeof(void*),polys->getNumGeometries());
    for( i = 0; i < polys->getNumGeometries(); i++ )
    {
        ahInGeoms[i] = polys->getGeometryRef(i)->exportToGEOS();
        if (!GEOSisValid(ahInGeoms[i])) ahInGeoms[i] = NULL;
    }
    poAreaReferenceLayer->ResetReading();
    while (OGRFeature *feature = poAreaReferenceLayer->GetNextFeatureRef())
    {
        GEOSGeom point = (GEOSGeom)feature->GetGeometryRef()->exportToGEOS();
        for (i = 0; i < polys->getNumGeometries(); i++ )
        {
            if (ahInGeoms[i] && GEOSWithin(point, ahInGeoms[i]))
            {
                OGRFeature* areaFeature = feature->Clone();
                areaFeature->SetGeometry( polys->getGeometryRef(i) );
                AddFeature(areaFeature);
                break;
            }
        }
        if (i == polys->getNumGeometries())
        {
            CPLDebug( "OGR_ILI", "Association between area and point failed.");
            feature->SetGeometry( &emptyPoly );
        }
        GEOSGeom_destroy( point );
    }
    for( i = 0; i < polys->getNumGeometries(); i++ )
        GEOSGeom_destroy( ahInGeoms[i] );
    CPLFree( ahInGeoms );
#endif
    poAreaReferenceLayer = 0;
    poAreaLineLayer = 0;
}

OGRGeometry* GML_BuildOGRGeometryFromList(const CPLXMLNode* const * papsGeometry,
                                          int bTryToMakeMultipolygons,
                                          int bInvertAxisOrderIfLatLong,
                                          const char* pszDefaultSRSName,
                                          int bConsiderEPSGAsURN,
                                          int bGetSecondaryGeometryOption,
                                          void* hCacheSRS)
{
    OGRGeometry* poGeom = NULL;
    int i;
    OGRGeometryCollection* poCollection = NULL;
    for(i=0;papsGeometry[i] != NULL;i++)
    {
        OGRGeometry* poSubGeom = GML2OGRGeometry_XMLNode( papsGeometry[i],
                                                          bGetSecondaryGeometryOption );
        if (poSubGeom)
        {
            if (poGeom == NULL)
                poGeom = poSubGeom;
            else
            {
                if (poCollection == NULL)
                {
                    if (bTryToMakeMultipolygons &&
                        wkbFlatten(poGeom->getGeometryType()) == wkbPolygon &&
                        wkbFlatten(poSubGeom->getGeometryType()) == wkbPolygon)
                    {
                        OGRGeometryCollection* poGeomColl = new OGRMultiPolygon();
                        poGeomColl->addGeometryDirectly(poGeom);
                        poGeomColl->addGeometryDirectly(poSubGeom);
                        poGeom = poGeomColl;
                    }
                    else if (bTryToMakeMultipolygons &&
                                wkbFlatten(poGeom->getGeometryType()) == wkbMultiPolygon &&
                                wkbFlatten(poSubGeom->getGeometryType()) == wkbPolygon)
                    {
                        OGRGeometryCollection* poGeomColl = (OGRGeometryCollection* )poGeom;
                        poGeomColl->addGeometryDirectly(poSubGeom);
                    }
                    else if (bTryToMakeMultipolygons &&
                                wkbFlatten(poGeom->getGeometryType()) == wkbMultiPolygon &&
                                wkbFlatten(poSubGeom->getGeometryType()) == wkbMultiPolygon)
                    {
                        OGRGeometryCollection* poGeomColl = (OGRGeometryCollection* )poGeom;
                        OGRGeometryCollection* poGeomColl2 = (OGRGeometryCollection* )poSubGeom;
                        int nCount = poGeomColl2->getNumGeometries();
                        int i;
                        for(i=0;i<nCount;i++)
                        {
                            poGeomColl->addGeometry(poGeomColl2->getGeometryRef(i));
                        }
                        delete poSubGeom;
                    }
                    else if (bTryToMakeMultipolygons &&
                                wkbFlatten(poGeom->getGeometryType()) == wkbMultiPolygon)
                    {
                        delete poGeom;
                        delete poSubGeom;
                        return GML_BuildOGRGeometryFromList(papsGeometry, FALSE,
                                                            bInvertAxisOrderIfLatLong,
                                                            pszDefaultSRSName,
                                                            bConsiderEPSGAsURN,
                                                            bGetSecondaryGeometryOption,
                                                            hCacheSRS);
                    }
                    else
                    {
                        poCollection = new OGRGeometryCollection();
                        poCollection->addGeometryDirectly(poGeom);
                        poGeom = poCollection;
                    }
                }
                if (poCollection != NULL)
                {
                    poCollection->addGeometryDirectly(poSubGeom);
                }
            }
        }
    }

    if ( poGeom != NULL && bInvertAxisOrderIfLatLong )
    {
        std::string osWork;
        const char* pszSRSName = GML_ExtractSrsNameFromGeometry(papsGeometry, osWork,
                                                          bConsiderEPSGAsURN);
        const char* pszNameLookup = pszSRSName ? pszSRSName : pszDefaultSRSName;
        if (pszNameLookup != NULL)
        {
            SRSCache* poSRSCache = (SRSCache*)hCacheSRS;
            int bSwap;
            if (strcmp(poSRSCache->osLastSRSName.c_str(), pszNameLookup) == 0)
            {
                bSwap = poSRSCache->bAxisInvertLastSRSName;
            }
            else
            {
                bSwap = GML_IsSRSLatLongOrder(pszNameLookup);
                poSRSCache->osLastSRSName = pszNameLookup;
                poSRSCache->bAxisInvertLastSRSName= bSwap;
            }
//.........这里部分代码省略.........

int OGRAVCBinLayer::FormPolygonGeometry( OGRFeature *poFeature, 
                                         AVCPal *psPAL )

{
/* -------------------------------------------------------------------- */
/*      Try to find the corresponding ARC layer if not already          */
/*      recorded.                                                       */
/* -------------------------------------------------------------------- */
    if( poArcLayer == NULL )
    {
        int	i;

        for( i = 0; i < poDS->GetLayerCount(); i++ )
        {
            OGRAVCBinLayer *poLayer = (OGRAVCBinLayer *) poDS->GetLayer(i);

            if( poLayer->eSectionType == AVCFileARC )
                poArcLayer = poLayer;
        }

        if( poArcLayer == NULL )
            return FALSE;
    }

/* -------------------------------------------------------------------- */
/*	Read all the arcs related to this polygon, making a working	*/
/*	copy of them since the one returned by AVC is temporary.	*/
/* -------------------------------------------------------------------- */
    OGRGeometryCollection oArcs;
    int		iArc;

    for( iArc = 0; iArc < psPAL->numArcs; iArc++ )
    {
        OGRFeature *poArc;

        if( psPAL->pasArcs[iArc].nArcId == 0 )
            continue;

        // If the other side of the line is the same polygon then this
        // arc is a "bridge" arc and can be discarded.  If we don't discard
        // it, then we should double it as bridge arcs seem to only appear
        // once.  But by discarding it we ensure a multi-ring polygon will be
        // properly formed. 
        if( psPAL->pasArcs[iArc].nAdjPoly == psPAL->nPolyId )
            continue;

        poArc = poArcLayer->GetFeature( ABS(psPAL->pasArcs[iArc].nArcId) );

        if( poArc == NULL )
            return FALSE;

        if( poArc->GetGeometryRef() == NULL )
            return FALSE;

        oArcs.addGeometry( poArc->GetGeometryRef() );
        OGRFeature::DestroyFeature( poArc );
    }

    OGRErr  eErr;
    OGRPolygon *poPolygon;

    poPolygon = (OGRPolygon *) 
        OGRBuildPolygonFromEdges( (OGRGeometryH) &oArcs, TRUE, FALSE,  
                                  0.0, &eErr );
    if( poPolygon != NULL )
        poFeature->SetGeometryDirectly( poPolygon );

    return eErr == OGRERR_NONE;
}

static OGRGeometryCollection* LoadGeometry( const char* pszDS,
                                            const char* pszSQL,
                                            const char* pszLyr,
                                            const char* pszWhere )
{
    OGRDataSource       *poDS;
    OGRLayer            *poLyr;
    OGRFeature          *poFeat;
    OGRGeometryCollection *poGeom = NULL;
        
    poDS = OGRSFDriverRegistrar::Open( pszDS, FALSE );
    if ( poDS == NULL )
        return NULL;

    if ( pszSQL != NULL )
        poLyr = poDS->ExecuteSQL( pszSQL, NULL, NULL ); 
    else if ( pszLyr != NULL )
        poLyr = poDS->GetLayerByName( pszLyr );
    else
        poLyr = poDS->GetLayer(0);
        
    if ( poLyr == NULL )
    {
        fprintf( stderr,
            "FAILURE: Failed to identify source layer from datasource.\n" );
        OGRDataSource::DestroyDataSource( poDS );
        return NULL;
    }
    
    if ( pszWhere )
        poLyr->SetAttributeFilter( pszWhere );
        
    while ( (poFeat = poLyr->GetNextFeature()) != NULL )
    {
        OGRGeometry* poSrcGeom = poFeat->GetGeometryRef();
        if ( poSrcGeom )
        {
            OGRwkbGeometryType eType =
                wkbFlatten( poSrcGeom->getGeometryType() );
            
            if ( poGeom == NULL )
                poGeom = new OGRMultiPolygon();

            if ( eType == wkbPolygon )
                poGeom->addGeometry( poSrcGeom );
            else if ( eType == wkbMultiPolygon )
            {
                int iGeom;
                int nGeomCount =
                    ((OGRMultiPolygon *)poSrcGeom)->getNumGeometries();

                for ( iGeom = 0; iGeom < nGeomCount; iGeom++ )
                {
                    poGeom->addGeometry(
                        ((OGRMultiPolygon *)poSrcGeom)->getGeometryRef(iGeom) );
                }
            }
            else
            {
                fprintf( stderr, "FAILURE: Geometry not of polygon type.\n" );
                OGRGeometryFactory::destroyGeometry( poGeom );
                OGRFeature::DestroyFeature( poFeat );
                if ( pszSQL != NULL )
                    poDS->ReleaseResultSet( poLyr );
                OGRDataSource::DestroyDataSource( poDS );
                return NULL;
            }
        }
    
        OGRFeature::DestroyFeature( poFeat );
    }
    
    if( pszSQL != NULL )
        poDS->ReleaseResultSet( poLyr );
    OGRDataSource::DestroyDataSource( poDS );
    
    return poGeom;
}

void OGRILI1Layer::PolygonizeAreaLayer( OGRILI1Layer* poAreaLineLayer, int nAreaFieldIndex, int nPointFieldIndex )
{
    //add all lines from poAreaLineLayer to collection
    OGRGeometryCollection *gc = new OGRGeometryCollection();
    poAreaLineLayer->ResetReading();
    while (OGRFeature *feature = poAreaLineLayer->GetNextFeatureRef())
        gc->addGeometry(feature->GetGeometryRef());

    //polygonize lines
    CPLDebug( "OGR_ILI", "Polygonizing layer %s with %d multilines", poAreaLineLayer->GetLayerDefn()->GetName(), gc->getNumGeometries());
    poAreaLineLayer = 0;
    OGRMultiPolygon* polys = Polygonize( gc , false);
    CPLDebug( "OGR_ILI", "Resulting polygons: %d", polys->getNumGeometries());
    if (polys->getNumGeometries() != GetFeatureCount())
    {
        CPLDebug( "OGR_ILI", "Feature count of layer %s: " CPL_FRMT_GIB, GetLayerDefn()->GetName(), GetFeatureCount());
        CPLDebug( "OGR_ILI", "Polygonizing again with crossing line fix");
        delete polys;
        polys = Polygonize( gc, true ); //try again with crossing line fix
        CPLDebug( "OGR_ILI", "Resulting polygons: %d", polys->getNumGeometries());
    }
    delete gc;

    //associate polygon feature with data row according to centroid
#if defined(HAVE_GEOS)
    int i;
    OGRPolygon emptyPoly;
    GEOSGeom *ahInGeoms = NULL;

    CPLDebug( "OGR_ILI", "Associating layer %s with area polygons", GetLayerDefn()->GetName());
    ahInGeoms = (GEOSGeom *) CPLCalloc(sizeof(void*), polys->getNumGeometries());
    GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();
    for( i = 0; i < polys->getNumGeometries(); i++ )
    {
        ahInGeoms[i] = polys->getGeometryRef(i)->exportToGEOS(hGEOSCtxt);
        if (!GEOSisValid_r(hGEOSCtxt, ahInGeoms[i])) ahInGeoms[i] = NULL;
    }
    for ( int nFidx = 0; nFidx < nFeatures; nFidx++)
    {
        OGRFeature *feature = papoFeatures[nFidx];
        OGRGeometry* geomRef = feature->GetGeomFieldRef(nPointFieldIndex);
        if( !geomRef )
        {
            continue;
        }
        GEOSGeom point = (GEOSGeom)(geomRef->exportToGEOS(hGEOSCtxt));
        for (i = 0; i < polys->getNumGeometries(); i++ )
        {
            if (ahInGeoms[i] && GEOSWithin_r(hGEOSCtxt, point, ahInGeoms[i]))
            {
                feature->SetGeomField(nAreaFieldIndex, polys->getGeometryRef(i));
                break;
            }
        }
        if (i == polys->getNumGeometries())
        {
            CPLDebug( "OGR_ILI", "Association between area and point failed.");
            feature->SetGeometry( &emptyPoly );
        }
        GEOSGeom_destroy_r( hGEOSCtxt, point );
    }
    for( i = 0; i < polys->getNumGeometries(); i++ )
        GEOSGeom_destroy_r( hGEOSCtxt, ahInGeoms[i] );
    CPLFree( ahInGeoms );
    OGRGeometry::freeGEOSContext( hGEOSCtxt );
#endif
    poAreaLineLayer = 0;
    delete polys;
}