C++ OGRGeometry::exportToGEOS方法代码示例

static CPLErr
BlendMaskGenerator(
#ifndef HAVE_GEOS
                    CPL_UNUSED int nXOff, CPL_UNUSED int nYOff,
                    CPL_UNUSED int nXSize, CPL_UNUSED int nYSize,
                    CPL_UNUSED GByte *pabyPolyMask,
                    CPL_UNUSED float *pafValidityMask,
                    CPL_UNUSED OGRGeometryH hPolygon,
                    CPL_UNUSED double dfBlendDist
#else
                    int nXOff, int nYOff, int nXSize, int nYSize,
                    GByte *pabyPolyMask, float *pafValidityMask,
                    OGRGeometryH hPolygon, double dfBlendDist
#endif
)
{
#ifndef HAVE_GEOS
    CPLError( CE_Failure, CPLE_AppDefined,
              "Blend distance support not available without the GEOS library.");
    return CE_Failure;

#else /* HAVE_GEOS */

/* -------------------------------------------------------------------- */
/*      Convert the polygon into a collection of lines so that we       */
/*      measure distance from the edge even on the inside.              */
/* -------------------------------------------------------------------- */
    OGRGeometry *poLines
        = OGRGeometryFactory::forceToMultiLineString(
            ((OGRGeometry *) hPolygon)->clone() );

/* -------------------------------------------------------------------- */
/*      Prepare a clipping polygon a bit bigger than the area of        */
/*      interest in the hopes of simplifying the cutline down to        */
/*      stuff that will be relavent for this area of interest.          */
/* -------------------------------------------------------------------- */
    CPLString osClipRectWKT;

    osClipRectWKT.Printf( "POLYGON((%g %g,%g %g,%g %g,%g %g,%g %g))", 
                          nXOff - (dfBlendDist+1), 
                          nYOff - (dfBlendDist+1), 
                          nXOff + nXSize + (dfBlendDist+1), 
                          nYOff - (dfBlendDist+1), 
                          nXOff + nXSize + (dfBlendDist+1), 
                          nYOff + nYSize + (dfBlendDist+1), 
                          nXOff - (dfBlendDist+1), 
                          nYOff + nYSize + (dfBlendDist+1), 
                          nXOff - (dfBlendDist+1), 
                          nYOff - (dfBlendDist+1) );
    
    OGRPolygon *poClipRect = NULL;
    char *pszWKT = (char *) osClipRectWKT.c_str();
    
    OGRGeometryFactory::createFromWkt( &pszWKT, NULL, 
                                       (OGRGeometry**) (&poClipRect) );

    if( poClipRect )
    {

        /***** if it doesnt intersect the polym zero the mask and return *****/

        if ( ! ((OGRGeometry *) hPolygon)->Intersects( poClipRect ) )
        {
            
            memset( pafValidityMask, 0, sizeof(float) * nXSize * nYSize );

            delete poLines;
            delete poClipRect;

            return CE_None;
        }

        /***** if it doesnt intersect the line at all just return *****/
         
        else if ( ! ((OGRGeometry *) poLines)->Intersects( poClipRect ) )
        {
            delete poLines;
            delete poClipRect;

            return CE_None;
        }

        OGRGeometry *poClippedLines = 
            poLines->Intersection( poClipRect );
        delete poLines;
        poLines = poClippedLines;
        delete poClipRect;
    }

/* -------------------------------------------------------------------- */
/*      Convert our polygon into GEOS format, and compute an            */
/*      envelope to accelerate later distance operations.               */
/* -------------------------------------------------------------------- */
    OGREnvelope sEnvelope;
    int iXMin, iYMin, iXMax, iYMax;
    GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();
    GEOSGeom poGEOSPoly;

    poGEOSPoly = poLines->exportToGEOS(hGEOSCtxt);
    OGR_G_GetEnvelope( hPolygon, &sEnvelope );
//.........这里部分代码省略.........

void lasclip(std::string &outfile, std::string &shapefile,
		std::string &layername, std::vector<std::string> &files,
		std::set<int> &classes, bool quiet) {

	if (outfile.empty())
		g_argerr("An output file is required.");
	if (shapefile.empty())
		g_argerr("A shape file is required.");
	if (files.size() == 0)
		g_argerr("At least one input file is required.");
	if (classes.size() == 0)
		g_warn("No classes specified, matching all classes.");

	/* Attempt to open and load geometries from the shape file. */
	OGRRegisterAll();
	OGRLayer *layer;
	OGRFeature *feat;
	OGRGeometry *og;
	OGRwkbGeometryType type;
	gg::GeometryCollection *geomColl;
	gg::Geometry *geom;

	OGRDataSource *ds = OGRSFDriverRegistrar::Open(shapefile.c_str(), FALSE);
	if (ds == nullptr)
		g_runerr("Couldn't open shapefile.");
	if (layername.empty()) {
		layer = ds->GetLayer(0);
	} else {
		layer = ds->GetLayerByName(layername.c_str());
	}
	if (layer == nullptr)
		g_runerr("Couldn't get layer.");

	type = layer->GetGeomType();
	if (type != wkbPolygon && type != wkbMultiPolygon)
		g_runerr("Geometry must be polygon or multipolygon.");

	const GEOSContextHandle_t gctx = OGRGeometry::createGEOSContext();
	const gg::GeometryFactory *gf = gg::GeometryFactory::getDefaultInstance();
	const gg::CoordinateSequenceFactory *cf =
			gf->getCoordinateSequenceFactory();
	std::vector<gg::Geometry *> geoms;

	while ((feat = layer->GetNextFeature()) != NULL) {
		og = feat->GetGeometryRef();
		geom = (gg::Geometry *) og->exportToGEOS(gctx);
		geoms.push_back(geom);
	}

	GDALClose(ds);

	if (geoms.size() == 0)
		g_runerr("No geometries were found.");

	/* The geometry collection is used for checking whether a las file intersects
	 the region of interest. */
	geomColl = gf->createGeometryCollection(geoms);
	const gg::Envelope *env = geomColl->getEnvelopeInternal();
	Bounds cbounds(env->getMinX(), env->getMinY(), env->getMaxX(),
			env->getMaxY());

	/* Loop over files and figure out which ones are relevant. */
	liblas::ReaderFactory rf;
	liblas::Header *dsth = nullptr;
	std::vector<unsigned int> indices;

	for (unsigned int i = 0; i < files.size(); ++i) {

		std::ifstream in(files[i].c_str(), std::ios::in | std::ios::binary);
		liblas::Reader r = rf.CreateWithStream(in);
		liblas::Header h = r.GetHeader();

		if (i == 0)
			dsth = new liblas::Header(h);

		std::vector<gg::Coordinate> coords;
		coords.push_back(gg::Coordinate(h.GetMinX(), h.GetMinY()));
		coords.push_back(gg::Coordinate(h.GetMaxX(), h.GetMinY()));
		coords.push_back(gg::Coordinate(h.GetMaxX(), h.GetMaxY()));
		coords.push_back(gg::Coordinate(h.GetMinX(), h.GetMaxY()));
		coords.push_back(gg::Coordinate(h.GetMinX(), h.GetMinY()));

		gg::CoordinateSequence *cs = cf->create(&coords);
		gg::LinearRing *lr = gf->createLinearRing(cs);
		gg::Polygon *bounds = gf->createPolygon(lr, NULL);

		if (bounds->intersects(geomColl))
			indices.push_back(i);

		in.close();
	}

	if (indices.size() == 0)
		g_runerr("No files matched the given bounds.");

	std::ofstream out(outfile, std::ios::out | std::ios::binary);
	liblas::WriterFactory wf;
	liblas::Writer w(out, *dsth);
	liblas::Header::RecordsByReturnArray recs;
	int count = 0;
//.........这里部分代码省略.........

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());
    poAreaLineLayer = 0;
    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())
    {
        OGRGeometry* geomRef = feature->GetGeometryRef();
        if( !geomRef )
        {
            continue;
        }
        GEOSGeom point = (GEOSGeom)(geomRef->exportToGEOS());
        for (i = 0; i < polys->getNumGeometries(); i++ )
        {
            if (ahInGeoms[i] && GEOSWithin(point, ahInGeoms[i]))
            {
                OGRFeature* areaFeature = new OGRFeature(poFeatureDefn);
                areaFeature->SetFrom(feature);
                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
    poAreaLineLayer = 0;
    delete polys;
}

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;
}