diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/c_api.cpp | 17 | ||||
| -rw-r--r-- | src/coordinateoperation.cpp | 185 | ||||
| -rw-r--r-- | src/coordinatesystem.cpp | 2 | ||||
| -rw-r--r-- | src/crs.cpp | 125 | ||||
| -rw-r--r-- | src/internal.cpp | 10 | ||||
| -rw-r--r-- | src/io.cpp | 394 | ||||
| -rw-r--r-- | src/pj_ellps.c | 2 |
7 files changed, 571 insertions, 164 deletions
diff --git a/src/c_api.cpp b/src/c_api.cpp index 1720a7a1..fbba0f51 100644 --- a/src/c_api.cpp +++ b/src/c_api.cpp @@ -839,7 +839,10 @@ static const char *getOptionValue(const char *option, * <li>MULTILINE=YES/NO. Defaults to YES, except for WKT1_ESRI</li> * <li>INDENTATION_WIDTH=number. Defauls to 4 (when multiline output is * on).</li> - * <li>OUTPUT_AXIS=YES/NO. Defaults to YES, except for WKT1_ESRI.</li> + * <li>OUTPUT_AXIS=AUTO/YES/NO. In AUTO mode, axis will be output for WKT2 + * variants, for WKT1_GDAL for ProjectedCRS with easting/northing ordering + * (otherwise stripped), but not for WKT1_ESRI. Setting to YES will output + * them unconditionaly, and to NO will omit them unconditionaly.</li> * </ul> * @return a string, or NULL in case of error. */ @@ -890,7 +893,12 @@ const char *proj_obj_as_wkt(PJ_OBJ *obj, PJ_WKT_TYPE type, } else if ((value = getOptionValue(*iter, "INDENTATION_WIDTH="))) { formatter->setIndentationWidth(std::atoi(value)); } else if ((value = getOptionValue(*iter, "OUTPUT_AXIS="))) { - formatter->setOutputAxis(ci_equal(value, "YES")); + if (!ci_equal(value, "AUTO")) { + formatter->setOutputAxis( + ci_equal(value, "YES") + ? WKTFormatter::OutputAxisRule::YES + : WKTFormatter::OutputAxisRule::NO); + } } else { std::string msg("Unknown option :"); msg += *iter; @@ -925,7 +933,8 @@ const char *proj_obj_as_wkt(PJ_OBJ *obj, PJ_WKT_TYPE type, * @param options NULL-terminated list of strings with "KEY=VALUE" format. or * NULL. * The currently recognized option is USE_ETMERC=YES to use - * +proj=etmerc instead of +proj=tmerc + * +proj=etmerc instead of +proj=tmerc (or USE_ETMERC=NO to disable implicit + * use of etmerc by utm conversions) * @return a string, or NULL in case of error. */ const char *proj_obj_as_proj_string(PJ_OBJ *obj, PJ_PROJ_STRING_TYPE type, @@ -960,6 +969,8 @@ const char *proj_obj_as_proj_string(PJ_OBJ *obj, PJ_PROJ_STRING_TYPE type, if (options != nullptr && options[0] != nullptr) { if (ci_equal(options[0], "USE_ETMERC=YES")) { formatter->setUseETMercForTMerc(true); + } else if (ci_equal(options[0], "USE_ETMERC=NO")) { + formatter->setUseETMercForTMerc(false); } } obj->lastPROJString = exportable->exportToPROJString(formatter.get()); diff --git a/src/coordinateoperation.cpp b/src/coordinateoperation.cpp index 33e6f422..42cc806e 100644 --- a/src/coordinateoperation.cpp +++ b/src/coordinateoperation.cpp @@ -143,7 +143,7 @@ static std::set<std::string> buildSetEquivalentParameters() { std::set<std::string> set; - const char *const listOfEquivalentParameterNames[][5] = { + const char *const listOfEquivalentParameterNames[][7] = { {"latitude_of_point_1", "Latitude_Of_1st_Point", nullptr}, {"longitude_of_point_1", "Longitude_Of_1st_Point", nullptr}, {"latitude_of_point_2", "Latitude_Of_2nd_Point", nullptr}, @@ -157,11 +157,13 @@ static std::set<std::string> buildSetEquivalentParameters() { EPSG_NAME_PARAMETER_NORTHING_FALSE_ORIGIN, EPSG_NAME_PARAMETER_NORTHING_PROJECTION_CENTRE, nullptr}, - {EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN, + {WKT1_LATITUDE_OF_ORIGIN, WKT1_LATITUDE_OF_CENTER, + EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN, EPSG_NAME_PARAMETER_LATITUDE_FALSE_ORIGIN, EPSG_NAME_PARAMETER_LATITUDE_PROJECTION_CENTRE, nullptr}, - {EPSG_NAME_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN, + {WKT1_CENTRAL_MERIDIAN, WKT1_LONGITUDE_OF_CENTER, + EPSG_NAME_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN, EPSG_NAME_PARAMETER_LONGITUDE_FALSE_ORIGIN, EPSG_NAME_PARAMETER_LONGITUDE_PROJECTION_CENTRE, EPSG_NAME_PARAMETER_LONGITUDE_OF_ORIGIN, nullptr}, @@ -179,8 +181,7 @@ static std::set<std::string> buildSetEquivalentParameters() { return set; } -static bool areEquivalentParameters(const std::string &a, - const std::string &b) { +bool areEquivalentParameters(const std::string &a, const std::string &b) { static const std::set<std::string> setEquivalentParameters = buildSetEquivalentParameters(); @@ -690,6 +691,7 @@ struct OperationMethod::Private { util::optional<std::string> formula_{}; util::optional<metadata::Citation> formulaCitation_{}; std::vector<GeneralOperationParameterNNPtr> parameters_{}; + std::string projMethodOverride_{}; }; //! @endcond @@ -767,6 +769,7 @@ OperationMethodNNPtr OperationMethod::create( method->assignSelf(method); method->setProperties(properties); method->d->parameters_ = parameters; + properties.getStringValue("proj_method", method->d->projMethodOverride_); return method; } @@ -823,12 +826,17 @@ void OperationMethod::_exportToWKT(io::WKTFormatter *formatter) const { if (mapping == nullptr) { l_name = replaceAll(l_name, " ", "_"); } else { - if (mapping->wkt1_name == nullptr) { - throw io::FormattingException( - std::string("Unsupported conversion method: ") + - mapping->wkt2_name); + if (l_name == + PROJ_WKT2_NAME_METHOD_GEOSTATIONARY_SATELLITE_SWEEP_X) { + l_name = "Geostationary_Satellite"; + } else { + if (mapping->wkt1_name == nullptr) { + throw io::FormattingException( + std::string("Unsupported conversion method: ") + + mapping->wkt2_name); + } + l_name = mapping->wkt1_name; } - l_name = mapping->wkt1_name; } } formatter->addQuotedString(l_name); @@ -2870,7 +2878,7 @@ ConversionNNPtr Conversion::createEquidistantCylindrical( const common::Angle &longitudeNatOrigin, const common::Length &falseEasting, const common::Length &falseNorthing) { return create(properties, EPSG_CODE_METHOD_EQUIDISTANT_CYLINDRICAL, - createParams(latitudeFirstParallel, longitudeNatOrigin, + createParams(latitudeFirstParallel, 0.0, longitudeNatOrigin, falseEasting, falseNorthing)); } @@ -2905,7 +2913,7 @@ ConversionNNPtr Conversion::createEquidistantCylindricalSpherical( const common::Length &falseNorthing) { return create(properties, EPSG_CODE_METHOD_EQUIDISTANT_CYLINDRICAL_SPHERICAL, - createParams(latitudeFirstParallel, longitudeNatOrigin, + createParams(latitudeFirstParallel, 0.0, longitudeNatOrigin, falseEasting, falseNorthing)); } @@ -4781,8 +4789,30 @@ void Conversion::_exportToWKT(io::WKTFormatter *formatter) const { const MethodMapping *mapping = !isWKT2 ? getMapping(l_method.get()) : nullptr; - for (const auto ¶mValue : parameterValues()) { - paramValue->_exportToWKT(formatter, mapping); + for (const auto &genOpParamvalue : parameterValues()) { + + // EPSG has normally no Latitude of natural origin for Equidistant + // Cylindrical but PROJ can handle it, so output the parameter if + // not zero + if ((methodEPSGCode == EPSG_CODE_METHOD_EQUIDISTANT_CYLINDRICAL || + methodEPSGCode == + EPSG_CODE_METHOD_EQUIDISTANT_CYLINDRICAL_SPHERICAL)) { + auto opParamvalue = + dynamic_cast<const OperationParameterValue *>( + genOpParamvalue.get()); + if (opParamvalue && + opParamvalue->parameter()->getEPSGCode() == + EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN) { + const auto ¶mValue = opParamvalue->parameterValue(); + if (paramValue->type() == ParameterValue::Type::MEASURE) { + const auto &measure = paramValue->value(); + if (measure.getSIValue() == 0) { + continue; + } + } + } + } + genOpParamvalue->_exportToWKT(formatter, mapping); } } @@ -4899,14 +4929,29 @@ createPROJExtensionFromCustomProj(const Conversion *conv, // --------------------------------------------------------------------------- -void Conversion::addWKTExtensionNode(io::WKTFormatter *formatter) const { +bool Conversion::addWKTExtensionNode(io::WKTFormatter *formatter) const { const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2; if (!isWKT2) { - const auto &methodName = method()->nameStr(); - const int methodEPSGCode = method()->getEPSGCode(); - if (methodEPSGCode == - EPSG_CODE_METHOD_POPULAR_VISUALISATION_PSEUDO_MERCATOR || - nameStr() == "Popular Visualisation Mercator") { + const auto &l_method = method(); + const auto &methodName = l_method->nameStr(); + const int methodEPSGCode = l_method->getEPSGCode(); + int zone = 0; + bool north = true; + if (l_method->getPrivate()->projMethodOverride_ == "etmerc" && + !isUTM(zone, north)) { + auto projFormatter = io::PROJStringFormatter::create( + io::PROJStringFormatter::Convention::PROJ_4); + projFormatter->setUseETMercForTMerc(true); + formatter->startNode(io::WKTConstants::EXTENSION, false); + formatter->addQuotedString("PROJ4"); + _exportToPROJString(projFormatter.get()); + projFormatter->addParam("no_defs"); + formatter->addQuotedString(projFormatter->toString()); + formatter->endNode(); + return true; + } else if (methodEPSGCode == + EPSG_CODE_METHOD_POPULAR_VISUALISATION_PSEUDO_MERCATOR || + nameStr() == "Popular Visualisation Mercator") { auto projFormatter = io::PROJStringFormatter::create( io::PROJStringFormatter::Convention::PROJ_4); @@ -4915,6 +4960,7 @@ void Conversion::addWKTExtensionNode(io::WKTFormatter *formatter) const { formatter->addQuotedString("PROJ4"); formatter->addQuotedString(projFormatter->toString()); formatter->endNode(); + return true; } } else if (starts_with(methodName, "PROJ ")) { auto projFormatter = io::PROJStringFormatter::create( @@ -4925,9 +4971,22 @@ void Conversion::addWKTExtensionNode(io::WKTFormatter *formatter) const { formatter->addQuotedString("PROJ4"); formatter->addQuotedString(projFormatter->toString()); formatter->endNode(); + return true; } + } else if (methodName == + PROJ_WKT2_NAME_METHOD_GEOSTATIONARY_SATELLITE_SWEEP_X) { + auto projFormatter = io::PROJStringFormatter::create( + io::PROJStringFormatter::Convention::PROJ_4); + formatter->startNode(io::WKTConstants::EXTENSION, false); + formatter->addQuotedString("PROJ4"); + _exportToPROJString(projFormatter.get()); + projFormatter->addParam("no_defs"); + formatter->addQuotedString(projFormatter->toString()); + formatter->endNode(); + return true; } } + return false; } // --------------------------------------------------------------------------- @@ -4981,15 +5040,16 @@ void Conversion::_exportToPROJString( // Check for UTM int zone = 0; bool north = true; - if (isUTM(zone, north)) { + bool etMercSettingSet = false; + useETMerc = formatter->getUseETMercForTMerc(etMercSettingSet) || + l_method->getPrivate()->projMethodOverride_ == "etmerc"; + if (isUTM(zone, north) && !(etMercSettingSet && !useETMerc)) { bConversionDone = true; formatter->addStep("utm"); formatter->addParam("zone", zone); if (!north) { formatter->addParam("south"); } - } else { - useETMerc = formatter->getUseETMercForTMerc(); } } else if (methodEPSGCode == EPSG_CODE_METHOD_HOTINE_OBLIQUE_MERCATOR_VARIANT_A) { @@ -5077,6 +5137,21 @@ void Conversion::_exportToPROJString( std::string("Unsupported value for ") + EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN); } + } else if (methodEPSGCode == EPSG_CODE_METHOD_MERCATOR_VARIANT_B) { + const auto &scaleFactor = parameterValueMeasure(WKT1_SCALE_FACTOR, 0); + if (scaleFactor.unit().type() != common::UnitOfMeasure::Type::UNKNOWN && + std::fabs(scaleFactor.getSIValue() - 1.0) > 1e-10) { + throw io::FormattingException( + "Unexpected presence of scale factor in Mercator (variant B)"); + } + double latitudeOrigin = parameterValueNumeric( + EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN, + common::UnitOfMeasure::DEGREE); + if (latitudeOrigin != 0) { + throw io::FormattingException( + std::string("Unsupported value for ") + + EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN); + } // PROJ.4 specific hack for webmercator } else if (formatter->convention() == io::PROJStringFormatter::Convention::PROJ_4 && @@ -5265,12 +5340,15 @@ bool Conversion::isUTM(int &zone, bool &north) const { const auto &measure = l_parameterValue->value(); if (epsg_code == EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN && - measure.value() == UTM_LATITUDE_OF_NATURAL_ORIGIN) { + std::fabs(measure.value() - + UTM_LATITUDE_OF_NATURAL_ORIGIN) < 1e-10) { bLatitudeNatOriginUTM = true; } else if ( epsg_code == EPSG_CODE_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN && - measure.unit() == common::UnitOfMeasure::DEGREE) { + measure.unit()._isEquivalentTo( + common::UnitOfMeasure::DEGREE, + util::IComparable::Criterion::EQUIVALENT)) { double dfZone = (measure.value() + 183.0) / 6.0; if (dfZone > 0.9 && dfZone < 60.1 && std::abs(dfZone - std::round(dfZone)) < 1e-10) { @@ -5279,20 +5357,29 @@ bool Conversion::isUTM(int &zone, bool &north) const { } else if ( epsg_code == EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN && - measure.value() == UTM_SCALE_FACTOR) { + measure.unit()._isEquivalentTo( + common::UnitOfMeasure::SCALE_UNITY, + util::IComparable::Criterion::EQUIVALENT) && + std::fabs(measure.value() - UTM_SCALE_FACTOR) < 1e-10) { bScaleFactorUTM = true; } else if (epsg_code == EPSG_CODE_PARAMETER_FALSE_EASTING && measure.value() == UTM_FALSE_EASTING && - measure.unit() == common::UnitOfMeasure::METRE) { + measure.unit()._isEquivalentTo( + common::UnitOfMeasure::METRE, + util::IComparable::Criterion::EQUIVALENT)) { bFalseEastingUTM = true; } else if (epsg_code == EPSG_CODE_PARAMETER_FALSE_NORTHING && - measure.unit() == common::UnitOfMeasure::METRE) { - if (measure.value() == UTM_NORTH_FALSE_NORTHING) { + measure.unit()._isEquivalentTo( + common::UnitOfMeasure::METRE, + util::IComparable::Criterion::EQUIVALENT)) { + if (std::fabs(measure.value() - + UTM_NORTH_FALSE_NORTHING) < 1e-10) { bFalseNorthingUTM = true; north = true; - } else if (measure.value() == - UTM_SOUTH_FALSE_NORTHING) { + } else if (std::fabs(measure.value() - + UTM_SOUTH_FALSE_NORTHING) < + 1e-10) { bFalseNorthingUTM = true; north = false; } @@ -6197,32 +6284,36 @@ TransformationNNPtr Transformation::createTOWGS84( "Invalid number of elements in TOWGS84Parameters"); } - crs::CRSPtr transformSourceCRS = sourceCRSIn->extractGeographicCRS(); + crs::CRSPtr transformSourceCRS = sourceCRSIn->extractGeodeticCRS(); if (!transformSourceCRS) { throw InvalidOperation( - "Cannot find GeographicCRS in sourceCRS of TOWGS84 transformation"); + "Cannot find GeodeticCRS in sourceCRS of TOWGS84 transformation"); } + util::PropertyMap properties; + properties.set(common::IdentifiedObject::NAME_KEY, + concat("Transformation from ", transformSourceCRS->nameStr(), + " to WGS84")); + + auto targetCRS = + dynamic_cast<const crs::GeographicCRS *>(transformSourceCRS.get()) + ? util::nn_static_pointer_cast<crs::CRS>( + crs::GeographicCRS::EPSG_4326) + : util::nn_static_pointer_cast<crs::CRS>( + crs::GeodeticCRS::EPSG_4978); + if (TOWGS84Parameters.size() == 3) { return createGeocentricTranslations( - util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, - concat("Transformation from ", - transformSourceCRS->nameStr(), - " to WGS84")), - NN_NO_CHECK(transformSourceCRS), crs::GeographicCRS::EPSG_4326, + properties, NN_NO_CHECK(transformSourceCRS), targetCRS, TOWGS84Parameters[0], TOWGS84Parameters[1], TOWGS84Parameters[2], {}); } - return createPositionVector( - util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, - concat("Transformation from ", - transformSourceCRS->nameStr(), - " to WGS84")), - NN_NO_CHECK(transformSourceCRS), crs::GeographicCRS::EPSG_4326, - TOWGS84Parameters[0], TOWGS84Parameters[1], TOWGS84Parameters[2], - TOWGS84Parameters[3], TOWGS84Parameters[4], TOWGS84Parameters[5], - TOWGS84Parameters[6], {}); + return createPositionVector(properties, NN_NO_CHECK(transformSourceCRS), + targetCRS, TOWGS84Parameters[0], + TOWGS84Parameters[1], TOWGS84Parameters[2], + TOWGS84Parameters[3], TOWGS84Parameters[4], + TOWGS84Parameters[5], TOWGS84Parameters[6], {}); } // --------------------------------------------------------------------------- diff --git a/src/coordinatesystem.cpp b/src/coordinatesystem.cpp index 1c84cf2c..f8d2d2b3 100644 --- a/src/coordinatesystem.cpp +++ b/src/coordinatesystem.cpp @@ -460,7 +460,7 @@ CoordinateSystem::axisList() PROJ_CONST_DEFN { void CoordinateSystem::_exportToWKT( io::WKTFormatter *formatter) const // throw(FormattingException) { - if (!formatter->outputAxis()) { + if (formatter->outputAxis() != io::WKTFormatter::OutputAxisRule::YES) { return; } const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2; diff --git a/src/crs.cpp b/src/crs.cpp index dab704b4..a204a037 100644 --- a/src/crs.cpp +++ b/src/crs.cpp @@ -87,6 +87,7 @@ namespace crs { //! @cond Doxygen_Suppress struct CRS::Private { BoundCRSPtr canonicalBoundCRS_{}; + std::string extensionProj4_{}; }; //! @endcond @@ -745,6 +746,8 @@ GeodeticCRS::create(const util::PropertyMap &properties, GeodeticCRS::nn_make_shared<GeodeticCRS>(datum, datumEnsemble, cs)); crs->assignSelf(crs); crs->setProperties(properties); + properties.getStringValue("EXTENSION_PROJ4", + crs->CRS::getPrivate()->extensionProj4_); return crs; } @@ -789,6 +792,8 @@ GeodeticCRS::create(const util::PropertyMap &properties, GeodeticCRS::nn_make_shared<GeodeticCRS>(datum, datumEnsemble, cs)); crs->assignSelf(crs); crs->setProperties(properties); + properties.getStringValue("EXTENSION_PROJ4", + crs->CRS::getPrivate()->extensionProj4_); return crs; } @@ -850,6 +855,17 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const { } cs->_exportToWKT(formatter); ObjectUsage::baseExportToWKT(formatter); + + if (!isWKT2 && !formatter->useESRIDialect()) { + const auto &extensionProj4 = CRS::getPrivate()->extensionProj4_; + if (!extensionProj4.empty()) { + formatter->startNode(io::WKTConstants::EXTENSION, false); + formatter->addQuotedString("PROJ4"); + formatter->addQuotedString(extensionProj4); + formatter->endNode(); + } + } + formatter->endNode(); } //! @endcond @@ -862,7 +878,8 @@ void GeodeticCRS::addGeocentricUnitConversionIntoPROJString( const auto &axisList = coordinateSystem()->axisList(); const auto &unit = axisList[0]->unit(); - if (unit != common::UnitOfMeasure::METRE) { + if (!unit._isEquivalentTo(common::UnitOfMeasure::METRE, + util::IComparable::Criterion::EQUIVALENT)) { if (formatter->convention() == io::PROJStringFormatter::Convention::PROJ_4) { io::FormattingException::Throw("GeodeticCRS::exportToPROJString(): " @@ -885,6 +902,9 @@ void GeodeticCRS::addGeocentricUnitConversionIntoPROJString( const auto &toSI = unit.conversionToSI(); formatter->addParam("xy_out", toSI); formatter->addParam("z_out", toSI); + } else if (formatter->convention() == + io::PROJStringFormatter::Convention::PROJ_4) { + formatter->addParam("units", "m"); } } //! @endcond @@ -895,6 +915,16 @@ void GeodeticCRS::addGeocentricUnitConversionIntoPROJString( void GeodeticCRS::_exportToPROJString( io::PROJStringFormatter *formatter) const // throw(io::FormattingException) { + if (formatter->convention() == + io::PROJStringFormatter::Convention::PROJ_4) { + const auto &extensionProj4 = CRS::getPrivate()->extensionProj4_; + if (!extensionProj4.empty()) { + formatter->ingestPROJString(extensionProj4); + formatter->addNoDefs(false); + return; + } + } + if (!isGeocentric()) { io::FormattingException::Throw( "GeodeticCRS::exportToPROJString() only " @@ -907,14 +937,7 @@ void GeodeticCRS::_exportToPROJString( } else { formatter->addStep("cart"); } - ellipsoid()->_exportToPROJString(formatter); - if (formatter->convention() == - io::PROJStringFormatter::Convention::PROJ_4) { - const auto &TOWGS84Params = formatter->getTOWGS84Parameters(); - if (TOWGS84Params.size() == 7) { - formatter->addParam("towgs84", TOWGS84Params); - } - } + addDatumInfoToPROJString(formatter); addGeocentricUnitConversionIntoPROJString(formatter); } //! @endcond @@ -1421,6 +1444,8 @@ GeographicCRS::create(const util::PropertyMap &properties, GeographicCRS::nn_make_shared<GeographicCRS>(datum, datumEnsemble, cs)); crs->assignSelf(crs); crs->setProperties(properties); + properties.getStringValue("EXTENSION_PROJ4", + crs->CRS::getPrivate()->extensionProj4_); return crs; } @@ -1626,6 +1651,16 @@ void GeographicCRS::addAngularUnitConvertAndAxisSwap( void GeographicCRS::_exportToPROJString( io::PROJStringFormatter *formatter) const // throw(io::FormattingException) { + if (formatter->convention() == + io::PROJStringFormatter::Convention::PROJ_4) { + const auto &extensionProj4 = CRS::getPrivate()->extensionProj4_; + if (!extensionProj4.empty()) { + formatter->ingestPROJString(extensionProj4); + formatter->addNoDefs(false); + return; + } + } + if (!formatter->omitProjLongLatIfPossible() || primeMeridian()->longitude().getSIValue() != 0.0 || !formatter->getTOWGS84Parameters().empty() || @@ -1763,7 +1798,15 @@ void VerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const { if (!isWKT2) { axisList[0]->unit()._exportToWKT(formatter); } + + const auto oldAxisOutputRule = formatter->outputAxis(); + if (oldAxisOutputRule == + io::WKTFormatter::OutputAxisRule::WKT1_GDAL_EPSG_STYLE) { + formatter->setOutputAxis(io::WKTFormatter::OutputAxisRule::YES); + } cs->_exportToWKT(formatter); + formatter->setOutputAxis(oldAxisOutputRule); + ObjectUsage::baseExportToWKT(formatter); formatter->endNode(); } @@ -2231,6 +2274,22 @@ const cs::CartesianCSNNPtr &ProjectedCRS::coordinateSystem() PROJ_CONST_DEFN { void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const { const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2; + const auto &l_coordinateSystem = d->coordinateSystem(); + const auto &axisList = l_coordinateSystem->axisList(); + + const auto exportAxis = [&l_coordinateSystem, &axisList, &formatter]() { + const auto oldAxisOutputRule = formatter->outputAxis(); + if (oldAxisOutputRule == + io::WKTFormatter::OutputAxisRule::WKT1_GDAL_EPSG_STYLE) { + if (&axisList[0]->direction() == &cs::AxisDirection::EAST && + &axisList[1]->direction() == &cs::AxisDirection::NORTH) { + formatter->setOutputAxis(io::WKTFormatter::OutputAxisRule::YES); + } + } + l_coordinateSystem->_exportToWKT(formatter); + formatter->setOutputAxis(oldAxisOutputRule); + }; + if (!isWKT2 && !formatter->useESRIDialect() && starts_with(nameStr(), "Popular Visualisation CRS / Mercator")) { formatter->startNode(io::WKTConstants::PROJCS, !identifiers().empty()); @@ -2263,9 +2322,8 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const { formatter->add(0.0); formatter->endNode(); - const auto &axisList = d->coordinateSystem()->axisList(); axisList[0]->unit()._exportToWKT(formatter); - d->coordinateSystem()->_exportToWKT(formatter); + exportAxis(); derivingConversionRef()->addWKTExtensionNode(formatter); ObjectUsage::baseExportToWKT(formatter); formatter->endNode(); @@ -2321,7 +2379,6 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const { l_baseCRS->_exportToWKT(formatter); } - const auto &axisList = d->coordinateSystem()->axisList(); formatter->pushAxisLinearUnit( common::UnitOfMeasure::create(axisList[0]->unit())); @@ -2338,10 +2395,18 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const { axisList[0]->unit()._exportToWKT(formatter); } - d->coordinateSystem()->_exportToWKT(formatter); + exportAxis(); if (!isWKT2 && !formatter->useESRIDialect()) { - derivingConversionRef()->addWKTExtensionNode(formatter); + const auto &extensionProj4 = CRS::getPrivate()->extensionProj4_; + if (!extensionProj4.empty()) { + formatter->startNode(io::WKTConstants::EXTENSION, false); + formatter->addQuotedString("PROJ4"); + formatter->addQuotedString(extensionProj4); + formatter->endNode(); + } else { + derivingConversionRef()->addWKTExtensionNode(formatter); + } } ObjectUsage::baseExportToWKT(formatter); @@ -2355,6 +2420,16 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const { void ProjectedCRS::_exportToPROJString( io::PROJStringFormatter *formatter) const // throw(io::FormattingException) { + if (formatter->convention() == + io::PROJStringFormatter::Convention::PROJ_4) { + const auto &extensionProj4 = CRS::getPrivate()->extensionProj4_; + if (!extensionProj4.empty()) { + formatter->ingestPROJString(extensionProj4); + formatter->addNoDefs(false); + return; + } + } + baseExportToPROJString(formatter); } @@ -2384,6 +2459,8 @@ ProjectedCRS::create(const util::PropertyMap &properties, crs->assignSelf(crs); crs->setProperties(properties); crs->setDerivingConversionCRS(); + properties.getStringValue("EXTENSION_PROJ4", + crs->CRS::getPrivate()->extensionProj4_); return crs; } @@ -2404,7 +2481,8 @@ void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter, bool axisSpecFound) const { const auto &axisList = d->coordinateSystem()->axisList(); const auto &unit = axisList[0]->unit(); - if (unit != common::UnitOfMeasure::METRE) { + if (!unit._isEquivalentTo(common::UnitOfMeasure::METRE, + util::IComparable::Criterion::EQUIVALENT)) { auto projUnit = unit.exportToPROJString(); const double toSI = unit.conversionToSI(); if (formatter->convention() == @@ -2426,6 +2504,12 @@ void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter, formatter->addParam("units", projUnit); } } + } else if (formatter->convention() == + io::PROJStringFormatter::Convention::PROJ_4) { + // could come from the hardcoded def of webmerc + if (!formatter->hasParam("units")) { + formatter->addParam("units", "m"); + } } if (formatter->convention() == @@ -3234,8 +3318,17 @@ BoundCRS::create(const CRSNNPtr &baseCRSIn, const CRSNNPtr &hubCRSIn, BoundCRSNNPtr BoundCRS::createFromTOWGS84(const CRSNNPtr &baseCRSIn, const std::vector<double> &TOWGS84Parameters) { + + auto geodCRS = baseCRSIn->extractGeodeticCRS(); + auto targetCRS = + geodCRS.get() == nullptr || + dynamic_cast<const crs::GeographicCRS *>(geodCRS.get()) + ? util::nn_static_pointer_cast<crs::CRS>( + crs::GeographicCRS::EPSG_4326) + : util::nn_static_pointer_cast<crs::CRS>( + crs::GeodeticCRS::EPSG_4978); return create( - baseCRSIn, GeographicCRS::EPSG_4326, + baseCRSIn, targetCRS, operation::Transformation::createTOWGS84(baseCRSIn, TOWGS84Parameters)); } diff --git a/src/internal.cpp b/src/internal.cpp index aa27192a..4bec1bf9 100644 --- a/src/internal.cpp +++ b/src/internal.cpp @@ -107,6 +107,16 @@ bool ci_equal(const char *a, const char *b) noexcept { // --------------------------------------------------------------------------- +bool ci_less(const std::string &a, const std::string &b) noexcept { +#ifdef _MSC_VER + return _stricmp(a.c_str(), b.c_str()) < 0; +#else + return strcasecmp(a.c_str(), b.c_str()) < 0; +#endif +} + +// --------------------------------------------------------------------------- + /** * Convert to lower case. */ @@ -128,7 +128,7 @@ struct WKTFormatter::Private { bool primeMeridianInDegree_ = false; bool use2018Keywords_ = false; bool useESRIDialect_ = false; - bool outputAxis_ = true; + OutputAxisRule outputAxis_ = WKTFormatter::OutputAxisRule::YES; }; Params params_{}; DatabaseContextPtr dbContext_{}; @@ -225,10 +225,9 @@ WKTFormatter &WKTFormatter::setIndentationWidth(int width) noexcept { // --------------------------------------------------------------------------- /** \brief Set whether AXIS nodes should be output. - * - * This can typically be set to false for some variants of WKT1_GDAL. */ -WKTFormatter &WKTFormatter::setOutputAxis(bool outputAxisIn) noexcept { +WKTFormatter & +WKTFormatter::setOutputAxis(OutputAxisRule outputAxisIn) noexcept { d->params_.outputAxis_ = outputAxisIn; return *this; } @@ -308,6 +307,8 @@ WKTFormatter::WKTFormatter(Convention convention) d->params_.outputAxisOrder_ = false; d->params_.forceUNITKeyword_ = true; d->params_.primeMeridianInDegree_ = true; + d->params_.outputAxis_ = + WKTFormatter::OutputAxisRule::WKT1_GDAL_EPSG_STYLE; break; case Convention::WKT1_ESRI: @@ -317,7 +318,7 @@ WKTFormatter::WKTFormatter(Convention convention) d->params_.primeMeridianInDegree_ = true; d->params_.useESRIDialect_ = true; d->params_.multiLine_ = false; - d->params_.outputAxis_ = false; + d->params_.outputAxis_ = WKTFormatter::OutputAxisRule::NO; break; default: @@ -568,7 +569,9 @@ const UnitOfMeasureNNPtr &WKTFormatter::axisAngularUnit() const { // --------------------------------------------------------------------------- -bool WKTFormatter::outputAxis() const { return d->params_.outputAxis_; } +WKTFormatter::OutputAxisRule WKTFormatter::outputAxis() const { + return d->params_.outputAxis_; +} // --------------------------------------------------------------------------- @@ -1185,6 +1188,9 @@ struct WKTParser::Private { const UnitOfMeasure &defaultLinearUnit, const UnitOfMeasure &defaultAngularUnit); + static void addExtensionProj4ToProp(const WKTNode::Private *nodeP, + PropertyMap &props); + ConversionNNPtr buildConversion(const WKTNodeNNPtr &node, const UnitOfMeasure &defaultLinearUnit, const UnitOfMeasure &defaultAngularUnit); @@ -1192,6 +1198,9 @@ struct WKTParser::Private { static bool hasWebMercPROJ4String(const WKTNodeNNPtr &projCRSNode, const WKTNodeNNPtr &projectionNode); + static bool projectionHasParameter(const WKTNodeNNPtr &projCRSNode, + const char *paramName); + ConversionNNPtr buildProjection(const WKTNodeNNPtr &projCRSNode, const WKTNodeNNPtr &projectionNode, const UnitOfMeasure &defaultLinearUnit, @@ -1940,6 +1949,7 @@ GeodeticReferenceFrameNNPtr WKTParser::Private::buildGeodeticReferenceFrame( throw ParsingException("Invalid TOWGS84 node"); } } + auto &extensionNode = nodeP->lookForChild(WKTConstants::EXTENSION); const auto &extensionChildren = extensionNode->GP()->children(); if (extensionChildren.size() == 2) { @@ -2403,6 +2413,19 @@ WKTParser::Private::buildCS(const WKTNodeNNPtr &node, /* maybe null */ // --------------------------------------------------------------------------- +void WKTParser::Private::addExtensionProj4ToProp(const WKTNode::Private *nodeP, + PropertyMap &props) { + auto &extensionNode = nodeP->lookForChild(WKTConstants::EXTENSION); + const auto &extensionChildren = extensionNode->GP()->children(); + if (extensionChildren.size() == 2) { + if (ci_equal(stripQuotes(extensionChildren[0]), "PROJ4")) { + props.set("EXTENSION_PROJ4", stripQuotes(extensionChildren[1])); + } + } +} + +// --------------------------------------------------------------------------- + GeodeticCRSNNPtr WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) { const auto *nodeP = node->GP(); @@ -2450,6 +2473,9 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) { angularUnit = primeMeridian->longitude().unit(); } + auto props = buildProperties(node); + addExtensionProj4ToProp(nodeP, props); + auto datum = !isNull(datumNode) ? buildGeodeticReferenceFrame(datumNode, primeMeridian, dynamicNode) @@ -2465,8 +2491,7 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) { if (ellipsoidalCS) { assert(!ci_equal(nodeName, WKTConstants::GEOCCS)); try { - return GeographicCRS::create(buildProperties(node), datum, - datumEnsemble, + return GeographicCRS::create(props, datum, datumEnsemble, NN_NO_CHECK(ellipsoidalCS)); } catch (const util::Exception &e) { throw ParsingException(std::string("buildGeodeticCRS: ") + @@ -2487,8 +2512,8 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) { "Cartesian CS for a GeodeticCRS should have 3 axis"); } try { - return GeodeticCRS::create(buildProperties(node), datum, - datumEnsemble, NN_NO_CHECK(cartesianCS)); + return GeodeticCRS::create(props, datum, datumEnsemble, + NN_NO_CHECK(cartesianCS)); } catch (const util::Exception &e) { throw ParsingException(std::string("buildGeodeticCRS: ") + e.what()); @@ -2498,8 +2523,8 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) { auto sphericalCS = nn_dynamic_pointer_cast<SphericalCS>(cs); if (sphericalCS) { try { - return GeodeticCRS::create(buildProperties(node), datum, - datumEnsemble, NN_NO_CHECK(sphericalCS)); + return GeodeticCRS::create(props, datum, datumEnsemble, + NN_NO_CHECK(sphericalCS)); } catch (const util::Exception &e) { throw ParsingException(std::string("buildGeodeticCRS: ") + e.what()); @@ -2582,7 +2607,8 @@ UnitOfMeasure WKTParser::Private::guessUnitForParameter( ci_find(paramName, "meridian") != std::string::npos || ci_find(paramName, "parallel") != std::string::npos || ci_find(paramName, "azimuth") != std::string::npos || - ci_find(paramName, "angle") != std::string::npos) { + ci_find(paramName, "angle") != std::string::npos || + ci_find(paramName, "heading") != std::string::npos) { unit = defaultAngularUnit; } else if (ci_find(paramName, "easting") != std::string::npos || ci_find(paramName, "northing") != std::string::npos || @@ -2839,8 +2865,15 @@ ConversionNNPtr WKTParser::Private::buildProjectionFromESRI( defaultLinearUnit, defaultAngularUnit); } + struct ci_less_struct { + bool operator()(const std::string &lhs, const std::string &rhs) const + noexcept { + return ci_less(lhs, rhs); + } + }; + // Build a map of present parameters - std::map<std::string, std::string> mapParamNameToValue; + std::map<std::string, std::string, ci_less_struct> mapParamNameToValue; for (const auto &childNode : projCRSNode->GP()->children()) { if (ci_equal(childNode->GP()->value(), WKTConstants::PARAMETER)) { const auto &childNodeChildren = childNode->GP()->children(); @@ -2888,7 +2921,7 @@ ConversionNNPtr WKTParser::Private::buildProjectionFromESRI( const auto *wkt2_mapping = getMapping(esriMapping->wkt2_name); assert(wkt2_mapping); - if (esriProjectionName == "Stereographic") { + if (ci_equal(esriProjectionName, "Stereographic")) { try { if (std::fabs(io::asDouble( mapParamNameToValue["Latitude_Of_Origin"])) == 90.0) { @@ -2937,9 +2970,20 @@ ConversionNNPtr WKTParser::Private::buildProjectionFromESRI( if (iter == mapWKT2NameToESRIName.end()) { continue; } - auto iter2 = mapParamNameToValue.find(iter->second); - if (iter2 == mapParamNameToValue.end()) { - continue; + const auto &esriParamName = iter->second; + auto iter2 = mapParamNameToValue.find(esriParamName); + auto mapParamNameToValueEnd = mapParamNameToValue.end(); + if (iter2 == mapParamNameToValueEnd) { + // In case we don't find a direct match, try the aliases + for (iter2 = mapParamNameToValue.begin(); + iter2 != mapParamNameToValueEnd; ++iter2) { + if (areEquivalentParameters(iter2->first, esriParamName)) { + break; + } + } + if (iter2 == mapParamNameToValueEnd) { + continue; + } } PropertyMap propertiesParameter; @@ -2987,13 +3031,31 @@ WKTParser::Private::buildProjection(const WKTNodeNNPtr &projCRSNode, return buildProjectionStandard(projCRSNode, projectionNode, defaultLinearUnit, defaultAngularUnit); } + +// --------------------------------------------------------------------------- + +bool WKTParser::Private::projectionHasParameter(const WKTNodeNNPtr &projCRSNode, + const char *paramName) { + for (const auto &childNode : projCRSNode->GP()->children()) { + if (ci_equal(childNode->GP()->value(), WKTConstants::PARAMETER)) { + const auto &childNodeChildren = childNode->GP()->children(); + if (childNodeChildren.size() == 2 && + metadata::Identifier::isEquivalentName( + stripQuotes(childNodeChildren[0]).c_str(), paramName)) { + return true; + } + } + } + return false; +} + // --------------------------------------------------------------------------- ConversionNNPtr WKTParser::Private::buildProjectionStandard( const WKTNodeNNPtr &projCRSNode, const WKTNodeNNPtr &projectionNode, const UnitOfMeasure &defaultLinearUnit, const UnitOfMeasure &defaultAngularUnit) { - const std::string wkt1ProjectionName = + std::string wkt1ProjectionName = stripQuotes(projectionNode->GP()->children()[0]); std::vector<OperationParameterNNPtr> parameters; @@ -3003,23 +3065,32 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard( auto &extensionNode = projCRSNode->lookForChild(WKTConstants::EXTENSION); const auto &extensionChildren = extensionNode->GP()->children(); + bool gdal_3026_hack = false; if (metadata::Identifier::isEquivalentName(wkt1ProjectionName.c_str(), "Mercator_1SP") && - projCRSNode->countChildrenOfName("center_latitude") == 0) { + !projectionHasParameter(projCRSNode, "center_latitude")) { - // The latitude of origin, which should always be zero, is - // missing - // in GDAL WKT1, but provisionned in the EPSG Mercator_1SP - // definition, - // so add it manually. - PropertyMap propertiesParameter; - propertiesParameter.set(IdentifiedObject::NAME_KEY, - "Latitude of natural origin"); - propertiesParameter.set(Identifier::CODE_KEY, 8801); - propertiesParameter.set(Identifier::CODESPACE_KEY, Identifier::EPSG); - parameters.push_back(OperationParameter::create(propertiesParameter)); - values.push_back( - ParameterValue::create(Measure(0, UnitOfMeasure::DEGREE))); + // Hack for https://trac.osgeo.org/gdal/ticket/3026 + if (projectionHasParameter(projCRSNode, "latitude_of_origin")) { + wkt1ProjectionName = "Mercator_2SP"; + gdal_3026_hack = true; + } else { + // The latitude of origin, which should always be zero, is + // missing + // in GDAL WKT1, but provisionned in the EPSG Mercator_1SP + // definition, + // so add it manually. + PropertyMap propertiesParameter; + propertiesParameter.set(IdentifiedObject::NAME_KEY, + "Latitude of natural origin"); + propertiesParameter.set(Identifier::CODE_KEY, 8801); + propertiesParameter.set(Identifier::CODESPACE_KEY, + Identifier::EPSG); + parameters.push_back( + OperationParameter::create(propertiesParameter)); + values.push_back( + ParameterValue::create(Measure(0, UnitOfMeasure::DEGREE))); + } } else if (metadata::Identifier::isEquivalentName( wkt1ProjectionName.c_str(), "Polar_Stereographic")) { @@ -3036,7 +3107,7 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard( defaultLinearUnit, defaultAngularUnit); mapParameters.insert(std::pair<std::string, Measure>( - wkt1ParameterName, Measure(val, unit))); + tolower(wkt1ParameterName), Measure(val, unit))); } catch (const std::exception &) { } } @@ -3051,25 +3122,28 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard( Measure falseEasting = mapParameters["false_easting"]; Measure falseNorthing = mapParameters["false_northing"]; if (latitudeOfOrigin.unit() != UnitOfMeasure::NONE && - std::fabs(std::fabs(latitudeOfOrigin.convertToUnit( - UnitOfMeasure::DEGREE)) - - 90.0) < 1e-10) { - return Conversion::createPolarStereographicVariantA( + scaleFactor.getSIValue() == 1.0) { + return Conversion::createPolarStereographicVariantB( PropertyMap().set(IdentifiedObject::NAME_KEY, "unnamed"), Angle(latitudeOfOrigin.value(), latitudeOfOrigin.unit()), Angle(centralMeridian.value(), centralMeridian.unit()), - Scale(scaleFactor.value(), scaleFactor.unit()), Length(falseEasting.value(), falseEasting.unit()), Length(falseNorthing.value(), falseNorthing.unit())); - } else if (latitudeOfOrigin.unit() != UnitOfMeasure::NONE && - scaleFactor.getSIValue() == 1.0) { - return Conversion::createPolarStereographicVariantB( + } + + if (latitudeOfOrigin.unit() != UnitOfMeasure::NONE && + std::fabs(std::fabs(latitudeOfOrigin.convertToUnit( + UnitOfMeasure::DEGREE)) - + 90.0) < 1e-10) { + return Conversion::createPolarStereographicVariantA( PropertyMap().set(IdentifiedObject::NAME_KEY, "unnamed"), Angle(latitudeOfOrigin.value(), latitudeOfOrigin.unit()), Angle(centralMeridian.value(), centralMeridian.unit()), + Scale(scaleFactor.value(), scaleFactor.unit()), Length(falseEasting.value(), falseEasting.unit()), Length(falseNorthing.value(), falseNorthing.unit())); } + tryToIdentifyWKT1Method = false; // Import GDAL PROJ4 extension nodes } else if (extensionChildren.size() == 2 && @@ -3124,13 +3198,32 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard( if (childNodeChildren.size() < 2) { ThrowNotEnoughChildren(WKTConstants::PARAMETER); } + const auto ¶mValue = childNodeChildren[1]->GP()->value(); + PropertyMap propertiesParameter; const std::string wkt1ParameterName( stripQuotes(childNodeChildren[0])); std::string parameterName(wkt1ParameterName); + if (gdal_3026_hack) { + if (ci_equal(parameterName, "latitude_of_origin")) { + parameterName = "standard_parallel_1"; + } else if (ci_equal(parameterName, "scale_factor") && + paramValue == "1") { + continue; + } + } const auto *paramMapping = mapping ? getMappingFromWKT1(mapping, parameterName) : nullptr; - if (paramMapping) { + if (mapping && + mapping->epsg_code == EPSG_CODE_METHOD_MERCATOR_VARIANT_B && + ci_equal(parameterName, "latitude_of_origin")) { + parameterName = EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN; + propertiesParameter.set( + Identifier::CODE_KEY, + EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN); + propertiesParameter.set(Identifier::CODESPACE_KEY, + Identifier::EPSG); + } else if (paramMapping) { parameterName = paramMapping->wkt2_name; if (paramMapping->epsg_code != 0) { propertiesParameter.set(Identifier::CODE_KEY, @@ -3142,7 +3235,6 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard( propertiesParameter.set(IdentifiedObject::NAME_KEY, parameterName); parameters.push_back( OperationParameter::create(propertiesParameter)); - const auto ¶mValue = childNodeChildren[1]->GP()->value(); try { double val = io::asDouble(paramValue); auto unit = guessUnitForParameter( @@ -3302,6 +3394,8 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) { } } + addExtensionProj4ToProp(nodeP, props); + return ProjectedCRS::create(props, baseGeodCRS, conversion, NN_NO_CHECK(cartesianCS)); } @@ -4226,6 +4320,14 @@ IPROJStringExportable::~IPROJStringExportable() = default; std::string IPROJStringExportable::exportToPROJString( PROJStringFormatter *formatter) const { _exportToPROJString(formatter); + if (formatter->getAddNoDefs() && + formatter->convention() == + io::PROJStringFormatter::Convention::PROJ_4 && + dynamic_cast<const crs::CRS *>(this)) { + if (!formatter->hasParam("no_defs")) { + formatter->addParam("no_defs"); + } + } return formatter->toString(); } //! @endcond @@ -4291,6 +4393,8 @@ struct PROJStringFormatter::Private { bool omitZUnitConversion_ = false; DatabaseContextPtr dbContext_{}; bool useETMercForTMerc_ = false; + bool useETMercForTMercSet_ = false; + bool addNoDefs_ = true; std::string result_{}; @@ -4346,6 +4450,7 @@ PROJStringFormatter::create(Convention conventionIn, * instead of tmerc */ void PROJStringFormatter::setUseETMercForTMerc(bool flag) { d->useETMercForTMerc_ = flag; + d->useETMercForTMercSet_ = true; } // --------------------------------------------------------------------------- @@ -4768,7 +4873,8 @@ PROJStringFormatter::Convention PROJStringFormatter::convention() const { // --------------------------------------------------------------------------- -bool PROJStringFormatter::getUseETMercForTMerc() const { +bool PROJStringFormatter::getUseETMercForTMerc(bool &settingSetOut) const { + settingSetOut = d->useETMercForTMercSet_; return d->useETMercForTMerc_; } @@ -4962,6 +5068,27 @@ void PROJStringFormatter::setCurrentStepInverted(bool inverted) { // --------------------------------------------------------------------------- +bool PROJStringFormatter::hasParam(const char *paramName) const { + if (!d->steps_.empty()) { + for (const auto ¶mValue : d->steps_.back().paramValues) { + if (paramValue.keyEquals(paramName)) { + return true; + } + } + } + return false; +} + +// --------------------------------------------------------------------------- + +void PROJStringFormatter::addNoDefs(bool b) { d->addNoDefs_ = b; } + +// --------------------------------------------------------------------------- + +bool PROJStringFormatter::getAddNoDefs() const { return d->addNoDefs_; } + +// --------------------------------------------------------------------------- + void PROJStringFormatter::addParam(const std::string ¶mName) { if (d->steps_.empty()) { d->addStep(); @@ -5136,6 +5263,8 @@ struct PROJStringParser::Private { DatabaseContextPtr dbContext_{}; std::vector<std::string> warningList_{}; + std::string projString_{}; + std::vector<Step> steps_{}; std::vector<Step::KeyValue> globalParamValues_{}; std::string title_{}; @@ -5173,6 +5302,16 @@ struct PROJStringParser::Private { } // cppcheck-suppress functionStatic + const std::string &getParamValueK(const Step &step) { + for (const auto &pair : step.paramValues) { + if (ci_equal(pair.key, "k") || ci_equal(pair.key, "k_0")) { + return pair.value; + } + } + return emptyString; + } + + // cppcheck-suppress functionStatic std::string guessBodyName(double a); PrimeMeridianNNPtr buildPrimeMeridian(const Step &step); @@ -5383,7 +5522,7 @@ static const struct DatumDesc { // --------------------------------------------------------------------------- -static bool isGeodeticStep(const std::string &name) { +static bool isGeographicStep(const std::string &name) { return name == "longlat" || name == "lonlat" || name == "latlong" || name == "latlon"; } @@ -5483,19 +5622,38 @@ PROJStringParser::Private::buildDatum(const Step &step, const auto &aStr = getParamValue(step, "a"); const auto &bStr = getParamValue(step, "b"); const auto &rfStr = getParamValue(step, "rf"); + const auto &fStr = getParamValue(step, "f"); const auto &RStr = getParamValue(step, "R"); const util::optional<std::string> optionalEmptyString{}; PrimeMeridianNNPtr pm(buildPrimeMeridian(step)); PropertyMap grfMap; + // It is arguable that we allow the prime meridian of a datum defined by + // its name to be overriden, but this is found at least in a regression test + // of GDAL. So let's keep the ellipsoid part of the datum in that case and + // use the specified prime meridian. + const auto overridePmIfNeeded = + [&pm](const GeodeticReferenceFrameNNPtr &grf) { + if (pm->_isEquivalentTo(PrimeMeridian::GREENWICH.get())) { + return grf; + } else { + return GeodeticReferenceFrame::create( + PropertyMap().set(IdentifiedObject::NAME_KEY, + "Unknown based on " + + grf->ellipsoid()->nameStr() + + " ellipsoid"), + grf->ellipsoid(), grf->anchorDefinition(), pm); + } + }; + if (!datumStr.empty()) { if (datumStr == "WGS84") { - return GeodeticReferenceFrame::EPSG_6326; + return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6326); } else if (datumStr == "NAD83") { - return GeodeticReferenceFrame::EPSG_6269; + return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6269); } else if (datumStr == "NAD27") { - return GeodeticReferenceFrame::EPSG_6267; + return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6267); } else { for (const auto &datumDesc : datumDescs) { @@ -5621,6 +5779,29 @@ PROJStringParser::Private::buildDatum(const Step &step, ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm)); } + else if (!aStr.empty() && !fStr.empty()) { + double a; + try { + a = c_locale_stod(aStr); + } catch (const std::invalid_argument &) { + throw ParsingException("Invalid a value"); + } + double f; + try { + f = c_locale_stod(fStr); + } catch (const std::invalid_argument &) { + throw ParsingException("Invalid f value"); + } + auto ellipsoid = Ellipsoid::createFlattenedSphere( + createMapWithUnknownName(), Length(a), + Scale(f != 0.0 ? 1.0 / f : 0.0), guessBodyName(a)) + ->identify(); + return GeodeticReferenceFrame::create( + grfMap.set(IdentifiedObject::NAME_KEY, + title.empty() ? "unknown" : title.c_str()), + ellipsoid, optionalEmptyString, fixupPrimeMeridan(ellipsoid, pm)); + } + else if (!RStr.empty()) { double R; try { @@ -5637,7 +5818,7 @@ PROJStringParser::Private::buildDatum(const Step &step, } if (!aStr.empty() && bStr.empty() && rfStr.empty()) { - throw ParsingException("a found, but b or rf missing"); + throw ParsingException("a found, but b, f or rf missing"); } if (!bStr.empty() && aStr.empty()) { @@ -5648,14 +5829,11 @@ PROJStringParser::Private::buildDatum(const Step &step, throw ParsingException("rf found, but a missing"); } - if (pm->_isEquivalentTo(PrimeMeridian::GREENWICH.get())) { - return GeodeticReferenceFrame::EPSG_6326; - } else { - return GeodeticReferenceFrame::create( - grfMap.set(IdentifiedObject::NAME_KEY, - "Unknown based on WGS84 ellipsoid"), - Ellipsoid::WGS84, optionalEmptyString, pm); + if (!fStr.empty() && aStr.empty()) { + throw ParsingException("f found, but a missing"); } + + return overridePmIfNeeded(GeodeticReferenceFrame::EPSG_6326); } // --------------------------------------------------------------------------- @@ -5830,15 +6008,20 @@ PROJStringParser::Private::buildGeographicCRS(int iStep, int iUnitConvert, bool ignorePROJAxis) { const auto &step = steps_[iStep]; - const auto &title = isGeodeticStep(step.name) ? title_ : emptyString; + const bool l_isGeographicStep = isGeographicStep(step.name); + const auto &title = l_isGeographicStep ? title_ : emptyString; auto datum = buildDatum(step, title); + auto props = PropertyMap().set(IdentifiedObject::NAME_KEY, + title.empty() ? "unknown" : title); + if (l_isGeographicStep && hasParamValue(step, "wktext")) { + props.set("EXTENSION_PROJ4", projString_); + } + return GeographicCRS::create( - PropertyMap().set(IdentifiedObject::NAME_KEY, - title.empty() ? "unknown" : title), - datum, buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, ignoreVUnits, - ignorePROJAxis)); + props, datum, buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, + ignoreVUnits, ignorePROJAxis)); } // --------------------------------------------------------------------------- @@ -5889,11 +6072,14 @@ PROJStringParser::Private::buildGeocentricCRS(int iStep, int iUnitConvert) { } } + auto props = PropertyMap().set(IdentifiedObject::NAME_KEY, + title.empty() ? "unknown" : title); + if (hasParamValue(step, "wktext")) { + props.set("EXTENSION_PROJ4", projString_); + } + auto cs = CartesianCS::createGeocentric(unit); - return GeodeticCRS::create( - PropertyMap().set(IdentifiedObject::NAME_KEY, - title.empty() ? "unknown" : title), - datum, cs); + return GeodeticCRS::create(props, datum, cs); } // --------------------------------------------------------------------------- @@ -5997,30 +6183,29 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( if (!buildPrimeMeridian(step)->longitude()._isEquivalentTo( geogCRS->primeMeridian()->longitude(), util::IComparable::Criterion::EQUIVALENT)) { - throw ParsingException("inconsistant pm values between projectedCRS " + throw ParsingException("inconsistent pm values between projectedCRS " "and its base geographicalCRS"); } auto axisType = AxisType::REGULAR; - if (step.name == "tmerc" && getParamValue(step, "axis") == "wsu") { + if (step.name == "tmerc" && + ((getParamValue(step, "axis") == "wsu" && iAxisSwap < 0) || + (iAxisSwap > 0 && + getParamValue(steps_[iAxisSwap], "order") == "-1,-2"))) { mapping = getMapping(EPSG_CODE_METHOD_TRANSVERSE_MERCATOR_SOUTH_ORIENTATED); } else if (step.name == "etmerc") { mapping = getMapping(EPSG_CODE_METHOD_TRANSVERSE_MERCATOR); - // TODO: we loose the link to the proj etmerc method here. Add some - // property to Conversion to keep it ? } else if (step.name == "lcc") { const auto &lat_0 = getParamValue(step, "lat_0"); const auto &lat_1 = getParamValue(step, "lat_1"); const auto &lat_2 = getParamValue(step, "lat_2"); - const auto &k_0 = getParamValue(step, "k_0"); - const auto &k = getParamValue(step, "k"); + const auto &k = getParamValueK(step); if (lat_2.empty() && !lat_0.empty() && !lat_1.empty() && getAngularValue(lat_0) == getAngularValue(lat_1)) { mapping = getMapping(EPSG_CODE_METHOD_LAMBERT_CONIC_CONFORMAL_1SP); - } else if ((!k_0.empty() && getNumericValue(k_0) != 1.0) || - (!k.empty() && getNumericValue(k) != 1.0)) { + } else if (!k.empty() && getNumericValue(k) != 1.0) { mapping = getMapping( EPSG_CODE_METHOD_LAMBERT_CONIC_CONFORMAL_2SP_MICHIGAN); } else { @@ -6058,15 +6243,17 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( step.paramValues.emplace_back(Step::KeyValue("gamma", "90")); step.paramValues.emplace_back( Step::KeyValue("lonc", getParamValue(step, "lon_0"))); - } else if (step.name == "krovak" && getParamValue(step, "axis") == "swu") { + } else if (step.name == "krovak" && + ((getParamValue(step, "axis") == "swu" && iAxisSwap < 0) || + (iAxisSwap > 0 && + getParamValue(steps_[iAxisSwap], "order") == "-2,-1"))) { mapping = getMapping(EPSG_CODE_METHOD_KROVAK); } else if (step.name == "merc") { if (hasParamValue(step, "a") && hasParamValue(step, "b") && getParamValue(step, "a") == getParamValue(step, "b") && (!hasParamValue(step, "lat_ts") || getAngularValue(getParamValue(step, "lat_ts")) == 0.0) && - (!hasParamValue(step, "k") || - getNumericValue(getParamValue(step, "k")) == 1.0) && + getNumericValue(getParamValueK(step)) == 1.0 && getParamValue(step, "nadgrids") == "@null") { mapping = getMapping( EPSG_CODE_METHOD_POPULAR_VISUALISATION_PSEUDO_MERCATOR); @@ -6091,7 +6278,16 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( } else { axisType = AxisType::SOUTH_POLE; } - if (hasParamValue(step, "lat_ts")) { + const auto &lat_ts = getParamValue(step, "lat_ts"); + const auto &k = getParamValueK(step); + if (!lat_ts.empty() && + std::fabs(getAngularValue(lat_ts) - lat_0) > 1e-10 && + !k.empty() && std::fabs(getNumericValue(k) - 1) > 1e-10) { + throw ParsingException("lat_ts != lat_0 and k != 1 not " + "supported for Polar Stereographic"); + } + if (!lat_ts.empty() && + (k.empty() || std::fabs(getNumericValue(k) - 1) < 1e-10)) { mapping = getMapping(EPSG_CODE_METHOD_POLAR_STEREOGRAPHIC_VARIANT_B); } else { @@ -6161,14 +6357,10 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( const auto *param = mapping->params[i]; std::string proj_name(param->proj_name ? param->proj_name : ""); const std::string *paramValue = - !proj_name.empty() ? &getParamValue(step, proj_name) - : &emptyString; - // k and k_0 may be used indifferently - if (paramValue->empty() && proj_name == "k") { - paramValue = &getParamValue(step, "k_0"); - } else if (paramValue->empty() && proj_name == "k_0") { - paramValue = &getParamValue(step, "k"); - } + (proj_name == "k" || proj_name == "k_0") + ? &getParamValueK(step) + : !proj_name.empty() ? &getParamValue(step, proj_name) + : &emptyString; double value = 0; if (!paramValue->empty()) { bool hasError = false; @@ -6233,6 +6425,10 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( : UnitOfMeasure::NONE))); } + if (step.name == "etmerc") { + methodMap.set("proj_method", "etmerc"); + } + conv = Conversion::create(mapWithUnknownName, methodMap, parameters, values) .as_nullable(); @@ -6296,10 +6492,14 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( auto cs = CartesianCS::create(emptyPropertyMap, axis[0], axis[1]); - CRSNNPtr crs = ProjectedCRS::create( - PropertyMap().set(IdentifiedObject::NAME_KEY, - title.empty() ? "unknown" : title), - geogCRS, NN_NO_CHECK(conv), cs); + auto props = PropertyMap().set(IdentifiedObject::NAME_KEY, + title.empty() ? "unknown" : title); + + if (hasParamValue(step, "wktext")) { + props.set("EXTENSION_PROJ4", projString_); + } + + CRSNNPtr crs = ProjectedCRS::create(props, geogCRS, NN_NO_CHECK(conv), cs); if (!hasParamValue(step, "geoidgrids") && (hasParamValue(step, "vunits") || hasParamValue(step, "vto_meter"))) { @@ -6322,7 +6522,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( static bool isDatumDefiningParam(const std::string ¶m) { return (param == "datum" || param == "ellps" || param == "a" || - param == "b" || param == "rf" || param == "R"); + param == "b" || param == "rf" || param == "f" || param == "R"); } // --------------------------------------------------------------------------- @@ -6585,6 +6785,7 @@ PROJStringParser::createFromPROJString(const std::string &projString) { std::string vunits; std::string vto_meter; + d->projString_ = projString; PROJStringSyntaxParser(projString, d->steps_, d->globalParamValues_, d->title_, vunits, vto_meter); @@ -6613,10 +6814,11 @@ PROJStringParser::createFromPROJString(const std::string &projString) { if ((d->steps_.size() == 1 || (d->steps_.size() == 2 && d->steps_[1].name == "unitconvert")) && isGeocentricStep(d->steps_[0].name)) { - return d->buildGeocentricCRS( - 0, (d->steps_.size() == 2 && d->steps_[1].name == "unitconvert") - ? 1 - : -1); + return d->buildBoundOrCompoundCRSIfNeeded( + 0, d->buildGeocentricCRS(0, (d->steps_.size() == 2 && + d->steps_[1].name == "unitconvert") + ? 1 + : -1)); } int iFirstGeogStep = -1; @@ -6633,7 +6835,7 @@ PROJStringParser::createFromPROJString(const std::string &projString) { bool unexpectedStructure = false; for (int i = 0; i < static_cast<int>(d->steps_.size()); i++) { const auto &stepName = d->steps_[i].name; - if (isGeodeticStep(stepName)) { + if (isGeographicStep(stepName)) { if (iFirstGeogStep < 0) { iFirstGeogStep = i; } else if (iSecondGeogStep < 0) { diff --git a/src/pj_ellps.c b/src/pj_ellps.c index 08b81c3f..8b3b8f0a 100644 --- a/src/pj_ellps.c +++ b/src/pj_ellps.c @@ -11,7 +11,7 @@ pj_ellps[] = { {"SGS85", "a=6378136.0", "rf=298.257", "Soviet Geodetic System 85"}, {"GRS80", "a=6378137.0", "rf=298.257222101", "GRS 1980(IUGG, 1980)"}, {"IAU76", "a=6378140.0", "rf=298.257", "IAU 1976"}, -{"airy", "a=6377563.396", "b=6356256.910", "Airy 1830"}, +{"airy", "a=6377563.396", "rf=299.3249646", "Airy 1830"}, {"APL4.9", "a=6378137.0", "rf=298.25", "Appl. Physics. 1965"}, {"NWL9D", "a=6378145.0", "rf=298.25", "Naval Weapons Lab., 1965"}, {"mod_airy", "a=6377340.189", "b=6356034.446", "Modified Airy"}, |