EVOLUTION-MANAGER

Edit File: ogr_fromepsg.cpp

/****************************************************************************** * * Project: OpenGIS Simple Features Reference Implementation * Purpose: Generate an OGRSpatialReference object based on an EPSG * PROJCS, or GEOGCS code. * Author: Frank Warmerdam, warmerdam@pobox.com * ****************************************************************************** * Copyright (c) 2000, Frank Warmerdam * Copyright (c) 2008-2013, Even Rouault <even dot rouault at mines-paris dot org> * * 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 "cpl_port.h" #include "ogr_srs_api.h" #include <ctype.h> #include <cmath> #include <cstddef> #include <cstdio> #include <cstring> #include <cstdlib> #include <map> #include <memory> #include <string> #include <vector> #include <limits> #include "cpl_conv.h" #include "cpl_csv.h" #include "cpl_error.h" #include "cpl_multiproc.h" #include "cpl_string.h" #include "ogr_core.h" #include "ogr_p.h" #include "ogr_spatialref.h" CPL_CVSID("$Id: ogr_fromepsg.cpp a5d470bdd8e89e3146e962ff5909fa30e2cb6cdf 2018-04-18 22:09:11 +0200 Even Rouault $") extern void OGRsnPrintDouble( char * pszStrBuf, size_t size, double dfValue ); int EPSGGetWGS84Transform( int nGeogCS, std::vector<CPLString>& asTransform ); void OGREPSGDatumNameMassage( char ** ppszDatum ); void CleanupFindMatchesCacheAndMutex(); static const char * const apszDatumEquiv[] = { "Militar_Geographische_Institut", "Militar_Geographische_Institute", "World_Geodetic_System_1984", "WGS_1984", "WGS_72_Transit_Broadcast_Ephemeris", "WGS_1972_Transit_Broadcast_Ephemeris", "World_Geodetic_System_1972", "WGS_1972", "European_Terrestrial_Reference_System_89", "European_Reference_System_1989", nullptr }; static CPLMutex* hFindMatchesMutex = nullptr; static std::vector<OGRSpatialReference*>* papoSRSCache_PROJCS = nullptr; static std::vector<OGRSpatialReference*>* papoSRSCache_GEOGCS = nullptr; static std::map<CPLString, int>* poMapESRIPROJCSNameToEPSGCode = nullptr; static std::map<CPLString, int>* poMapESRIGEOGCSNameToEPSGCode = nullptr; /************************************************************************/ /* OGREPSGDatumNameMassage() */ /* */ /* Massage an EPSG datum name into WMT format. Also transform */ /* specific exception cases into WKT versions. */ /************************************************************************/ void OGREPSGDatumNameMassage( char ** ppszDatum ) { char *pszDatum = *ppszDatum; if( pszDatum[0] == '\0' ) return; /* -------------------------------------------------------------------- */ /* Translate non-alphanumeric values to underscores. */ /* -------------------------------------------------------------------- */ for( int i = 0; pszDatum[i] != '\0'; i++ ) { if( pszDatum[i] != '+' && !(pszDatum[i] >= 'A' && pszDatum[i] <= 'Z') && !(pszDatum[i] >= 'a' && pszDatum[i] <= 'z') && !(pszDatum[i] >= '0' && pszDatum[i] <= '9') ) { pszDatum[i] = '_'; } } /* -------------------------------------------------------------------- */ /* Remove repeated and trailing underscores. */ /* -------------------------------------------------------------------- */ int j = 0; // Used after for loop. for( int i = 1; pszDatum[i] != '\0'; i++ ) { if( pszDatum[j] == '_' && pszDatum[i] == '_' ) continue; pszDatum[++j] = pszDatum[i]; } if( pszDatum[j] == '_' ) pszDatum[j] = '\0'; else pszDatum[j+1] = '\0'; /* -------------------------------------------------------------------- */ /* Search for datum equivelences. Specific massaged names get */ /* mapped to OpenGIS specified names. */ /* -------------------------------------------------------------------- */ for( int i = 0; apszDatumEquiv[i] != nullptr; i += 2 ) { if( EQUAL(*ppszDatum, apszDatumEquiv[i]) ) { CPLFree( *ppszDatum ); *ppszDatum = CPLStrdup( apszDatumEquiv[i+1] ); break; } } } /************************************************************************/ /* EPSGAngleStringToDD() */ /* */ /* Convert an angle in the specified units to decimal degrees. */ /************************************************************************/ static double EPSGAngleStringToDD( const char * pszAngle, int nUOMAngle ) { double dfAngle = 0.0; if( nUOMAngle == 9110 ) // DDD.MMSSsss { dfAngle = std::abs(atoi(pszAngle)); const char *pszDecimal = strchr(pszAngle, '.'); if( pszDecimal != nullptr && strlen(pszDecimal) > 1 ) { char szMinutes[3] = { '\0', '\0', '\0' }; szMinutes[0] = pszDecimal[1]; if( pszDecimal[2] >= '0' && pszDecimal[2] <= '9' ) szMinutes[1] = pszDecimal[2]; else szMinutes[1] = '0'; dfAngle += atoi(szMinutes) / 60.0; if( strlen(pszDecimal) > 3 ) { char szSeconds[64] = { '\0' }; szSeconds[0] = pszDecimal[3]; if( pszDecimal[4] >= '0' && pszDecimal[4] <= '9' ) { szSeconds[1] = pszDecimal[4]; szSeconds[2] = '.'; strncpy( szSeconds+3, pszDecimal + 5, sizeof(szSeconds)-3-1 ); szSeconds[sizeof(szSeconds)-1] = 0; } else { szSeconds[1] = '0'; szSeconds[2] = '\0'; } dfAngle += CPLAtof(szSeconds) / 3600.0; } } if( pszAngle[0] == '-' ) dfAngle *= -1; } else if( nUOMAngle == 9105 || nUOMAngle == 9106 ) // Grad. { dfAngle = 180.0 * (CPLAtof(pszAngle ) / 200.0); } else if( nUOMAngle == 9101 ) // Radians. { dfAngle = 180.0 * (CPLAtof(pszAngle ) / M_PI); } else if( nUOMAngle == 9103 ) // Arc-minute. { dfAngle = CPLAtof(pszAngle) / 60.0; } else if( nUOMAngle == 9104 ) // Arc-second. { dfAngle = CPLAtof(pszAngle) / 3600.0; } else // Decimal degrees. Some cases missing, but seemingly never used. { CPLAssert( nUOMAngle == 9102 || nUOMAngle == 0 ); dfAngle = CPLAtof(pszAngle ); } return dfAngle; } /************************************************************************/ /* EPSGGetUOMAngleInfo() */ /************************************************************************/ static bool EPSGGetUOMAngleInfo( int nUOMAngleCode, char **ppszUOMName, double * pdfInDegrees ) { // We do a special override of some of the DMS formats name // This will also solve accuracy problems when computing // the dfInDegree value from the CSV values (#3643). if( nUOMAngleCode == 9102 || nUOMAngleCode == 9107 || nUOMAngleCode == 9108 || nUOMAngleCode == 9110 || nUOMAngleCode == 9122 ) { if( ppszUOMName != nullptr ) *ppszUOMName = CPLStrdup("degree"); if( pdfInDegrees != nullptr ) *pdfInDegrees = 1.0; return true; } const char *pszFilename = CSVFilename( "unit_of_measure.csv" ); char szSearchKey[24] = { '\0' }; snprintf( szSearchKey, sizeof(szSearchKey), "%d", nUOMAngleCode ); const char *pszUOMName = CSVGetField( pszFilename, "UOM_CODE", szSearchKey, CC_Integer, "UNIT_OF_MEAS_NAME" ); /* -------------------------------------------------------------------- */ /* If the file is found, read from there. Note that FactorC is */ /* an empty field for any of the DMS style formats, and in this */ /* case we really want to return the default InDegrees value */ /* (1.0) from above. */ /* -------------------------------------------------------------------- */ double dfInDegrees = 1.0; if( !EQUAL( pszUOMName, "" ) ) { const double dfFactorB = CPLAtof(CSVGetField( pszFilename, "UOM_CODE", szSearchKey, CC_Integer, "FACTOR_B" )); const double dfFactorC = CPLAtof(CSVGetField( pszFilename, "UOM_CODE", szSearchKey, CC_Integer, "FACTOR_C" )); if( dfFactorC != 0.0 ) dfInDegrees = (dfFactorB / dfFactorC) * (180.0 / M_PI); // For some reason, (FactorB) is not very precise in EPSG, use // a more exact form for grads. if( nUOMAngleCode == 9105 ) dfInDegrees = 180.0 / 200.0; } /* -------------------------------------------------------------------- */ /* Otherwise handle a few well known units directly. */ /* -------------------------------------------------------------------- */ else { switch( nUOMAngleCode ) { case 9101: pszUOMName = "radian"; dfInDegrees = 180.0 / M_PI; break; case 9103: pszUOMName = "arc-minute"; dfInDegrees = 1 / 60.0; break; case 9104: pszUOMName = "arc-second"; dfInDegrees = 1 / 3600.0; break; case 9105: pszUOMName = "grad"; dfInDegrees = 180.0 / 200.0; break; case 9106: pszUOMName = "gon"; dfInDegrees = 180.0 / 200.0; break; case 9109: pszUOMName = "microradian"; dfInDegrees = 180.0 / (M_PI * 1000000.0); break; default: return false; } } /* -------------------------------------------------------------------- */ /* Return to caller. */ /* -------------------------------------------------------------------- */ if( ppszUOMName != nullptr ) *ppszUOMName = CPLStrdup( pszUOMName ); if( pdfInDegrees != nullptr ) *pdfInDegrees = dfInDegrees; return true; } /************************************************************************/ /* EPSGGetUOMLengthInfo() */ /* */ /* Note: This function should eventually also know how to */ /* lookup length aliases in the UOM_LE_ALIAS table. */ /************************************************************************/ static bool EPSGGetUOMLengthInfo( int nUOMLengthCode, char **ppszUOMName, double * pdfInMeters ) { /* -------------------------------------------------------------------- */ /* We short cut meter to save work in the most common case. */ /* -------------------------------------------------------------------- */ if( nUOMLengthCode == 9001 ) { if( ppszUOMName != nullptr ) *ppszUOMName = CPLStrdup( "metre" ); if( pdfInMeters != nullptr ) *pdfInMeters = 1.0; return true; } /* -------------------------------------------------------------------- */ /* Search the units database for this unit. If we don't find */ /* it return failure. */ /* -------------------------------------------------------------------- */ const char *uom_filename = CSVFilename( "unit_of_measure.csv" ); char szSearchKey[24] = { '\0' }; snprintf( szSearchKey, sizeof(szSearchKey), "%d", nUOMLengthCode ); char **papszUnitsRecord = CSVScanFileByName( uom_filename, "UOM_CODE", szSearchKey, CC_Integer ); if( papszUnitsRecord == nullptr ) return false; /* -------------------------------------------------------------------- */ /* Get the name, if requested. */ /* -------------------------------------------------------------------- */ if( ppszUOMName != nullptr ) { const int iNameField = CSVGetFileFieldId( uom_filename, "UNIT_OF_MEAS_NAME" ); *ppszUOMName = CPLStrdup( CSLGetField(papszUnitsRecord, iNameField) ); } /* -------------------------------------------------------------------- */ /* Get the A and B factor fields, and create the multiplicative */ /* factor. */ /* -------------------------------------------------------------------- */ if( pdfInMeters != nullptr ) { const int iBFactorField = CSVGetFileFieldId( uom_filename, "FACTOR_B" ); const int iCFactorField = CSVGetFileFieldId( uom_filename, "FACTOR_C" ); if( CPLAtof(CSLGetField(papszUnitsRecord, iCFactorField)) > 0.0 ) *pdfInMeters = CPLAtof(CSLGetField(papszUnitsRecord, iBFactorField)) / CPLAtof(CSLGetField(papszUnitsRecord, iCFactorField)); else *pdfInMeters = 0.0; } return true; } /************************************************************************/ /* EPSGNegateString() */ /************************************************************************/ static void EPSGNegateString(CPLString& osValue) { if( osValue.compare("0") == 0 ) return; if( osValue[0] == '-' ) { osValue = osValue.substr(1); return; } if( osValue[0] == '+' ) { osValue[0] = '-'; return; } osValue = "-" + osValue; } /************************************************************************/ /* EPSGGetWGS84Transform() */ /* */ /* The following code attempts to find a bursa-wolf */ /* transformation from this GeogCS to WGS84 (4326). */ /* */ /* Faults: */ /* o I think there are codes other than 9603 and 9607 that */ /* return compatible, or easily transformed parameters. */ /* o Only the first path from the given GeogCS is checked due */ /* to limitations in the CSV API. */ /************************************************************************/ int EPSGGetWGS84Transform( int nGeogCS, std::vector<CPLString>& asTransform ) { /* -------------------------------------------------------------------- */ /* Fetch the line from the GCS table. */ /* -------------------------------------------------------------------- */ const char *pszFilename = CSVFilename("gcs.override.csv"); char szCode[32] = { '\0' }; snprintf( szCode, sizeof(szCode), "%d", nGeogCS ); char **papszLine = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE", szCode, CC_Integer ); if( papszLine == nullptr ) { pszFilename = CSVFilename("gcs.csv"); snprintf( szCode, sizeof(szCode), "%d", nGeogCS ); papszLine = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE", szCode, CC_Integer ); } if( papszLine == nullptr ) return FALSE; /* -------------------------------------------------------------------- */ /* Verify that the method code is one of our accepted ones. */ /* -------------------------------------------------------------------- */ const int nMethodCode = atoi(CSLGetField( papszLine, CSVGetFileFieldId(pszFilename, "COORD_OP_METHOD_CODE"))); if( nMethodCode != 9603 && nMethodCode != 9607 && nMethodCode != 9606 ) return FALSE; /* -------------------------------------------------------------------- */ /* Fetch the transformation parameters. */ /* -------------------------------------------------------------------- */ const int iDXField = CSVGetFileFieldId(pszFilename, "DX"); if( iDXField < 0 || CSLCount(papszLine) < iDXField + 7 ) return FALSE; asTransform.resize(0); for( int iField = 0; iField < 7; iField++ ) { const char* pszValue = papszLine[iDXField+iField]; if( pszValue[0] ) asTransform.push_back(pszValue); else asTransform.push_back("0"); } /* -------------------------------------------------------------------- */ /* 9607 - coordinate frame rotation has reverse signs on the */ /* rotational coefficients. Fix up now since we internal */ /* operate according to method 9606 (position vector 7-parameter). */ /* -------------------------------------------------------------------- */ if( nMethodCode == 9607 ) { EPSGNegateString(asTransform[3]); EPSGNegateString(asTransform[4]); EPSGNegateString(asTransform[5]); } return TRUE; } /************************************************************************/ /* EPSGGetPMInfo() */ /* */ /* Get the offset between a given prime meridian and Greenwich */ /* in degrees. */ /************************************************************************/ static bool EPSGGetPMInfo( int nPMCode, char ** ppszName, double *pdfOffset ) { /* -------------------------------------------------------------------- */ /* Use a special short cut for Greenwich, since it is so common. */ /* -------------------------------------------------------------------- */ // FIXME? Where does 7022 come from ? Let's keep it just in case // 8901 is the official current code for Greenwich. if( nPMCode == 7022 /* PM_Greenwich */ || nPMCode == 8901 ) { if( pdfOffset != nullptr ) *pdfOffset = 0.0; if( ppszName != nullptr ) *ppszName = CPLStrdup( "Greenwich" ); return true; } /* -------------------------------------------------------------------- */ /* Search the database for the corresponding datum code. */ /* -------------------------------------------------------------------- */ char szSearchKey[24] = { '\0' }; snprintf( szSearchKey, sizeof(szSearchKey), "%d", nPMCode ); const char *PM_FILENAME = CSVFilename("prime_meridian.csv"); const int nUOMAngle = atoi(CSVGetField( PM_FILENAME, "PRIME_MERIDIAN_CODE", szSearchKey, CC_Integer, "UOM_CODE" ) ); if( nUOMAngle < 1 ) return false; /* -------------------------------------------------------------------- */ /* Get the PM offset. */ /* -------------------------------------------------------------------- */ if( pdfOffset != nullptr ) { *pdfOffset = EPSGAngleStringToDD( CSVGetField( PM_FILENAME, "PRIME_MERIDIAN_CODE", szSearchKey, CC_Integer, "GREENWICH_LONGITUDE" ), nUOMAngle ); } /* -------------------------------------------------------------------- */ /* Get the name, if requested. */ /* -------------------------------------------------------------------- */ if( ppszName != nullptr ) *ppszName = CPLStrdup( CSVGetField( PM_FILENAME, "PRIME_MERIDIAN_CODE", szSearchKey, CC_Integer, "PRIME_MERIDIAN_NAME" )); return true; } /************************************************************************/ /* EPSGGetGCSInfo() */ /* */ /* Fetch the datum, and prime meridian related to a particular */ /* GCS. */ /************************************************************************/ static bool EPSGGetGCSInfo( int nGCSCode, char ** ppszName, int * pnDatum, char **ppszDatumName, int * pnPM, int *pnEllipsoid, int *pnUOMAngle, int * pnCoordSysCode ) { /* -------------------------------------------------------------------- */ /* Search the database for the corresponding datum code. */ /* -------------------------------------------------------------------- */ const char *pszFilename = CSVFilename("gcs.override.csv"); char szSearchKey[24] = { '\0' }; snprintf( szSearchKey, sizeof(szSearchKey), "%d", nGCSCode ); int nDatum = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer, "DATUM_CODE" ) ); if( nDatum < 1 ) { pszFilename = CSVFilename("gcs.csv"); snprintf( szSearchKey, sizeof(szSearchKey), "%d", nGCSCode ); nDatum = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer, "DATUM_CODE" ) ); } if( nDatum < 1 ) return false; if( pnDatum != nullptr ) *pnDatum = nDatum; /* -------------------------------------------------------------------- */ /* Get the PM. */ /* -------------------------------------------------------------------- */ const int nPM = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer, "PRIME_MERIDIAN_CODE" ) ); if( nPM < 1 ) return false; if( pnPM != nullptr ) *pnPM = nPM; /* -------------------------------------------------------------------- */ /* Get the Ellipsoid. */ /* -------------------------------------------------------------------- */ const int nEllipsoid = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer, "ELLIPSOID_CODE" ) ); if( nEllipsoid < 1 ) return false; if( pnEllipsoid != nullptr ) *pnEllipsoid = nEllipsoid; /* -------------------------------------------------------------------- */ /* Get the angular units. */ /* -------------------------------------------------------------------- */ const int nUOMAngle = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer, "UOM_CODE" ) ); if( nUOMAngle < 1 ) return false; if( pnUOMAngle != nullptr ) *pnUOMAngle = nUOMAngle; /* -------------------------------------------------------------------- */ /* Get the name, if requested. */ /* -------------------------------------------------------------------- */ if( ppszName != nullptr ) *ppszName = CPLStrdup(CSVGetField( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer, "COORD_REF_SYS_NAME" )); /* -------------------------------------------------------------------- */ /* Get the datum name, if requested. */ /* -------------------------------------------------------------------- */ if( ppszDatumName != nullptr ) *ppszDatumName = CPLStrdup(CSVGetField( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer, "DATUM_NAME" )); /* -------------------------------------------------------------------- */ /* Get the CoordSysCode */ /* -------------------------------------------------------------------- */ const int nCSC = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer, "COORD_SYS_CODE" ) ); if( pnCoordSysCode != nullptr ) *pnCoordSysCode = nCSC; return TRUE; } /************************************************************************/ /* OSRGetEllipsoidInfo() */ /************************************************************************/ /** * Fetch info about an ellipsoid. * * This helper function will return ellipsoid parameters corresponding to EPSG * code provided. Axes are always returned in meters. Semi major computed * based on inverse flattening where that is provided. * * @param nCode EPSG code of the requested ellipsoid * * @param ppszName pointer to string where ellipsoid name will be returned. It * is caller responsibility to free this string after using with CPLFree(). * * @param pdfSemiMajor pointer to variable where semi major axis will be * returned. * * @param pdfInvFlattening pointer to variable where inverse flattening will * be returned. * * @return OGRERR_NONE on success or an error code in case of failure. **/ OGRErr OSRGetEllipsoidInfo( int nCode, char ** ppszName, double * pdfSemiMajor, double * pdfInvFlattening ) { /* -------------------------------------------------------------------- */ /* Get the semi major axis. */ /* -------------------------------------------------------------------- */ char szSearchKey[24] = { '\0' }; snprintf( szSearchKey, sizeof(szSearchKey), "%d", nCode ); szSearchKey[sizeof(szSearchKey) - 1] = '\n'; double dfSemiMajor = CPLAtof(CSVGetField( CSVFilename("ellipsoid.csv" ), "ELLIPSOID_CODE", szSearchKey, CC_Integer, "SEMI_MAJOR_AXIS" ) ); if( dfSemiMajor == 0.0 ) return OGRERR_UNSUPPORTED_SRS; /* -------------------------------------------------------------------- */ /* Get the translation factor into meters. */ /* -------------------------------------------------------------------- */ const int nUOMLength = atoi(CSVGetField( CSVFilename("ellipsoid.csv" ), "ELLIPSOID_CODE", szSearchKey, CC_Integer, "UOM_CODE" )); double dfToMeters = 1.0; if( !EPSGGetUOMLengthInfo( nUOMLength, nullptr, &dfToMeters ) ) { dfToMeters = 1.0; } dfSemiMajor *= dfToMeters; if( pdfSemiMajor != nullptr ) *pdfSemiMajor = dfSemiMajor; /* -------------------------------------------------------------------- */ /* Get the semi-minor if requested. If the Semi-minor axis */ /* isn't available, compute it based on the inverse flattening. */ /* -------------------------------------------------------------------- */ if( pdfInvFlattening != nullptr ) { *pdfInvFlattening = CPLAtof(CSVGetField( CSVFilename("ellipsoid.csv" ), "ELLIPSOID_CODE", szSearchKey, CC_Integer, "INV_FLATTENING" )); if( *pdfInvFlattening == 0.0 ) { const double dfSemiMinor = CPLAtof(CSVGetField( CSVFilename("ellipsoid.csv" ), "ELLIPSOID_CODE", szSearchKey, CC_Integer, "SEMI_MINOR_AXIS" )) * dfToMeters; if( dfSemiMajor == 0.0 ) *pdfInvFlattening = 0.0; else *pdfInvFlattening = OSRCalcInvFlattening(dfSemiMajor, dfSemiMinor); } } /* -------------------------------------------------------------------- */ /* Get the name, if requested. */ /* -------------------------------------------------------------------- */ if( ppszName != nullptr ) *ppszName = CPLStrdup(CSVGetField( CSVFilename("ellipsoid.csv" ), "ELLIPSOID_CODE", szSearchKey, CC_Integer, "ELLIPSOID_NAME" )); return OGRERR_NONE; } constexpr int CoLatConeAxis = 1036; // See #4223. constexpr int NatOriginLat = 8801; constexpr int NatOriginLong = 8802; constexpr int NatOriginScaleFactor = 8805; constexpr int FalseEasting = 8806; constexpr int FalseNorthing = 8807; constexpr int ProjCenterLat = 8811; constexpr int ProjCenterLong = 8812; constexpr int Azimuth = 8813; constexpr int AngleRectifiedToSkewedGrid = 8814; constexpr int InitialLineScaleFactor = 8815; constexpr int ProjCenterEasting = 8816; constexpr int ProjCenterNorthing = 8817; constexpr int PseudoStdParallelLat = 8818; constexpr int PseudoStdParallelScaleFactor = 8819; constexpr int FalseOriginLat = 8821; constexpr int FalseOriginLong = 8822; constexpr int StdParallel1Lat = 8823; constexpr int StdParallel2Lat = 8824; constexpr int FalseOriginEasting = 8826; constexpr int FalseOriginNorthing = 8827; constexpr int SphericalOriginLat = 8828; constexpr int SphericalOriginLong = 8829; #if 0 constexpr int InitialLongitude = 8830; constexpr int ZoneWidth = 8831; #endif constexpr int PolarLatStdParallel = 8832; constexpr int PolarLongOrigin = 8833; /************************************************************************/ /* EPSGGetProjTRFInfo() */ /* */ /* Transform a PROJECTION_TRF_CODE into a projection method, */ /* and a set of parameters. The parameters identify will */ /* depend on the returned method, but they will all have been */ /* normalized into degrees and meters. */ /************************************************************************/ static bool EPSGGetProjTRFInfo( int nPCS, int * pnProjMethod, int *panParmIds, double * padfProjParms ) { /* -------------------------------------------------------------------- */ /* Get the proj method. If this fails to return a meaningful */ /* number, then the whole function fails. */ /* -------------------------------------------------------------------- */ CPLString osFilename = CSVFilename( "pcs.override.csv" ); char szTRFCode[16] = { '\0' }; snprintf( szTRFCode, sizeof(szTRFCode), "%d", nPCS ); int nProjMethod = atoi( CSVGetField( osFilename, "COORD_REF_SYS_CODE", szTRFCode, CC_Integer, "COORD_OP_METHOD_CODE" ) ); if( nProjMethod == 0 ) { osFilename = CSVFilename( "pcs.csv" ); snprintf( szTRFCode, sizeof(szTRFCode), "%d", nPCS ); nProjMethod = atoi( CSVGetField( osFilename, "COORD_REF_SYS_CODE", szTRFCode, CC_Integer, "COORD_OP_METHOD_CODE" ) ); if( nProjMethod == 0 ) return false; } /* -------------------------------------------------------------------- */ /* Get the parameters for this projection. */ /* -------------------------------------------------------------------- */ double adfProjParms[7] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; for( int i = 0; i < 7; i++ ) { if( panParmIds == nullptr ) { CPLError( CE_Failure, CPLE_AppDefined, "panParmIds cannot be NULL." ); return false; } char szParamUOMID[32] = { '\0' }; char szParamValueID[32] = { '\0' }; char szParamCodeID[32] = { '\0' }; snprintf( szParamCodeID, sizeof(szParamCodeID), "PARAMETER_CODE_%d", i+1 ); snprintf( szParamUOMID, sizeof(szParamUOMID), "PARAMETER_UOM_%d", i+1 ); snprintf( szParamValueID, sizeof(szParamValueID), "PARAMETER_VALUE_%d", i+1 ); panParmIds[i] = atoi(CSVGetField( osFilename, "COORD_REF_SYS_CODE", szTRFCode, CC_Integer, szParamCodeID )); int nUOM = atoi(CSVGetField( osFilename, "COORD_REF_SYS_CODE", szTRFCode, CC_Integer, szParamUOMID )); char *pszValue = CPLStrdup( CSVGetField( osFilename, "COORD_REF_SYS_CODE", szTRFCode, CC_Integer, szParamValueID )); // there is a bug in the EPSG 6.2.2 database for PCS 2935 and 2936 // such that they have foot units for the scale factor. Avoid this. if( (panParmIds[i] == NatOriginScaleFactor || panParmIds[i] == InitialLineScaleFactor || panParmIds[i] == PseudoStdParallelScaleFactor) && nUOM < 9200 ) nUOM = 9201; if( nUOM >= 9100 && nUOM < 9200 ) { adfProjParms[i] = EPSGAngleStringToDD( pszValue, nUOM ); } else if( nUOM > 9000 && nUOM < 9100 ) { double dfInMeters = 0.0; if( !EPSGGetUOMLengthInfo( nUOM, nullptr, &dfInMeters ) ) dfInMeters = 1.0; adfProjParms[i] = CPLAtof(pszValue) * dfInMeters; } else if( EQUAL(pszValue, "") ) // Null field. { adfProjParms[i] = 0.0; } else // Really, should consider looking up other scaling factors. { if( nUOM != 9201 ) CPLDebug( "OGR", "Non-unity scale factor units! (UOM=%d, PCS=%d)", nUOM, nPCS ); adfProjParms[i] = CPLAtof(pszValue); } CPLFree( pszValue ); } /* -------------------------------------------------------------------- */ /* Transfer requested data into passed variables. */ /* -------------------------------------------------------------------- */ if( pnProjMethod != nullptr ) *pnProjMethod = nProjMethod; if( padfProjParms != nullptr ) { for( int i = 0; i < 7; i++ ) padfProjParms[i] = adfProjParms[i]; } return true; } /************************************************************************/ /* EPSGGetPCSInfo() */ /************************************************************************/ static bool EPSGGetPCSInfo( int nPCSCode, char **ppszEPSGName, int *pnUOMLengthCode, int *pnUOMAngleCode, int *pnGeogCS, int *pnTRFCode, int *pnCoordSysCode, double* adfTOWGS84 ) { /* -------------------------------------------------------------------- */ /* Search the units database for this unit. If we don't find */ /* it return failure. */ /* -------------------------------------------------------------------- */ const char *pszFilename = CSVFilename( "pcs.override.csv" ); char szSearchKey[24] = { '\0' }; snprintf( szSearchKey, sizeof(szSearchKey), "%d", nPCSCode ); char **papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer ); if( papszRecord == nullptr ) { pszFilename = CSVFilename( "pcs.csv" ); snprintf( szSearchKey, sizeof(szSearchKey), "%d", nPCSCode ); papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer ); } if( papszRecord == nullptr ) return false; /* -------------------------------------------------------------------- */ /* Get the name, if requested. */ /* -------------------------------------------------------------------- */ if( ppszEPSGName != nullptr ) { CPLString osPCSName = CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "COORD_REF_SYS_NAME")); const char *pszDeprecated = CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "DEPRECATED") ); if( pszDeprecated != nullptr && *pszDeprecated == '1' ) osPCSName += " (deprecated)"; *ppszEPSGName = CPLStrdup(osPCSName); } /* -------------------------------------------------------------------- */ /* Get the UOM Length code, if requested. */ /* -------------------------------------------------------------------- */ if( pnUOMLengthCode != nullptr ) { const char *pszValue = CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "UOM_CODE")); if( atoi(pszValue) > 0 ) *pnUOMLengthCode = atoi(pszValue); else *pnUOMLengthCode = 0; } /* -------------------------------------------------------------------- */ /* Get the UOM Angle code, if requested. */ /* -------------------------------------------------------------------- */ if( pnUOMAngleCode != nullptr ) { const char *pszValue = CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "UOM_ANGLE_CODE") ); if( atoi(pszValue) > 0 ) *pnUOMAngleCode = atoi(pszValue); else *pnUOMAngleCode = 0; } /* -------------------------------------------------------------------- */ /* Get the GeogCS (Datum with PM) code, if requested. */ /* -------------------------------------------------------------------- */ if( pnGeogCS != nullptr ) { const char *pszValue = CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "SOURCE_GEOGCRS_CODE")); if( atoi(pszValue) > 0 ) *pnGeogCS = atoi(pszValue); else *pnGeogCS = 0; } /* -------------------------------------------------------------------- */ /* Get the GeogCS (Datum with PM) code, if requested. */ /* -------------------------------------------------------------------- */ if( pnTRFCode != nullptr ) { const char *pszValue = CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "COORD_OP_CODE")); if( atoi(pszValue) > 0 ) *pnTRFCode = atoi(pszValue); else *pnTRFCode = 0; } /* -------------------------------------------------------------------- */ /* Get the CoordSysCode */ /* -------------------------------------------------------------------- */ const int nCSC = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer, "COORD_SYS_CODE" ) ); if( pnCoordSysCode != nullptr ) *pnCoordSysCode = nCSC; /* -------------------------------------------------------------------- */ /* Get the TOWGS84 (override) parameters */ /* -------------------------------------------------------------------- */ const char *pszDX = CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "DX")); if( pszDX[0] != '\0' ) { adfTOWGS84[0] = CPLAtof(pszDX); adfTOWGS84[1] = CPLAtof(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "DY"))); adfTOWGS84[2] = CPLAtof(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "DZ"))); adfTOWGS84[3] = CPLAtof(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "RX"))); adfTOWGS84[4] = CPLAtof(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "RY"))); adfTOWGS84[5] = CPLAtof(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "RZ"))); adfTOWGS84[6] = CPLAtof(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "DS"))); } return true; } /************************************************************************/ /* SetEPSGAxisInfo() */ /************************************************************************/ static OGRErr SetEPSGAxisInfo( OGRSpatialReference *poSRS, const char *pszTargetKey, int nCoordSysCode ) { /* -------------------------------------------------------------------- */ /* Special cases for well known and common values. We short */ /* circuit these to save time doing file lookups. */ /* -------------------------------------------------------------------- */ // Conventional and common Easting/Northing values. if( nCoordSysCode >= 4400 && nCoordSysCode <= 4410 ) { return poSRS->SetAxes( pszTargetKey, "Easting", OAO_East, "Northing", OAO_North ); } // Conventional and common Easting/Northing values. if( nCoordSysCode >= 6400 && nCoordSysCode <= 6423 ) { return poSRS->SetAxes( pszTargetKey, "Latitude", OAO_North, "Longitude", OAO_East ); } /* -------------------------------------------------------------------- */ /* Get the definition from the coordinate_axis.csv file. */ /* -------------------------------------------------------------------- */ char szSearchKey[24] = { '\0' }; const char *pszFilename = CSVFilename( "coordinate_axis.csv" ); snprintf( szSearchKey, sizeof(szSearchKey), "%d", nCoordSysCode ); char **papszRecord = CSVScanFileByName( pszFilename, "COORD_SYS_CODE", szSearchKey, CC_Integer ); char **papszAxis1 = nullptr; char **papszAxis2 = nullptr; if( papszRecord != nullptr ) { papszAxis1 = CSLDuplicate( papszRecord ); papszRecord = CSVGetNextLine( pszFilename ); if( CSLCount(papszRecord) > 0 && EQUAL(papszRecord[0], papszAxis1[0]) ) papszAxis2 = CSLDuplicate( papszRecord ); } if( papszAxis2 == nullptr ) { CSLDestroy( papszAxis1 ); CPLError( CE_Failure, CPLE_AppDefined, "Failed to find entries for COORD_SYS_CODE %d " "in coordinate_axis.csv", nCoordSysCode ); return OGRERR_FAILURE; } /* -------------------------------------------------------------------- */ /* Confirm the records are complete, and work out which columns */ /* are which. */ /* -------------------------------------------------------------------- */ const int iAxisOrientationField = CSVGetFileFieldId( pszFilename, "coord_axis_orientation" ); const int iAxisAbbrevField = CSVGetFileFieldId( pszFilename, "coord_axis_abbreviation" ); const int iAxisOrderField = CSVGetFileFieldId( pszFilename, "coord_axis_order" ); const int iAxisNameCodeField = CSVGetFileFieldId( pszFilename, "coord_axis_name_code" ); // Check that all fields are available and that the axis_order field // is the one with highest index. if( !(iAxisOrientationField >= 0 && iAxisOrientationField < iAxisOrderField && iAxisAbbrevField >= 0 && iAxisAbbrevField < iAxisOrderField && iAxisOrderField >= 0 && iAxisNameCodeField >= 0 && iAxisNameCodeField < iAxisOrderField) ) { CSLDestroy( papszAxis1 ); CSLDestroy( papszAxis2 ); CPLError( CE_Failure, CPLE_AppDefined, "coordinate_axis.csv corrupted" ); return OGRERR_FAILURE; } if( CSLCount(papszAxis1) < iAxisOrderField+1 || CSLCount(papszAxis2) < iAxisOrderField+1 ) { CSLDestroy( papszAxis1 ); CSLDestroy( papszAxis2 ); CPLError( CE_Failure, CPLE_AppDefined, "Axis records appear incomplete for COORD_SYS_CODE %d " "in coordinate_axis.csv", nCoordSysCode ); return OGRERR_FAILURE; } /* -------------------------------------------------------------------- */ /* Do we need to switch the axes around? */ /* -------------------------------------------------------------------- */ if( atoi(papszAxis2[iAxisOrderField]) < atoi(papszAxis1[iAxisOrderField]) ) { papszRecord = papszAxis1; papszAxis1 = papszAxis2; papszAxis2 = papszRecord; } /* -------------------------------------------------------------------- */ /* Work out axis enumeration values. */ /* -------------------------------------------------------------------- */ OGRAxisOrientation eOAxis1 = OAO_Other; OGRAxisOrientation eOAxis2 = OAO_Other; constexpr int anCodes[7] = { -1, 9907, 9909, 9906, 9908, -1, -1 }; for( int iAO = 0; iAO < 7; iAO++ ) { const OGRAxisOrientation eAO = static_cast<OGRAxisOrientation>(iAO); if( EQUAL(papszAxis1[iAxisOrientationField], OSRAxisEnumToName(eAO)) ) eOAxis1 = eAO; if( EQUAL(papszAxis2[iAxisOrientationField], OSRAxisEnumToName(eAO)) ) eOAxis2 = eAO; if( eOAxis1 == OAO_Other && anCodes[iAO] == atoi(papszAxis1[iAxisNameCodeField]) ) eOAxis1 = eAO; if( eOAxis2 == OAO_Other && anCodes[iAO] == atoi(papszAxis2[iAxisNameCodeField]) ) eOAxis2 = eAO; } /* -------------------------------------------------------------------- */ /* Work out the axis name. We try to expand the abbreviation */ /* to a longer name. */ /* -------------------------------------------------------------------- */ const char *apszAxisName[2] = { papszAxis1[iAxisAbbrevField], papszAxis2[iAxisAbbrevField] }; for( int iAO = 0; iAO < 2; iAO++ ) { if( EQUAL(apszAxisName[iAO], "N") ) apszAxisName[iAO] = "Northing"; else if( EQUAL(apszAxisName[iAO], "E") ) apszAxisName[iAO] = "Easting"; else if( EQUAL(apszAxisName[iAO], "S") ) apszAxisName[iAO] = "Southing"; else if( EQUAL(apszAxisName[iAO], "W") ) apszAxisName[iAO] = "Westing"; } /* -------------------------------------------------------------------- */ /* Set the axes. */ /* -------------------------------------------------------------------- */ const OGRErr eResult = poSRS->SetAxes( pszTargetKey, apszAxisName[0], eOAxis1, apszAxisName[1], eOAxis2 ); CSLDestroy( papszAxis1 ); CSLDestroy( papszAxis2 ); return eResult; } /************************************************************************/ /* SetEPSGGeogCS() */ /* */ /* FLAWS: */ /* o Units are all hardcoded. */ /************************************************************************/ static OGRErr SetEPSGGeogCS( OGRSpatialReference * poSRS, int nGeogCS ) { int nDatumCode = 0; int nPMCode = 0; int nUOMAngle = 0; int nEllipsoidCode = 0; int nCSC = 0; char *pszGeogCSName = nullptr; char *pszDatumName = nullptr; char *pszAngleName = nullptr; if( !EPSGGetGCSInfo( nGeogCS, &pszGeogCSName, &nDatumCode, &pszDatumName, &nPMCode, &nEllipsoidCode, &nUOMAngle, &nCSC ) ) return OGRERR_UNSUPPORTED_SRS; char *pszPMName = nullptr; double dfPMOffset = 0.0; if( !EPSGGetPMInfo( nPMCode, &pszPMName, &dfPMOffset ) ) { CPLFree( pszDatumName ); CPLFree( pszGeogCSName ); return OGRERR_UNSUPPORTED_SRS; } OGREPSGDatumNameMassage( &pszDatumName ); char *pszEllipsoidName = nullptr; double dfSemiMajor = 0.0; double dfInvFlattening = 0.0; if( OSRGetEllipsoidInfo( nEllipsoidCode, &pszEllipsoidName, &dfSemiMajor, &dfInvFlattening ) != OGRERR_NONE ) { CPLFree( pszDatumName ); CPLFree( pszGeogCSName ); CPLFree( pszPMName ); return OGRERR_UNSUPPORTED_SRS; } double dfAngleInDegrees = 0.0; double dfAngleInRadians = 0.0; if( !EPSGGetUOMAngleInfo( nUOMAngle, &pszAngleName, &dfAngleInDegrees ) ) { pszAngleName = CPLStrdup("degree"); dfAngleInDegrees = 1.0; nUOMAngle = -1; } if( dfAngleInDegrees == 1.0 ) dfAngleInRadians = CPLAtof(SRS_UA_DEGREE_CONV); else dfAngleInRadians = CPLAtof(SRS_UA_DEGREE_CONV) * dfAngleInDegrees; poSRS->SetGeogCS( pszGeogCSName, pszDatumName, pszEllipsoidName, dfSemiMajor, dfInvFlattening, pszPMName, dfPMOffset, pszAngleName, dfAngleInRadians ); std::vector<CPLString> asBursaTransform; if( EPSGGetWGS84Transform( nGeogCS, asBursaTransform ) ) { OGR_SRSNode *poWGS84 = new OGR_SRSNode( "TOWGS84" ); for( int iCoeff = 0; iCoeff < 7; iCoeff++ ) { poWGS84->AddChild( new OGR_SRSNode( asBursaTransform[iCoeff].c_str() ) ); } poSRS->GetAttrNode( "DATUM" )->AddChild( poWGS84 ); } poSRS->SetAuthority( "GEOGCS", "EPSG", nGeogCS ); poSRS->SetAuthority( "DATUM", "EPSG", nDatumCode ); poSRS->SetAuthority( "SPHEROID", "EPSG", nEllipsoidCode ); poSRS->SetAuthority( "PRIMEM", "EPSG", nPMCode ); if( nUOMAngle > 0 ) poSRS->SetAuthority( "GEOGCS|UNIT", "EPSG", nUOMAngle ); CPLFree( pszAngleName ); CPLFree( pszDatumName ); CPLFree( pszEllipsoidName ); CPLFree( pszGeogCSName ); CPLFree( pszPMName ); /* -------------------------------------------------------------------- */ /* Set axes */ /* -------------------------------------------------------------------- */ if( nCSC > 0 ) { SetEPSGAxisInfo( poSRS, "GEOGCS", nCSC ); CPLErrorReset(); } return OGRERR_NONE; } /************************************************************************/ /* OGR_FetchParm() */ /* */ /* Fetch a parameter from the parm list, based on its EPSG */ /* parameter code. */ /************************************************************************/ static double OGR_FetchParm( double *padfProjParms, int *panParmIds, int nTargetId, double /* dfFromGreenwich */) { /* -------------------------------------------------------------------- */ /* Set default in meters/degrees. */ /* -------------------------------------------------------------------- */ double dfResult = 0.0; switch( nTargetId ) { case NatOriginScaleFactor: case InitialLineScaleFactor: case PseudoStdParallelScaleFactor: dfResult = 1.0; break; case AngleRectifiedToSkewedGrid: dfResult = 90.0; break; default: dfResult = 0.0; } /* -------------------------------------------------------------------- */ /* Try to find actual value in parameter list. */ /* -------------------------------------------------------------------- */ for( int i = 0; i < 7; i++ ) { if( panParmIds[i] == nTargetId ) { dfResult = padfProjParms[i]; break; } } /* -------------------------------------------------------------------- */ /* EPSG longitudes are relative to greenwich. The follow code */ /* could be used to make them relative to the prime meridian of */ /* the associated GCS if that was appropriate. However, the */ /* SetNormProjParm() method expects longitudes relative to */ /* greenwich, so there is nothing for us to do. */ /* -------------------------------------------------------------------- */ #if 0 switch( nTargetId ) { case NatOriginLong: case ProjCenterLong: case FalseOriginLong: case SphericalOriginLong: case InitialLongitude: // Note that the EPSG values are already relative to greenwich. // This shift is really making it relative to the provided prime // meridian, so that when SetTM() and company the correction back // ends up back relative to greenwich. dfResult = dfResult + dfFromGreenwich; break; default: ; } #endif return dfResult; } #define OGR_FP(x) OGR_FetchParm(adfProjParms, anParmIds, (x), dfFromGreenwich) /************************************************************************/ /* SetEPSGProjCS() */ /************************************************************************/ static OGRErr SetEPSGProjCS( OGRSpatialReference * poSRS, int nPCSCode ) { int nGCSCode = 0; int nUOMAngleCode = 0; int nUOMLength = 0; int nTRFCode = 0; int nCSC = 0; char *pszPCSName = nullptr; double adfTOWGS84[7] = { std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity() }; if( !EPSGGetPCSInfo( nPCSCode, &pszPCSName, &nUOMLength, &nUOMAngleCode, &nGCSCode, &nTRFCode, &nCSC, adfTOWGS84 ) ) { CPLFree(pszPCSName); return OGRERR_UNSUPPORTED_SRS; } poSRS->SetNode( "PROJCS", pszPCSName ); /* -------------------------------------------------------------------- */ /* Set GEOGCS. */ /* -------------------------------------------------------------------- */ const OGRErr nErr = SetEPSGGeogCS( poSRS, nGCSCode ); if( nErr != OGRERR_NONE ) { CPLFree(pszPCSName); return nErr; } /* -------------------------------------------------------------------- */ /* Set overridden TOWGS84 parameters */ /* -------------------------------------------------------------------- */ if( adfTOWGS84[0] != std::numeric_limits<double>::infinity() ) { poSRS->SetTOWGS84( adfTOWGS84[0], adfTOWGS84[1], adfTOWGS84[2], adfTOWGS84[3], adfTOWGS84[4], adfTOWGS84[5], adfTOWGS84[6] ); } // Used by OGR_FP macro const double dfFromGreenwich = poSRS->GetPrimeMeridian(); /* -------------------------------------------------------------------- */ /* Set linear units. */ /* -------------------------------------------------------------------- */ char *pszUOMLengthName = nullptr; double dfInMeters = 0.0; if( !EPSGGetUOMLengthInfo( nUOMLength, &pszUOMLengthName, &dfInMeters ) ) { CPLFree(pszPCSName); return OGRERR_UNSUPPORTED_SRS; } poSRS->SetLinearUnits( pszUOMLengthName, dfInMeters ); poSRS->SetAuthority( "PROJCS|UNIT", "EPSG", nUOMLength ); CPLFree( pszUOMLengthName ); CPLFree( pszPCSName ); /* -------------------------------------------------------------------- */ /* Set projection and parameters. */ /* -------------------------------------------------------------------- */ int nProjMethod = 0; int anParmIds[7] = {}; double adfProjParms[7] = {}; if( !EPSGGetProjTRFInfo( nPCSCode, &nProjMethod, anParmIds, adfProjParms )) return OGRERR_UNSUPPORTED_SRS; switch( nProjMethod ) { case 9801: case 9817: // Really LCC near conformal. poSRS->SetLCC1SP( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ), OGR_FP( NatOriginScaleFactor ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); break; case 9802: poSRS->SetLCC( OGR_FP( StdParallel1Lat ), OGR_FP( StdParallel2Lat ), OGR_FP( FalseOriginLat ), OGR_FP( FalseOriginLong ), OGR_FP( FalseOriginEasting ), OGR_FP( FalseOriginNorthing )); break; case 9803: poSRS->SetLCCB( OGR_FP( StdParallel1Lat ), OGR_FP( StdParallel2Lat ), OGR_FP( FalseOriginLat ), OGR_FP( FalseOriginLong ), OGR_FP( FalseOriginEasting ), OGR_FP( FalseOriginNorthing )); break; case 9805: poSRS->SetMercator2SP( OGR_FP( StdParallel1Lat ), OGR_FP( NatOriginLat ), OGR_FP(NatOriginLong), OGR_FP( FalseEasting ), OGR_FP(FalseNorthing) ); break; case 9804: case 9841: // Mercator 1SP (Spherical). case 1024: // Google Mercator. poSRS->SetMercator( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ), OGR_FP( NatOriginScaleFactor ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); // override hack for google mercator. if( nProjMethod == 1024 || nProjMethod == 9841 ) { poSRS->SetExtension( "PROJCS", "PROJ4", "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 " "+x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null " "+wktext +no_defs" ); } break; case 9806: poSRS->SetCS( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); break; case 9807: poSRS->SetTM( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ), OGR_FP( NatOriginScaleFactor ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); break; case 9808: poSRS->SetTMSO( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ), OGR_FP( NatOriginScaleFactor ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); break; case 9809: poSRS->SetOS( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ), OGR_FP( NatOriginScaleFactor ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); break; case 9810: poSRS->SetPS( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ), OGR_FP( NatOriginScaleFactor ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); break; case 9811: poSRS->SetNZMG( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); break; case 9812: case 9813: { poSRS->SetHOM( OGR_FP( ProjCenterLat ), OGR_FP( ProjCenterLong ), OGR_FP( Azimuth ), OGR_FP( AngleRectifiedToSkewedGrid ), OGR_FP( InitialLineScaleFactor ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); OGR_SRSNode *poNode = poSRS->GetAttrNode( "PROJECTION" )->GetChild( 0 ); if( nProjMethod == 9813 ) poNode->SetValue( SRS_PT_LABORDE_OBLIQUE_MERCATOR ); break; } case 9814: // NOTE: This is no longer used. Swiss Oblique Mercator gets // implemented using 9815 instead. poSRS->SetSOC( OGR_FP( ProjCenterLat ), OGR_FP( ProjCenterLong ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); break; case 9815: poSRS->SetHOMAC( OGR_FP( ProjCenterLat ), OGR_FP( ProjCenterLong ), OGR_FP( Azimuth ), OGR_FP( AngleRectifiedToSkewedGrid ), OGR_FP( InitialLineScaleFactor ), OGR_FP( ProjCenterEasting ), OGR_FP( ProjCenterNorthing ) ); break; case 9816: poSRS->SetTMG( OGR_FP( FalseOriginLat ), OGR_FP( FalseOriginLong ), OGR_FP( FalseOriginEasting ), OGR_FP( FalseOriginNorthing ) ); break; case 9818: poSRS->SetPolyconic( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); break; case 1041: // Used by EPSG:5514. case 9819: { double dfCenterLong = OGR_FP( ProjCenterLong ); if( dfCenterLong == 0.0 ) // See ticket #2559. dfCenterLong = OGR_FP( PolarLongOrigin ); double dfAzimuth = OGR_FP( CoLatConeAxis ); // See ticket #4223. if( dfAzimuth == 0.0 ) dfAzimuth = OGR_FP( Azimuth ); poSRS->SetKrovak( OGR_FP( ProjCenterLat ), dfCenterLong, dfAzimuth, OGR_FP( PseudoStdParallelLat ), OGR_FP( PseudoStdParallelScaleFactor ), OGR_FP( ProjCenterEasting ), OGR_FP( ProjCenterNorthing ) ); } break; case 9820: case 1027: // Used by EPSG:2163, 3408, 3409, 3973 and 3974. poSRS->SetLAEA( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); break; case 9821: // DEPRECATED : this is the spherical form, and really needs // different equations which give different results but PROJ.4 // doesn't seem to support the spherical form. poSRS->SetLAEA( OGR_FP( SphericalOriginLat ), OGR_FP( SphericalOriginLong ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); break; case 9822: // Albers (Conic) Equal Area. poSRS->SetACEA( OGR_FP( StdParallel1Lat ), OGR_FP( StdParallel2Lat ), OGR_FP( FalseOriginLat ), OGR_FP( FalseOriginLong ), OGR_FP( FalseOriginEasting ), OGR_FP( FalseOriginNorthing ) ); break; case 9823: // Equidistant Cylindrical / Plate Carre / Equirectangular. case 9842: case 1028: case 1029: poSRS->SetEquirectangular( OGR_FP( NatOriginLat ), OGR_FP( NatOriginLong ), 0.0, 0.0 ); break; case 9829: // Polar Stereographic (Variant B). poSRS->SetPS( OGR_FP( PolarLatStdParallel ), OGR_FP(PolarLongOrigin), 1.0, OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); break; case 9834: // Lambert Cylindrical Equal Area (Spherical) bug #2659. case 9835: // Lambert Cylindrical Equal Area (Ellipsoidal). poSRS->SetCEA( OGR_FP( StdParallel1Lat ), OGR_FP( NatOriginLong ), OGR_FP( FalseEasting ), OGR_FP( FalseNorthing ) ); break; default: CPLDebug( "EPSG", "No WKT support for projection method %d.", nProjMethod ); return OGRERR_UNSUPPORTED_SRS; } /* -------------------------------------------------------------------- */ /* Set overall PCS authority code. */ /* -------------------------------------------------------------------- */ poSRS->SetAuthority( "PROJCS", "EPSG", nPCSCode ); /* -------------------------------------------------------------------- */ /* Set axes */ /* -------------------------------------------------------------------- */ if( nCSC > 0 ) { SetEPSGAxisInfo( poSRS, "PROJCS", nCSC ); CPLErrorReset(); } return OGRERR_NONE; } /************************************************************************/ /* SetEPSGVertCS() */ /************************************************************************/ static OGRErr SetEPSGVertCS( OGRSpatialReference * poSRS, int nVertCSCode ) { /* -------------------------------------------------------------------- */ /* Fetch record from the vertcs.csv or override file. */ /* -------------------------------------------------------------------- */ const char *pszFilename = CSVFilename( "vertcs.override.csv" ); char szSearchKey[24] = {}; snprintf( szSearchKey, sizeof(szSearchKey), "%d", nVertCSCode ); char **papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer ); if( papszRecord == nullptr ) { pszFilename = CSVFilename( "vertcs.csv" ); papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer ); } if( papszRecord == nullptr ) return OGRERR_UNSUPPORTED_SRS; /* -------------------------------------------------------------------- */ /* Setup the basic VERT_CS. */ /* -------------------------------------------------------------------- */ poSRS->SetVertCS( CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "COORD_REF_SYS_NAME")), CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "DATUM_NAME")) ); /* -------------------------------------------------------------------- */ /* Should we add a geoidgrids extension node? */ /* -------------------------------------------------------------------- */ const char *pszMethod = CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "COORD_OP_METHOD_CODE_1")); if( pszMethod && EQUAL(pszMethod, "9665") ) { const char *pszParm11 = CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "PARM_1_1")); poSRS->SetExtension( "VERT_CS|VERT_DATUM", "PROJ4_GRIDS", pszParm11 ); } /* -------------------------------------------------------------------- */ /* Setup the VERT_DATUM node. */ /* -------------------------------------------------------------------- */ poSRS->SetAuthority( "VERT_CS|VERT_DATUM", "EPSG", atoi(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "DATUM_CODE"))) ); /* -------------------------------------------------------------------- */ /* Set linear units. */ /* -------------------------------------------------------------------- */ const int nUOM_CODE = atoi(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "UOM_CODE"))); char *pszUOMLengthName = nullptr; double dfInMeters = 0.0; if( !EPSGGetUOMLengthInfo( nUOM_CODE, &pszUOMLengthName, &dfInMeters ) ) { CPLError( CE_Failure, CPLE_AppDefined, "Failed to lookup UOM CODE %d", nUOM_CODE ); } else { poSRS->SetTargetLinearUnits( "VERT_CS", pszUOMLengthName, dfInMeters ); poSRS->SetAuthority( "VERT_CS|UNIT", "EPSG", nUOM_CODE ); CPLFree( pszUOMLengthName ); } /* -------------------------------------------------------------------- */ /* Set overall authority code. */ /* -------------------------------------------------------------------- */ poSRS->SetAuthority( "VERT_CS", "EPSG", nVertCSCode ); return OGRERR_NONE; } /************************************************************************/ /* SetEPSGCompdCS() */ /************************************************************************/ static OGRErr SetEPSGCompdCS( OGRSpatialReference * poSRS, int nCCSCode ) { /* -------------------------------------------------------------------- */ /* Fetch record from the compdcs.csv or override file. */ /* -------------------------------------------------------------------- */ char szSearchKey[24] = {}; snprintf( szSearchKey, sizeof(szSearchKey), "%d", nCCSCode ); // So far no override file needed. // pszFilename = CSVFilename( "compdcs.override.csv" ); // papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE", // szSearchKey, CC_Integer ); //if( papszRecord == NULL ) const char *pszFilename = CSVFilename( "compdcs.csv" ); char **papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer ); if( papszRecord == nullptr ) return OGRERR_UNSUPPORTED_SRS; /* -------------------------------------------------------------------- */ /* Fetch subinformation now before anything messes with the */ /* last loaded record. */ /* -------------------------------------------------------------------- */ const int nPCSCode = atoi(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "CMPD_HORIZCRS_CODE"))); const int nVertCSCode = atoi(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "CMPD_VERTCRS_CODE"))); /* -------------------------------------------------------------------- */ /* Set the COMPD_CS node with a name. */ /* -------------------------------------------------------------------- */ poSRS->SetNode( "COMPD_CS", CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "COORD_REF_SYS_NAME")) ); /* -------------------------------------------------------------------- */ /* Lookup the projected coordinate system. Can the */ /* horizontal CRS be a GCS? */ /* -------------------------------------------------------------------- */ OGRSpatialReference oPCS; OGRErr eErr = SetEPSGProjCS( &oPCS, nPCSCode ); if( eErr != OGRERR_NONE ) { // perhaps it is a GCS? eErr = SetEPSGGeogCS( &oPCS, nPCSCode ); } if( eErr != OGRERR_NONE ) { return eErr; } poSRS->GetRoot()->AddChild( oPCS.GetRoot()->Clone() ); /* -------------------------------------------------------------------- */ /* Lookup the VertCS. */ /* -------------------------------------------------------------------- */ OGRSpatialReference oVertCS; eErr = SetEPSGVertCS( &oVertCS, nVertCSCode ); if( eErr != OGRERR_NONE ) return eErr; poSRS->GetRoot()->AddChild( oVertCS.GetRoot()->Clone() ); /* -------------------------------------------------------------------- */ /* Set overall authority code. */ /* -------------------------------------------------------------------- */ poSRS->SetAuthority( "COMPD_CS", "EPSG", nCCSCode ); return OGRERR_NONE; } /************************************************************************/ /* SetEPSGGeocCS() */ /************************************************************************/ static OGRErr SetEPSGGeocCS( OGRSpatialReference * poSRS, int nGCSCode ) { /* -------------------------------------------------------------------- */ /* Fetch record from the geoccs.csv or override file. */ /* -------------------------------------------------------------------- */ char szSearchKey[24] = {}; snprintf( szSearchKey, sizeof(szSearchKey), "%d", nGCSCode ); // So far no override file needed. // pszFilename = CSVFilename( "compdcs.override.csv" ); // papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE", // szSearchKey, CC_Integer ); // if( papszRecord == NULL ) const char *pszFilename = CSVFilename( "geoccs.csv" ); char **papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE", szSearchKey, CC_Integer ); if( papszRecord == nullptr ) return OGRERR_UNSUPPORTED_SRS; /* -------------------------------------------------------------------- */ /* Set the GEOCCS node with a name. */ /* -------------------------------------------------------------------- */ poSRS->Clear(); poSRS->SetGeocCS( CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "COORD_REF_SYS_NAME")) ); /* -------------------------------------------------------------------- */ /* Get datum related information. */ /* -------------------------------------------------------------------- */ const int nDatumCode = atoi(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "DATUM_CODE"))); char *pszDatumName = CPLStrdup( CSLGetField(papszRecord, CSVGetFileFieldId(pszFilename, "DATUM_NAME")) ); OGREPSGDatumNameMassage( &pszDatumName ); const int nEllipsoidCode = atoi(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "ELLIPSOID_CODE"))); const int nPMCode = atoi(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "PRIME_MERIDIAN_CODE"))); /* -------------------------------------------------------------------- */ /* Get prime meridian information. */ /* -------------------------------------------------------------------- */ char *pszPMName = nullptr; double dfPMOffset = 0.0; if( !EPSGGetPMInfo( nPMCode, &pszPMName, &dfPMOffset ) ) { CPLFree( pszDatumName ); return OGRERR_UNSUPPORTED_SRS; } /* -------------------------------------------------------------------- */ /* Get the ellipsoid information. */ /* -------------------------------------------------------------------- */ char *pszEllipsoidName = nullptr; double dfSemiMajor = 0.0; double dfInvFlattening = 0.0; if( OSRGetEllipsoidInfo( nEllipsoidCode, &pszEllipsoidName, &dfSemiMajor, &dfInvFlattening ) != OGRERR_NONE ) { CPLFree( pszDatumName ); CPLFree( pszPMName ); return OGRERR_UNSUPPORTED_SRS; } /* -------------------------------------------------------------------- */ /* Setup the spheroid. */ /* -------------------------------------------------------------------- */ OGR_SRSNode *poSpheroid = new OGR_SRSNode( "SPHEROID" ); poSpheroid->AddChild( new OGR_SRSNode( pszEllipsoidName ) ); char szValue[128] = {}; OGRsnPrintDouble( szValue, sizeof(szValue), dfSemiMajor ); poSpheroid->AddChild( new OGR_SRSNode(szValue) ); OGRsnPrintDouble( szValue, sizeof(szValue), dfInvFlattening ); poSpheroid->AddChild( new OGR_SRSNode(szValue) ); CPLFree( pszEllipsoidName ); /* -------------------------------------------------------------------- */ /* Setup the Datum. */ /* -------------------------------------------------------------------- */ OGR_SRSNode *poDatum = new OGR_SRSNode( "DATUM" ); poDatum->AddChild( new OGR_SRSNode(pszDatumName) ); poDatum->AddChild( poSpheroid ); poSRS->GetRoot()->AddChild( poDatum ); CPLFree( pszDatumName ); /* -------------------------------------------------------------------- */ /* Setup the prime meridian. */ /* -------------------------------------------------------------------- */ if( dfPMOffset == 0.0 ) strcpy( szValue, "0" ); else OGRsnPrintDouble( szValue, sizeof(szValue), dfPMOffset ); OGR_SRSNode *poPM = new OGR_SRSNode( "PRIMEM" ); poPM->AddChild( new OGR_SRSNode( pszPMName ) ); poPM->AddChild( new OGR_SRSNode( szValue ) ); poSRS->GetRoot()->AddChild( poPM ); CPLFree( pszPMName ); /* -------------------------------------------------------------------- */ /* Should we try to lookup a datum transform? */ /* -------------------------------------------------------------------- */ #if 0 if( EPSGGetWGS84Transform( nGeogCS, adfBursaTransform ) ) { char szValue[100] = {}; OGR_SRSNode *poWGS84 = new OGR_SRSNode( "TOWGS84" ); for( int iCoeff = 0; iCoeff < 7; iCoeff++ ) { CPLsnprintf( szValue, sizeof(szValue), "%g", adfBursaTransform[iCoeff] ); poWGS84->AddChild( new OGR_SRSNode( szValue ) ); } poSRS->GetAttrNode( "DATUM" )->AddChild( poWGS84 ); } #endif /* -------------------------------------------------------------------- */ /* Set linear units. */ /* -------------------------------------------------------------------- */ int nUOMLength = atoi(CSLGetField( papszRecord, CSVGetFileFieldId(pszFilename, "UOM_CODE"))); double dfInMeters = 1.0; char *pszUOMLengthName = nullptr; if( !EPSGGetUOMLengthInfo( nUOMLength, &pszUOMLengthName, &dfInMeters ) ) { return OGRERR_UNSUPPORTED_SRS; } poSRS->SetLinearUnits( pszUOMLengthName, dfInMeters ); poSRS->SetAuthority( "GEOCCS|UNIT", "EPSG", nUOMLength ); CPLFree( pszUOMLengthName ); /* -------------------------------------------------------------------- */ /* Set axes */ /* -------------------------------------------------------------------- */ OGR_SRSNode *poAxis = new OGR_SRSNode( "AXIS" ); poAxis->AddChild( new OGR_SRSNode( "Geocentric X" ) ); poAxis->AddChild( new OGR_SRSNode( OSRAxisEnumToName(OAO_Other) ) ); poSRS->GetRoot()->AddChild( poAxis ); poAxis = new OGR_SRSNode( "AXIS" ); poAxis->AddChild( new OGR_SRSNode( "Geocentric Y" ) ); poAxis->AddChild( new OGR_SRSNode( OSRAxisEnumToName(OAO_Other) ) ); poSRS->GetRoot()->AddChild( poAxis ); poAxis = new OGR_SRSNode( "AXIS" ); poAxis->AddChild( new OGR_SRSNode( "Geocentric Z" ) ); poAxis->AddChild( new OGR_SRSNode( OSRAxisEnumToName(OAO_North) ) ); poSRS->GetRoot()->AddChild( poAxis ); /* -------------------------------------------------------------------- */ /* Set the authority codes. */ /* -------------------------------------------------------------------- */ poSRS->SetAuthority( "DATUM", "EPSG", nDatumCode ); poSRS->SetAuthority( "SPHEROID", "EPSG", nEllipsoidCode ); poSRS->SetAuthority( "PRIMEM", "EPSG", nPMCode ); poSRS->SetAuthority( "GEOCCS", "EPSG", nGCSCode ); return OGRERR_NONE; } /************************************************************************/ /* importFromEPSG() */ /************************************************************************/ /** * \brief Initialize SRS based on EPSG GCS or PCS code. * * This method will initialize the spatial reference based on the * passed in EPSG GCS or PCS code. The coordinate system definitions * are normally read from the EPSG derived support files such as * pcs.csv, gcs.csv, pcs.override.csv, gcs.override.csv and falling * back to search for a PROJ.4 epsg init file or a definition in epsg.wkt. * * These support files are normally searched for in /usr/local/share/gdal * or in the directory identified by the GDAL_DATA configuration option. * See CPLFindFile() for details. * * This method is relatively expensive, and generally involves quite a bit * of text file scanning. Reasonable efforts should be made to avoid calling * it many times for the same coordinate system. * * This method is similar to importFromEPSGA() except that EPSG preferred * axis ordering will *not* be applied for geographic coordinate systems. * EPSG normally defines geographic coordinate systems to use lat/long * contrary to typical GIS use). Since OGR 1.10.0, EPSG preferred * axis ordering will also *not* be applied for projected coordinate systems * that use northing/easting order. * * This method is the same as the C function OSRImportFromEPSG(). * * @param nCode a GCS or PCS code from the horizontal coordinate system table. * * @return OGRERR_NONE on success, or an error code on failure. */ OGRErr OGRSpatialReference::importFromEPSG( int nCode ) { OGRErr eErr = importFromEPSGA( nCode ); // Strip any GCS axis settings found. if( eErr == OGRERR_NONE ) { OGR_SRSNode *poGEOGCS = GetAttrNode( "GEOGCS" ); if( poGEOGCS != nullptr ) poGEOGCS->StripNodes( "AXIS" ); OGR_SRSNode *poPROJCS = GetAttrNode( "PROJCS" ); if( poPROJCS != nullptr && EPSGTreatsAsNorthingEasting() ) poPROJCS->StripNodes( "AXIS" ); } return eErr; } /************************************************************************/ /* OSRImportFromEPSG() */ /************************************************************************/ /** * \brief Initialize SRS based on EPSG GCS or PCS code. * * This function is the same as OGRSpatialReference::importFromEPSG(). */ OGRErr CPL_STDCALL OSRImportFromEPSG( OGRSpatialReferenceH hSRS, int nCode ) { VALIDATE_POINTER1( hSRS, "OSRImportFromEPSG", OGRERR_FAILURE ); return reinterpret_cast<OGRSpatialReference *>(hSRS)-> importFromEPSG( nCode ); } /************************************************************************/ /* importFromEPSGA() */ /************************************************************************/ /** * \brief Initialize SRS based on EPSG GCS or PCS code. * * This method will initialize the spatial reference based on the * passed in EPSG GCS or PCS code. * * This method is similar to importFromEPSG() except that EPSG preferred axis * ordering *will* be applied for geographic and projected coordinate systems. * EPSG normally defines geographic coordinate systems to use lat/long, and also * there are also a few projected coordinate systems that use northing/easting * order contrary to typical GIS use). See * OGRSpatialReference::importFromEPSG() for more details on operation of this * method. * * This method is the same as the C function OSRImportFromEPSGA(). * * @param nCode a GCS or PCS code from the horizontal coordinate system table. * * @return OGRERR_NONE on success, or an error code on failure. */ OGRErr OGRSpatialReference::importFromEPSGA( int nCode ) { return importFromEPSGAInternal(nCode, nullptr); } /************************************************************************/ /* importFromEPSGAInternal() */ /************************************************************************/ OGRErr OGRSpatialReference::importFromEPSGAInternal( int nCode, const char* pszSRSType ) { const int nCodeIn = nCode; // HACK to support 3D WGS84 if( nCode == 4979 ) nCode = 4326; bNormInfoSet = FALSE; /* -------------------------------------------------------------------- */ /* Clear any existing definition. */ /* -------------------------------------------------------------------- */ if( GetRoot() != nullptr ) { delete poRoot; poRoot = nullptr; } /* -------------------------------------------------------------------- */ /* Verify that we can find the required filename(s). */ /* -------------------------------------------------------------------- */ if( CSVScanFileByName( CSVFilename( "gcs.csv" ), "COORD_REF_SYS_CODE", "4269", CC_Integer ) == nullptr ) { CPLError( CE_Failure, CPLE_OpenFailed, "Unable to open EPSG support file %s. " "Try setting the GDAL_DATA environment variable to point to " "the directory containing EPSG csv files.", CSVFilename( "gcs.csv" ) ); return OGRERR_FAILURE; } /* -------------------------------------------------------------------- */ /* Try this as various sorts of objects till one works. */ /* -------------------------------------------------------------------- */ OGRErr eErr; if( pszSRSType && EQUAL(pszSRSType, "GEOGCS") ) { eErr = SetEPSGGeogCS( this, nCode ); if( eErr != OGRERR_NONE ) return eErr; } else if( pszSRSType && EQUAL(pszSRSType, "PROJCS") ) { eErr = SetEPSGProjCS( this, nCode ); if( eErr != OGRERR_NONE ) return eErr; } else { eErr = SetEPSGGeogCS( this, nCode ); } if( eErr == OGRERR_UNSUPPORTED_SRS ) eErr = SetEPSGProjCS( this, nCode ); if( eErr == OGRERR_UNSUPPORTED_SRS ) eErr = SetEPSGVertCS( this, nCode ); if( eErr == OGRERR_UNSUPPORTED_SRS ) eErr = SetEPSGCompdCS( this, nCode ); if( eErr == OGRERR_UNSUPPORTED_SRS ) eErr = SetEPSGGeocCS( this, nCode ); /* -------------------------------------------------------------------- */ /* If we get it as an unsupported code, try looking it up in */ /* the epsg.wkt coordinate system dictionary. */ /* -------------------------------------------------------------------- */ if( eErr == OGRERR_UNSUPPORTED_SRS ) { char szCode[32] = { 0 }; snprintf( szCode, sizeof(szCode), "%d", nCode ); eErr = importFromDict( "epsg.wkt", szCode ); } /* -------------------------------------------------------------------- */ /* If we get it as an unsupported code, try looking it up in */ /* the PROJ.4 support file(s). */ /* -------------------------------------------------------------------- */ if( eErr == OGRERR_UNSUPPORTED_SRS ) { char szWrkDefn[100] = {}; snprintf( szWrkDefn, sizeof(szWrkDefn), "+init=epsg:%d", nCode ); char *pszNormalized = OCTProj4Normalize( szWrkDefn ); if( strstr(pszNormalized, "proj=") != nullptr ) eErr = importFromProj4( pszNormalized ); CPLFree( pszNormalized ); } /* -------------------------------------------------------------------- */ /* Push in authority information if we were successful, and it */ /* is not already present. */ /* -------------------------------------------------------------------- */ const char *pszAuthName = nullptr; if( IsProjected() ) pszAuthName = GetAuthorityName( "PROJCS" ); else pszAuthName = GetAuthorityName( "GEOGCS" ); if( eErr == OGRERR_NONE && (pszAuthName == nullptr || nCode != nCodeIn) ) { if( IsProjected() ) SetAuthority( "PROJCS", "EPSG", nCodeIn ); else if( IsGeographic() ) SetAuthority( "GEOGCS", "EPSG", nCodeIn ); } /* -------------------------------------------------------------------- */ /* Otherwise officially issue an error message. */ /* -------------------------------------------------------------------- */ if( eErr == OGRERR_UNSUPPORTED_SRS ) { CPLError( CE_Failure, CPLE_NotSupported, "EPSG PCS/GCS code %d not found in EPSG support files. " "Is this a valid EPSG coordinate system?", nCode ); } /* -------------------------------------------------------------------- */ /* To the extent possible, we want to return the results in as */ /* close to standard OGC format as possible, so we fixup the */ /* ordering. */ /* -------------------------------------------------------------------- */ if( eErr == OGRERR_NONE ) { eErr = FixupOrdering(); } return eErr; } /************************************************************************/ /* OSRImportFromEPSGA() */ /************************************************************************/ /** * \brief Initialize SRS based on EPSG GCS or PCS code. * * This function is the same as OGRSpatialReference::importFromEPSGA(). */ OGRErr CPL_STDCALL OSRImportFromEPSGA( OGRSpatialReferenceH hSRS, int nCode ) { VALIDATE_POINTER1( hSRS, "OSRImportFromEPSGA", OGRERR_FAILURE ); return reinterpret_cast<OGRSpatialReference *>(hSRS)-> importFromEPSGA( nCode ); } /************************************************************************/ /* SetStatePlane() */ /************************************************************************/ /** * \brief Set State Plane projection definition. * * This will attempt to generate a complete definition of a state plane * zone based on generating the entire SRS from the EPSG tables. If the * EPSG tables are unavailable, it will produce a stubbed LOCAL_CS definition * and return OGRERR_FAILURE. * * This method is the same as the C function OSRSetStatePlaneWithUnits(). * * @param nZone State plane zone number, in the USGS numbering scheme (as * distinct from the Arc/Info and Erdas numbering scheme. * * @param bNAD83 TRUE if the NAD83 zone definition should be used or FALSE * if the NAD27 zone definition should be used. * * @param pszOverrideUnitName Linear unit name to apply overriding the * legal definition for this zone. * * @param dfOverrideUnit Linear unit conversion factor to apply overriding * the legal definition for this zone. * * @return OGRERR_NONE on success, or OGRERR_FAILURE on failure, mostly likely * due to the EPSG tables not being accessible. */ OGRErr OGRSpatialReference::SetStatePlane( int nZone, int bNAD83, const char *pszOverrideUnitName, double dfOverrideUnit ) { /* -------------------------------------------------------------------- */ /* Get the index id from stateplane.csv. */ /* -------------------------------------------------------------------- */ if( !bNAD83 && nZone > INT_MAX - 10000 ) return OGRERR_FAILURE; const int nAdjustedId = bNAD83 ? nZone : nZone + 10000; /* -------------------------------------------------------------------- */ /* Turn this into a PCS code. We assume there will only be one */ /* PCS corresponding to each Proj_ code since the proj code */ /* already effectively indicates NAD27 or NAD83. */ /* -------------------------------------------------------------------- */ char szID[32] = {}; snprintf( szID, sizeof(szID), "%d", nAdjustedId ); const int nPCSCode = atoi( CSVGetField( CSVFilename( "stateplane.csv" ), "ID", szID, CC_Integer, "EPSG_PCS_CODE" ) ); if( nPCSCode < 1 ) { static bool bFailureReported = false; if( !bFailureReported ) { bFailureReported = true; CPLError( CE_Warning, CPLE_OpenFailed, "Unable to find state plane zone in stateplane.csv, " "likely because the GDAL data files cannot be found. " "Using incomplete definition of state plane zone." ); } Clear(); if( bNAD83 ) { char szName[128] = {}; snprintf( szName, sizeof(szName), "State Plane Zone %d / NAD83", nZone ); SetLocalCS( szName ); SetLinearUnits( SRS_UL_METER, 1.0 ); } else { char szName[128] = {}; snprintf( szName, sizeof(szName), "State Plane Zone %d / NAD27", nZone ); SetLocalCS( szName ); SetLinearUnits( SRS_UL_US_FOOT, CPLAtof(SRS_UL_US_FOOT_CONV) ); } return OGRERR_FAILURE; } /* -------------------------------------------------------------------- */ /* Define based on a full EPSG definition of the zone. */ /* -------------------------------------------------------------------- */ OGRErr eErr = importFromEPSG( nPCSCode ); if( eErr != OGRERR_NONE ) return eErr; /* -------------------------------------------------------------------- */ /* Apply units override if required. */ /* */ /* We will need to adjust the linear projection parameter to */ /* match the provided units, and clear the authority code. */ /* -------------------------------------------------------------------- */ if( dfOverrideUnit != 0.0 && fabs(dfOverrideUnit - GetLinearUnits()) > 0.0000000001 ) { const double dfFalseEasting = GetNormProjParm( SRS_PP_FALSE_EASTING ); const double dfFalseNorthing = GetNormProjParm( SRS_PP_FALSE_NORTHING); SetLinearUnits( pszOverrideUnitName, dfOverrideUnit ); SetNormProjParm( SRS_PP_FALSE_EASTING, dfFalseEasting ); SetNormProjParm( SRS_PP_FALSE_NORTHING, dfFalseNorthing ); OGR_SRSNode * const poPROJCS = GetAttrNode( "PROJCS" ); if( poPROJCS != nullptr && poPROJCS->FindChild( "AUTHORITY" ) != -1 ) { poPROJCS->DestroyChild( poPROJCS->FindChild( "AUTHORITY" ) ); } } return OGRERR_NONE; } /************************************************************************/ /* OSRSetStatePlane() */ /************************************************************************/ /** * \brief Set State Plane projection definition. * * This function is the same as OGRSpatialReference::SetStatePlane(). */ OGRErr OSRSetStatePlane( OGRSpatialReferenceH hSRS, int nZone, int bNAD83 ) { VALIDATE_POINTER1( hSRS, "OSRSetStatePlane", OGRERR_FAILURE ); return reinterpret_cast<OGRSpatialReference *>(hSRS)-> SetStatePlane( nZone, bNAD83 ); } /************************************************************************/ /* OSRSetStatePlaneWithUnits() */ /************************************************************************/ /** * \brief Set State Plane projection definition. * * This function is the same as OGRSpatialReference::SetStatePlane(). */ OGRErr OSRSetStatePlaneWithUnits( OGRSpatialReferenceH hSRS, int nZone, int bNAD83, const char *pszOverrideUnitName, double dfOverrideUnit ) { VALIDATE_POINTER1( hSRS, "OSRSetStatePlaneWithUnits", OGRERR_FAILURE ); return reinterpret_cast<OGRSpatialReference *>(hSRS)-> SetStatePlane( nZone, bNAD83, pszOverrideUnitName, dfOverrideUnit ); } /************************************************************************/ /* GetEPSGGeogCS() */ /************************************************************************/ /** Try to establish what the EPSG code for this coordinate systems * GEOGCS might be. Returns -1 if no reasonable guess can be made. * * @return EPSG code */ // TODO: We really need to do some name lookups. int OGRSpatialReference::GetEPSGGeogCS() const { const char *pszAuthName = GetAuthorityName( "GEOGCS" ); /* -------------------------------------------------------------------- */ /* Do we already have it? */ /* -------------------------------------------------------------------- */ if( pszAuthName != nullptr && EQUAL(pszAuthName, "epsg") ) return atoi(GetAuthorityCode( "GEOGCS" )); /* -------------------------------------------------------------------- */ /* Get the datum and geogcs names. */ /* -------------------------------------------------------------------- */ const char *pszGEOGCS = GetAttrValue( "GEOGCS" ); const char *pszDatum = GetAttrValue( "DATUM" ); // We can only operate on coordinate systems with a geogcs. if( pszGEOGCS == nullptr || pszDatum == nullptr ) return -1; /* -------------------------------------------------------------------- */ /* Is this a "well known" geographic coordinate system? */ /* -------------------------------------------------------------------- */ const bool bWGS = strstr(pszGEOGCS, "WGS") != nullptr || strstr(pszDatum, "WGS") || strstr(pszGEOGCS, "World Geodetic System") || strstr(pszGEOGCS, "World_Geodetic_System") || strstr(pszDatum, "World Geodetic System") || strstr(pszDatum, "World_Geodetic_System"); const bool bNAD = strstr(pszGEOGCS, "NAD") != nullptr || strstr(pszDatum, "NAD") || strstr(pszGEOGCS, "North American") || strstr(pszGEOGCS, "North_American") || strstr(pszDatum, "North American") || strstr(pszDatum, "North_American"); if( bWGS && (strstr(pszGEOGCS, "84") || strstr(pszDatum, "84")) ) return 4326; if( bWGS && (strstr(pszGEOGCS, "72") || strstr(pszDatum, "72")) ) return 4322; if( bNAD && (strstr(pszGEOGCS, "83") || strstr(pszDatum, "83")) ) return 4269; if( bNAD && (strstr(pszGEOGCS, "27") || strstr(pszDatum, "27")) ) return 4267; /* -------------------------------------------------------------------- */ /* If we know the datum, associate the most likely GCS with */ /* it. */ /* -------------------------------------------------------------------- */ pszAuthName = GetAuthorityName( "GEOGCS|DATUM" ); if( pszAuthName != nullptr && EQUAL(pszAuthName, "epsg") && GetPrimeMeridian() == 0.0 ) { const int nDatum = atoi(GetAuthorityCode("GEOGCS|DATUM")); if( nDatum >= 6000 && nDatum <= 6999 ) return nDatum - 2000; } return -1; } /************************************************************************/ /* AutoIdentifyEPSG() */ /************************************************************************/ /** * \brief Set EPSG authority info if possible. * * This method inspects a WKT definition, and adds EPSG authority nodes * where an aspect of the coordinate system can be easily and safely * corresponded with an EPSG identifier. In practice, this method will * evolve over time. In theory it can add authority nodes for any object * (i.e. spheroid, datum, GEOGCS, units, and PROJCS) that could have an * authority node. Mostly this is useful to inserting appropriate * PROJCS codes for common formulations (like UTM n WGS84). * * If it success the OGRSpatialReference is updated in place, and the * method return OGRERR_NONE. If the method fails to identify the * general coordinate system OGRERR_UNSUPPORTED_SRS is returned but no * error message is posted via CPLError(). * * This method is the same as the C function OSRAutoIdentifyEPSG(). * * Since GDAL 2.3, the FindMatches() method can also be used for improved * matching by researching the EPSG catalog. * * @return OGRERR_NONE or OGRERR_UNSUPPORTED_SRS. */ OGRErr OGRSpatialReference::AutoIdentifyEPSG() { /* -------------------------------------------------------------------- */ /* Do we have a GEOGCS node, but no authority? If so, try */ /* guessing it. */ /* -------------------------------------------------------------------- */ if( (IsProjected() || IsGeographic()) && GetAuthorityCode( "GEOGCS" ) == nullptr ) { const int nGCS = GetEPSGGeogCS(); if( nGCS != -1 ) SetAuthority( "GEOGCS", "EPSG", nGCS ); } if( IsProjected() && GetAuthorityCode( "PROJCS") == nullptr ) { const char *pszProjection = GetAttrValue( "PROJECTION" ); /* -------------------------------------------------------------------- */ /* Is this a UTM coordinate system with a common GEOGCS? */ /* -------------------------------------------------------------------- */ int nZone = 0; int bNorth = FALSE; if( (nZone = GetUTMZone( &bNorth )) != 0 ) { const char *pszAuthName = GetAuthorityName( "PROJCS|GEOGCS" ); const char *pszAuthCode = GetAuthorityCode( "PROJCS|GEOGCS" ); if( pszAuthName == nullptr || pszAuthCode == nullptr ) { // Don't exactly recognise datum. } else if( EQUAL(pszAuthName, "EPSG") && atoi(pszAuthCode) == 4326 ) { // WGS84 if( bNorth ) SetAuthority( "PROJCS", "EPSG", 32600 + nZone ); else SetAuthority( "PROJCS", "EPSG", 32700 + nZone ); } else if( EQUAL(pszAuthName, "EPSG") && atoi(pszAuthCode) == 4267 && nZone >= 3 && nZone <= 22 && bNorth ) { SetAuthority( "PROJCS", "EPSG", 26700 + nZone ); // NAD27 } else if( EQUAL(pszAuthName, "EPSG") && atoi(pszAuthCode) == 4269 && nZone >= 3 && nZone <= 23 && bNorth ) { SetAuthority( "PROJCS", "EPSG", 26900 + nZone ); // NAD83 } else if( EQUAL(pszAuthName, "EPSG") && atoi(pszAuthCode) == 4322 ) { // WGS72 if( bNorth ) SetAuthority( "PROJCS", "EPSG", 32200 + nZone ); else SetAuthority( "PROJCS", "EPSG", 32300 + nZone ); } } /* -------------------------------------------------------------------- */ /* Is this a Polar Stereographic system on WGS 84 ? */ /* -------------------------------------------------------------------- */ else if ( pszProjection != nullptr && EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC) ) { const char *pszAuthName = GetAuthorityName( "PROJCS|GEOGCS" ); const char *pszAuthCode = GetAuthorityCode( "PROJCS|GEOGCS" ); const double dfLatOrigin = GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 ); if( pszAuthName != nullptr && EQUAL(pszAuthName, "EPSG") && pszAuthCode != nullptr && atoi(pszAuthCode) == 4326 && fabs( fabs(dfLatOrigin ) - 71.0 ) < 1e-15 && fabs(GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 )) < 1e-15 && fabs(GetProjParm( SRS_PP_SCALE_FACTOR, 1.0 ) - 1.0) < 1e-15 && fabs(GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 )) < 1e-15 && fabs(GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 )) < 1e-15 && fabs(GetLinearUnits() - 1.0) < 1e-15 ) { if( dfLatOrigin > 0 ) // Arctic Polar Stereographic SetAuthority( "PROJCS", "EPSG", 3995 ); else // Antarctic Polar Stereographic SetAuthority( "PROJCS", "EPSG", 3031 ); } } } /* -------------------------------------------------------------------- */ /* Return. */ /* -------------------------------------------------------------------- */ if( IsProjected() && GetAuthorityCode("PROJCS") != nullptr ) return OGRERR_NONE; if( IsGeographic() && GetAuthorityCode("GEOGCS") != nullptr ) return OGRERR_NONE; return OGRERR_UNSUPPORTED_SRS; } /************************************************************************/ /* OSRAutoIdentifyEPSG() */ /************************************************************************/ /** * \brief Set EPSG authority info if possible. * * This function is the same as OGRSpatialReference::AutoIdentifyEPSG(). * * Since GDAL 2.3, the OSRFindMatches() function can also be used for improved * matching by researching the EPSG catalog. * */ OGRErr OSRAutoIdentifyEPSG( OGRSpatialReferenceH hSRS ) { VALIDATE_POINTER1( hSRS, "OSRAutoIdentifyEPSG", OGRERR_FAILURE ); return reinterpret_cast<OGRSpatialReference *>(hSRS)->AutoIdentifyEPSG(); } /************************************************************************/ /* EPSGTreatsAsLatLong() */ /************************************************************************/ /** * \brief This method returns TRUE if EPSG feels this geographic coordinate * system should be treated as having lat/long coordinate ordering. * * Currently this returns TRUE for all geographic coordinate systems * with an EPSG code set, and AXIS values set defining it as lat, long. * Note that coordinate systems with an EPSG code and no axis settings * will be assumed to not be lat/long. * * FALSE will be returned for all coordinate systems that are not geographic, * or that do not have an EPSG code set. * * This method is the same as the C function OSREPSGTreatsAsLatLong(). * * @return TRUE or FALSE. */ int OGRSpatialReference::EPSGTreatsAsLatLong() const { if( !IsGeographic() ) return FALSE; const char *pszAuth = GetAuthorityName( "GEOGCS" ); if( pszAuth == nullptr || !EQUAL(pszAuth, "EPSG") ) return FALSE; const OGR_SRSNode * const poFirstAxis = GetAttrNode( "GEOGCS|AXIS" ); if( poFirstAxis == nullptr ) return FALSE; if( poFirstAxis->GetChildCount() >= 2 && EQUAL(poFirstAxis->GetChild(1)->GetValue(), "NORTH") ) return TRUE; return FALSE; } /************************************************************************/ /* OSREPSGTreatsAsLatLong() */ /************************************************************************/ /** * \brief This function returns TRUE if EPSG feels this geographic coordinate * system should be treated as having lat/long coordinate ordering. * * This function is the same as OGRSpatialReference::OSREPSGTreatsAsLatLong(). */ int OSREPSGTreatsAsLatLong( OGRSpatialReferenceH hSRS ) { VALIDATE_POINTER1( hSRS, "OSREPSGTreatsAsLatLong", OGRERR_FAILURE ); return reinterpret_cast<OGRSpatialReference *>(hSRS)->EPSGTreatsAsLatLong(); } /************************************************************************/ /* EPSGTreatsAsNorthingEasting() */ /************************************************************************/ /** * \brief This method returns TRUE if EPSG feels this projected coordinate * system should be treated as having northing/easting coordinate ordering. * * Currently this returns TRUE for all projected coordinate systems * with an EPSG code set, and AXIS values set defining it as northing, easting. * * FALSE will be returned for all coordinate systems that are not projected, * or that do not have an EPSG code set. * * This method is the same as the C function EPSGTreatsAsNorthingEasting(). * * @return TRUE or FALSE. * * @since OGR 1.10.0 */ int OGRSpatialReference::EPSGTreatsAsNorthingEasting() const { if( !IsProjected() ) return FALSE; const char *pszAuth = GetAuthorityName( "PROJCS" ); if( pszAuth == nullptr || !EQUAL(pszAuth, "EPSG") ) return FALSE; const OGR_SRSNode * const poFirstAxis = GetAttrNode( "PROJCS|AXIS" ); if( poFirstAxis == nullptr ) return FALSE; if( poFirstAxis->GetChildCount() >= 2 && EQUAL(poFirstAxis->GetChild(1)->GetValue(), "NORTH") ) return TRUE; return FALSE; } /************************************************************************/ /* OSREPSGTreatsAsNorthingEasting() */ /************************************************************************/ /** * \brief This function returns TRUE if EPSG feels this geographic coordinate * system should be treated as having northing/easting coordinate ordering. * * This function is the same as * OGRSpatialReference::EPSGTreatsAsNorthingEasting(). * * @since OGR 1.10.0 */ int OSREPSGTreatsAsNorthingEasting( OGRSpatialReferenceH hSRS ) { VALIDATE_POINTER1( hSRS, "OSREPSGTreatsAsNorthingEasting", OGRERR_FAILURE ); return reinterpret_cast<OGRSpatialReference *>(hSRS)-> EPSGTreatsAsNorthingEasting(); } /************************************************************************/ /* CleanupFindMatchesCacheAndMutex() */ /************************************************************************/ void CleanupFindMatchesCacheAndMutex() { if( hFindMatchesMutex != nullptr ) { CPLDestroyMutex(hFindMatchesMutex); hFindMatchesMutex = nullptr; } if( papoSRSCache_GEOGCS ) { for( size_t i = 0; i < papoSRSCache_GEOGCS->size(); i++ ) delete (*papoSRSCache_GEOGCS)[i]; delete papoSRSCache_GEOGCS; papoSRSCache_GEOGCS = nullptr; } if( papoSRSCache_PROJCS ) { for( size_t i = 0; i < papoSRSCache_PROJCS->size(); i++ ) delete (*papoSRSCache_PROJCS)[i]; delete papoSRSCache_PROJCS; papoSRSCache_PROJCS = nullptr; } delete poMapESRIPROJCSNameToEPSGCode; poMapESRIPROJCSNameToEPSGCode = nullptr; delete poMapESRIGEOGCSNameToEPSGCode; poMapESRIGEOGCSNameToEPSGCode = nullptr; } /************************************************************************/ /* MassageSRSName() */ /************************************************************************/ /* Transform a SRS name typically coming from EPSG or ESRI into a simplified * form that can be compared. */ static CPLString MassageSRSName(const char* pszInput, bool bExtraMassaging) { CPLString osRet; bool bLastWasSep = false; for( int i = 0; pszInput[i] != '\0'; ++ i) { if( isdigit( pszInput[i] ) ) { if( (i > 0 && isalpha( pszInput[i-1] )) || bLastWasSep ) osRet += "_"; bLastWasSep = false; /* Abbreviate 19xx as xx */ if( pszInput[i] == '1' && pszInput[i+1] == '9' && pszInput[i+2] != '\0' && isdigit(pszInput[i+2]) && (i == 0 || !isdigit(pszInput[i-1])) ) { i ++; continue; } osRet += pszInput[i]; } else if( isalpha( pszInput[i] ) ) { if( bLastWasSep ) osRet += "_"; osRet += pszInput[i]; bLastWasSep = false; } else { bLastWasSep = true; } } osRet.tolower(); osRet.replaceAll("gauss_kruger", "gk"); // EPSG -> ESRI osRet.replaceAll("rt_90_25", "rt_90_2_5"); // ESRI -> EPSG osRet.replaceAll("rt_38_25", "rt_38_2_5"); // ESRI -> EPSG osRet.replaceAll("_zone_", "_"); // EPSG -> ESRI osRet.replaceAll("_stateplane_", "_"); // ESRI -> EPSG osRet.replaceAll("_nsidc_", "_"); // EPSG -> ESRI osRet.replaceAll("_I_", "_1_"); // ESRI -> EPSG osRet.replaceAll("_II_", "_2_"); // ESRI -> EPSG osRet.replaceAll("_III_", "_3_"); // ESRI -> EPSG osRet.replaceAll("_IV_", "_4_"); // ESRI -> EPSG osRet.replaceAll("_V_", "_5_"); // ESRI -> EPSG osRet.replaceAll("pulkovo_42_adj_83_", "pulkovo_42_83_"); // ESRI -> EPSG osRet.replaceAll("_old_fips", "_deprecated_fips"); if( bExtraMassaging ) osRet.replaceAll("_deprecated", ""); // EPSG -> ESRI // _FIPS_XXXX_Feet --> _ftUS ESRI -> EPSG // _FIPS_XXXX_Ft_US --> _ftUS ESRI -> EPSG // _FIPS_XXXX --> "" ESRI -> EPSG size_t nPos = osRet.find("_fips_"); if( nPos != std::string::npos ) { size_t nPos2 = osRet.find("_feet", nPos + strlen("_fips_")); if( nPos2 != std::string::npos && nPos2 + strlen("_feet") == osRet.size() ) { osRet.resize(nPos); osRet += "_ftus"; } else { nPos2 = osRet.find("_ft_us", nPos + strlen("_fips_")); if( nPos2 != std::string::npos && nPos2 + strlen("_ft_us") == osRet.size() ) { osRet.resize(nPos); osRet += "_ftus"; } else if( osRet.find('_', nPos + strlen("_fips_")) == std::string::npos ) { osRet.resize(nPos); } } } return osRet; } /************************************************************************/ /* IngestDict() */ /************************************************************************/ static void IngestDict(const char* pszDictFile, const char* pszSRSType, std::vector<OGRSpatialReference*>* papoSRSCache, VSILFILE* fpOut) { /* -------------------------------------------------------------------- */ /* Find and open file. */ /* -------------------------------------------------------------------- */ const char *pszFilename = CPLFindFile( "gdal", pszDictFile ); if( pszFilename == nullptr ) return; VSILFILE *fp = VSIFOpenL( pszFilename, "rb" ); if( fp == nullptr ) return; if( fpOut ) { VSIFPrintfL(fpOut, "# From %s\n", pszDictFile); } /* -------------------------------------------------------------------- */ /* Process lines. */ /* -------------------------------------------------------------------- */ const char *pszLine = nullptr; while( (pszLine = CPLReadLineL(fp)) != nullptr ) { if( pszLine[0] == '#' ) continue; const char* pszComma = strchr(pszLine, ','); if( pszComma ) { const char* pszWKT = pszComma + 1; if( STARTS_WITH(pszWKT, pszSRSType) ) { OGRSpatialReference* poSRS = new OGRSpatialReference(); if( poSRS->SetFromUserInput(pszWKT) == OGRERR_NONE ) { const char *pszProjection = poSRS->GetAttrValue( "PROJECTION" ); if( pszProjection && EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP) ) { // Remove duplicate Standard_Parallel_1 double dfLatOrigin = poSRS->GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN); double dfStdParallel1 = poSRS->GetProjParm(SRS_PP_STANDARD_PARALLEL_1); if( dfLatOrigin == dfStdParallel1 ) { OGR_SRSNode *poPROJCS = poSRS->GetAttrNode( "PROJCS" ); if( poPROJCS ) { const int iChild = poSRS->FindProjParm( SRS_PP_STANDARD_PARALLEL_1, poPROJCS ); if( iChild != -1 ) poPROJCS->DestroyChild( iChild); } } } poSRS->morphFromESRI(); papoSRSCache->push_back(poSRS); if( fpOut ) { VSIFPrintfL(fpOut, "%s\n", pszWKT); } } else { delete poSRS; } } } } /* -------------------------------------------------------------------- */ /* Cleanup */ /* -------------------------------------------------------------------- */ VSIFCloseL( fp ); } /************************************************************************/ /* BuildESRICSNameCache() */ /************************************************************************/ static void BuildESRICSNameCache(const char* pszSRSType, std::map<CPLString, int>* poMapCSNameToCode, VSILFILE* fpOut) { const char *pszFilename = CPLFindFile( "gdal", "esri_epsg.wkt" ); if( pszFilename == nullptr ) return; VSILFILE *fp = VSIFOpenL( pszFilename, "rb" ); if( fp == nullptr ) return; /* -------------------------------------------------------------------- */ /* Process lines. */ /* -------------------------------------------------------------------- */ const char *pszLine = nullptr; while( (pszLine = CPLReadLineL(fp)) != nullptr ) { if( pszLine[0] == '#' ) continue; const char* pszComma = strchr(pszLine, ','); if( pszComma ) { const char* pszWKT = pszComma + 1; if( STARTS_WITH(pszWKT, pszSRSType) ) { OGRSpatialReference oSRS; if( oSRS.SetFromUserInput(pszWKT) == OGRERR_NONE ) { const char* pszSRSName = oSRS.GetAttrValue(pszSRSType); const char* pszAuthCode = oSRS.GetAuthorityCode(nullptr); if( pszSRSName && pszAuthCode ) { (*poMapCSNameToCode)[pszSRSName] = atoi(pszAuthCode); VSIFPrintfL(fpOut, "%s,%s\n", pszSRSName, pszAuthCode); } } } } } VSIFCloseL( fp ); } /************************************************************************/ /* GetSRSCache() */ /************************************************************************/ const std::vector<OGRSpatialReference*>* OGRSpatialReference::GetSRSCache( const char* pszSRSType, const std::map<CPLString, int>*& poMapCSNameToCodeOut) { if( pszSRSType == nullptr ) return nullptr; std::vector<OGRSpatialReference*>* papoSRSCache = nullptr; std::map<CPLString, int>* poMapCSNameToCode = nullptr; CPLMutexHolderD(&hFindMatchesMutex); CPLString osFilename; const char* pszFilename = ""; if( EQUAL(pszSRSType, "PROJCS") ) { pszFilename = "pcs.csv"; if( papoSRSCache_PROJCS != nullptr ) return papoSRSCache_PROJCS; papoSRSCache_PROJCS = new std::vector<OGRSpatialReference*>(); papoSRSCache = papoSRSCache_PROJCS; poMapESRIPROJCSNameToEPSGCode = new std::map<CPLString, int>(); poMapCSNameToCode = poMapESRIPROJCSNameToEPSGCode; } else if( EQUAL(pszSRSType, "GEOGCS") ) { pszFilename = "gcs.csv" ; if( papoSRSCache_GEOGCS != nullptr ) return papoSRSCache_GEOGCS; papoSRSCache_GEOGCS = new std::vector<OGRSpatialReference*>(); papoSRSCache = papoSRSCache_GEOGCS; poMapESRIGEOGCSNameToEPSGCode = new std::map<CPLString, int>(); poMapCSNameToCode = poMapESRIGEOGCSNameToEPSGCode; } else { return nullptr; } poMapCSNameToCodeOut = poMapCSNameToCode; // First try to look an already built SRS cache in ~/.gdal/X.Y/srs_cache const char* pszHome = CPLGetHomeDir(); const char* pszCSVFilename = CSVFilename(pszFilename); CPLString osCacheFilename; CPLString osCacheDirectory = CPLGetConfigOption("OSR_SRS_CACHE_DIRECTORY", ""); if( (pszHome != nullptr || !osCacheDirectory.empty()) && CPLTestBool(CPLGetConfigOption("OSR_SRS_CACHE", "YES")) ) { if( osCacheDirectory.empty() ) { osCacheDirectory = CPLFormFilename(pszHome, ".gdal", nullptr); // Version this, because might be sensitive to GDAL / // EPSG versions osCacheDirectory = CPLFormFilename(osCacheDirectory, CPLSPrintf("%d.%d", GDAL_VERSION_MAJOR, GDAL_VERSION_MINOR), nullptr); osCacheDirectory = CPLFormFilename(osCacheDirectory, "srs_cache", nullptr); } osCacheFilename = CPLFormFilename(osCacheDirectory, CPLResetExtension(pszFilename, "wkt"), nullptr); VSIStatBufL sStatCache; VSIStatBufL sStatCSV; VSILFILE* fp = nullptr; VSILFILE* fpESRINames = nullptr; if( VSIStatL((osCacheFilename + ".gz").c_str(), &sStatCache) == 0 && VSIStatL(pszCSVFilename, &sStatCSV) == 0 && sStatCache.st_mtime >= sStatCSV.st_mtime && (VSIStatL(CPLResetExtension( pszCSVFilename, "override.csv"), &sStatCSV) != 0 || sStatCache.st_mtime >= sStatCSV.st_mtime) ) { fp = VSIFOpenL(("/vsigzip/" + osCacheFilename + ".gz").c_str(), "rb"); if( fp ) { fpESRINames = VSIFOpenL( (CPLString("/vsigzip/") + CPLResetExtension(osCacheFilename, "esri.gz")).c_str(), "rb"); } } if( fp ) { CPLDebug("OSR", "Using %s cache", osCacheFilename.c_str()); const char* pszLine; while( (pszLine = CPLReadLineL(fp)) != nullptr ) { OGRSpatialReference* poSRS = new OGRSpatialReference(); poSRS->SetFromUserInput(pszLine); papoSRSCache->push_back(poSRS); const char* pszSRSName = poSRS->GetAttrValue(pszSRSType); const char* pszAuthCode = poSRS->GetAuthorityCode(nullptr); if( pszSRSName && pszAuthCode ) { (*poMapCSNameToCode)[pszSRSName] = atoi(pszAuthCode); } } VSIFCloseL(fp); if( fpESRINames ) { while( (pszLine = CPLReadLineL(fpESRINames)) != nullptr ) { const char* pszComma = strchr(pszLine, ','); if( pszComma ) { CPLString osName(pszLine); osName.resize(pszComma - pszLine); (*poMapCSNameToCode)[osName] = atoi(pszComma + 1); } } VSIFCloseL(fpESRINames); } return papoSRSCache; } } // If no already built cache, ingest the EPSG database and write the // cache CPLDebug("OSR", "Building %s cache", pszSRSType); VSILFILE* fp = VSIFOpenL(pszCSVFilename, "rb"); if( fp == nullptr ) { return nullptr; } VSILFILE* fpOut = nullptr; if( !osCacheFilename.empty() ) { CPLString osDirname(CPLGetDirname(osCacheFilename)); CPLString osDirnameParent(CPLGetDirname(osDirname)); CPLString osDirnameGrantParent( CPLGetDirname(osDirnameParent)); VSIMkdir( osDirnameGrantParent, 0755 ); VSIMkdir( osDirnameParent, 0755 ); VSIMkdir( osDirname, 0755 ); fpOut = VSIFOpenL( ("/vsigzip/" + osCacheFilename + ".gz").c_str(), "wb"); if( fpOut != nullptr ) { VSIFPrintfL(fpOut, "# From %s\n", pszFilename); } } const char* pszLine; OGRSpatialReference* poSRS = nullptr; CPLPushErrorHandler(CPLQuietErrorHandler); while( (pszLine = CPLReadLineL(fp)) != nullptr ) { if( poSRS == nullptr ) poSRS = new OGRSpatialReference(); int nCode = atoi(pszLine); if( nCode > 0 && poSRS->importFromEPSGAInternal(nCode, pszSRSType) == OGRERR_NONE ) { // Strip AXIS like in importFromEPSG() OGR_SRSNode *poGEOGCS = poSRS->GetAttrNode( "GEOGCS" ); if( poGEOGCS != nullptr ) poGEOGCS->StripNodes( "AXIS" ); OGR_SRSNode *poPROJCS = poSRS->GetAttrNode( "PROJCS" ); if( poPROJCS != nullptr && poSRS->EPSGTreatsAsNorthingEasting() ) poPROJCS->StripNodes( "AXIS" ); if( fpOut ) { char* pszWKT = nullptr; poSRS->exportToWkt(&pszWKT); if( pszWKT ) { VSIFPrintfL(fpOut, "%s\n", pszWKT); CPLFree(pszWKT); } } papoSRSCache->push_back(poSRS); const char* pszSRSName = poSRS->GetAttrValue(pszSRSType); if( pszSRSName ) { (*poMapCSNameToCode)[pszSRSName] = nCode; } poSRS = nullptr; } } CPLPopErrorHandler(); delete poSRS; VSIFCloseL(fp); IngestDict("esri_extra.wkt", pszSRSType, papoSRSCache, fpOut); if( fpOut ) VSIFCloseL(fpOut); fpOut = VSIFOpenL( (CPLString("/vsigzip/") + CPLResetExtension(osCacheFilename, "esri.gz")).c_str(), "wb"); if( fpOut != nullptr ) { BuildESRICSNameCache(pszSRSType, poMapCSNameToCode, fpOut); VSIFCloseL(fpOut); } return papoSRSCache; } /************************************************************************/ /* FindMatches() */ /************************************************************************/ /** * \brief Try to identify a match between the passed SRS and a related SRS * in a catalog (currently EPSG only) * * Matching may be partial, or may fail. * Returned entries will be sorted by decreasing match confidence (first * entry has the highest match confidence). * * The exact way matching is done may change in future versions. * * The current algorithm is: * - try first AutoIdentifyEPSG(). If it succeeds, return the corresponding SRS * - otherwise iterate over all SRS from the EPSG catalog (as found in GDAL * pcs.csv and gcs.csv files+esri_extra.wkt), and find those that match the * input SRS using the IsSame() function (ignoring TOWGS84 clauses) * - if there is a single match using IsSame() or one of the matches has the * same SRS name, return it with 100% confidence * - if a SRS has the same SRS name, but does not pass the IsSame() criteria, * return it with 50% confidence. * - otherwise return all candidate SRS that pass the IsSame() criteria with a * 90% confidence. * * A pre-built SRS cache in ~/.gdal/X.Y/srs_cache will be used if existing, * otherwise it will be built at the first run of this function. * * This method is the same as OSRFindMatches(). * * @param papszOptions NULL terminated list of options or NULL * @param pnEntries Output parameter. Number of values in the returned array. * @param ppanMatchConfidence Output parameter (or NULL). *ppanMatchConfidence * will be allocated to an array of *pnEntries whose values between 0 and 100 * indicate the confidence in the match. 100 is the highest confidence level. * The array must be freed with CPLFree(). * * @return an array of SRS that match the passed SRS, or NULL. Must be freed with * OSRFreeSRSArray() * * @since GDAL 2.3 */ OGRSpatialReferenceH* OGRSpatialReference::FindMatches( char** papszOptions, int* pnEntries, int** ppanMatchConfidence ) const { CPL_IGNORE_RET_VAL(papszOptions); if( pnEntries ) *pnEntries = 0; if( ppanMatchConfidence ) *ppanMatchConfidence = nullptr; OGRSpatialReference oSRSClone(*this); if( oSRSClone.AutoIdentifyEPSG() == OGRERR_NONE ) { const char* pszCode = oSRSClone.GetAuthorityCode(nullptr); if( pszCode ) oSRSClone.importFromEPSG(atoi(pszCode)); OGRSpatialReferenceH* pahRet = static_cast<OGRSpatialReferenceH*>( CPLCalloc(sizeof(OGRSpatialReferenceH), 2)); pahRet[0] = reinterpret_cast<OGRSpatialReferenceH>(oSRSClone.Clone()); if( pnEntries ) *pnEntries = 1; if( ppanMatchConfidence ) { *ppanMatchConfidence = static_cast<int*>(CPLMalloc(sizeof(int))); (*ppanMatchConfidence)[0] = 100; } return pahRet; } const char* pszSRSType = ""; if( IsProjected() ) { pszSRSType = "PROJCS"; } else if( IsGeographic() ) { pszSRSType = "GEOGCS"; } else { return nullptr; } const char*pszSRSName = GetAttrValue(pszSRSType); if( pszSRSName == nullptr ) return nullptr; const std::map<CPLString, int>* poMapCSNameToCode = nullptr; const std::vector<OGRSpatialReference*>* papoSRSCache = GetSRSCache(pszSRSType, poMapCSNameToCode); if( papoSRSCache == nullptr ) return nullptr; // If we have an exact match with a coordinate system name coming from // EPSG entries (either ours or ESRI), and the SRS are equivalent, then // use that exact match const char* apszOptions[] = { "TOWGS84=ONLY_IF_IN_BOTH", nullptr }; if( poMapCSNameToCode ) { auto oIter = poMapCSNameToCode->find(pszSRSName); if( oIter != poMapCSNameToCode->end() ) { OGRSpatialReference oSRS; if( oSRS.importFromEPSG(oIter->second) == OGRERR_NONE ) { if( IsSame(&oSRS, apszOptions) ) { OGRSpatialReferenceH* pahRet = static_cast<OGRSpatialReferenceH*>( CPLCalloc(sizeof(OGRSpatialReferenceH), 2)); pahRet[0] = reinterpret_cast<OGRSpatialReferenceH>( oSRS.Clone()); if( pnEntries ) *pnEntries = 1; if( ppanMatchConfidence ) { *ppanMatchConfidence = static_cast<int*>( CPLMalloc(sizeof(int))); (*ppanMatchConfidence)[0] = 100; } return pahRet; } } } } std::vector< OGRSpatialReference* > apoSameSRS; CPLString osSRSName(MassageSRSName(pszSRSName, false)); CPLString osSRSNameExtra(MassageSRSName(osSRSName, true)); std::vector<size_t> anMatchingSRSNameIndices; for(size_t i = 0; i < papoSRSCache->size(); i++ ) { const OGRSpatialReference* poOtherSRS = (*papoSRSCache)[i]; #ifdef notdef if( EQUAL(poOtherSRS->GetAuthorityCode(nullptr), "2765") ) { printf("brkpt\n"); /* ok */ } #endif const char* pszOtherSRSName = poOtherSRS->GetAttrValue(pszSRSType); if( pszOtherSRSName == nullptr ) continue; CPLString osOtherSRSName( MassageSRSName(pszOtherSRSName, false)); if( EQUAL( osSRSName, osOtherSRSName ) ) { anMatchingSRSNameIndices.push_back(i); } if( IsSame(poOtherSRS, apszOptions) ) { apoSameSRS.push_back( poOtherSRS->Clone() ); } } const size_t nSameCount = apoSameSRS.size(); if( nSameCount == 1 ) { OGRSpatialReferenceH* pahRet = static_cast<OGRSpatialReferenceH*>( CPLCalloc(sizeof(OGRSpatialReferenceH), 2)); pahRet[0] = reinterpret_cast<OGRSpatialReferenceH>(apoSameSRS[0]); if( pnEntries ) *pnEntries = 1; if( ppanMatchConfidence ) { *ppanMatchConfidence = static_cast<int*>(CPLMalloc(sizeof(int))); (*ppanMatchConfidence)[0] = 100; } return pahRet; } int nCountExtraMatches = 0; size_t iExtraMatch = 0; for(size_t i=0; i<nSameCount; i++) { CPLString osOtherSRSName( MassageSRSName(apoSameSRS[i]->GetAttrValue(pszSRSType), false)); CPLString osOtherSRSNameExtra(MassageSRSName(osOtherSRSName, true)); if( EQUAL(osSRSName, osOtherSRSName) ) { OGRSpatialReferenceH* pahRet = static_cast<OGRSpatialReferenceH*>( CPLCalloc(sizeof(OGRSpatialReferenceH), 2)); pahRet[0] = reinterpret_cast<OGRSpatialReferenceH>(apoSameSRS[i]); if( pnEntries ) *pnEntries = 1; if( ppanMatchConfidence ) { *ppanMatchConfidence = static_cast<int*>(CPLMalloc(sizeof(int))); (*ppanMatchConfidence)[0] = 100; } for(size_t j=0; j<nSameCount; j++) { if( i != j ) delete apoSameSRS[j]; } return pahRet; } else if ( EQUAL(osSRSNameExtra, osOtherSRSNameExtra) ) { nCountExtraMatches ++; iExtraMatch = i; } } if( nCountExtraMatches == 1 ) { OGRSpatialReferenceH* pahRet = static_cast<OGRSpatialReferenceH*>( CPLCalloc(sizeof(OGRSpatialReferenceH), 2)); pahRet[0] = reinterpret_cast<OGRSpatialReferenceH>(apoSameSRS[iExtraMatch]); if( pnEntries ) *pnEntries = 1; if( ppanMatchConfidence ) { *ppanMatchConfidence = static_cast<int*>(CPLMalloc(sizeof(int))); (*ppanMatchConfidence)[0] = 100; } for(size_t j=0; j<nSameCount; j++) { if( iExtraMatch != j ) delete apoSameSRS[j]; } return pahRet; } if( nSameCount == 0 && anMatchingSRSNameIndices.size() == 1 ) { const OGRSpatialReference* poOtherSRS = (*papoSRSCache)[anMatchingSRSNameIndices[0]]; OGRSpatialReferenceH* pahRet = static_cast<OGRSpatialReferenceH*>( CPLCalloc(sizeof(OGRSpatialReferenceH), 2)); pahRet[0] = reinterpret_cast<OGRSpatialReferenceH>(poOtherSRS->Clone()); if( pnEntries ) *pnEntries = 1; if( ppanMatchConfidence ) { *ppanMatchConfidence = static_cast<int*>(CPLMalloc(sizeof(int))); (*ppanMatchConfidence)[0] = 50; } return pahRet; } if( nSameCount == 0 ) return nullptr; if( pnEntries ) *pnEntries = static_cast<int>(nSameCount); OGRSpatialReferenceH* pahRet = static_cast<OGRSpatialReferenceH*>( CPLCalloc(sizeof(OGRSpatialReferenceH), nSameCount + 1)); if( ppanMatchConfidence ) { *ppanMatchConfidence = static_cast<int*>( CPLMalloc(sizeof(int) * (nSameCount + 1))); } for(size_t i=0; i<nSameCount; i++) { pahRet[i] = reinterpret_cast<OGRSpatialReferenceH>(apoSameSRS[i]); if( ppanMatchConfidence ) (*ppanMatchConfidence)[i] = 90; // Arbitrary... } pahRet[ nSameCount ] = nullptr; return pahRet; }