[Gdal-dev] Tool to convert shapefile from [-180, 180] to [0, 360]
range
Daniel Morissette
dmorissette at mapgears.com
Wed Jun 14 16:59:25 EDT 2006
OOpps.. the file I just sent contained a typo in the bNeedToShift test.
Works fine for points still since minx == maxx, but would not have
worked too well for some polygons or linestrings.
Here is what it should have been:
bNeedToShift = ((oEnv.MinX+oEnv.MaxX)/2.0<0.0) ? TRUE : FALSE;
Daniel
Daniel Morissette wrote:
> Mateusz Loskot wrote:
>
>>
>> I created small utility on basis of ogr2ogr according your needs, just
>> for fun and practice :-)
>> I attached it as a 180to360.cpp file, so it's a C++ application.
>> I hope this attachement will not be cutted off by the gdal-dev mailer.
>> It works only with shapefiles.
>>
>> This utility just adds 180 to every x coordinate.
>> So, it works well only if input coordinates are in range frm -180 to
>> 180, but it does not check it.
>
>
> Hi Mateusz,
>
> Thanks for the little program, that's what I was thinking of doing too.
> In case anyone is interested, I have attached a new version that I
> modified it a little bit to:
>
> 1- Shift by +360 instead of +180. i.e. -50 degrees = +310 degrees
> 2- Shift coordinates only if the centroid of the source feature is in
> the [-180,0] range
>
> One last thing to do would be support for geometries other than points
> (at least linestrings and polygons).
>
> Daniel
>
>
> ------------------------------------------------------------------------
>
> // Converts longtitude format from -180;180 to 0;360 degrees
> // This program was created on basis of ogr2ogr
>
> #include "ogrsf_frmts.h"
> #include "ogr_p.h"
> #include "cpl_conv.h"
> #include "cpl_string.h"
> #include "ogr_api.h"
>
> static void Usage();
>
> static int TranslateLayer( 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;
>
> /************************************************************************/
> /* 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 != 3 )
> {
> Usage();
> exit( -nArgc );
> }
>
> pszDataSource = papszArgv[1];
> pszDestDataSource = papszArgv[2];
>
> /* -------------------------------------------------------------------- */
> /* 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 );
> exit( 1 );
> }
>
> /* -------------------------------------------------------------------- */
> /* Try opening the output datasource as an existing, writable */
> /* -------------------------------------------------------------------- */
> OGRDataSource *poODS;
>
> OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar();
> OGRSFDriver *poDriver = NULL;
> int iDriver;
>
> poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(pszFormat);
> if( poDriver == NULL )
> {
> printf( "Unable to use driver `%s'.\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 );
> }
>
> /* -------------------------------------------------------------------- */
> /* 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( !TranslateLayer( 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. */
> /* -------------------------------------------------------------------- */
> delete poODS;
> delete poDS;
>
> CSLDestroy(papszSelFields);
> CSLDestroy( papszArgv );
> CSLDestroy( papszLayers );
>
> OGRCleanupAll();
>
> #ifdef DBMALLOC
> malloc_dump(1);
> #endif
>
> return 0;
> }
>
> /************************************************************************/
> /* Usage() */
> /************************************************************************/
>
> static void Usage()
>
> {
> OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar();
>
> printf( "Usage: ogr2ogr <input dst> <output dst>\n"
> );
>
> exit( 1 );
> }
>
> /************************************************************************/
> /* TranslateLayer() */
> /************************************************************************/
>
> static int TranslateLayer( 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();
>
> /* -------------------------------------------------------------------- */
> /* 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 layer does not exist, then create it. */
> /* -------------------------------------------------------------------- */
> if( poDstLayer == NULL )
> {
> if( eGType == -2 )
> eGType = poFDefn->GetGeomType();
>
> CPLAssert( poDstDS->TestCapability( ODsCCreateLayer ) );
> CPLErrorReset();
>
> poDstLayer = poDstDS->CreateLayer( pszNewLayerName, poOutputSRS,
> (OGRwkbGeometryType) eGType,
> papszLCO );
>
> if( poDstLayer == NULL )
> return FALSE;
>
> bAppend = 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;
>
> for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
> {
> poDstLayer->CreateField( poFDefn->GetFieldDefn(iField) );
> }
>
> /* -------------------------------------------------------------------- */
> /* Transfer features. */
> /* -------------------------------------------------------------------- */
> OGRFeature *poFeature;
>
> poSrcLayer->ResetReading();
>
> while( TRUE )
> {
> OGRFeature *poDstFeature = NULL;
>
> poFeature = poSrcLayer->GetNextFeature();
>
> if( poFeature == NULL )
> break;
>
> CPLErrorReset();
> poDstFeature = OGRFeature::CreateFeature( poDstLayer->GetLayerDefn() );
>
> if( poDstFeature->SetFrom( poFeature, TRUE ) != OGRERR_NONE )
> {
> 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() );
>
> OGRGeometry* poGeom = poDstFeature->GetGeometryRef();
>
> // We will shift coordinates by +360 only if object centroid is in the
> // [-180, 0] range, otherwise we leave the coordinates untouched
> GBool bNeedToShift = FALSE;
> if (poGeom)
> {
> OGREnvelope oEnv;
> poGeom->getEnvelope(&oEnv);
> bNeedToShift = ((oEnv.MinX+oEnv.MinX)/2<0.0) ? TRUE : FALSE;
> }
>
> if( bNeedToShift )
> {
> OGRwkbGeometryType wkbType = poGeom->getGeometryType();
> if ( wkbPoint == wkbType )
> {
> OGRPoint* poPoint = static_cast<OGRPoint*>(poGeom);
> double x = poPoint->getX();
> if (bNeedToShift)
> poPoint->setX(x + 360);
> }
> }
>
> OGRFeature::DestroyFeature( poFeature );
>
> CPLErrorReset();
> if( poDstLayer->CreateFeature( poDstFeature ) != OGRERR_NONE
> && !bSkipFailures )
> {
> OGRFeature::DestroyFeature( poDstFeature );
> return FALSE;
> }
>
> OGRFeature::DestroyFeature( poDstFeature );
> }
>
> return TRUE;
> }
>
--
Daniel Morissette
http://www.mapgears.com/
More information about the Gdal-dev
mailing list