diff --git a/gdal/apps/gdalsrsinfo.cpp b/gdal/apps/gdalsrsinfo.cpp index 00b602d..0265e75 100644 --- a/gdal/apps/gdalsrsinfo.cpp +++ b/gdal/apps/gdalsrsinfo.cpp @@ -43,8 +43,6 @@ CPLErr PrintSRS( const OGRSpatialReference &oSRS, int bPretty, int bPrintSep ); void PrintSRSOutputTypes( const OGRSpatialReference &oSRS, const char ** papszOutputTypes ); -int FindEPSG( const OGRSpatialReference &oSRS ); -int SearchCSVForWKT( const char *pszFileCSV, const char *pszTarget ); /************************************************************************/ /* Usage() */ @@ -176,7 +174,7 @@ int main( int argc, char ** argv ) if ( bFindEPSG ) { CPLError( CE_Warning, CPLE_AppDefined, "EPSG detection is experimental and requires new data files (see bug #4345)" ); - nEPSGCode = FindEPSG( oSRS ); + nEPSGCode = oSRS.FindEPSG( ); /* If found, replace oSRS based on EPSG code */ if(nEPSGCode != -1) { CPLDebug( "gdalsrsinfo", @@ -252,13 +250,13 @@ int FindSRS( const char *pszInput, OGRSpatialReference &oSRS, int bDebug ) { int bGotSRS = FALSE; - VSILFILE *fp = NULL; + // VSILFILE *fp = NULL; GDALDataset *poGDALDS = NULL; OGRDataSource *poOGRDS = NULL; OGRLayer *poLayer = NULL; char *pszProjection = NULL; CPLErrorHandler oErrorHandler = NULL; - int bIsFile = FALSE; + // int bIsFile = FALSE; OGRErr eErr = CE_None; /* temporarily supress error messages we may get from xOpen() */ @@ -266,11 +264,15 @@ int FindSRS( const char *pszInput, OGRSpatialReference &oSRS, int bDebug ) oErrorHandler = CPLSetErrorHandler ( CPLQuietErrorHandler ); /* If argument is a file, try to open it with GDAL and OGROpen() */ - fp = VSIFOpenL( pszInput, "r" ); - if ( fp ) { + /* remove this requirement, to support NETCDF:file.nc:file syntax */ + /* TODO - support OGR layer(s) */ + // fp = VSIFOpenL( pszInput, "r" ); + // if ( fp ) { - bIsFile = TRUE; - VSIFCloseL( fp ); + // CPLDebug( "gdalsrsinfo", "argument %s is a file", pszInput ); + + // bIsFile = TRUE; + // VSIFCloseL( fp ); /* try to open with GDAL */ CPLDebug( "gdalsrsinfo", "trying to open with GDAL" ); @@ -309,7 +311,7 @@ int FindSRS( const char *pszInput, OGRSpatialReference &oSRS, int bDebug ) CPLDebug( "gdalsrsinfo", "did not open with OGR" ); } - } + // } /* Try ESRI file */ if ( ! bGotSRS && (strstr(pszInput,".prj") != NULL) ) { @@ -376,9 +378,6 @@ CPLErr PrintSRS( const OGRSpatialReference &oSRS, if ( ! pszOutputType || EQUAL(pszOutputType,"")) return CE_None; - CPLDebug( "gdalsrsinfo", "PrintSRS( oSRS, %s, %d, %d )\n", - pszOutputType, bPretty, bPrintSep ); - char *pszOutput = NULL; if ( EQUAL("proj4", pszOutputType ) ) { @@ -466,173 +465,3 @@ void PrintSRSOutputTypes( const OGRSpatialReference &oSRS, printf( "\n" ); } } - -/************************************************************************/ -/* SearchCSVForWKT() */ -/* */ -/* Search CSV file for target WKT, return EPSG code (or -1). */ -/* For saving space, the file can be compressed (gz) */ -/* If CSV file is absent are absent the function silently exits */ -/************************************************************************/ -int SearchCSVForWKT( const char *pszFileCSV, const char *pszTarget ) -{ - const char *pszFilename = NULL; - const char *pszWKT = NULL; - char szTemp[1024]; - int nPos = 0; - const char *pszTemp = NULL; - - VSILFILE *fp = NULL; - OGRSpatialReference oSRS; - int nCode = 0; - int nFound = -1; - - CPLDebug( "gdalsrsinfo", - "SearchCSVForWKT()\nfile=%s\nWKT=%s\n", - pszFileCSV, pszTarget); - -/* -------------------------------------------------------------------- */ -/* Find and open file. */ -/* -------------------------------------------------------------------- */ - // pszFilename = pszFileCSV; - pszFilename = CPLFindFile( "gdal", pszFileCSV ); - if( pszFilename == NULL ) - { - CPLDebug( "gdalsrsinfo", "could not find support file %s", - pszFileCSV ); - // return OGRERR_UNSUPPORTED_SRS; - return -1; - } - - /* support gzipped file */ - if ( strstr( pszFileCSV,".gz") != NULL ) - sprintf( szTemp, "/vsigzip/%s", pszFilename); - else - sprintf( szTemp, "%s", pszFilename); - - CPLDebug( "gdalsrsinfo", "SearchCSVForWKT() using file %s", - szTemp ); - - fp = VSIFOpenL( szTemp, "r" ); - if( fp == NULL ) - { - CPLDebug( "gdalsrsinfo", "could not open support file %s", - pszFilename ); - - // return OGRERR_UNSUPPORTED_SRS; - return -1; - } - -/* -------------------------------------------------------------------- */ -/* Process lines. */ -/* -------------------------------------------------------------------- */ - const char *pszLine; - - while( (pszLine = CPLReadLine2L(fp,-1,NULL)) != NULL ) - - { - // CPLDebug( "gdalsrsinfo", "read line %s", pszLine ); - - if( pszLine[0] == '#' ) - continue; - /* do nothing */; - - // else if( EQUALN(pszLine,"include ",8) ) - // { - // eErr = importFromDict( pszLine + 8, pszCode ); - // if( eErr != OGRERR_UNSUPPORTED_SRS ) - // break; - // } - - // else if( strstr(pszLine,",") == NULL ) - // /* do nothing */; - - pszTemp = strstr(pszLine,","); - if (pszTemp) - { - nPos = pszTemp - pszLine; - - if ( nPos == 0 ) - continue; - - strncpy( szTemp, pszLine, nPos ); - szTemp[nPos] = '\0'; - nCode = atoi(szTemp); - - pszWKT = (char *) pszLine + nPos +1; - - // CPLDebug( "gdalsrsinfo", - // "code=%d\nWKT=\n[%s]\ntarget=\n[%s]\n", - // nCode,pszWKT, pszTarget ); - - if ( EQUAL(pszTarget,pszWKT) ) - { - nFound = nCode; - CPLDebug( "gdalsrsinfo", "found EPSG:%d\n" - "current=%s\ntarget= %s\n", - nCode, pszWKT, pszTarget ); - break; - } - } - } - - VSIFCloseL( fp ); - - return nFound; - -} - -/* TODO - - search for well-known values (AutoIdentifyEPSG()) - - - should we search .override.csv files? - - - fix precision differences (namely in degree: 17 vs 15) so we can use epsg_ogc_simple -target: - orig: GEOGCS["SAD69",DATUM["D_South_American_1969",SPHEROID["GRS_1967_Modified",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] - ESRI: GEOGCS["SAD69",DATUM["D_South_American_1969",SPHEROID["GRS_1967_Truncated",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] - OGC: GEOGCS["SAD69",DATUM["South_American_Datum_1969",SPHEROID["GRS_1967_Modified",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] -database: - ESRI: GEOGCS["SAD69",DATUM["D_South_American_1969",SPHEROID["GRS_1967_Truncated",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] - OGC: GEOGCS["SAD69",DATUM["South_American_Datum_1969",SPHEROID["GRS 1967 Modified",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]] -*/ - -/************************************************************************/ -/* FindEPSG() */ -/* */ -/* Return EPSG code corresponding to spatial reference (or -1) */ -/************************************************************************/ -int FindEPSG( const OGRSpatialReference &oSRS ) -{ - char *pszWKT = NULL; - char *pszESRI = NULL; - int nFound = -1; - OGRSpatialReference *poSRS = NULL; - - poSRS = oSRS.Clone(); - poSRS->StripCTParms(); - poSRS->exportToWkt( &pszWKT ); - OGRSpatialReference::DestroySpatialReference( poSRS ); - - poSRS = oSRS.Clone(); - poSRS->morphToESRI( ); - poSRS->exportToWkt( &pszESRI ); - OGRSpatialReference::DestroySpatialReference( poSRS ); - - CPLDebug( "gdalsrsinfo", "FindEPSG()\nWKT (OGC)= %s\nWKT (ESRI)=%s", - pszWKT,pszESRI ); - - /* search for EPSG code in epsg_*.wkt.gz files */ - /* using ESRI WKT for now, as it seems to work best */ - nFound = SearchCSVForWKT( "epsg_esri.wkt.gz", pszESRI ); - if ( nFound == -1 ) - nFound = SearchCSVForWKT( "epsg_ogc_simple.wkt.gz", pszESRI ); - if ( nFound == -1 ) - nFound = SearchCSVForWKT( "epsg_ogc.wkt.gz", pszESRI ); - - - CPLFree( pszWKT ); - CPLFree( pszESRI ); - - return nFound; -} diff --git a/gdal/ogr/ogr_fromepsg.cpp b/gdal/ogr/ogr_fromepsg.cpp index b17be29..5983da5 100644 --- a/gdal/ogr/ogr_fromepsg.cpp +++ b/gdal/ogr/ogr_fromepsg.cpp @@ -2692,4 +2692,93 @@ int OSREPSGTreatsAsLatLong( OGRSpatialReferenceH hSRS ) return ((OGRSpatialReference *) hSRS)->EPSGTreatsAsLatLong(); } + +/************************************************************************/ +/* FindEPSG() */ +/************************************************************************/ +/** + * \brief This method returns the EPSG code corresponding to this + * spatial reference, or -1 if none is found. +* + * This function looks through a number of dictionnary files containing + * WKT definitions of EPSG codes known to OGR + * (by calling OGRSpatialReference::SearchDictForWKT()). + * If these files are not installed no EPSG code will be found. + */ + +int OGRSpatialReference::FindEPSG( ) +{ +/* TODO + - search for well-known values (AutoIdentifyEPSG()) + + - should we search .override.csv files? + + - fix precision differences (namely in degree: 17 vs 15) so we can use epsg_ogc_simple +target: + orig: GEOGCS["SAD69",DATUM["D_South_American_1969",SPHEROID["GRS_1967_Modified",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] + ESRI: GEOGCS["SAD69",DATUM["D_South_American_1969",SPHEROID["GRS_1967_Truncated",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] + OGC: GEOGCS["SAD69",DATUM["South_American_Datum_1969",SPHEROID["GRS_1967_Modified",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] +database: + ESRI: GEOGCS["SAD69",DATUM["D_South_American_1969",SPHEROID["GRS_1967_Truncated",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] + OGC: GEOGCS["SAD69",DATUM["South_American_Datum_1969",SPHEROID["GRS 1967 Modified",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]] +*/ + char *pszWKT = NULL; + char *pszWKTNoCT = NULL; + char *pszWKTESRI = NULL; + int nFound = -1; + OGRSpatialReference *poSRS = NULL; + OGRSpatialReference *poSRSNoCT = NULL; + OGRSpatialReference *poSRSESRI = NULL; + + poSRS = this->Clone(); + poSRS->exportToWkt( &pszWKT ); + + poSRSNoCT = this->Clone(); + poSRSNoCT->StripCTParms(); + poSRSNoCT->exportToWkt( &pszWKTNoCT ); + + poSRSESRI = this->Clone(); + poSRSESRI->morphToESRI( ); + poSRSESRI->exportToWkt( &pszWKTESRI ); + + CPLDebug( "OGRSpatialReference", + "FindEPSG()\n WKT (OGC)= %s\n WKT (no CT)=%s\n WKT (ESRI)=%s", + pszWKT, pszWKTNoCT, pszWKTESRI ); + + /* search for EPSG code in epsg_*.wkt.gz files */ + /* using ESRI WKT first for now, as it seems to work best */ + /* TODO should we test IsSame() ? */ + nFound = OGRSpatialReference::SearchDictForWKT( "epsg_esri.wkt.gz", pszWKTESRI ); + // nFound = OGRSpatialReference::SearchDictForWKT( "epsg_esri.wkt", pszWKTESRI ); + if ( nFound == -1 ) + nFound = SearchDictForWKT( "epsg_ogc_simple.wkt.gz", pszWKTNoCT ); + if ( nFound == -1 ) + nFound = SearchDictForWKT( "epsg_ogc.wkt.gz", pszWKT ); + + OGRSpatialReference::DestroySpatialReference( poSRS ); + OGRSpatialReference::DestroySpatialReference( poSRSNoCT ); + OGRSpatialReference::DestroySpatialReference( poSRSESRI ); + CPLFree( pszWKT ); + CPLFree( pszWKTNoCT ); + CPLFree( pszWKTESRI ); + + return nFound; +} +/************************************************************************/ +/* OSRFindEPSG() */ +/************************************************************************/ +/** + * \brief This method returns the EPSG code corresponding to this + * spatial reference, or -1 if none is found. + * + * This function is the same as OGRSpatialReference::FindEPSG(). + */ + +int OSRFindEPSG( OGRSpatialReferenceH hSRS ) + +{ + VALIDATE_POINTER1( hSRS, "OSRFindEPSG", CE_Failure ); + + return ((OGRSpatialReference *) hSRS)->FindEPSG(); +} diff --git a/gdal/ogr/ogr_spatialref.h b/gdal/ogr/ogr_spatialref.h index 921a4d2..53d62a3 100644 --- a/gdal/ogr/ogr_spatialref.h +++ b/gdal/ogr/ogr_spatialref.h @@ -528,8 +528,13 @@ class CPL_DLL OGRSpatialReference OGRErr ImportFromESRIWisconsinWKT( const char* pszPrjName, double dfCentralMeridian, double dfLatOfOrigin, const char* pszUnitsName, const char* pszCSName = 0 ); + + /** Find EPSG */ + int FindEPSG( ); + static int SearchDictForWKT( const char *pszDictFile, const char *pszWKT ); }; + /************************************************************************/ /* OGRCoordinateTransformation */ /* */ diff --git a/gdal/ogr/ogr_srs_api.h b/gdal/ogr/ogr_srs_api.h index 778468e..9c4a6cc 100644 --- a/gdal/ogr/ogr_srs_api.h +++ b/gdal/ogr/ogr_srs_api.h @@ -649,6 +649,10 @@ OGRErr CPL_DLL OSRSetWagner( OGRSpatialReferenceH hSRS, int nVariation, void CPL_DLL OSRCleanup( void ); +/** Find EPSG */ +int OSRFindEPSG( OGRSpatialReferenceH hSRS ); +int OSRSearchDictForWKT( const char *pszDictFile, const char *pszWKT ); + /* -------------------------------------------------------------------- */ /* OGRCoordinateTransform C API. */ /* -------------------------------------------------------------------- */ diff --git a/gdal/ogr/ogr_srs_dict.cpp b/gdal/ogr/ogr_srs_dict.cpp index 89b6590..b127543 100644 --- a/gdal/ogr/ogr_srs_dict.cpp +++ b/gdal/ogr/ogr_srs_dict.cpp @@ -132,3 +132,121 @@ OGRErr OSRImportFromDict( OGRSpatialReferenceH hSRS, return ((OGRSpatialReference *) hSRS)->importFromDict( pszDictFile, pszCode ); } + + + +/************************************************************************/ +/* SearchDictForWKT() */ +/************************************************************************/ +/** + * \brief Search SRS from WKT dictionary and return corresponding EPSG code. + * + * This method will attempt to find the indicated coordinate system identity + * in the indicated dictionary file. If found, the EPSG code corresponding + * to the WKT representation is returned. If not found, -1 is returned. + * + * @param pszDictFile the name of the dictionary file to load. + * + * @param pszWKT the WKT to lookup in the dictionary. + * + * @return EPSG code on success, or -1 if the code isn't found + */ + +int OGRSpatialReference::SearchDictForWKT( const char *pszDictFile, + const char *pszWKT ) +{ + const char *pszFilename = NULL; + const char *pszTemp = NULL; + char szTemp[1024]; + int nPos = 0; + // const char *pszTemp = NULL; + + VSILFILE *fp = NULL; + int nCode = 0; + int nFound = -1; + + CPLDebug( "OSRSearchDictForWKT()", "file=%s\n WKT=%s", + pszDictFile, pszWKT); + +/* -------------------------------------------------------------------- */ +/* Find and open file. */ +/* -------------------------------------------------------------------- */ + pszFilename = CPLDefaultFindFile( "gdal", pszDictFile ); + if( pszFilename == NULL ) + { + CPLDebug( "OSRSearchDictForWKT()", "could not find support file %s", + pszDictFile ); + return -1; + } + + /* For saving space, the file can be compressed (gz) */ + if ( strstr( pszDictFile,".gz") != NULL ) + sprintf( szTemp, "/vsigzip/%s", pszFilename); + else + sprintf( szTemp, "%s", pszFilename); + + fp = VSIFOpenL( szTemp, "rb" ); + + if( fp == NULL ) + { + CPLDebug( "OSRSearchDictForWKT()", "could not open support file %s", + pszFilename ); + return -1; + } + +/* -------------------------------------------------------------------- */ +/* Process lines. */ +/* -------------------------------------------------------------------- */ + const char *pszLine; + + while( (pszLine = CPLReadLine2L(fp,-1,NULL)) != NULL ) + { + /* do nothing */; + if( pszLine[0] == '#' ) + continue; + + pszTemp = strstr(pszLine,","); + if (pszTemp) + { + nPos = pszTemp - pszLine; + + /* do nothing */; + if ( nPos == 0 ) + continue; + + strncpy( szTemp, pszLine, nPos ); + szTemp[nPos] = '\0'; + nCode = atoi(szTemp); + + pszTemp = (char *) pszLine + nPos +1; + + if ( EQUAL(pszWKT,pszTemp) ) + { + nFound = nCode; + CPLDebug( "OSRSearchDictForWKT()", "found EPSG:%d WKT=%s", + nCode, pszWKT ); + break; + } + } + } + + if ( fp ) VSIFCloseL( fp ); + + return nFound; + +} + +/************************************************************************/ +/* OSRSearchDictForWKT() */ +/************************************************************************/ +/** + * This function is the same as the C++ method + * OGRSpatialReference::SearchDictForWKT() + */ + +int OSRSearchDictForWKT( const char *pszDictFile, const char *pszWKT ) + +{ + return OGRSpatialReference::SearchDictForWKT( pszDictFile, + pszWKT ); +} diff --git a/gdal/ogr/ogr_srs_esri.cpp b/gdal/ogr/ogr_srs_esri.cpp index 5655668..1e6e6db 100644 --- a/gdal/ogr/ogr_srs_esri.cpp +++ b/gdal/ogr/ogr_srs_esri.cpp @@ -657,6 +657,46 @@ OGRErr OGRSpatialReference::importFromESRI( char **papszPrj ) if( eErr == OGRERR_NONE ) eErr = morphFromESRI(); + if( eErr == OGRERR_NONE ) + { + if ( EQUAL(CPLGetConfigOption("OGR_ESRI_FIND_EPSG","NO"),"YES" ) ) + { + int nCode = -1; + nCode = this->FindEPSG(); + if ( nCode != -1 ) + { + OGRSpatialReference oSRSTemp; + char *pszWKT3; + this->exportToWkt( &pszWKT3 ); + + /* first try import and morph on oSRSTemp, if works do it on this */ + if( oSRSTemp.importFromEPSG( nCode ) == OGRERR_NONE ) + { + if( oSRSTemp.morphFromESRI() == OGRERR_NONE ) + { + eErr = this->morphFromESRI(); + if( eErr == OGRERR_NONE ) + eErr = this->importFromEPSG( nCode ); + /* if something went wrong, restore original */ + if ( eErr != OGRERR_NONE ) + { + char *pszWKT4 = pszWKT3; + eErr = importFromWkt( &pszWKT4 ); + } + else + { + CPLDebug( "OGR_ESRI", + "Updated spatial reference from EPSG code %d", + nCode ); + } + } + } + + CPLFree( pszWKT3 ); + } + } + } + return eErr; }