C++ OGRMultiPolygon类代码示例

std::vector<Polygon> Polygon::polygons() const
{
    std::vector<Polygon> polys;

    OGRwkbGeometryType t = m_geom->getGeometryType();

    if (t == wkbPolygon || t == wkbPolygon25D)
        polys.emplace_back(*this);
    else if (t == wkbMultiPolygon || t == wkbMultiPolygon25D)
    {
        // Not until GDAL 2.3
        /**
        OGRMultiPolygon *mPoly = m_geom->toMultiPolygon();
        for (auto it = mPoly->begin(); it != mPoly->end(); ++it)
        {
            Polygon p;
            p.m_geom.reset((*it)->clone());
            polys.push_back(p);
        }
        **/
        OGRMultiPolygon *mPoly = static_cast<OGRMultiPolygon *>(m_geom.get());
        for (int i = 0; i < mPoly->getNumGeometries(); ++i)
        {
            Polygon p;
            p.m_geom.reset(mPoly->getGeometryRef(i)->clone());
            polys.push_back(p);
        }
    }
    return polys;
}

 OGRMultiPolygon* make()
 {
     OGRMultiPolygon* poCollection = new OGRMultiPolygon();
     
     poCollection->addGeometryDirectly(make<OGRPolygon>());
     
     return poCollection;
 }

 OGRFeature* create_area_feature(const shared_ptr<Osmium::OSM::Area const>& area)
 {
     OGRFeature* feature = OGRFeature::CreateFeature(m_layer_polygon->GetLayerDefn());
     Osmium::Geometry::MultiPolygon mp(*area);
     OGRMultiPolygon* ogrgeom = Osmium::Geometry::create_ogr_geometry(mp);
     ogrgeom->transform(m_transformation);
     feature->SetGeometryDirectly(ogrgeom);
     sprintf(longint, "%ld", area->from_way() ? area->orig_id() : -area->orig_id());
     feature->SetField("osm_id", longint);
     feature->SetField("z_order", calculate_z_order(area.get()));
     feature->SetField("way_area", ogrgeom->get_Area());
     return feature;
 }

OGRGeometry *OGRGeometryFactory::forceToMultiPolygon( OGRGeometry *poGeom )

{
    if( poGeom == NULL )
        return NULL;

/* -------------------------------------------------------------------- */
/*      Check for the case of a geometrycollection that can be          */
/*      promoted to MultiPolygon.                                       */
/* -------------------------------------------------------------------- */
    if( wkbFlatten(poGeom->getGeometryType()) == wkbGeometryCollection )
    {
        int iGeom;
        int bAllPoly = TRUE;
        OGRGeometryCollection *poGC = (OGRGeometryCollection *) poGeom;

        for( iGeom = 0; iGeom < poGC->getNumGeometries(); iGeom++ )
        {
            if( wkbFlatten(poGC->getGeometryRef(iGeom)->getGeometryType())
                != wkbPolygon )
                bAllPoly = FALSE;
        }

        if( !bAllPoly )
            return poGeom;
        
        OGRMultiPolygon *poMP = new OGRMultiPolygon();

        while( poGC->getNumGeometries() > 0 )
        {
            poMP->addGeometryDirectly( poGC->getGeometryRef(0) );
            poGC->removeGeometry( 0, FALSE );
        }

        delete poGC;

        return poMP;
    }

/* -------------------------------------------------------------------- */
/*      Eventually we should try to split the polygon into component    */
/*      island polygons.  But thats alot of work and can be put off.    */
/* -------------------------------------------------------------------- */
    if( wkbFlatten(poGeom->getGeometryType()) != wkbPolygon )
        return poGeom;

    OGRMultiPolygon *poMP = new OGRMultiPolygon();
    poMP->addGeometryDirectly( poGeom );

    return poMP;
}

void GDALMergeFaces::run()
{
	OGRGeometry * geo;
	leadingView.resetReading();

	OGRFeature * f;
	std::set<int> indizes;
	while (f = leadingView.getNextFeature()) {
		indizes.insert(f->GetFieldAsInteger(this->attriubteName.c_str()));
	}

	//int cluster_id = 1;
	int counter = 1;


	foreach (int cluster_id, indizes) {
	//while(geo = joinCluster(cluster_id)) {
		//cluster_id++;
		geo = joinCluster(cluster_id);

		if (counter % 100 == 0) {
		DM::Logger(DM::Standard) << "merged " << counter << "/" << indizes.size();
		}
		counter++;
		if (!geo)
			continue;



		if (wkbMultiPolygon == geo->getGeometryType()){
			geo = geo->UnionCascaded();
			OGRMultiPolygon * mgeo = (OGRMultiPolygon*) geo;
			if (mgeo->getNumGeometries() == 0) {
				continue;
			}
			geo = mgeo->getGeometryRef(0);

			int n = mgeo->getNumGeometries();
			for (int i = 0; i < n; i++) {
				OGRFeature * f = combinedView.createFeature();

				f->SetGeometry(mgeo->getGeometryRef(i));
				f->SetField("test_id", counter);
			}
			continue;
		}
		OGRFeature * f = combinedView.createFeature();
		f->SetGeometry(geo);
	}

int main() {
    // Read in raster data for night time lights
    int band_number = 1; // only one band, starts with one
    Raster* raster = import_raster("raster.tif", band_number);

    // Read in shapefile data containing municipality administrative regions
    int layer_number = 0; // only one layer, starts with zero
    OGRLayer* shapelayer = import_shapefile("MEX_adm2.shp", layer_number);

    shapelayer->SetAttributeFilter("ID_1 = 1834");       // Filter for Yucatan
    const int idx_of_number_field = 5;             // Column number of municipality number
    const int idx_of_name_field = 6;               // Column number of municipality name

    OGRFeature* poFeature;
    int feature_ctr = 0;
    while( (poFeature = shapelayer->GetNextFeature()) != NULL ) {
        cerr << "Feature: " << feature_ctr++ << "\t";
        int feature_num = poFeature->GetFieldAsInteger(idx_of_number_field);
        string feature_name = poFeature->GetFieldAsString(idx_of_name_field);

        OGRFeatureDefn *poFDefn = shapelayer->GetLayerDefn();
        for( int iField = 0; iField < poFDefn->GetFieldCount(); iField++ ) {
            OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
            if( poFieldDefn->GetType() == OFTString )  cerr << poFeature->GetFieldAsString(iField) << ",";
        }

        OGRGeometry* poGeometry = poFeature->GetGeometryRef();
        if( poGeometry != NULL) {
            // For contiguous regions
            if ( wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon ) {
                cerr << " polygon" << endl;
                report_raster_data_within_polygon(raster, (OGRPolygon *) poGeometry, feature_num, feature_name);
            // For disjoint regions
            } else if ( wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon ) {
                cerr << " multipolygon" << endl;
                OGRMultiPolygon *multipolygon = (OGRMultiPolygon *) poGeometry;
                for (int i = 0; i<multipolygon->getNumGeometries(); i++) {
                    report_raster_data_within_polygon(raster, (OGRPolygon*) multipolygon->getGeometryRef(i), feature_num, feature_name);
                }
            // Is this really the right shapefile?
            } else {
                cerr << "No polygon or multipolygon geometry for this feature: " << poGeometry->getGeometryName() << endl;
            }
        } else {
            cerr << "No geometry for this feature" << endl;
        }
    }
    OGRFeature::DestroyFeature( poFeature );
}

OGRGeometry *OGRMultiPolygon::clone() const

{
    OGRMultiPolygon     *poNewGC;

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

    for( int i = 0; i < getNumGeometries(); i++ )
    {
        poNewGC->addGeometry( getGeometryRef(i) );
    }

    return poNewGC;
}

void Polygon::simplify(double distance_tolerance, double area_tolerance)
{
    throwNoGeos();

    auto deleteSmallRings = [area_tolerance](OGRGeometry *geom)
    {
// Missing until GDAL 2.3.
//        OGRPolygon *poly = geom->toPolygon();
        OGRPolygon *poly = static_cast<OGRPolygon *>(geom);

        std::vector<int> deleteRings;
        for (int i = 0; i < poly->getNumInteriorRings(); ++i)
        {
            OGRLinearRing *lr = poly->getInteriorRing(i);
            if (lr->get_Area() < area_tolerance)
                deleteRings.push_back(i + 1);
        }
        // Note that interior rings are in a list with the exterior ring,
        // which is why the ring numbers are offset by one when used in
        // this context (what a mess).
        for (auto i : deleteRings)
// Missing until 2.3
//            poly->removeRing(i, true);
            OGR_G_RemoveGeometry(gdal::toHandle(poly), i, true);
    };

    OGRGeometry *g = m_geom->SimplifyPreserveTopology(distance_tolerance);
    m_geom.reset(g);

    OGRwkbGeometryType t = m_geom->getGeometryType();
    if (t == wkbPolygon || t == wkbPolygon25D)
        deleteSmallRings(m_geom.get());
    else if (t == wkbMultiPolygon || t == wkbMultiPolygon25D)
    {
// Missing until 2.3
/**
        OGRMultiPolygon *mpoly = m_geom->toMultiPolygon();
        for (auto it = mpoly->begin(); it != mpoly->end(); ++it)
            deleteSmallRings(*it);
**/
        OGRMultiPolygon *mpoly = static_cast<OGRMultiPolygon *>(m_geom.get());
        for (int i = 0; i < mpoly->getNumGeometries(); ++i)
            deleteSmallRings(mpoly->getGeometryRef(i));
    }
}

//.........这里部分代码省略.........
                    Position pos((SUMOReal) cgeom->getX(j), (SUMOReal) cgeom->getY(j));
                    if (!geoConvHelper.x2cartesian(pos)) {
                        WRITE_ERROR("Unable to project coordinates for polygon '" + id + "'.");
                    }
                    shape.push_back_noDoublePos(pos);
                }
                Polygon* poly = new Polygon(id, type, color, shape, false, (SUMOReal)layer);
                if (!toFill.insert(id, poly, layer)) {
                    WRITE_ERROR("Polygon '" + id + "' could not be added.");
                    delete poly;
                }
            }
            break;
            case wkbPolygon: {
                OGRLinearRing* cgeom = ((OGRPolygon*) poGeometry)->getExteriorRing();
                PositionVector shape;
                for (int j = 0; j < cgeom->getNumPoints(); j++) {
                    Position pos((SUMOReal) cgeom->getX(j), (SUMOReal) cgeom->getY(j));
                    if (!geoConvHelper.x2cartesian(pos)) {
                        WRITE_ERROR("Unable to project coordinates for polygon '" + id + "'.");
                    }
                    shape.push_back_noDoublePos(pos);
                }
                Polygon* poly = new Polygon(id, type, color, shape, true, (SUMOReal)layer);
                if (!toFill.insert(id, poly, layer)) {
                    WRITE_ERROR("Polygon '" + id + "' could not be added.");
                    delete poly;
                }
            }
            break;
            case wkbMultiPoint: {
                OGRMultiPoint* cgeom = (OGRMultiPoint*) poGeometry;
                for (int i = 0; i < cgeom->getNumGeometries(); ++i) {
                    OGRPoint* cgeom2 = (OGRPoint*) cgeom->getGeometryRef(i);
                    Position pos((SUMOReal) cgeom2->getX(), (SUMOReal) cgeom2->getY());
                    std::string tid = id + "#" + toString(i);
                    if (!geoConvHelper.x2cartesian(pos)) {
                        WRITE_ERROR("Unable to project coordinates for POI '" + tid + "'.");
                    }
                    PointOfInterest* poi = new PointOfInterest(tid, type, color, pos, (SUMOReal)layer);
                    if (!toFill.insert(tid, poi, layer)) {
                        WRITE_ERROR("POI '" + tid + "' could not be added.");
                        delete poi;
                    }
                }
            }
            break;
            case wkbMultiLineString: {
                OGRMultiLineString* cgeom = (OGRMultiLineString*) poGeometry;
                for (int i = 0; i < cgeom->getNumGeometries(); ++i) {
                    OGRLineString* cgeom2 = (OGRLineString*) cgeom->getGeometryRef(i);
                    PositionVector shape;
                    std::string tid = id + "#" + toString(i);
                    for (int j = 0; j < cgeom2->getNumPoints(); j++) {
                        Position pos((SUMOReal) cgeom2->getX(j), (SUMOReal) cgeom2->getY(j));
                        if (!geoConvHelper.x2cartesian(pos)) {
                            WRITE_ERROR("Unable to project coordinates for polygon '" + tid + "'.");
                        }
                        shape.push_back_noDoublePos(pos);
                    }
                    Polygon* poly = new Polygon(tid, type, color, shape, false, (SUMOReal)layer);
                    if (!toFill.insert(tid, poly, layer)) {
                        WRITE_ERROR("Polygon '" + tid + "' could not be added.");
                        delete poly;
                    }
                }
            }
            break;
            case wkbMultiPolygon: {
                OGRMultiPolygon* cgeom = (OGRMultiPolygon*) poGeometry;
                for (int i = 0; i < cgeom->getNumGeometries(); ++i) {
                    OGRLinearRing* cgeom2 = ((OGRPolygon*) cgeom->getGeometryRef(i))->getExteriorRing();
                    PositionVector shape;
                    std::string tid = id + "#" + toString(i);
                    for (int j = 0; j < cgeom2->getNumPoints(); j++) {
                        Position pos((SUMOReal) cgeom2->getX(j), (SUMOReal) cgeom2->getY(j));
                        if (!geoConvHelper.x2cartesian(pos)) {
                            WRITE_ERROR("Unable to project coordinates for polygon '" + tid + "'.");
                        }
                        shape.push_back_noDoublePos(pos);
                    }
                    Polygon* poly = new Polygon(tid, type, color, shape, true, (SUMOReal)layer);
                    if (!toFill.insert(tid, poly, layer)) {
                        WRITE_ERROR("Polygon '" + tid + "' could not be added.");
                        delete poly;
                    }
                }
            }
            break;
            default:
                WRITE_WARNING("Unsupported shape type occured (id='" + id + "').");
                break;
        }
        OGRFeature::DestroyFeature(poFeature);
    }
    PROGRESS_DONE_MESSAGE();
#else
    WRITE_ERROR("SUMO was compiled without GDAL support.");
#endif
}

static CPLErr OGRPolygonContourWriter( double dfLevelMin, double dfLevelMax,
                                const OGRMultiPolygon& multipoly,
                                void *pInfo )

{
    OGRContourWriterInfo *poInfo = static_cast<OGRContourWriterInfo *>(pInfo);

    OGRFeatureDefnH hFDefn =
        OGR_L_GetLayerDefn( static_cast<OGRLayerH>(poInfo->hLayer) );

    OGRFeatureH hFeat = OGR_F_Create( hFDefn );

    if( poInfo->nIDField != -1 )
        OGR_F_SetFieldInteger( hFeat, poInfo->nIDField, poInfo->nNextID++ );

    if( poInfo->nElevFieldMin != -1 )
        OGR_F_SetFieldDouble( hFeat, poInfo->nElevFieldMin, dfLevelMin );

    if( poInfo->nElevFieldMax != -1 )
        OGR_F_SetFieldDouble( hFeat, poInfo->nElevFieldMax, dfLevelMax );

    const bool bHasZ = wkbHasZ(OGR_FD_GetGeomType(hFDefn));
    OGRGeometryH hGeom = OGR_G_CreateGeometry(
        bHasZ ? wkbMultiPolygon25D : wkbMultiPolygon );

    for ( int iPart = 0; iPart < multipoly.getNumGeometries(); iPart++ )
    {
        OGRPolygon* poNewPoly = new OGRPolygon();
        const OGRPolygon* poPolygon = static_cast<const OGRPolygon*>(multipoly.getGeometryRef(iPart));

        for ( int iRing = 0; iRing < poPolygon->getNumInteriorRings() + 1; iRing++ )
        {
            const OGRLinearRing* poRing = iRing == 0 ?
                poPolygon->getExteriorRing()
                : poPolygon->getInteriorRing(iRing - 1);

            OGRLinearRing* poNewRing = new OGRLinearRing();
            for ( int iPoint = 0; iPoint < poRing->getNumPoints(); iPoint++ )
            {
                const double dfX = poInfo->adfGeoTransform[0]
                    + poInfo->adfGeoTransform[1] * poRing->getX(iPoint)
                    + poInfo->adfGeoTransform[2] * poRing->getY(iPoint);
                const double dfY = poInfo->adfGeoTransform[3]
                    + poInfo->adfGeoTransform[4] * poRing->getX(iPoint)
                    + poInfo->adfGeoTransform[5] * poRing->getY(iPoint);
                if( bHasZ )
                    OGR_G_SetPoint( OGRGeometry::ToHandle( poNewRing ), iPoint, dfX, dfY, dfLevelMax );
                else
                    OGR_G_SetPoint_2D( OGRGeometry::ToHandle( poNewRing ), iPoint, dfX, dfY );
            }
            poNewPoly->addRingDirectly( poNewRing );
        }
        OGR_G_AddGeometryDirectly( hGeom, OGRGeometry::ToHandle( poNewPoly ) );
    }

    OGR_F_SetGeometryDirectly( hFeat, hGeom );

    const OGRErr eErr =
        OGR_L_CreateFeature(static_cast<OGRLayerH>(poInfo->hLayer), hFeat);
    OGR_F_Destroy( hFeat );

    return eErr == OGRERR_NONE ? CE_None : CE_Failure;
}

//.........这里部分代码省略.........

        if( !ParseGMLCoordinates( psNode, &oPoints ) )
            return NULL;

        if( oPoints.getNumPoints() < 2 )
            return NULL;

        OGRLinearRing *poBoxRing = new OGRLinearRing();
        OGRPolygon *poBoxPoly = new OGRPolygon();

        poBoxRing->setNumPoints( 5 );
        poBoxRing->setPoint( 
            0, oPoints.getX(0), oPoints.getY(0), oPoints.getZ(0) );
        poBoxRing->setPoint( 
            1, oPoints.getX(1), oPoints.getY(0), oPoints.getZ(0) );
        poBoxRing->setPoint( 
            2, oPoints.getX(1), oPoints.getY(1), oPoints.getZ(1) );
        poBoxRing->setPoint( 
            3, oPoints.getX(0), oPoints.getY(1), oPoints.getZ(0) );
        poBoxRing->setPoint( 
            4, oPoints.getX(0), oPoints.getY(0), oPoints.getZ(0) );

        poBoxPoly->addRingDirectly( poBoxRing );

        return poBoxPoly;
    }

/* -------------------------------------------------------------------- */
/*      MultiPolygon                                                    */
/* -------------------------------------------------------------------- */
    if( EQUAL(pszBaseGeometry,"MultiPolygon") )
    {
        CPLXMLNode *psChild;
        OGRMultiPolygon *poMPoly = new OGRMultiPolygon();

        // Find all inner rings 
        for( psChild = psNode->psChild; 
             psChild != NULL;
             psChild = psChild->psNext ) 
        {
            if( psChild->eType == CXT_Element
                && EQUAL(BareGMLElement(psChild->pszValue),"polygonMember") )
            {
                OGRPolygon *poPolygon;

                poPolygon = (OGRPolygon *) 
                    GML2OGRGeometry_XMLNode( psChild->psChild );

                if( poPolygon == NULL )
                {
                    delete poMPoly;
                    return NULL;
                }

                if( !EQUAL(poPolygon->getGeometryName(),"POLYGON") )
                {
                    CPLError( CE_Failure, CPLE_AppDefined, 
                              "Got %.500s geometry as polygonMember instead of MULTIPOLYGON.",
                              poPolygon->getGeometryName() );
                    delete poPolygon;
                    delete poMPoly;
                    return NULL;
                }

                poMPoly->addGeometryDirectly( poPolygon );
            }

//.........这里部分代码省略.........
					double x,y;
					x= poLinearRing->getX( currentPoint );
					y= poLinearRing->getY( currentPoint );
					if(!poTransform->Transform(1, &x, &y))
					{
						Log::Inst().Warning("(Warning) Failed to project vector map.");
						OGRDataset->DestroyDataSource(OGRDataset);
						return false;
					}
				
					// project and store the points
					geoVector->pointX[currentPoint] = x;
					geoVector->pointY[currentPoint] = y;

					if(x < minX) minX=x;
					if(y < minY) minY=y;
					if(x > maxX) maxX=x;
					if(y > maxY) maxY=y;
				}
				else 			
				{
					geoVector->pointX[currentPoint] = poLinearRing->getX( currentPoint );
					geoVector->pointY[currentPoint] = poLinearRing->getY( currentPoint );

				}
			}

			vectorMapController->GetVectorMapModel()->AddGeoVector( geoVector );			
		}

		// Handle MULTIPOLYGONS
		else if ( wkbFlatten( poGeometry->getGeometryType() ) == wkbMultiPolygon )
		{
			OGRMultiPolygon *poMultiPolygon = (OGRMultiPolygon *) poGeometry;

			for ( int currGeometry = 0; currGeometry < poMultiPolygon->getNumGeometries(); currGeometry++ )
			{
				GeoVector* geoVector = new GeoVector();

				// OGRPolygon http://www.gdal.org/ogr/classOGRPolygon.html
				OGRPolygon *poPolygon = ( OGRPolygon* )poMultiPolygon->getGeometryRef( currGeometry );

				// Retrieve the EXTERNAL ring of the multipolygon
				OGRLinearRing *poLinearRing = poPolygon->getExteriorRing();

				geoVector->SetGeometryType( wkbLinearRing );
				geoVector->SetNumberOfPoints( poLinearRing->getNumPoints() );

				for ( int currentPoint = 0; currentPoint < poLinearRing->getNumPoints(); currentPoint++ )
				{

					 if(needProjection)
					{
						double x,y;
						x= poLinearRing->getX( currentPoint );
						y= poLinearRing->getY( currentPoint );
						if(!poTransform->Transform(1, &x, &y))
						{
							Log::Inst().Warning("(Warning) Failed to project vector map.");
							OGRDataset->DestroyDataSource(OGRDataset);
							return false;
						}
				
					// project and store the points
						geoVector->pointX[currentPoint] = x;
						geoVector->pointY[currentPoint] = y;

//.........这里部分代码省略.........
                }
                WriteFeatureAttributes(fp, poFeature);
                VSIFPrintfL( fp, "%d", nBNAPoints);
                n = ring->getNumPoints();
                int nbPair = 0;
                for(i=0;i<n;i++)
                {
                    VSIFPrintfL( fp, "%s", ((nbPair % nbPairPerLine) == 0) ? partialEol : " ");
                    WriteCoord(fp, ring->getX(i), ring->getY(i));
                    nbPair++;
                }
                for(i=0;i<nInteriorRings;i++)
                {
                    ring = polygon->getInteriorRing(i);
                    n = ring->getNumPoints();
                    for(j=0;j<n;j++)
                    {
                        VSIFPrintfL( fp, "%s", ((nbPair % nbPairPerLine) == 0) ? partialEol : " ");
                        WriteCoord(fp, ring->getX(j), ring->getY(j));
                        nbPair++;
                    }
                    VSIFPrintfL( fp, "%s", ((nbPair % nbPairPerLine) == 0) ? partialEol : " ");
                    WriteCoord(fp, firstX, firstY);
                    nbPair++;
                }
                VSIFPrintfL( fp, "%s", eol);
            }
            break;
        }

        case wkbMultiPolygon:
        case wkbMultiPolygon25D:
        {
            OGRMultiPolygon* multipolygon = (OGRMultiPolygon*)poGeom;
            int N = multipolygon->getNumGeometries();
            int nBNAPoints = 0;
            double firstX = 0, firstY = 0; 
            for(i=0;i<N;i++)
            {
                OGRPolygon* polygon = (OGRPolygon*)multipolygon->getGeometryRef(i);
                OGRLinearRing* ring = polygon->getExteriorRing();
                if (ring == NULL)
                    continue;

                if (nBNAPoints)
                    nBNAPoints ++;
                else
                {
                    firstX = ring->getX(0);
                    firstY = ring->getY(0);
                }
                nBNAPoints += ring->getNumPoints();
                int nInteriorRings = polygon->getNumInteriorRings();
                for(j=0;j<nInteriorRings;j++)
                {
                    nBNAPoints += polygon->getInteriorRing(j)->getNumPoints() + 1;
                }
            }
            if (nBNAPoints <= 3)
            {
                CPLError( CE_Failure, CPLE_AppDefined, "Invalid geometry" );
                return OGRERR_FAILURE;
            }
            WriteFeatureAttributes(fp, poFeature);
            VSIFPrintfL( fp, "%d", nBNAPoints);
            int nbPair = 0;

//.........这里部分代码省略.........
					x/=scaleFactor; y/=scaleFactor;

					if(i==0)
						path = WPainterPath( WPointF(x, y));
					else
						path.lineTo( x , y);
				}

				painter.drawPath(path); 
				poGeometry->Centroid(centroid); 
			}

			else if( (poGeometry != NULL)  &&  wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon) 
			{
				OGRPolygon *poPoint = (OGRPolygon *) poGeometry;
				OGRLinearRing *extring = poPoint->getExteriorRing();

				int n = extring->getNumPoints();
				double x, y;
				for(int i=0;i<n;i++)
				{
					x = extring->getX(i);   y = extring->getY(i);
					x/=scaleFactor; y/=scaleFactor;
					if(i==0)
						path = WPainterPath( WPointF(x , y));

					else
						path.lineTo( x , y);
				}
				painter.drawPath(path);
				poGeometry->Centroid(centroid);

			}
			else if( (poGeometry != NULL)  &&  wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon) 
			{
				double x, y;
				OGRMultiPolygon *poPoint = (OGRMultiPolygon *) poGeometry;
				int p = poPoint->getNumGeometries();
				for(int k=0;k<p;k++) {
					OGRGeometry* geom = poPoint->getGeometryRef(k);
					OGRPolygon *poly = (OGRPolygon *) geom;
					OGRLinearRing  *ring = poly->getExteriorRing();

					for(int i=0;i<ring->getNumPoints();i++)
					{
						x = ring->getX(i);   y = ring->getY(i);
						x/=scaleFactor; y/=scaleFactor;
						if(i==0)
							path = WPainterPath( WPointF(x , y));
						else
							path.lineTo( x , y);

					}
					painter.drawPath(path); 
					poGeometry->Centroid(centroid);


				}

			}
			if(labelindex>0 && !centroid->IsEmpty()){
				LABELS l={centroid->getX(),centroid->getY(),label};
				label_list.push_back(l);
			}
		}
		painter.restore(); 


		//labelling the contents
if(increase_width<100 && increase_height<100){
		painter.setWindow(0.0, 0.0, (paintDevice->width()).value(),(paintDevice->height()).value());
		font= new WFont(WFont::SansSerif);
		font->setSize(WLength(10*labelpercentage));
		painter.setFont(*font);
		pen.setColor(labelcolor);
		painter.setPen(pen);
		std::vector<LABELS>::iterator the_iterator = label_list.begin();
		double x, y, minx=(xMin+(-x_pos_shift+increase_width/2)/100*gWidth* widthFactor),miny=(yMax+(y_pos_shift-increase_height/2)/100*gHeight*widthFactor);
		
		double multx=(paintDevice->width().value())/(gWidth* widthFactor*(100-increase_width)/100);
		double multy=(paintDevice->height().value())/(gHeight*widthFactor*(-1+increase_height/100));

		while( the_iterator != label_list.end() ) {
		
		x=((*the_iterator).x/scaleFactor-minx)*multx;
		y=((*the_iterator).y/scaleFactor-miny)*multy;
	
		painter.drawText(WRectF( x-(*the_iterator).label.size()*5*labelpercentage, y-5*labelpercentage, (*the_iterator).label.size()	*10*labelpercentage,10*labelpercentage),AlignCenter,(*the_iterator).label);
			++the_iterator;
		}
		pen.setColor(red);
		painter.setPen(pen);
		painter.setFont(*font);
		//painter.drawText(WRectF(paintDevice->width().value()-dfile.size()*10*labelpercentage,paintDevice->height().value()-10*labelpercentage*(paintDevice->height()).value(), dfile.size()*10*labelpercentage,10*labelpercentage ),AlignCenter,dfile); //this text is not seen in the picture when painted.
}
	}					



}

int main(int argc, char **argv)
{
  // Get data from ogr
  OGRRegisterAll();
  std::cout << "Opening: " << argv[1] << std::endl;

  OGRDataSource *shp = OGRSFDriverRegistrar::Open(argv[1], FALSE);
  IsValid(shp, "Error opening file.");

  std::cout << "Shape contains " << shp->GetLayerCount() << " layers." << std::endl;
  OGRLayer *layer = shp->GetLayerByName(argv[2]);
  IsValid(layer, "Couldn't grab layer");

  OGRSpatialReference *srcSRS = NULL;
  srcSRS = layer->GetSpatialRef();

  // Set up writing
  const char *kDriverName = "ESRI Shapefile";
  OGRSFDriver *shpDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(kDriverName);
  IsValid(shpDriver, "Couldn't grab the shapefile driver.");

  IsValid(argv[3], "Please provide a output shp.");
  std::cout << "Writing to: " << argv[3] << std::endl;
  OGRDataSource *shpOut = shpDriver->CreateDataSource(argv[3], NULL);
  IsValid(shpOut, "Couldn't open output file");

  OGRLayer *outLayer = shpOut->CreateLayer(layer->GetName(), srcSRS, wkbMultiLineString, NULL);
  IsValid(outLayer, "Couldn't create an output layer");

  // copy over the fields from the source file
  OGRFeatureDefn *source = layer->GetLayerDefn();
  for(int i=0; i < source->GetFieldCount(); i++){
    OGRFieldDefn *field = source->GetFieldDefn(i);
    if(outLayer->CreateField(field) != OGRERR_NONE) {
      std::cout << "Couldn't make layer" << std::endl; exit(1);
    };
  }

  // Loop through features and grab the hull and put it into CGAL then
  // skeletonize the points
  OGRFeature *feature;
  int count = 0;
  while((feature = layer->GetNextFeature()) != NULL)
  {
    OGRMultiPolygon *geometry = dynamic_cast<OGRMultiPolygon *>(OGRGeometryFactory::forceToMultiPolygon(feature->GetGeometryRef()));
    IsValid(geometry, "No geometry.");

    OGRFeature *outFeature = OGRFeature::CreateFeature(outLayer->GetLayerDefn());
    IsValid(outFeature, "Couldn't make a feature.");

    for(int i=0; i < source->GetFieldCount(); i++){
      OGRField *field = feature->GetRawFieldRef(i);
      outFeature->SetField(i, field);
    }

    OGRGeometry* line = NULL;
    for(int i=0; i < geometry->getNumGeometries(); i++){
      OGRGeometry* segment = BuildMultiLine(geometry->getGeometryRef(i));
      if(segment != NULL){
        if(line == NULL) { line = new OGRLineString; }
        OGRGeometry* tmp = line->Union(segment);
        if(tmp != NULL){
          delete line;
          line = tmp;
        }
        delete segment;
      }
    }
    outFeature->SetGeometry(line);
    if(outLayer->CreateFeature(outFeature) != OGRERR_NONE){
      std::cout << "Couldn't create feature." << std::endl;
      exit(1);
    }

    // clean up
    OGRFeature::DestroyFeature(outFeature);
    std::cout << std::endl << ++count << std::endl;
  }

  // cleanup
  OGRDataSource::DestroyDataSource(shp);
  OGRDataSource::DestroyDataSource(shpOut);
  return 0;
}