EVOLUTION-MANAGER
Edit File: ogrdissolve.cpp
/****************************************************************************** * * Project: OpenGIS Simple Features Reference Implementation * Purpose: Allow a user to dissolve geometries based on an attribute. * Author: Howard Butler, hobu.inc@gmail.com * ****************************************************************************** * Copyright (c) 2007, Howard Butler * * Permission is hereby granted, free of charge, to any person obtaining a * copy of this software and associated documentation files (the "Software"), * to deal in the Software without restriction, including without limitation * the rights to use, copy, modify, merge, publish, distribute, sublicense, * and/or sell copies of the Software, and to permit persons to whom the * Software is furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included * in all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER * DEALINGS IN THE SOFTWARE. ****************************************************************************/ #include "ogrsf_frmts.h" #include "ogr_p.h" #include "cpl_conv.h" #include "cpl_string.h" #include "ogr_api.h" #include <map> #include <list> CPL_CVSID("$Id: ogrdissolve.cpp 35931 2016-10-25 16:20:52Z goatbar $"); static void Usage(); static int DissolveLayer( OGRDataSource *poSrcDS, OGRLayer * poSrcLayer, OGRDataSource *poDstDS, char ** papszLSCO, const char *pszNewLayerName, int bTransform, OGRSpatialReference *poOutputSRS, OGRSpatialReference *poSourceSRS, char **papszSelFields, int bAppend, int eGType, int bOverwrite ); static int bSkipFailures = FALSE; static int nGroupTransactions = 200; static int bPreserveFID = FALSE; static int nFIDToFetch = OGRNullFID; typedef std::multimap<CPLString, OGRGeometry*> StringGeometryMMap; typedef std::map<CPLString, OGRGeometryCollection*> StringGeometryColMap; typedef std::map<CPLString, OGRGeometry*> StringGeometryMap; typedef std::list<OGRGeometry*> GeometriesList; /************************************************************************/ /* main() */ /************************************************************************/ int main( int nArgc, char ** papszArgv ) { const char *pszFormat = "ESRI Shapefile"; const char *pszDataSource = NULL; const char *pszDestDataSource = NULL; char **papszLayers = NULL; char **papszDSCO = NULL, **papszLCO = NULL; int bTransform = FALSE; int bAppend = FALSE, bUpdate = FALSE, bOverwrite = FALSE; const char *pszOutputSRSDef = NULL; const char *pszSourceSRSDef = NULL; OGRSpatialReference *poOutputSRS = NULL; OGRSpatialReference *poSourceSRS = NULL; const char *pszNewLayerName = NULL; const char *pszWHERE = NULL; OGRGeometry *poSpatialFilter = NULL; const char *pszSelect; char **papszSelFields = NULL; const char *pszSQLStatement = NULL; int eGType = -2; /* -------------------------------------------------------------------- */ /* Register format(s). */ /* -------------------------------------------------------------------- */ OGRRegisterAll(); /* -------------------------------------------------------------------- */ /* Processing command line arguments. */ /* -------------------------------------------------------------------- */ nArgc = OGRGeneralCmdLineProcessor( nArgc, &papszArgv, 0 ); if( nArgc < 1 ) exit( -nArgc ); for( int iArg = 1; iArg < nArgc; iArg++ ) { if( EQUAL(papszArgv[iArg],"-f") && iArg < nArgc-1 ) { pszFormat = papszArgv[++iArg]; } else if( EQUAL(papszArgv[iArg],"-dsco") && iArg < nArgc-1 ) { papszDSCO = CSLAddString(papszDSCO, papszArgv[++iArg] ); } else if( EQUAL(papszArgv[iArg],"-lco") && iArg < nArgc-1 ) { papszLCO = CSLAddString(papszLCO, papszArgv[++iArg] ); } else if( EQUAL(papszArgv[iArg],"-preserve_fid") ) { bPreserveFID = TRUE; } else if( STARTS_WITH_CI(papszArgv[iArg], "-skip") ) { bSkipFailures = TRUE; } else if( EQUAL(papszArgv[iArg],"-append") ) { bAppend = TRUE; } else if( EQUAL(papszArgv[iArg],"-overwrite") ) { bOverwrite = TRUE; } else if( EQUAL(papszArgv[iArg],"-update") ) { bUpdate = TRUE; } else if( EQUAL(papszArgv[iArg],"-fid") && papszArgv[iArg+1] != NULL ) { nFIDToFetch = atoi(papszArgv[++iArg]); } else if( EQUAL(papszArgv[iArg],"-sql") && papszArgv[iArg+1] != NULL ) { pszSQLStatement = papszArgv[++iArg]; } else if( EQUAL(papszArgv[iArg],"-nln") && iArg < nArgc-1 ) { pszNewLayerName = papszArgv[++iArg]; } else if( EQUAL(papszArgv[iArg],"-nlt") && iArg < nArgc-1 ) { int bIs3D = FALSE; CPLString osGeomName = papszArgv[iArg+1]; if (strlen(papszArgv[iArg+1]) > 3 && STARTS_WITH_CI(papszArgv[iArg+1] + strlen(papszArgv[iArg+1]) - 3, "25D")) { bIs3D = TRUE; osGeomName.resize(osGeomName.size() - 3); } else if (strlen(papszArgv[iArg+1]) > 1 && STARTS_WITH_CI(papszArgv[iArg+1] + strlen(papszArgv[iArg+1]) - 1, "Z")) { bIs3D = TRUE; osGeomName.resize(osGeomName.size() - 1); } if( EQUAL(osGeomName,"NONE") ) eGType = wkbNone; else if( EQUAL(osGeomName,"GEOMETRY") ) eGType = wkbUnknown; else { eGType = OGRFromOGCGeomType(osGeomName); if (eGType == wkbUnknown) { fprintf( stderr, "-nlt %s: type not recognised.\n", papszArgv[iArg+1] ); exit( 1 ); } } if (eGType != wkbNone && bIs3D) eGType = wkbSetZ((OGRwkbGeometryType)eGType); iArg++; } else if( EQUAL(papszArgv[iArg],"-tg") && iArg < nArgc-1 ) { nGroupTransactions = atoi(papszArgv[++iArg]); } else if( EQUAL(papszArgv[iArg],"-s_srs") && iArg < nArgc-1 ) { pszSourceSRSDef = papszArgv[++iArg]; } else if( EQUAL(papszArgv[iArg],"-a_srs") && iArg < nArgc-1 ) { pszOutputSRSDef = papszArgv[++iArg]; } else if( EQUAL(papszArgv[iArg],"-t_srs") && iArg < nArgc-1 ) { pszOutputSRSDef = papszArgv[++iArg]; bTransform = TRUE; } else if( EQUAL(papszArgv[iArg],"-spat") && papszArgv[iArg+1] != NULL && papszArgv[iArg+2] != NULL && papszArgv[iArg+3] != NULL && papszArgv[iArg+4] != NULL ) { OGRLinearRing oRing; oRing.addPoint( CPLAtof(papszArgv[iArg+1]), CPLAtof(papszArgv[iArg+2]) ); oRing.addPoint( CPLAtof(papszArgv[iArg+1]), CPLAtof(papszArgv[iArg+4]) ); oRing.addPoint( CPLAtof(papszArgv[iArg+3]), CPLAtof(papszArgv[iArg+4]) ); oRing.addPoint( CPLAtof(papszArgv[iArg+3]), CPLAtof(papszArgv[iArg+2]) ); oRing.addPoint( CPLAtof(papszArgv[iArg+1]), CPLAtof(papszArgv[iArg+2]) ); poSpatialFilter = new OGRPolygon(); ((OGRPolygon *) poSpatialFilter)->addRing( &oRing ); iArg += 4; } else if( EQUAL(papszArgv[iArg],"-where") && papszArgv[iArg+1] != NULL ) { pszWHERE = papszArgv[++iArg]; } else if( EQUAL(papszArgv[iArg],"-select") && papszArgv[iArg+1] != NULL) { pszSelect = papszArgv[++iArg]; papszSelFields = CSLTokenizeStringComplex(pszSelect, " ,", FALSE, FALSE ); } else if( papszArgv[iArg][0] == '-' ) { Usage(); } else if( pszDestDataSource == NULL ) pszDestDataSource = papszArgv[iArg]; else if( pszDataSource == NULL ) pszDataSource = papszArgv[iArg]; else papszLayers = CSLAddString( papszLayers, papszArgv[iArg] ); } if( pszDataSource == NULL ) Usage(); /* -------------------------------------------------------------------- */ /* Open data source. */ /* -------------------------------------------------------------------- */ OGRDataSource *poDS; poDS = OGRSFDriverRegistrar::Open( pszDataSource, FALSE ); /* -------------------------------------------------------------------- */ /* Report failure */ /* -------------------------------------------------------------------- */ if( poDS == NULL ) { OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar(); printf( "FAILURE:\n" "Unable to open datasource `%s' with the following drivers.\n", pszDataSource ); for( int iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ ) { printf( " -> %s\n", poR->GetDriver(iDriver)->GetName() ); } exit( 1 ); } /* -------------------------------------------------------------------- */ /* Try opening the output datasource as an existing, writable */ /* -------------------------------------------------------------------- */ OGRDataSource *poODS; if( bUpdate ) { poODS = OGRSFDriverRegistrar::Open( pszDestDataSource, TRUE ); if( poODS == NULL ) { printf( "FAILURE:\n" "Unable to open existing output datasource `%s'.\n", pszDestDataSource ); exit( 1 ); } } /* -------------------------------------------------------------------- */ /* Find the output driver. */ /* -------------------------------------------------------------------- */ else { OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar(); OGRSFDriver *poDriver = NULL; int iDriver; for( iDriver = 0; iDriver < poR->GetDriverCount() && poDriver == NULL; iDriver++ ) { if( EQUAL(poR->GetDriver(iDriver)->GetName(),pszFormat) ) { poDriver = poR->GetDriver(iDriver); } } if( poDriver == NULL ) { printf( "Unable to find driver `%s'.\n", pszFormat ); printf( "The following drivers are available:\n" ); for( iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ ) { printf( " -> `%s'\n", poR->GetDriver(iDriver)->GetName() ); } exit( 1 ); } if( !poDriver->TestCapability( ODrCCreateDataSource ) ) { printf( "%s driver does not support data source creation.\n", pszFormat ); exit( 1 ); } /* -------------------------------------------------------------------- */ /* Create the output data source. */ /* -------------------------------------------------------------------- */ poODS = poDriver->CreateDataSource( pszDestDataSource, papszDSCO ); if( poODS == NULL ) { printf( "%s driver failed to create %s\n", pszFormat, pszDestDataSource ); exit( 1 ); } } /* -------------------------------------------------------------------- */ /* Parse the output SRS definition if possible. */ /* -------------------------------------------------------------------- */ if( pszOutputSRSDef != NULL ) { poOutputSRS = new OGRSpatialReference(); if( poOutputSRS->SetFromUserInput( pszOutputSRSDef ) != OGRERR_NONE ) { printf( "Failed to process SRS definition: %s\n", pszOutputSRSDef ); exit( 1 ); } } /* -------------------------------------------------------------------- */ /* Parse the source SRS definition if possible. */ /* -------------------------------------------------------------------- */ if( pszSourceSRSDef != NULL ) { poSourceSRS = new OGRSpatialReference(); if( poSourceSRS->SetFromUserInput( pszSourceSRSDef ) != OGRERR_NONE ) { printf( "Failed to process SRS definition: %s\n", pszSourceSRSDef ); exit( 1 ); } } /* -------------------------------------------------------------------- */ /* Special case for -sql clause. No source layers required. */ /* -------------------------------------------------------------------- */ if( pszSQLStatement != NULL ) { OGRLayer *poResultSet; if( pszWHERE != NULL ) printf( "-where clause ignored in combination with -sql.\n" ); if( CSLCount(papszLayers) > 0 ) printf( "layer names ignored in combination with -sql.\n" ); poResultSet = poDS->ExecuteSQL( pszSQLStatement, poSpatialFilter, NULL ); if( poResultSet != NULL ) { if( !DissolveLayer( poDS, poResultSet, poODS, papszLCO, pszNewLayerName, bTransform, poOutputSRS, poSourceSRS, papszSelFields, bAppend, eGType, bOverwrite )) { CPLError( CE_Failure, CPLE_AppDefined, "Terminating translation prematurely after failed\n" "translation from sql statement." ); exit( 1 ); } poDS->ReleaseResultSet( poResultSet ); } } /* -------------------------------------------------------------------- */ /* Process each data source layer. */ /* -------------------------------------------------------------------- */ for( int iLayer = 0; pszSQLStatement == NULL && iLayer < poDS->GetLayerCount(); iLayer++ ) { OGRLayer *poLayer = poDS->GetLayer(iLayer); if( poLayer == NULL ) { printf( "FAILURE: Couldn't fetch advertised layer %d!\n", iLayer ); exit( 1 ); } if( CSLCount(papszLayers) == 0 || CSLFindString( papszLayers, poLayer->GetLayerDefn()->GetName() ) != -1 ) { if( pszWHERE != NULL ) poLayer->SetAttributeFilter( pszWHERE ); if( poSpatialFilter != NULL ) poLayer->SetSpatialFilter( poSpatialFilter ); if( !DissolveLayer( poDS, poLayer, poODS, papszLCO, pszNewLayerName, bTransform, poOutputSRS, poSourceSRS, papszSelFields, bAppend, eGType, bOverwrite ) && !bSkipFailures ) { CPLError( CE_Failure, CPLE_AppDefined, "Terminating translation prematurely after failed\n" "translation of layer %s\n", poLayer->GetLayerDefn()->GetName() ); exit( 1 ); } } } /* -------------------------------------------------------------------- */ /* Close down. */ /* -------------------------------------------------------------------- */ OGRSpatialReference::DestroySpatialReference( poOutputSRS ); OGRSpatialReference::DestroySpatialReference( poSourceSRS ) ; OGRDataSource::DestroyDataSource( poODS ); OGRDataSource::DestroyDataSource( poDS ); CSLDestroy(papszSelFields); CSLDestroy( papszArgv ); CSLDestroy( papszLayers ); CSLDestroy( papszDSCO ); CSLDestroy( papszLCO ); OGRCleanupAll(); #ifdef DBMALLOC malloc_dump(1); #endif return 0; } /************************************************************************/ /* Usage() */ /************************************************************************/ static void Usage() { OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar(); printf( "Usage: ogr2ogr [--help-general] [-skipfailures] [-append] [-update]\n" " [-select field_list] [-where restricted_where] \n" " [-sql <sql statement>] \n" " [-spat xmin ymin xmax ymax] [-preserve_fid] [-fid FID]\n" " [-a_srs srs_def] [-t_srs srs_def] [-s_srs srs_def]\n" " [-f format_name] [-overwrite] [[-dsco NAME=VALUE] ...]\n" " dst_datasource_name src_datasource_name\n" " [-lco NAME=VALUE] [-nln name] [-nlt type] [layer [layer ...]]\n" "\n" " -f format_name: output file format name, possible values are:\n"); for( int iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ ) { OGRSFDriver *poDriver = poR->GetDriver(iDriver); if( poDriver->TestCapability( ODrCCreateDataSource ) ) printf( " -f \"%s\"\n", poDriver->GetName() ); } printf( " -append: Append to existing layer instead of creating new if it exists\n" " -overwrite: delete the output layer and recreate it empty\n" " -update: Open existing output datasource in update mode\n" " -select field_list: Comma-delimited list of fields from input layer to\n" " copy to the new layer (defaults to all)\n" " -where restricted_where: Attribute query (like SQL WHERE)\n" " -sql statement: Execute given SQL statement and save result.\n" " -skipfailures: skip features or layers that fail to convert\n" " -spat xmin ymin xmax ymax: spatial query extents\n" " -dsco NAME=VALUE: Dataset creation option (format specific)\n" " -lco NAME=VALUE: Layer creation option (format specific)\n" " -nln name: Assign an alternate name to the new layer\n" " -nlt type: Force a geometry type for new layer. One of NONE, GEOMETRY,\n" " POINT, LINESTRING, POLYGON, GEOMETRYCOLLECTION, MULTIPOINT, MULTILINE,\n" " MULTIPOLYGON, or MULTILINESTRING. Add \"25D\" for 3D layers.\n" " Default is type of source layer.\n" ); printf(" -a_srs srs_def: Assign an output SRS\n" " -t_srs srs_def: Reproject/transform to this SRS on output\n" " -s_srs srs_def: Override source SRS\n" "\n" " Srs_def can be a full WKT definition (hard to escape properly),\n" " or a well known definition (i.e. EPSG:4326) or a file with a WKT\n" " definition.\n" ); exit( 1 ); } StringGeometryMap* CollectGeometries( OGRLayer* poSrcLayer, const char** papszFields) { /* -------------------------------------------------------------------- */ /* CollectGeometries returns a dictionary where the keys are the */ /* values in the fields that the user has selected and the values */ /* are a GeometryCollection of all of the geometries for records */ /* with that value. */ /* -------------------------------------------------------------------- */ StringGeometryMMap poGeometriesMap; OGRFeature *poFeature; poSrcLayer->ResetReading(); int iField; /* -------------------------------------------------------------------- */ /* Get all of the features and put them in a multi map. This may */ /* include values for which the selected fields is NULL. */ /* -------------------------------------------------------------------- */ while( true ) { poFeature = poSrcLayer->GetNextFeature(); if( poFeature == NULL ) break; CPLString poKey(""); for( iField=0; papszFields[iField] != NULL; iField++) { int nField = poFeature->GetFieldIndex(papszFields[iField]); poKey = poKey + poFeature->GetFieldAsString(nField); } if (poFeature->GetGeometryRef()->IsValid()) { poGeometriesMap.insert(std::make_pair( CPLString( poKey), poFeature->GetGeometryRef() ) ); } else { CPLDebug("CollectGeometries", "Geometry was invalid not adding!!!!"); } } /* -------------------------------------------------------------------- */ /* Loop through our features map and get a unique list of field */ /* values. This could be done using something other than a map */ /* of course, but it was convenient. */ /* -------------------------------------------------------------------- */ typedef std::map<CPLString, int> StringIntMap; StringIntMap::const_iterator ipos; StringIntMap poFieldsmap; StringGeometryMMap::const_iterator pos; for (pos = poGeometriesMap.begin(); pos != poGeometriesMap.end(); ++pos) { /* we currently throw out any null field values at this time */ // if (!(pos->first.empty())) { poFieldsmap[CPLString(pos->first.c_str())] = 1; // } } /* -------------------------------------------------------------------- */ /* Make a new map of GeometryCollection for each value in the */ /* poFieldsmap. This is a 1:1 relationship, and all of the */ /* geometries for a given field are all put into the same */ /* GeometryCollection. After we build the poCollections, we will */ /* use GEOS' buffer(0) trick to have GEOS perform the segmentation */ /* -------------------------------------------------------------------- */ StringGeometryColMap poCollections; CPLDebug("CollectGeometries", "Field map size: %d", poFieldsmap.size()); for (ipos = poFieldsmap.begin(); ipos != poFieldsmap.end(); ++ipos) { CPLString fid = ipos->first; CPLDebug ("CollectGeometries", "First %s Second %d", ipos->first.c_str(), ipos->second); OGRGeometryCollection* geom = new OGRGeometryCollection; for (pos = poGeometriesMap.lower_bound(fid); pos != poGeometriesMap.upper_bound(fid); ++pos) { geom->addGeometry(pos->second); } poCollections.insert(std::make_pair(fid, geom)); } CPLDebug("CollectGeometries", "Geo map size: %d", poCollections.size()); /* -------------------------------------------------------------------- */ /* Loop through our poCollections map and buffer(0) each */ /* GeometryCollection. GEOS will collapse the geometries down */ /* -------------------------------------------------------------------- */ StringGeometryMap* buffers = new StringGeometryMap; StringGeometryColMap::const_iterator collections_i; for ( collections_i = poCollections.begin(); collections_i != poCollections.end(); ++collections_i){ CPLDebug( "CollectGeometries", "poCollections Geometry size %d", collections_i->second->getNumGeometries()); OGRGeometry* buffer = collections_i->second->Buffer(0); buffers->insert(std::make_pair(collections_i->first, buffer)); } for (collections_i = poCollections.begin(); collections_i != poCollections.end(); ++collections_i) { delete collections_i->second; } return buffers; } GeometriesList* FlattenGeometries(GeometriesList* input) { GeometriesList::const_iterator geometry_i; GeometriesList* output = new GeometriesList; CPLDebug("CollectGeometries", "Input geometries in FlattenGeometries size: %d", input->size()); for ( geometry_i = input->begin(); geometry_i != input->end(); ++geometry_i){ OGRGeometry* buffer = (*geometry_i); // int nGeometries = buffer->getNumGeometries(); OGRwkbGeometryType iGType = buffer->getGeometryType(); if (iGType == wkbPolygon) { OGRPolygon* geom = (OGRPolygon*)buffer; output->push_back((OGRGeometry*)geom); CPLDebug( "CollectGeometries", "Collapsing wkbPolygon geometries......" ); } if (iGType == wkbMultiPolygon) { OGRMultiPolygon* geom = (OGRMultiPolygon*)buffer; for (int i=0; i< geom->getNumGeometries(); i++) { OGRPolygon* g = (OGRPolygon*)geom->getGeometryRef(i); output->push_back((OGRGeometry*)g); } CPLDebug( "CollectGeometries", "Collapsing wkbMultiPolygon geometries......" ); } if (iGType == wkbGeometryCollection) { OGRGeometryCollection* geom = (OGRGeometryCollection*)buffer; GeometriesList* collection = new GeometriesList; GeometriesList::const_iterator g_i; for (int i=0; i< geom->getNumGeometries(); i++) { OGRGeometry* g = (OGRGeometry*)geom->getGeometryRef(i); collection->push_back(g); } GeometriesList* collapsed = FlattenGeometries(collection); for( g_i = collapsed->begin(); g_i != collapsed->end(); g_i++ ) { output->push_back((OGRGeometry*)(*g_i)); CPLDebug( "CollectGeometries", "Collapsing wkbGeometryCollection geometries." ); } } // CPLDebug( "CollectGeometries", // "Buffered Geometry size %d", // nGeometries); } return output; } /************************************************************************/ /* DissolveLayer() */ /************************************************************************/ static int DissolveLayer( OGRDataSource *poSrcDS, OGRLayer * poSrcLayer, OGRDataSource *poDstDS, char **papszLCO, const char *pszNewLayerName, int bTransform, OGRSpatialReference *poOutputSRS, OGRSpatialReference *poSourceSRS, char **papszSelFields, int bAppend, int eGType, int bOverwrite ) { OGRLayer *poDstLayer; OGRFeatureDefn *poFDefn; OGRErr eErr; int bForceToPolygon = FALSE; int bForceToMultiPolygon = FALSE; if( pszNewLayerName == NULL ) pszNewLayerName = poSrcLayer->GetLayerDefn()->GetName(); if( wkbFlatten(eGType) == wkbPolygon ) bForceToPolygon = TRUE; else if( wkbFlatten(eGType) == wkbMultiPolygon ) bForceToMultiPolygon = TRUE; /* -------------------------------------------------------------------- */ /* Setup coordinate transformation if we need it. */ /* -------------------------------------------------------------------- */ OGRCoordinateTransformation *poCT = NULL; if( bTransform ) { if( poSourceSRS == NULL ) poSourceSRS = poSrcLayer->GetSpatialRef(); if( poSourceSRS == NULL ) { printf( "Can't transform coordinates, source layer has no\n" "coordinate system. Use -s_srs to set one.\n" ); exit( 1 ); } CPLAssert( NULL != poSourceSRS ); CPLAssert( NULL != poOutputSRS ); poCT = OGRCreateCoordinateTransformation( poSourceSRS, poOutputSRS ); if( poCT == NULL ) { char *pszWKT = NULL; printf("Failed to create coordinate transformation between the\n" "following coordinate systems. This may be because they\n" "are not transformable, or because projection services\n" "(PROJ.4 DLL/.so) could not be loaded.\n" ); poSourceSRS->exportToPrettyWkt( &pszWKT, FALSE ); printf( "Source:\n%s\n", pszWKT ); poOutputSRS->exportToPrettyWkt( &pszWKT, FALSE ); printf( "Target:\n%s\n", pszWKT ); exit( 1 ); } } /* -------------------------------------------------------------------- */ /* Get other info. */ /* -------------------------------------------------------------------- */ poFDefn = poSrcLayer->GetLayerDefn(); if( poOutputSRS == NULL ) poOutputSRS = poSrcLayer->GetSpatialRef(); /* -------------------------------------------------------------------- */ /* Find the layer. */ /* -------------------------------------------------------------------- */ int iLayer = -1; poDstLayer = NULL; for( iLayer = 0; iLayer < poDstDS->GetLayerCount(); iLayer++ ) { OGRLayer *poLayer = poDstDS->GetLayer(iLayer); if( poLayer != NULL && EQUAL(poLayer->GetLayerDefn()->GetName(),pszNewLayerName) ) { poDstLayer = poLayer; break; } } /* -------------------------------------------------------------------- */ /* If the user requested overwrite, and we have the layer in */ /* question we need to delete it now so it will get recreated */ /* (overwritten). */ /* -------------------------------------------------------------------- */ if( poDstLayer != NULL && bOverwrite ) { if( poDstDS->DeleteLayer( iLayer ) != OGRERR_NONE ) { fprintf( stderr, "DeleteLayer() failed when overwrite requested.\n" ); return FALSE; } poDstLayer = NULL; } /* -------------------------------------------------------------------- */ /* If the layer does not exist, then create it. */ /* -------------------------------------------------------------------- */ if( poDstLayer == NULL ) { if( eGType == -2 ) eGType = poFDefn->GetGeomType(); if( !poDstDS->TestCapability( ODsCCreateLayer ) ) { fprintf( stderr, "Layer %s not found, and CreateLayer not supported by driver.", pszNewLayerName ); return FALSE; } CPLErrorReset(); poDstLayer = poDstDS->CreateLayer( pszNewLayerName, poOutputSRS, (OGRwkbGeometryType) eGType, papszLCO ); if( poDstLayer == NULL ) return FALSE; bAppend = FALSE; } /* -------------------------------------------------------------------- */ /* Otherwise we will append to it, if append was requested. */ /* -------------------------------------------------------------------- */ else if( !bAppend ) { printf( "FAILED: Layer %s already exists, and -append not specified.\n" " Consider using -append, or -overwrite.\n", pszNewLayerName ); return FALSE; } /* -------------------------------------------------------------------- */ /* Add fields. Default to copy all field. */ /* If only a subset of all fields requested, then output only */ /* the selected fields, and in the order that they were */ /* selected. */ /* -------------------------------------------------------------------- */ int iField; if (papszSelFields && !bAppend ) { for( iField=0; papszSelFields[iField] != NULL; iField++) { int iSrcField = poFDefn->GetFieldIndex(papszSelFields[iField]); if (iSrcField >= 0) poDstLayer->CreateField( poFDefn->GetFieldDefn(iSrcField) ); else { printf( "Field '%s' not found in source layer.\n", papszSelFields[iField] ); if( !bSkipFailures ) return FALSE; } } } else if( !bAppend ) { for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ ) poDstLayer->CreateField( poFDefn->GetFieldDefn(iField) ); } /* -------------------------------------------------------------------- */ /* Transfer features. */ /* -------------------------------------------------------------------- */ OGRFeature *poFeature; int nFeaturesInTransaction = 0; poSrcLayer->ResetReading(); if( nGroupTransactions ) poDstLayer->StartTransaction(); StringGeometryMap* buffers = CollectGeometries(poSrcLayer, (const char**)papszSelFields); StringGeometryMap::const_iterator buffers_i; GeometriesList* input = new GeometriesList; CPLDebug("CollectGeometries", "Buffers size: %d", buffers->size()); for ( buffers_i = buffers->begin(); buffers_i != buffers->end(); ++buffers_i){ input->push_back(buffers_i->second); } GeometriesList* geometries = FlattenGeometries(input); GeometriesList::const_iterator g_i; for (g_i=geometries->begin(); g_i!=geometries->end(); g_i++) { OGRFeature* feature = new OGRFeature(poFDefn); feature->SetGeometry((*g_i)); feature->SetField("TAXDIST","fid"); poDstLayer->CreateFeature(feature); } if( nGroupTransactions ) poDstLayer->CommitTransaction(); // getGeometryType // if( pszNewLayerName == NULL ) // pszNewLayerName = poSrcLayer->GetLayerDefn()->GetName(); // // if( wkbFlatten(eGType) == wkbPolygon ) // bForceToPolygon = TRUE; // else if( wkbFlatten(eGType) == wkbMultiPolygon ) // bForceToMultiPolygon = TRUE; // // /* -------------------------------------------------------------------- */ // /* Setup coordinate transformation if we need it. */ // /* -------------------------------------------------------------------- */ // OGRCoordinateTransformation *poCT = NULL; // // if( bTransform ) // { // if( poSourceSRS == NULL ) // poSourceSRS = poSrcLayer->GetSpatialRef(); // // if( poSourceSRS == NULL ) // { // printf( "Can't transform coordinates, source layer has no\n" // "coordinate system. Use -s_srs to set one.\n" ); // exit( 1 ); // } // // CPLAssert( NULL != poSourceSRS ); // CPLAssert( NULL != poOutputSRS ); // // poCT = OGRCreateCoordinateTransformation( poSourceSRS, poOutputSRS ); // if( poCT == NULL ) // { // char *pszWKT = NULL; // // printf("Failed to create coordinate transformation between the\n" // "following coordinate systems. This may be because they\n" // "are not transformable, or because projection services\n" // "(PROJ.4 DLL/.so) could not be loaded.\n" ); // // poSourceSRS->exportToPrettyWkt( &pszWKT, FALSE ); // printf( "Source:\n%s\n", pszWKT ); // // poOutputSRS->exportToPrettyWkt( &pszWKT, FALSE ); // printf( "Target:\n%s\n", pszWKT ); // exit( 1 ); // } // } // // /* -------------------------------------------------------------------- */ // /* Get other info. */ // /* -------------------------------------------------------------------- */ // poFDefn = poSrcLayer->GetLayerDefn(); // // if( poOutputSRS == NULL ) // poOutputSRS = poSrcLayer->GetSpatialRef(); // // /* -------------------------------------------------------------------- */ // /* Find the layer. */ // /* -------------------------------------------------------------------- */ // int iLayer = -1; // poDstLayer = NULL; // // for( iLayer = 0; iLayer < poDstDS->GetLayerCount(); iLayer++ ) // { // OGRLayer *poLayer = poDstDS->GetLayer(iLayer); // // if( poLayer != NULL // && EQUAL(poLayer->GetLayerDefn()->GetName(),pszNewLayerName) ) // { // poDstLayer = poLayer; // break; // } // } // // /* -------------------------------------------------------------------- */ // /* If the user requested overwrite, and we have the layer in */ // /* question we need to delete it now so it will get recreated */ // /* (overwritten). */ // /* -------------------------------------------------------------------- */ // if( poDstLayer != NULL && bOverwrite ) // { // if( poDstDS->DeleteLayer( iLayer ) != OGRERR_NONE ) // { // fprintf( stderr, // "DeleteLayer() failed when overwrite requested.\n" ); // return FALSE; // } // poDstLayer = NULL; // } // // /* -------------------------------------------------------------------- */ // /* If the layer does not exist, then create it. */ // /* -------------------------------------------------------------------- */ // if( poDstLayer == NULL ) // { // if( eGType == -2 ) // eGType = poFDefn->GetGeomType(); // // if( !poDstDS->TestCapability( ODsCCreateLayer ) ) // { // fprintf( stderr, // "Layer %s not found, and CreateLayer not supported by driver.", // pszNewLayerName ); // return FALSE; // } // // CPLErrorReset(); // // poDstLayer = poDstDS->CreateLayer( pszNewLayerName, poOutputSRS, // (OGRwkbGeometryType) eGType, // papszLCO ); // // if( poDstLayer == NULL ) // return FALSE; // // bAppend = FALSE; // } // // /* -------------------------------------------------------------------- */ // /* Otherwise we will append to it, if append was requested. */ // /* -------------------------------------------------------------------- */ // else if( !bAppend ) // { // printf( "FAILED: Layer %s already exists, and -append not specified.\n" // " Consider using -append, or -overwrite.\n", // pszNewLayerName ); // return FALSE; // } // // /* -------------------------------------------------------------------- */ // /* Add fields. Default to copy all field. */ // /* If only a subset of all fields requested, then output only */ // /* the selected fields, and in the order that they were */ // /* selected. */ // /* -------------------------------------------------------------------- */ // int iField; // // if (papszSelFields && !bAppend ) // { // for( iField=0; papszSelFields[iField] != NULL; iField++) // { // int iSrcField = poFDefn->GetFieldIndex(papszSelFields[iField]); // if (iSrcField >= 0) // poDstLayer->CreateField( poFDefn->GetFieldDefn(iSrcField) ); // else // { // printf( "Field '%s' not found in source layer.\n", // papszSelFields[iField] ); // if( !bSkipFailures ) // return FALSE; // } // } // } // else if( !bAppend ) // { // for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ ) // poDstLayer->CreateField( poFDefn->GetFieldDefn(iField) ); // } // // /* -------------------------------------------------------------------- */ // /* Transfer features. */ // /* -------------------------------------------------------------------- */ // OGRFeature *poFeature; // int nFeaturesInTransaction = 0; // // poSrcLayer->ResetReading(); // // if( nGroupTransactions ) // poDstLayer->StartTransaction(); // // while( true ) // { // OGRFeature *poDstFeature = NULL; // // if( nFIDToFetch != OGRNullFID ) // { // // Only fetch feature on first pass. // if( nFeaturesInTransaction == 0 ) // poFeature = poSrcLayer->GetFeature(nFIDToFetch); // else // poFeature = NULL; // } // else // poFeature = poSrcLayer->GetNextFeature(); // // if( poFeature == NULL ) // break; // // if( ++nFeaturesInTransaction == nGroupTransactions ) // { // poDstLayer->CommitTransaction(); // poDstLayer->StartTransaction(); // nFeaturesInTransaction = 0; // } // // CPLErrorReset(); // poFDefn = poSrcLayer->GetLayerDefn(); // // if( poDstFeature->SetFrom( poFeature, TRUE ) != OGRERR_NONE ) // { // if( nGroupTransactions ) // poDstLayer->CommitTransaction(); // // CPLError( CE_Failure, CPLE_AppDefined, // "Unable to translate feature %d from layer %s.\n", // poFeature->GetFID(), poFDefn->GetName() ); // // OGRFeature::DestroyFeature( poFeature ); // OGRFeature::DestroyFeature( poDstFeature ); // return FALSE; // } // // if( bPreserveFID ) // poDstFeature->SetFID( poFeature->GetFID() ); // // if( poCT && poDstFeature->GetGeometryRef() != NULL ) // { // eErr = poDstFeature->GetGeometryRef()->transform( poCT ); // if( eErr != OGRERR_NONE ) // { // if( nGroupTransactions ) // poDstLayer->CommitTransaction(); // // printf( "Failed to transform feature %d.\n", // (int) poFeature->GetFID() ); // if( !bSkipFailures ) // { // OGRFeature::DestroyFeature( poFeature ); // OGRFeature::DestroyFeature( poDstFeature ); // return FALSE; // } // } // } // // if( poDstFeature->GetGeometryRef() != NULL && bForceToPolygon ) // { // poDstFeature->SetGeometryDirectly( // OGRGeometryFactory::forceToPolygon( // poDstFeature->StealGeometry() ) ); // } // // if( poDstFeature->GetGeometryRef() != NULL && bForceToMultiPolygon ) // { // poDstFeature->SetGeometryDirectly( // OGRGeometryFactory::forceToMultiPolygon( // poDstFeature->StealGeometry() ) ); // } // // OGRFeature::DestroyFeature( poFeature ); // // CPLErrorReset(); // if( poDstLayer->CreateFeature( poDstFeature ) != OGRERR_NONE // && !bSkipFailures ) // { // if( nGroupTransactions ) // poDstLayer->RollbackTransaction(); // // OGRFeature::DestroyFeature( poDstFeature ); // return FALSE; // } // // OGRFeature::DestroyFeature( poDstFeature ); // } // // if( nGroupTransactions ) // poDstLayer->CommitTransaction(); // // /* -------------------------------------------------------------------- */ // /* Cleaning */ // /* -------------------------------------------------------------------- */ // delete poCT; return TRUE; }