diff --git a/include/proj/coordinateoperation.hpp b/include/proj/coordinateoperation.hpp index 5e3ecf0b9f..293ba480c7 100644 --- a/include/proj/coordinateoperation.hpp +++ b/include/proj/coordinateoperation.hpp @@ -497,6 +497,7 @@ class PROJ_GCC_DLL OperationMethod : public common::IdentifiedObject { PROJ_INTERNAL OperationMethod(); PROJ_INTERNAL OperationMethod(const OperationMethod &other); INLINED_MAKE_SHARED + friend class Conversion; private: PROJ_OPAQUE_PRIVATE_DATA @@ -1296,7 +1297,7 @@ class PROJ_GCC_DLL Conversion : public SingleOperation { INLINED_MAKE_SHARED PROJ_FRIEND(crs::ProjectedCRS); - PROJ_INTERNAL void addWKTExtensionNode(io::WKTFormatter *formatter) const; + PROJ_INTERNAL bool addWKTExtensionNode(io::WKTFormatter *formatter) const; private: PROJ_OPAQUE_PRIVATE_DATA diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp index 252c90e3da..82b2bd4999 100644 --- a/include/proj/crs.hpp +++ b/include/proj/crs.hpp @@ -472,7 +472,7 @@ class PROJ_GCC_DLL DerivedCRS : virtual public SingleCRS { PROJ_PRIVATE : //! @cond Doxygen_Suppress - const operation::ConversionNNPtr & + PROJ_INTERNAL const operation::ConversionNNPtr & derivingConversionRef() PROJ_CONST_DECL; //! @endcond diff --git a/include/proj/internal/coordinateoperation_constants.hpp b/include/proj/internal/coordinateoperation_constants.hpp index e999b3a64e..d54951938c 100644 --- a/include/proj/internal/coordinateoperation_constants.hpp +++ b/include/proj/internal/coordinateoperation_constants.hpp @@ -88,6 +88,11 @@ static const ParamMapping paramScaleFactor = { EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN, WKT1_SCALE_FACTOR, common::UnitOfMeasure::Type::SCALE, k_0}; +static const ParamMapping paramScaleFactorK = { + EPSG_NAME_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN, + EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN, WKT1_SCALE_FACTOR, + common::UnitOfMeasure::Type::SCALE, k}; + static const ParamMapping paramFalseEasting = { EPSG_NAME_PARAMETER_FALSE_EASTING, EPSG_CODE_PARAMETER_FALSE_EASTING, WKT1_FALSE_EASTING, common::UnitOfMeasure::Type::LINEAR, x_0}; @@ -130,6 +135,10 @@ static const ParamMapping *const paramsNatOriginScale[] = { ¶mLatitudeNatOrigin, ¶mLongitudeNatOrigin, ¶mScaleFactor, ¶mFalseEasting, ¶mFalseNorthing, nullptr}; +static const ParamMapping *const paramsNatOriginScaleK[] = { + ¶mLatitudeNatOrigin, ¶mLongitudeNatOrigin, ¶mScaleFactorK, + ¶mFalseEasting, ¶mFalseNorthing, nullptr}; + static const ParamMapping paramLatFirstPoint = { "Latitude of 1st point", 0, "Latitude_Of_1st_Point", common::UnitOfMeasure::Type::ANGULAR, lat_1}; @@ -243,9 +252,11 @@ static const ParamMapping *const paramsEQDC[] = {¶mLatNatLatCenter, static const ParamMapping *const paramsLonNatOrigin[] = { ¶mLongitudeNatOrigin, ¶mFalseEasting, ¶mFalseNorthing, nullptr}; -static const ParamMapping *const paramsEqc[] = // same as paramsCEA - {¶mLat1stParallelLatTs, ¶mLongitudeNatOrigin, ¶mFalseEasting, - ¶mFalseNorthing, nullptr}; +static const ParamMapping *const paramsEqc[] = { + ¶mLat1stParallelLatTs, + ¶mLatitudeNatOrigin, // extension of EPSG, but used by GDAL / PROJ + ¶mLongitudeNatOrigin, ¶mFalseEasting, + ¶mFalseNorthing, nullptr}; static const ParamMapping paramSatelliteHeight = { "Satellite Height", 0, "satellite_height", @@ -331,7 +342,7 @@ static const ParamMapping paramLonPoint2 = { common::UnitOfMeasure::Type::ANGULAR, lon_2}; static const ParamMapping *const paramsHomTwoPoint[] = { - ¶mLatCentreLatOrigin, + ¶mLatCentreLatCenter, ¶mLatPoint1, ¶mLonPoint1, ¶mLatPoint2, @@ -342,9 +353,8 @@ static const ParamMapping *const paramsHomTwoPoint[] = { nullptr}; static const ParamMapping *const paramsIMWP[] = { - ¶mLongitudeNatOrigin, ¶mLatitude1stStdParallel, - ¶mLatitude2ndStdParallel, ¶mFalseEasting, - ¶mFalseNorthing, nullptr}; + ¶mLongitudeNatOrigin, ¶mLatFirstPoint, ¶mLatSecondPoint, + ¶mFalseEasting, ¶mFalseNorthing, nullptr}; static const ParamMapping paramLonCentreLonCenter = { EPSG_NAME_PARAMETER_LONGITUDE_OF_ORIGIN, @@ -355,7 +365,7 @@ static const ParamMapping paramColatitudeConeAxis = { EPSG_NAME_PARAMETER_COLATITUDE_CONE_AXIS, EPSG_CODE_PARAMETER_COLATITUDE_CONE_AXIS, WKT1_AZIMUTH, common::UnitOfMeasure::Type::ANGULAR, - nullptr}; /* ignored by PROJ currently */ + "alpha"}; /* ignored by PROJ currently */ static const ParamMapping paramLatitudePseudoStdParallel = { EPSG_NAME_PARAMETER_LATITUDE_PSEUDO_STANDARD_PARALLEL, @@ -393,11 +403,6 @@ static const ParamMapping paramLatMerc1SP = { common::UnitOfMeasure::Type::ANGULAR, nullptr}; // always set to zero, not to be exported in PROJ strings -static const ParamMapping paramScaleFactorK = { - EPSG_NAME_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN, - EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN, WKT1_SCALE_FACTOR, - common::UnitOfMeasure::Type::SCALE, k}; - static const ParamMapping *const paramsMerc1SP[] = { ¶mLatMerc1SP, ¶mLongitudeNatOrigin, ¶mScaleFactorK, ¶mFalseEasting, ¶mFalseNorthing, nullptr}; @@ -492,12 +497,12 @@ static const ParamMapping *const paramsLabordeObliqueMercator[] = { static const MethodMapping methodMappings[] = { {EPSG_NAME_METHOD_TRANSVERSE_MERCATOR, EPSG_CODE_METHOD_TRANSVERSE_MERCATOR, - "Transverse_Mercator", "tmerc", nullptr, paramsNatOriginScale}, + "Transverse_Mercator", "tmerc", nullptr, paramsNatOriginScaleK}, {EPSG_NAME_METHOD_TRANSVERSE_MERCATOR_SOUTH_ORIENTATED, EPSG_CODE_METHOD_TRANSVERSE_MERCATOR_SOUTH_ORIENTATED, "Transverse_Mercator_South_Orientated", "tmerc", "axis=wsu", - paramsNatOriginScale}, + paramsNatOriginScaleK}, {PROJ_WKT2_NAME_METHOD_TWO_POINT_EQUIDISTANT, 0, "Two_Point_Equidistant", "tpeqd", nullptr, paramsTPEQD}, @@ -605,7 +610,7 @@ static const MethodMapping methodMappings[] = { {PROJ_WKT2_NAME_METHOD_INTERRUPTED_GOODE_HOMOLOSINE, 0, "Interrupted_Goode_Homolosine", "igh", nullptr, paramsLonNatOrigin}, - // No WKT1 representation fr sweep=x + // No proper WKT1 representation fr sweep=x {PROJ_WKT2_NAME_METHOD_GEOSTATIONARY_SATELLITE_SWEEP_X, 0, nullptr, "geos", "sweep=x", paramsGeos}, diff --git a/include/proj/internal/coordinateoperation_internal.hpp b/include/proj/internal/coordinateoperation_internal.hpp index 83e1e9871e..bafc3a3bf6 100644 --- a/include/proj/internal/coordinateoperation_internal.hpp +++ b/include/proj/internal/coordinateoperation_internal.hpp @@ -70,6 +70,8 @@ const ParamMapping *getMapping(const MethodMapping *mapping, const OperationParameterValue *param); const ParamMapping *getMappingFromWKT1(const MethodMapping *mapping, const std::string &wkt1_name); +bool areEquivalentParameters(const std::string &a, const std::string &b); + // --------------------------------------------------------------------------- struct ESRIParamMapping { diff --git a/include/proj/internal/internal.hpp b/include/proj/internal/internal.hpp index cbf6e2592b..85dd5ac3c1 100644 --- a/include/proj/internal/internal.hpp +++ b/include/proj/internal/internal.hpp @@ -125,6 +125,8 @@ inline bool starts_with(const std::string &str, const char *prefix) noexcept { return std::memcmp(str.c_str(), prefix, prefixSize) == 0; } +bool ci_less(const std::string &a, const std::string &b) noexcept; + bool ci_starts_with(const char *str, const char *prefix) noexcept; bool ci_starts_with(const std::string &str, const std::string &prefix) noexcept; diff --git a/include/proj/io.hpp b/include/proj/io.hpp index 689fede1a4..c649fa9f08 100644 --- a/include/proj/io.hpp +++ b/include/proj/io.hpp @@ -206,7 +206,19 @@ class PROJ_GCC_DLL WKTFormatter { PROJ_DLL WKTFormatter &setMultiLine(bool multiLine) noexcept; PROJ_DLL WKTFormatter &setIndentationWidth(int width) noexcept; - PROJ_DLL WKTFormatter &setOutputAxis(bool outputAxis) noexcept; + + /** Rule for output AXIS nodes */ + enum class OutputAxisRule { + /** Always include AXIS nodes */ + YES, + /** Never include AXIS nodes */ + NO, + /** Includes them only on PROJCS node if it uses Easting/Northing + *ordering. Typically used for WKT1_GDAL */ + WKT1_GDAL_EPSG_STYLE, + }; + + PROJ_DLL WKTFormatter &setOutputAxis(OutputAxisRule outputAxis) noexcept; PROJ_DLL WKTFormatter &setStrict(bool strict) noexcept; PROJ_DLL bool isStrict() const noexcept; @@ -273,7 +285,7 @@ class PROJ_GCC_DLL WKTFormatter { PROJ_INTERNAL bool isInverted() const; #endif - PROJ_INTERNAL bool outputAxis() const; + PROJ_INTERNAL OutputAxisRule outputAxis() const; PROJ_INTERNAL bool outputAxisOrder() const; PROJ_INTERNAL bool primeMeridianOmittedIfGreenwich() const; PROJ_INTERNAL bool ellipsoidUnitOmittedIfMetre() const; @@ -352,7 +364,7 @@ class PROJ_GCC_DLL PROJStringFormatter { startInversion(); PROJ_DLL void stopInversion(); PROJ_INTERNAL bool isInverted() const; - PROJ_INTERNAL bool getUseETMercForTMerc() const; + PROJ_INTERNAL bool getUseETMercForTMerc(bool &settingSetOut) const; PROJ_DLL void ingestPROJString(const std::string &str); // throw ParsingException @@ -373,6 +385,11 @@ class PROJ_GCC_DLL PROJStringFormatter { PROJ_DLL void addParam(const char *paramName, const std::vector &vals); + PROJ_INTERNAL bool hasParam(const char *paramName) const; + + PROJ_INTERNAL void addNoDefs(bool b); + PROJ_INTERNAL bool getAddNoDefs() const; + PROJ_INTERNAL std::set getUsedGridNames() const; PROJ_INTERNAL void setTOWGS84Parameters(const std::vector ¶ms); diff --git a/src/c_api.cpp b/src/c_api.cpp index 1720a7a10f..fbba0f5189 100644 --- a/src/c_api.cpp +++ b/src/c_api.cpp @@ -839,7 +839,10 @@ static const char *getOptionValue(const char *option, *
  • MULTILINE=YES/NO. Defaults to YES, except for WKT1_ESRI
    • *
  • INDENTATION_WIDTH=number. Defauls to 4 (when multiline output is * on).
    • - *
  • OUTPUT_AXIS=YES/NO. Defaults to YES, except for WKT1_ESRI.
    • + *
  • 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.
    • * * @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 33e6f422d7..42cc806e5d 100644 --- a/src/coordinateoperation.cpp +++ b/src/coordinateoperation.cpp @@ -143,7 +143,7 @@ static std::set buildSetEquivalentParameters() { std::set 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 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 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 setEquivalentParameters = buildSetEquivalentParameters(); @@ -690,6 +691,7 @@ struct OperationMethod::Private { util::optional formula_{}; util::optional formulaCitation_{}; std::vector 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( + 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(transformSourceCRS.get()) + ? util::nn_static_pointer_cast( + crs::GeographicCRS::EPSG_4326) + : util::nn_static_pointer_cast( + 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 1c84cf2cc1..f8d2d2b32d 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 dab704b4d4..a204a037f4 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(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(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(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 &TOWGS84Parameters) { + + auto geodCRS = baseCRSIn->extractGeodeticCRS(); + auto targetCRS = + geodCRS.get() == nullptr || + dynamic_cast(geodCRS.get()) + ? util::nn_static_pointer_cast( + crs::GeographicCRS::EPSG_4326) + : util::nn_static_pointer_cast( + 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 aa27192a51..4bec1bf9d1 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. */ diff --git a/src/io.cpp b/src/io.cpp index 834507451b..87caddd07c 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -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(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 mapParamNameToValue; + std::map 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 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( - 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(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 warningList_{}; + std::string projString_{}; + std::vector steps_{}; std::vector globalParamValues_{}; std::string title_{}; @@ -5172,6 +5301,16 @@ struct PROJStringParser::Private { return emptyString; } + // 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); @@ -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 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(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 08b81c3ffc..8b3b8f0a35 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"}, diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index 207a8cd189..8c9f114bd2 100644 --- a/test/unit/test_c_api.cpp +++ b/test/unit/test_c_api.cpp @@ -263,6 +263,22 @@ TEST_F(CApi, proj_obj_as_wkt) { EXPECT_TRUE(std::string(wkt).find("AXIS") == std::string::npos) << wkt; } + // OUTPUT_AXIS=AUTO + { + const char *const options[] = {"OUTPUT_AXIS=AUTO", nullptr}; + auto wkt = proj_obj_as_wkt(obj, PJ_WKT1_GDAL, options); + ASSERT_NE(wkt, nullptr); + EXPECT_TRUE(std::string(wkt).find("AXIS") == std::string::npos) << wkt; + } + + // OUTPUT_AXIS=YES + { + const char *const options[] = {"OUTPUT_AXIS=YES", nullptr}; + auto wkt = proj_obj_as_wkt(obj, PJ_WKT1_GDAL, options); + ASSERT_NE(wkt, nullptr); + EXPECT_TRUE(std::string(wkt).find("AXIS") != std::string::npos) << wkt; + } + // unsupported option { const char *const options[] = {"unsupported=yes", nullptr}; @@ -307,7 +323,7 @@ TEST_F(CApi, proj_obj_as_proj_string) { { auto proj_4 = proj_obj_as_proj_string(obj, PJ_PROJ_4, nullptr); ASSERT_NE(proj_4, nullptr); - EXPECT_EQ(std::string(proj_4), "+proj=longlat +datum=WGS84"); + EXPECT_EQ(std::string(proj_4), "+proj=longlat +datum=WGS84 +no_defs"); } } @@ -327,7 +343,7 @@ TEST_F(CApi, proj_obj_as_proj_string_incompatible_WKT1) { // --------------------------------------------------------------------------- -TEST_F(CApi, proj_obj_as_proj_string_etmerc_option) { +TEST_F(CApi, proj_obj_as_proj_string_etmerc_option_yes) { auto obj = proj_obj_create_from_proj_string(m_ctxt, "+proj=tmerc", nullptr); ObjectKeeper keeper(obj); ASSERT_NE(obj, nullptr); @@ -335,8 +351,24 @@ TEST_F(CApi, proj_obj_as_proj_string_etmerc_option) { const char *options[] = {"USE_ETMERC=YES", nullptr}; auto str = proj_obj_as_proj_string(obj, PJ_PROJ_4, options); ASSERT_NE(str, nullptr); - EXPECT_EQ(str, std::string("+proj=etmerc +lat_0=0 +lon_0=0 +k_0=1 +x_0=0 " - "+y_0=0 +datum=WGS84")); + EXPECT_EQ(str, std::string("+proj=etmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 " + "+y_0=0 +datum=WGS84 +units=m +no_defs")); +} + +// --------------------------------------------------------------------------- + +TEST_F(CApi, proj_obj_as_proj_string_etmerc_option_no) { + auto obj = + proj_obj_create_from_proj_string(m_ctxt, "+proj=utm +zone=31", nullptr); + ObjectKeeper keeper(obj); + ASSERT_NE(obj, nullptr); + + const char *options[] = {"USE_ETMERC=NO", nullptr}; + auto str = proj_obj_as_proj_string(obj, PJ_PROJ_4, options); + ASSERT_NE(str, nullptr); + EXPECT_EQ(str, std::string("+proj=tmerc +lat_0=0 +lon_0=3 +k=0.9996 " + "+x_0=500000 +y_0=0 +datum=WGS84 +units=m " + "+no_defs")); } // --------------------------------------------------------------------------- @@ -356,7 +388,8 @@ TEST_F(CApi, proj_obj_crs_create_bound_crs_to_WGS84) { EXPECT_EQ(std::string(proj_4), "+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 " "+y_0=500000 +ellps=krass " - "+towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69"); + "+towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69 " + "+units=m +no_defs"); } // --------------------------------------------------------------------------- diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp index 436982527c..84be1b22d7 100644 --- a/test/unit/test_crs.cpp +++ b/test/unit/test_crs.cpp @@ -302,22 +302,39 @@ TEST(crs, EPSG_4326_as_WKT2_2018_SIMPLIFIED) { TEST(crs, EPSG_4326_as_WKT1_GDAL) { auto crs = GeographicCRS::EPSG_4326; - WKTFormatterNNPtr f( - WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL)); - crs->exportToWKT(f.get()); - EXPECT_EQ(f->toString(), - "GEOGCS[\"WGS 84\",\n" - " DATUM[\"WGS_1984\",\n" - " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" - " AUTHORITY[\"EPSG\",\"7030\"]],\n" - " AUTHORITY[\"EPSG\",\"6326\"]],\n" - " PRIMEM[\"Greenwich\",0,\n" - " AUTHORITY[\"EPSG\",\"8901\"]],\n" - " UNIT[\"degree\",0.0174532925199433,\n" - " AUTHORITY[\"EPSG\",\"9122\"]],\n" - " AXIS[\"Latitude\",NORTH],\n" - " AXIS[\"Longitude\",EAST],\n" - " AUTHORITY[\"EPSG\",\"4326\"]]"); + auto wkt = crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()); + EXPECT_EQ(wkt, "GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4326\"]]"); +} + +// --------------------------------------------------------------------------- + +TEST(crs, EPSG_4326_as_WKT1_GDAL_with_axis) { + auto crs = GeographicCRS::EPSG_4326; + auto wkt = crs->exportToWKT( + &(WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL) + ->setOutputAxis(WKTFormatter::OutputAxisRule::YES))); + EXPECT_EQ(wkt, "GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AXIS[\"Latitude\",NORTH],\n" + " AXIS[\"Longitude\",EAST],\n" + " AUTHORITY[\"EPSG\",\"4326\"]]"); } // --------------------------------------------------------------------------- @@ -356,7 +373,7 @@ TEST(crs, EPSG_4326_as_PROJ_string) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=longlat +datum=WGS84"); + "+proj=longlat +datum=WGS84 +no_defs"); } // --------------------------------------------------------------------------- @@ -403,28 +420,45 @@ TEST(crs, EPSG_4979_as_WKT2_2018_SIMPLIFIED) { // --------------------------------------------------------------------------- -TEST(crs, EPSG_4979_as_WKT1_GDAL) { +TEST(crs, EPSG_4979_as_WKT1_GDAL_with_axis) { auto crs = GeographicCRS::EPSG_4979; - WKTFormatterNNPtr f( - WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL)); - crs->exportToWKT(f.get()); - // FIXME? WKT1 only supports 2 axis for GEOGCS. So this is an extension of + auto wkt = crs->exportToWKT( + &(WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL) + ->setOutputAxis(WKTFormatter::OutputAxisRule::YES))); + // WKT1 only supports 2 axis for GEOGCS. So this is an extension of // WKT1 as it // and GDAL doesn't really export such as beast, although it can import it - EXPECT_EQ(f->toString(), - "GEOGCS[\"WGS 84\",\n" - " DATUM[\"WGS_1984\",\n" - " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" - " AUTHORITY[\"EPSG\",\"7030\"]],\n" - " AUTHORITY[\"EPSG\",\"6326\"]],\n" - " PRIMEM[\"Greenwich\",0,\n" - " AUTHORITY[\"EPSG\",\"8901\"]],\n" - " UNIT[\"degree\",0.0174532925199433,\n" - " AUTHORITY[\"EPSG\",\"9122\"]],\n" - " AXIS[\"Latitude\",NORTH],\n" - " AXIS[\"Longitude\",EAST],\n" - " AXIS[\"Ellipsoidal height\",UP],\n" - " AUTHORITY[\"EPSG\",\"4979\"]]"); + EXPECT_EQ(wkt, "GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AXIS[\"Latitude\",NORTH],\n" + " AXIS[\"Longitude\",EAST],\n" + " AXIS[\"Ellipsoidal height\",UP],\n" + " AUTHORITY[\"EPSG\",\"4979\"]]"); +} + +// --------------------------------------------------------------------------- + +TEST(crs, EPSG_4979_as_WKT1_GDAL) { + auto crs = GeographicCRS::EPSG_4979; + auto wkt = crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()); + EXPECT_EQ(wkt, "GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4979\"]]"); } // --------------------------------------------------------------------------- @@ -512,8 +546,6 @@ TEST(crs, EPSG_4807_as_WKT1_GDAL) { " AUTHORITY[\"EPSG\",\"8903\"]],\n" " UNIT[\"grad\",0.015707963267949,\n" " AUTHORITY[\"EPSG\",\"9105\"]],\n" - " AXIS[\"Latitude\",NORTH],\n" - " AXIS[\"Longitude\",EAST],\n" " AUTHORITY[\"EPSG\",\"4807\"]]"); } @@ -554,7 +586,7 @@ TEST(crs, EPSG_4807_as_PROJ_string) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=longlat +ellps=clrk80ign +pm=paris"); + "+proj=longlat +ellps=clrk80ign +pm=paris +no_defs"); } // --------------------------------------------------------------------------- @@ -580,7 +612,7 @@ TEST(crs, EPSG_4267) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=longlat +datum=NAD27"); + "+proj=longlat +datum=NAD27 +no_defs"); } // --------------------------------------------------------------------------- @@ -619,7 +651,7 @@ TEST(crs, EPSG_4269) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=longlat +datum=NAD83"); + "+proj=longlat +datum=NAD83 +no_defs"); } // --------------------------------------------------------------------------- @@ -687,7 +719,7 @@ TEST(crs, EPSG_27561_projected_with_geodetic_in_grad_as_PROJ_string_and_WKT1) { PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), "+proj=lcc +lat_1=49.5 +lat_0=49.5 +lon_0=0 +k_0=0.999877341 " - "+x_0=600000 +y_0=200000 +ellps=clrk80ign +pm=paris"); + "+x_0=600000 +y_0=200000 +ellps=clrk80ign +pm=paris +units=m +no_defs"); auto nn_crs = NN_CHECK_ASSERT(crs); EXPECT_TRUE(nn_crs->isEquivalentTo(nn_crs.get())); @@ -704,9 +736,7 @@ TEST(crs, EPSG_27561_projected_with_geodetic_in_grad_as_PROJ_string_and_WKT1) { " DATUM[\"Nouvelle_Triangulation_Francaise_Paris\",\n" " SPHEROID[\"Clarke 1880 (IGN)\",6378249.2,293.4660213]],\n" " PRIMEM[\"Paris\",2.33722917000759],\n" - " UNIT[\"grad\",0.015707963268],\n" - " AXIS[\"Latitude\",NORTH],\n" - " AXIS[\"Longitude\",EAST]],\n" + " UNIT[\"grad\",0.015707963268]],\n" " PROJECTION[\"Lambert_Conformal_Conic_1SP\"],\n" " PARAMETER[\"latitude_of_origin\",55],\n" " PARAMETER[\"central_meridian\",0],\n" @@ -799,15 +829,15 @@ TEST(crs, EPSG_2222_projected_unit_foot_as_PROJ_string_and_WKT1) { EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline +step +proj=axisswap +order=2,1 +step " "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=tmerc " - "+lat_0=31 +lon_0=-110.166666666667 +k_0=0.9999 +x_0=213360 " + "+lat_0=31 +lon_0=-110.166666666667 +k=0.9999 +x_0=213360 " "+y_0=0 +ellps=GRS80 +step +proj=unitconvert +xy_in=m +z_in=m " "+xy_out=ft +z_out=ft"); EXPECT_EQ( crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=tmerc +lat_0=31 +lon_0=-110.166666666667 +k_0=0.9999 " - "+x_0=213360 +y_0=0 +datum=NAD83 +units=ft"); + "+proj=tmerc +lat_0=31 +lon_0=-110.166666666667 +k=0.9999 " + "+x_0=213360 +y_0=0 +datum=NAD83 +units=ft +no_defs"); auto wkt1 = crs->exportToWKT( WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()); @@ -819,9 +849,7 @@ TEST(crs, EPSG_2222_projected_unit_foot_as_PROJ_string_and_WKT1) { " PRIMEM[\"Greenwich\",0,\n" " AUTHORITY[\"EPSG\",\"8901\"]],\n" " UNIT[\"degree\",0.0174532925199433,\n" - " AUTHORITY[\"EPSG\",\"9122\"]],\n" - " AXIS[\"Latitude\",NORTH],\n" - " AXIS[\"Longitude\",EAST]],\n" + " AUTHORITY[\"EPSG\",\"9122\"]]],\n" " PROJECTION[\"Transverse_Mercator\"],\n" " PARAMETER[\"latitude_of_origin\",31],\n" " PARAMETER[\"central_meridian\",-110.166666666667],\n" @@ -867,9 +895,7 @@ TEST(crs, projected_with_parameter_unit_different_than_cs_unit_as_WKT1) { " SPHEROID[\"GRS 1980\",6378137,298.257222101]],\n" " PRIMEM[\"Greenwich\",0],\n" " UNIT[\"degree\",0.0174532925199433,\n" - " AUTHORITY[\"EPSG\",\"9122\"]],\n" - " AXIS[\"Latitude\",NORTH],\n" - " AXIS[\"Longitude\",EAST]],\n" + " AUTHORITY[\"EPSG\",\"9122\"]]],\n" " PROJECTION[\"Transverse_Mercator\"],\n" " PARAMETER[\"latitude_of_origin\",0],\n" " PARAMETER[\"central_meridian\",9],\n" @@ -1081,9 +1107,6 @@ TEST(crs, geocentricCRS_as_WKT1_GDAL) { " AUTHORITY[\"EPSG\",\"8901\"]],\n" " UNIT[\"metre\",1,\n" " AUTHORITY[\"EPSG\",\"9001\"]],\n" - " AXIS[\"Geocentric X\",OTHER],\n" - " AXIS[\"Geocentric Y\",OTHER],\n" - " AXIS[\"Geocentric Z\",NORTH],\n" " AUTHORITY[\"EPSG\",\"4328\"]]"); } @@ -1097,7 +1120,7 @@ TEST(crs, geocentricCRS_as_PROJ_string) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=geocent +ellps=WGS84"); + "+proj=geocent +datum=WGS84 +units=m +no_defs"); } // --------------------------------------------------------------------------- @@ -1597,8 +1620,6 @@ TEST(crs, projectedCRS_as_WKT1_GDAL) { " AUTHORITY[\"EPSG\",\"8901\"]],\n" " UNIT[\"degree\",0.0174532925199433,\n" " AUTHORITY[\"EPSG\",\"9122\"]],\n" - " AXIS[\"Latitude\",NORTH],\n" - " AXIS[\"Longitude\",EAST],\n" " AUTHORITY[\"EPSG\",\"4326\"]],\n" " PROJECTION[\"Transverse_Mercator\"],\n" " PARAMETER[\"latitude_of_origin\",0],\n" @@ -1654,7 +1675,7 @@ TEST(crs, projectedCRS_as_PROJ_string) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=utm +zone=31 +datum=WGS84"); + "+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs"); } // --------------------------------------------------------------------------- @@ -2870,8 +2891,6 @@ TEST(crs, compoundCRS_as_WKT1_GDAL) { " AUTHORITY[\"EPSG\",\"8901\"]],\n" " UNIT[\"degree\",0.0174532925199433,\n" " AUTHORITY[\"EPSG\",\"9122\"]],\n" - " AXIS[\"Latitude\",NORTH],\n" - " AXIS[\"Longitude\",EAST],\n" " AUTHORITY[\"EPSG\",\"4326\"]],\n" " PROJECTION[\"Transverse_Mercator\"],\n" " PARAMETER[\"latitude_of_origin\",0],\n" @@ -2913,7 +2932,7 @@ TEST(crs, compoundCRS_as_PROJ_string) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=utm +zone=31 +datum=WGS84 +vunits=m"); + "+proj=utm +zone=31 +datum=WGS84 +units=m +vunits=m +no_defs"); } // --------------------------------------------------------------------------- @@ -3301,9 +3320,7 @@ TEST(crs, boundCRS_to_WKT1) { " PRIMEM[\"Greenwich\",0,\n" " AUTHORITY[\"EPSG\",\"8901\"]],\n" " UNIT[\"degree\",0.0174532925199433,\n" - " AUTHORITY[\"EPSG\",\"9122\"]],\n" - " AXIS[\"Latitude\",NORTH],\n" - " AXIS[\"Longitude\",EAST]],\n" + " AUTHORITY[\"EPSG\",\"9122\"]]],\n" " PROJECTION[\"Transverse_Mercator\"],\n" " PARAMETER[\"latitude_of_origin\",0],\n" " PARAMETER[\"central_meridian\",3],\n" @@ -3339,7 +3356,7 @@ TEST(crs, boundCRS_geographicCRS_to_PROJ_string) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=longlat +ellps=WGS84 +towgs84=1,2,3,4,5,6,7"); + "+proj=longlat +ellps=WGS84 +towgs84=1,2,3,4,5,6,7 +no_defs"); } // --------------------------------------------------------------------------- @@ -3362,7 +3379,8 @@ TEST(crs, boundCRS_projectedCRS_to_PROJ_string) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=utm +zone=31 +ellps=WGS84 +towgs84=1,2,3,4,5,6,7"); + "+proj=utm +zone=31 +ellps=WGS84 +towgs84=1,2,3,4,5,6,7 +units=m " + "+no_defs"); } // --------------------------------------------------------------------------- @@ -3435,9 +3453,7 @@ TEST(crs, WKT1_DATUM_EXTENSION_to_WKT1_and_PROJ_string) { " SPHEROID[\"intl\",6378388,297],\n" " EXTENSION[\"PROJ4_GRIDS\",\"nzgd2kgrid0005.gsb\"]],\n" " PRIMEM[\"Greenwich\",0],\n" - " UNIT[\"degree\",0.0174532925199433],\n" - " AXIS[\"Longitude\",EAST],\n" - " AXIS[\"Latitude\",NORTH]],\n" + " UNIT[\"degree\",0.0174532925199433]],\n" " PROJECTION[\"New_Zealand_Map_Grid\"],\n" " PARAMETER[\"latitude_of_origin\",-41],\n" " PARAMETER[\"central_meridian\",173],\n" @@ -3461,7 +3477,7 @@ TEST(crs, WKT1_DATUM_EXTENSION_to_WKT1_and_PROJ_string) { PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), "+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 " - "+ellps=intl +nadgrids=nzgd2kgrid0005.gsb +units=m"); + "+ellps=intl +nadgrids=nzgd2kgrid0005.gsb +units=m +no_defs"); } // --------------------------------------------------------------------------- @@ -3473,7 +3489,7 @@ TEST(crs, WKT1_VERT_DATUM_EXTENSION_to_WKT1) { " AUTHORITY[\"EPSG\",\"1027\"]],\n" " UNIT[\"metre\",1,\n" " AUTHORITY[\"EPSG\",\"9001\"]],\n" - " AXIS[\"Up\",UP],\n" + " AXIS[\"Gravity-related height\",UP],\n" " AUTHORITY[\"EPSG\",\"3855\"]]"; auto obj = WKTParser().createFromWKT(wkt); @@ -3561,7 +3577,7 @@ TEST(crs, WKT1_VERT_DATUM_EXTENSION_to_PROJ_string) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+geoidgrids=egm08_25.gtx +vunits=m"); + "+geoidgrids=egm08_25.gtx +vunits=m +no_defs"); } // --------------------------------------------------------------------------- @@ -4197,9 +4213,7 @@ TEST(crs, engineeringCRS_WKT2) { TEST(crs, engineeringCRS_WKT1) { auto expected = "LOCAL_CS[\"Engineering CRS\",\n" - " LOCAL_DATUM[\"Engineering datum\",32767],\n" - " AXIS[\"Easting\",EAST],\n" - " AXIS[\"Northing\",NORTH]]"; + " LOCAL_DATUM[\"Engineering datum\",32767]]"; EXPECT_EQ( createEngineeringCRS()->exportToWKT( WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()), @@ -4522,7 +4536,8 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) { .get()), "+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 " "+y_0=500000 +ellps=krass " - "+towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69"); + "+towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69 " + "+units=m +no_defs"); } { // Pulkovo 42 Poland @@ -4538,7 +4553,8 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) { .get()), "+proj=sterea +lat_0=50.625 +lon_0=21.0833333333333 " "+k=0.9998 +x_0=4637000 +y_0=5647000 +ellps=krass " - "+towgs84=33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84"); + "+towgs84=33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84 " + "+units=m +no_defs"); } { // NTF (Paris) @@ -4553,7 +4569,7 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) { PROJStringFormatter::Convention::PROJ_4) .get()), "+proj=longlat +ellps=clrk80ign +pm=paris " - "+towgs84=-168,-60,320,0,0,0,0"); + "+towgs84=-168,-60,320,0,0,0,0 +no_defs"); } { // NTF (Paris) / Lambert zone II + NGF-IGN69 height @@ -4569,7 +4585,8 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) { .get()), "+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 " "+x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris " - "+towgs84=-168,-60,320,0,0,0,0 +vunits=m"); + "+towgs84=-168,-60,320,0,0,0,0 +units=m " + "+vunits=m +no_defs"); } { auto crs = createVerticalCRS(); @@ -4589,7 +4606,7 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) { .get()), "+proj=stere +lat_0=-90 +lat_ts=-67 +lon_0=140 +x_0=300000 " "+y_0=-2299363.482 +ellps=intl " - "+towgs84=324.912,153.282,172.026,0,0,0,0"); + "+towgs84=324.912,153.282,172.026,0,0,0,0 +units=m +no_defs"); } { auto factoryIGNF = @@ -4604,7 +4621,8 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) { PROJStringFormatter::Convention::PROJ_4) .get()), "+proj=geocent +ellps=intl " - "+towgs84=109.753,-528.133,-362.244,0,0,0,0"); + "+towgs84=109.753,-528.133,-362.244,0,0,0,0 +units=m " + "+no_defs"); } } diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 7d26d82acf..26597f097c 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -453,6 +453,69 @@ static std::string contentWKT2_EPSG_4326( // --------------------------------------------------------------------------- +TEST(wkt_parse, wkt1_geographic_with_PROJ4_extension) { + auto wkt = "GEOGCS[\"WGS 84\",\n" + " DATUM[\"unknown\",\n" + " SPHEROID[\"WGS84\",6378137,298.257223563]],\n" + " PRIMEM[\"Greenwich\",0],\n" + " UNIT[\"degree\",0.0174532925199433],\n" + " EXTENSION[\"PROJ4\",\"+proj=longlat +foo=bar +wktext\"]]"; + auto obj = WKTParser().createFromWKT(wkt); + + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ( + crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()), + wkt); + + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + "+proj=longlat +foo=bar +wktext"); + + EXPECT_TRUE( + crs->exportToWKT(WKTFormatter::create().get()).find("EXTENSION") == + std::string::npos); + + EXPECT_TRUE( + crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_ESRI).get()) + .find("EXTENSION") == std::string::npos); +} + +// --------------------------------------------------------------------------- + +TEST(wkt_parse, wkt1_geocentric_with_PROJ4_extension) { + auto wkt = "GEOCCS[\"WGS 84\",\n" + " DATUM[\"unknown\",\n" + " SPHEROID[\"WGS84\",6378137,298.257223563]],\n" + " PRIMEM[\"Greenwich\",0],\n" + " UNIT[\"Meter\",1],\n" + " EXTENSION[\"PROJ4\",\"+proj=geocent +foo=bar +wktext\"]]"; + auto obj = WKTParser().createFromWKT(wkt); + + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ( + crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()), + wkt); + + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + "+proj=geocent +foo=bar +wktext"); + + EXPECT_TRUE( + crs->exportToWKT(WKTFormatter::create().get()).find("EXTENSION") == + std::string::npos); +} + +// --------------------------------------------------------------------------- + TEST(wkt_parse, wkt2_GEODCRS_EPSG_4326) { auto obj = WKTParser().createFromWKT("GEODCRS" + contentWKT2_EPSG_4326); auto crs = nn_dynamic_pointer_cast(obj); @@ -860,6 +923,49 @@ TEST(wkt_parse, wkt1_projected_no_axis) { // --------------------------------------------------------------------------- +TEST(wkt_parse, wkt1_projected_with_PROJ4_extension) { + auto wkt = "PROJCS[\"unnamed\",\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"unknown\",\n" + " SPHEROID[\"WGS84\",6378137,298.257223563]],\n" + " PRIMEM[\"Greenwich\",0],\n" + " UNIT[\"degree\",0.0174532925199433]],\n" + " PROJECTION[\"Mercator_1SP\"],\n" + " PARAMETER[\"central_meridian\",0],\n" + " PARAMETER[\"scale_factor\",1],\n" + " PARAMETER[\"false_easting\",0],\n" + " PARAMETER[\"false_northing\",0],\n" + " UNIT[\"Meter\",1],\n" + " AXIS[\"Easting\",EAST],\n" + " AXIS[\"Northing\",NORTH],\n" + " EXTENSION[\"PROJ4\",\"+proj=merc +wktext\"]]"; + auto obj = WKTParser().createFromWKT(wkt); + + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ( + crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()), + wkt); + + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + "+proj=merc +wktext"); + + EXPECT_TRUE( + crs->exportToWKT(WKTFormatter::create().get()).find("EXTENSION") == + std::string::npos); + + EXPECT_TRUE( + crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_ESRI).get()) + .find("EXTENSION") == std::string::npos); +} + +// --------------------------------------------------------------------------- + TEST(wkt_parse, wkt1_krovak_south_west) { auto wkt = "PROJCS[\"S-JTSK / Krovak\"," @@ -942,13 +1048,15 @@ TEST(wkt_parse, wkt1_krovak_south_west) { auto expectedPROJString = "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad " "+step +proj=krovak +axis=swu +lat_0=49.5 " - "+lon_0=24.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel"; + "+lon_0=24.8333333333333 +alpha=30.2881397222222 " + "+k=0.9999 +x_0=0 +y_0=0 +ellps=bessel"; EXPECT_EQ(projString, expectedPROJString); obj = PROJStringParser().createFromPROJString(projString); crs = nn_dynamic_pointer_cast(obj); ASSERT_TRUE(crs != nullptr); auto wkt2 = crs->exportToWKT(WKTFormatter::create().get()); + EXPECT_TRUE(wkt2.find("METHOD[\"Krovak\"") != std::string::npos) << wkt2; EXPECT_TRUE( wkt2.find("PARAMETER[\"Latitude of pseudo standard parallel\",78.5,") != std::string::npos) @@ -959,6 +1067,17 @@ TEST(wkt_parse, wkt1_krovak_south_west) { << wkt2; EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()), expectedPROJString); + + obj = PROJStringParser().createFromPROJString( + "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=krovak +lat_0=49.5 " + "+lon_0=24.8333333333333 +alpha=30.2881397222222 " + "+k=0.9999 +x_0=0 +y_0=0 +ellps=bessel " + "+step +proj=axisswap +order=-2,-1"); + crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + wkt2 = crs->exportToWKT(WKTFormatter::create().get()); + EXPECT_TRUE(wkt2.find("METHOD[\"Krovak\"") != std::string::npos) << wkt2; } // --------------------------------------------------------------------------- @@ -1043,7 +1162,174 @@ TEST(wkt_parse, wkt1_krovak_north_oriented) { EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad " "+step +proj=krovak +lat_0=49.5 +lon_0=24.8333333333333 " - "+k=0.9999 +x_0=0 +y_0=0 +ellps=bessel"); + "+alpha=30.2881397222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel"); +} + +// --------------------------------------------------------------------------- + +TEST(wkt_parse, wkt1_polar_stereographic_latitude_of_origin_70) { + auto wkt = "PROJCS[\"unknown\",\n" + " GEOGCS[\"unknown\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]]],\n" + " PROJECTION[\"Polar_Stereographic\"],\n" + " PARAMETER[\"latitude_of_origin\",70],\n" + " PARAMETER[\"central_meridian\",2],\n" + " PARAMETER[\"false_easting\",3],\n" + " PARAMETER[\"false_northing\",4],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]]]"; + + auto obj = WKTParser().createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + + auto projString = crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()); + auto expectedPROJString = "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=2 " + "+x_0=3 +y_0=4 +datum=WGS84 +units=m +no_defs"; + EXPECT_EQ(projString, expectedPROJString); +} + +// --------------------------------------------------------------------------- + +TEST(wkt_parse, wkt1_polar_stereographic_latitude_of_origin_90) { + auto wkt = "PROJCS[\"unknown\",\n" + " GEOGCS[\"unknown\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]]],\n" + " PROJECTION[\"Polar_Stereographic\"],\n" + " PARAMETER[\"latitude_of_origin\",90],\n" + " PARAMETER[\"central_meridian\",2],\n" + " PARAMETER[\"false_easting\",3],\n" + " PARAMETER[\"false_northing\",4],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]]]"; + + auto obj = WKTParser().createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + + auto projString = crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()); + auto expectedPROJString = "+proj=stere +lat_0=90 +lat_ts=90 +lon_0=2 " + "+x_0=3 +y_0=4 +datum=WGS84 +units=m +no_defs"; + EXPECT_EQ(projString, expectedPROJString); +} + +// --------------------------------------------------------------------------- + +TEST(wkt_parse, wkt1_polar_stereographic_latitude_of_origin_90_scale_factor_1) { + auto wkt = "PROJCS[\"unknown\",\n" + " GEOGCS[\"unknown\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]]],\n" + " PROJECTION[\"Polar_Stereographic\"],\n" + " PARAMETER[\"latitude_of_origin\",90],\n" + " PARAMETER[\"central_meridian\",2],\n" + " PARAMETER[\"scale_factor\",1],\n" + " PARAMETER[\"false_easting\",3],\n" + " PARAMETER[\"false_northing\",4],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]]]"; + + auto obj = WKTParser().createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + + auto projString = crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()); + auto expectedPROJString = "+proj=stere +lat_0=90 +lat_ts=90 +lon_0=2 " + "+x_0=3 +y_0=4 +datum=WGS84 +units=m +no_defs"; + EXPECT_EQ(projString, expectedPROJString); +} + +// --------------------------------------------------------------------------- + +TEST(wkt_parse, wkt1_polar_stereographic_scale_factor) { + auto wkt = "PROJCS[\"unknown\",\n" + " GEOGCS[\"unknown\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]]],\n" + " PROJECTION[\"Polar_Stereographic\"],\n" + " PARAMETER[\"latitude_of_origin\",90],\n" + " PARAMETER[\"central_meridian\",2],\n" + " PARAMETER[\"scale_factor\",0.99],\n" + " PARAMETER[\"false_easting\",3],\n" + " PARAMETER[\"false_northing\",4],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]]]"; + + auto obj = WKTParser().createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + + auto projString = crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()); + auto expectedPROJString = "+proj=stere +lat_0=90 +lon_0=2 +k=0.99 +x_0=3 " + "+y_0=4 +datum=WGS84 +units=m +no_defs"; + EXPECT_EQ(projString, expectedPROJString); +} + +// --------------------------------------------------------------------------- + +TEST(wkt_parse, wkt1_Spherical_Cross_Track_Height) { + auto wkt = "PROJCS[\"unknown\",\n" + " GEOGCS[\"unknown\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]]],\n" + " PROJECTION[\"Spherical_Cross_Track_Height\"],\n" + " PARAMETER[\"peg_point_latitude\",1],\n" + " PARAMETER[\"peg_point_longitude\",2],\n" + " PARAMETER[\"peg_point_heading\",3],\n" + " PARAMETER[\"peg_point_height\",4],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]]]"; + + auto obj = WKTParser().createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + + auto projString = crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()); + auto expectedPROJString = "+proj=sch +plat_0=1 +plon_0=2 +phdg_0=3 +h_0=4 " + "+datum=WGS84 +units=m +no_defs"; + EXPECT_EQ(projString, expectedPROJString); } // --------------------------------------------------------------------------- @@ -3508,6 +3794,56 @@ TEST(wkt_parse, esri_projcs) { // --------------------------------------------------------------------------- +TEST(wkt_parse, wkt1_esri_case_insensitive_names) { + auto wkt = "PROJCS[\"WGS_1984_UTM_Zone_31N\",GEOGCS[\"GCS_WGS_1984\"," + "DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0," + "298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\"," + "0.0174532925199433]],PROJECTION[\"transverse_mercator\"]," + "PARAMETER[\"false_easting\",500000.0]," + "PARAMETER[\"false_northing\",0.0]," + "PARAMETER[\"central_meridian\",3.0]," + "PARAMETER[\"scale_factor\",0.9996]," + "PARAMETER[\"latitude_of_origin\",0.0],UNIT[\"Meter\",1.0]]"; + + auto obj = WKTParser().createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + + int zone = 0; + bool north = false; + EXPECT_TRUE(crs->derivingConversion()->isUTM(zone, north)); + EXPECT_EQ(zone, 31); + EXPECT_TRUE(north); +} + +// --------------------------------------------------------------------------- + +TEST(wkt_parse, wkt1_esri_non_expected_param_name) { + // We try to be lax on parameter names. + auto wkt = + "PROJCS[\"WGS_1984_UTM_Zone_31N\",GEOGCS[\"GCS_WGS_1984\"," + "DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0," + "298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\"," + "0.0174532925199433]],PROJECTION[\"transverse_mercator\"]," + "PARAMETER[\"false_easting\",500000.0]," + "PARAMETER[\"false_northing\",0.0]," + "PARAMETER[\"longitude_of_center\",3.0]," // should be Central_Meridian + "PARAMETER[\"scale_factor\",0.9996]," + "PARAMETER[\"latitude_of_origin\",0.0],UNIT[\"Meter\",1.0]]"; + + auto obj = WKTParser().createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + + int zone = 0; + bool north = false; + EXPECT_TRUE(crs->derivingConversion()->isUTM(zone, north)); + EXPECT_EQ(zone, 31); + EXPECT_TRUE(north); +} + +// --------------------------------------------------------------------------- + TEST(wkt_parse, wkt1_esri_krovak_south_west) { auto wkt = "PROJCS[\"S-JTSK_Krovak\",GEOGCS[\"GCS_S_JTSK\"," "DATUM[\"D_S_JTSK\"," @@ -5305,6 +5641,58 @@ TEST(io, projparse_longlat_a_rf) { // --------------------------------------------------------------------------- +TEST(io, projparse_longlat_a_f_zero) { + auto obj = + PROJStringParser().createFromPROJString("+proj=longlat +a=2 +f=0"); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + WKTFormatterNNPtr f(WKTFormatter::create()); + f->simulCurNodeHasId(); + crs->exportToWKT(f.get()); + auto expected = "GEODCRS[\"unknown\",\n" + " DATUM[\"unknown\",\n" + " ELLIPSOID[\"unknown\",2,0,\n" + " LENGTHUNIT[\"metre\",1]]],\n" + " PRIMEM[\"Reference meridian\",0,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " CS[ellipsoidal,2],\n" + " AXIS[\"longitude\",east,\n" + " ORDER[1],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " AXIS[\"latitude\",north,\n" + " ORDER[2],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]]]"; + EXPECT_EQ(f->toString(), expected); +} + +// --------------------------------------------------------------------------- + +TEST(io, projparse_longlat_a_f_non_zero) { + auto obj = + PROJStringParser().createFromPROJString("+proj=longlat +a=2 +f=0.5"); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + WKTFormatterNNPtr f(WKTFormatter::create()); + f->simulCurNodeHasId(); + crs->exportToWKT(f.get()); + auto expected = "GEODCRS[\"unknown\",\n" + " DATUM[\"unknown\",\n" + " ELLIPSOID[\"unknown\",2,2,\n" + " LENGTHUNIT[\"metre\",1]]],\n" + " PRIMEM[\"Reference meridian\",0,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " CS[ellipsoidal,2],\n" + " AXIS[\"longitude\",east,\n" + " ORDER[1],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " AXIS[\"latitude\",north,\n" + " ORDER[2],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]]]"; + EXPECT_EQ(f->toString(), expected); +} + +// --------------------------------------------------------------------------- + TEST(io, projparse_longlat_R) { auto obj = PROJStringParser().createFromPROJString("+proj=longlat +R=2"); auto crs = nn_dynamic_pointer_cast(obj); @@ -5339,7 +5727,7 @@ TEST(io, projparse_longlat_pm_paris) { f->simulCurNodeHasId(); crs->exportToWKT(f.get()); auto expected = "GEODCRS[\"unknown\",\n" - " DATUM[\"Unknown based on WGS84 ellipsoid\",\n" + " DATUM[\"Unknown based on WGS 84 ellipsoid\",\n" " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n" " LENGTHUNIT[\"metre\",1]]],\n" " PRIMEM[\"Paris\",2.5969213,\n" @@ -5393,7 +5781,7 @@ TEST(io, projparse_longlat_pm_numeric) { f->simulCurNodeHasId(); crs->exportToWKT(f.get()); auto expected = "GEODCRS[\"unknown\",\n" - " DATUM[\"Unknown based on WGS84 ellipsoid\",\n" + " DATUM[\"Unknown based on WGS 84 ellipsoid\",\n" " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n" " LENGTHUNIT[\"metre\",1]]],\n" " PRIMEM[\"unknown\",2.5,\n" @@ -5411,6 +5799,21 @@ TEST(io, projparse_longlat_pm_numeric) { // --------------------------------------------------------------------------- +TEST(io, projparse_longlat_pm_overriding_datum) { + // 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. + auto obj = PROJStringParser().createFromPROJString( + "+proj=longlat +datum=WGS84 +pm=ferro"); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ(crs->datum()->nameStr(), "Unknown based on WGS 84 ellipsoid"); + EXPECT_EQ(crs->datum()->primeMeridian()->nameStr(), "Ferro"); +} + +// --------------------------------------------------------------------------- + TEST(io, projparse_longlat_complex) { std::string input = "+proj=pipeline +step +proj=longlat +ellps=clrk80ign " @@ -5455,7 +5858,7 @@ TEST(io, projparse_longlat_towgs84_3_terms) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=longlat +ellps=GRS80 +towgs84=1.2,2,3,0,0,0,0"); + "+proj=longlat +ellps=GRS80 +towgs84=1.2,2,3,0,0,0,0 +no_defs"); } // --------------------------------------------------------------------------- @@ -5490,7 +5893,7 @@ TEST(io, projparse_longlat_towgs84_7_terms) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=longlat +ellps=GRS80 +towgs84=1.2,2,3,4,5,6,7"); + "+proj=longlat +ellps=GRS80 +towgs84=1.2,2,3,4,5,6,7 +no_defs"); } // --------------------------------------------------------------------------- @@ -5514,7 +5917,7 @@ TEST(io, projparse_longlat_nadgrids) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=longlat +ellps=GRS80 +nadgrids=foo.gsb"); + "+proj=longlat +ellps=GRS80 +nadgrids=foo.gsb +no_defs"); } // --------------------------------------------------------------------------- @@ -5542,7 +5945,7 @@ TEST(io, projparse_longlat_geoidgrids) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=longlat +ellps=GRS80 +geoidgrids=foo.gtx +vunits=m"); + "+proj=longlat +ellps=GRS80 +geoidgrids=foo.gtx +vunits=m +no_defs"); } // --------------------------------------------------------------------------- @@ -5839,6 +6242,13 @@ TEST(io, projparse_tmerc_south_oriented) { " LENGTHUNIT[\"metre\",1]]]"; EXPECT_EQ(f->toString(), expected); + + obj = PROJStringParser().createFromPROJString( + "+proj=pipeline +step +proj=tmerc +x_0=1 +lat_0=1 +k_0=2 +step " + "+proj=axisswap +order=-1,-2"); + crs = nn_dynamic_pointer_cast(obj); + EXPECT_EQ(crs->derivingConversion()->method()->nameStr(), + "Transverse Mercator (South Orientated)"); } // --------------------------------------------------------------------------- @@ -6093,14 +6503,26 @@ TEST(io, projparse_etmerc) { auto obj = PROJStringParser().createFromPROJString("+proj=etmerc"); auto crs = nn_dynamic_pointer_cast(obj); ASSERT_TRUE(crs != nullptr); - WKTFormatterNNPtr f(WKTFormatter::create()); - f->simulCurNodeHasId(); - f->setMultiLine(false); - crs->exportToWKT(f.get()); - auto wkt = f->toString(); - EXPECT_TRUE(wkt.find("METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]]") != - std::string::npos) - << wkt; + auto wkt2 = crs->exportToWKT( + &WKTFormatter::create()->simulCurNodeHasId().setMultiLine(false)); + EXPECT_TRUE( + wkt2.find("METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]]") != + std::string::npos) + << wkt2; + + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + "+proj=etmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 " + "+datum=WGS84 +units=m +no_defs"); + + auto wkt1 = crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()); + EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\",\"+proj=etmerc +lat_0=0 " + "+lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m " + "+no_defs\"]") != std::string::npos) + << wkt1; } // --------------------------------------------------------------------------- @@ -6181,6 +6603,39 @@ TEST(io, projparse_merc_stere_polar_variant_A) { // --------------------------------------------------------------------------- +TEST(io, projparse_merc_stere_polar_k_and_lat_ts) { + auto obj = PROJStringParser().createFromPROJString( + "+proj=stere +lat_0=90 +lat_ts=90 +k=1"); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + auto wkt = crs->exportToWKT( + &(WKTFormatter::create()->simulCurNodeHasId().setMultiLine(false))); + EXPECT_TRUE( + wkt.find( + "METHOD[\"Polar Stereographic (variant B)\",ID[\"EPSG\",9829]]") != + std::string::npos) + << wkt; + EXPECT_TRUE(wkt.find("PARAMETER[\"Latitude of standard parallel\",90") != + std::string::npos) + << wkt; + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + "+proj=stere +lat_0=90 +lat_ts=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 " + "+units=m +no_defs"); +} + +// --------------------------------------------------------------------------- + +TEST(io, projparse_merc_stere_polar_k_and_lat_ts_incompatible) { + EXPECT_THROW(PROJStringParser().createFromPROJString( + "+proj=stere +lat_0=90 +lat_ts=70 +k=0.994"), + ParsingException); +} + +// --------------------------------------------------------------------------- + TEST(io, projparse_merc_stere) { auto obj = PROJStringParser().createFromPROJString("+proj=stere +lat_0=30"); auto crs = nn_dynamic_pointer_cast(obj); @@ -6239,7 +6694,8 @@ TEST(io, projparse_utm_south) { // --------------------------------------------------------------------------- TEST(io, projparse_non_earth_ellipsoid) { - std::string input("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=1"); + std::string input("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=1 +units=m " + "+no_defs"); auto obj = PROJStringParser().createFromPROJString(input); auto crs = nn_dynamic_pointer_cast(obj); ASSERT_TRUE(crs != nullptr); @@ -6435,7 +6891,7 @@ TEST(io, projparse_projected_unknown) { "84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[" "\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\"," "\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\"," - "\"9122\"]],AXIS[\"Longitude\",EAST],AXIS[\"Latitude\",NORTH]]," + "\"9122\"]]]," "PROJECTION[\"custom_proj4\"],UNIT[\"metre\",1,AUTHORITY[\"EPSG\"," "\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],EXTENSION[" "\"PROJ4\",\"+proj=mbt_s +datum=WGS84 +unused_flag +lat_0=45 " @@ -6456,7 +6912,7 @@ TEST(io, projparse_projected_unknown) { PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), "+proj=mbt_s +unused_flag +lat_0=45 +lon_0=0 +k=1 +x_0=10 " - "+y_0=0 +datum=WGS84"); + "+y_0=0 +datum=WGS84 +units=m +no_defs"); EXPECT_EQ( crs->exportToPROJString( @@ -6512,6 +6968,39 @@ TEST(io, projparse_geocent) { // --------------------------------------------------------------------------- +TEST(io, projparse_geocent_towgs84) { + auto obj = PROJStringParser().createFromPROJString( + "+proj=geocent +ellps=WGS84 +towgs84=1,2,3"); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + WKTFormatterNNPtr f(WKTFormatter::create()); + f->simulCurNodeHasId(); + f->setMultiLine(false); + crs->exportToWKT(f.get()); + auto wkt = f->toString(); + EXPECT_TRUE( + wkt.find("METHOD[\"Geocentric translations (geocentric domain)") != + std::string::npos) + << wkt; + EXPECT_TRUE(wkt.find("PARAMETER[\"X-axis translation\",1") != + std::string::npos) + << wkt; + EXPECT_TRUE(wkt.find("PARAMETER[\"Y-axis translation\",2") != + std::string::npos) + << wkt; + EXPECT_TRUE(wkt.find("PARAMETER[\"Z-axis translation\",3") != + std::string::npos) + << wkt; + + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + "+proj=geocent +ellps=WGS84 +towgs84=1,2,3,0,0,0,0 +units=m +no_defs"); +} + +// --------------------------------------------------------------------------- + TEST(io, projparse_cart_unit) { std::string input( "+proj=pipeline +step +proj=cart +ellps=WGS84 +step " @@ -6544,6 +7033,48 @@ TEST(io, projparse_cart_unit_numeric) { // --------------------------------------------------------------------------- +TEST(io, projparse_longlat_wktext) { + std::string input("+proj=longlat +foo +wktext"); + auto obj = PROJStringParser().createFromPROJString(input); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + input); +} + +// --------------------------------------------------------------------------- + +TEST(io, projparse_geocent_wktext) { + std::string input("+proj=geocent +foo +wktext"); + auto obj = PROJStringParser().createFromPROJString(input); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + input); +} + +// --------------------------------------------------------------------------- + +TEST(io, projparse_projected_wktext) { + std::string input("+proj=merc +foo +wktext"); + auto obj = PROJStringParser().createFromPROJString(input); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + input); +} + +// --------------------------------------------------------------------------- + TEST(io, projparse_ob_tran_longlat) { std::string input( "+proj=pipeline +step +proj=axisswap +order=2,1 +step " @@ -6872,7 +7403,8 @@ TEST(io, projparse_longlat_title) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=longlat +ellps=intl +towgs84=109.753,-528.133,-362.244,0,0,0,0"); + "+proj=longlat +ellps=intl +towgs84=109.753,-528.133,-362.244,0,0,0,0 " + "+no_defs"); } // --------------------------------------------------------------------------- @@ -6895,7 +7427,7 @@ TEST(io, projparse_projected_title) { PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), "+proj=utm +zone=43 +south +ellps=intl " - "+towgs84=109.753,-528.133,-362.244,0,0,0,0"); + "+towgs84=109.753,-528.133,-362.244,0,0,0,0 +units=m +no_defs"); } // --------------------------------------------------------------------------- @@ -6992,6 +7524,10 @@ TEST(io, projparse_longlat_errors) { "+proj=longlat +a=1 +rf=invalid"), ParsingException); + EXPECT_THROW(PROJStringParser().createFromPROJString( + "+proj=longlat +a=1 +f=invalid"), + ParsingException); + EXPECT_THROW( PROJStringParser().createFromPROJString("+proj=longlat +R=invalid"), ParsingException); @@ -7005,6 +7541,9 @@ TEST(io, projparse_longlat_errors) { EXPECT_THROW(PROJStringParser().createFromPROJString("+proj=longlat +rf=1"), ParsingException); + EXPECT_THROW(PROJStringParser().createFromPROJString("+proj=longlat +f=0"), + ParsingException); + EXPECT_THROW( PROJStringParser().createFromPROJString("+proj=longlat +pm=unknown"), ParsingException); diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 2d2688a8a5..e62eddf907 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -1185,13 +1185,13 @@ TEST(operation, tmerc_export) { auto conv = Conversion::createTransverseMercator( PropertyMap(), Angle(1), Angle(2), Scale(3), Length(4), Length(5)); EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=tmerc +lat_0=1 +lon_0=2 +k_0=3 +x_0=4 +y_0=5"); + "+proj=tmerc +lat_0=1 +lon_0=2 +k=3 +x_0=4 +y_0=5"); { auto formatter = PROJStringFormatter::create(); formatter->setUseETMercForTMerc(true); EXPECT_EQ(conv->exportToPROJString(formatter.get()), - "+proj=etmerc +lat_0=1 +lon_0=2 +k_0=3 +x_0=4 +y_0=5"); + "+proj=etmerc +lat_0=1 +lon_0=2 +k=3 +x_0=4 +y_0=5"); } EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()), @@ -1270,7 +1270,7 @@ TEST(operation, tmerc_south_oriented_export) { PropertyMap(), Angle(1), Angle(2), Scale(3), Length(4), Length(5)); EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=tmerc +axis=wsu +lat_0=1 +lon_0=2 +k_0=3 +x_0=4 +y_0=5"); + "+proj=tmerc +axis=wsu +lat_0=1 +lon_0=2 +k=3 +x_0=4 +y_0=5"); EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()), "CONVERSION[\"Transverse Mercator (South Orientated)\",\n" @@ -1329,7 +1329,7 @@ TEST(operation, tmerc_south_oriented_export) { EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline +step +proj=axisswap +order=2,1 +step " "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=tmerc " - "+axis=wsu +lat_0=0 +lon_0=29 +k_0=1 +x_0=0 +y_0=0 +ellps=WGS84"); + "+axis=wsu +lat_0=0 +lon_0=29 +k=1 +x_0=0 +y_0=0 +ellps=WGS84"); } // --------------------------------------------------------------------------- @@ -1986,7 +1986,7 @@ TEST(operation, createEquidistantCylindrical) { PropertyMap(), Angle(1), Angle(2), Length(3), Length(4)); EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=eqc +lat_ts=1 +lon_0=2 +x_0=3 +y_0=4"); + "+proj=eqc +lat_ts=1 +lat_0=0 +lon_0=2 +x_0=3 +y_0=4"); EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()), "CONVERSION[\"Equidistant Cylindrical\",\n" @@ -2022,7 +2022,7 @@ TEST(operation, createEquidistantCylindricalSpherical) { PropertyMap(), Angle(1), Angle(2), Length(3), Length(4)); EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=eqc +lat_ts=1 +lon_0=2 +x_0=3 +y_0=4"); + "+proj=eqc +lat_ts=1 +lat_0=0 +lon_0=2 +x_0=3 +y_0=4"); EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()), "CONVERSION[\"Equidistant Cylindrical (Spherical)\",\n" @@ -2053,6 +2053,28 @@ TEST(operation, createEquidistantCylindricalSpherical) { // --------------------------------------------------------------------------- +TEST(operation, equidistant_cylindrical_lat_0) { + + auto obj = PROJStringParser().createFromPROJString( + "+proj=eqc +ellps=sphere +lat_0=-2 +lat_ts=1 +lon_0=-10"); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + + auto wkt = crs->exportToWKT(WKTFormatter::create().get()); + EXPECT_TRUE(wkt.find("PARAMETER[\"Latitude of natural origin\",-2") != + std::string::npos) + << wkt; + + auto projString = crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()); + EXPECT_EQ(projString, + "+proj=eqc +lat_ts=1 +lat_0=-2 +lon_0=-10 +x_0=0 +y_0=0 " + "+ellps=sphere +units=m +no_defs"); +} + +// --------------------------------------------------------------------------- + TEST(operation, gall_export) { auto conv = @@ -2173,10 +2195,18 @@ TEST(operation, geostationary_satellite_sweep_x_export) { " LENGTHUNIT[\"metre\",1],\n" " ID[\"EPSG\",8807]]]"); - EXPECT_THROW( - conv->exportToWKT( - WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()), - FormattingException); + auto crs = ProjectedCRS::create( + PropertyMap(), GeographicCRS::EPSG_4326, conv, + CartesianCS::createEastingNorthing(UnitOfMeasure::METRE)); + auto wkt1 = crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()); + EXPECT_TRUE(wkt1.find("PROJECTION[\"Geostationary_Satellite\"]") != + std::string::npos) + << wkt1; + EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\",\"+proj=geos +sweep=x +lon_0=1 " + "+h=2 +x_0=3 +y_0=4 +datum=WGS84 +units=m " + "+no_defs\"]]") != std::string::npos) + << wkt1; } // --------------------------------------------------------------------------- @@ -2410,7 +2440,7 @@ TEST(operation, hotine_oblique_mercator_two_point_natural_origin_export) { conv->exportToWKT( WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()), "PROJECTION[\"Hotine_Oblique_Mercator_Two_Point_Natural_Origin\"],\n" - "PARAMETER[\"latitude_of_origin\",1],\n" + "PARAMETER[\"latitude_of_center\",1],\n" "PARAMETER[\"latitude_of_point_1\",2],\n" "PARAMETER[\"longitude_of_point_1\",3],\n" "PARAMETER[\"latitude_of_point_2\",4],\n" @@ -2435,12 +2465,12 @@ TEST(operation, imw_polyconic_export) { " PARAMETER[\"Longitude of natural origin\",1,\n" " ANGLEUNIT[\"degree\",0.0174532925199433],\n" " ID[\"EPSG\",8802]],\n" - " PARAMETER[\"Latitude of 1st standard parallel\",3,\n" - " ANGLEUNIT[\"degree\",0.0174532925199433],\n" - " ID[\"EPSG\",8823]],\n" - " PARAMETER[\"Latitude of 2nd standard parallel\",4,\n" - " ANGLEUNIT[\"degree\",0.0174532925199433],\n" - " ID[\"EPSG\",8824]],\n" + " PARAMETER[\"Latitude of 1st point\",3,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433,\n" + " ID[\"EPSG\",9122]]],\n" + " PARAMETER[\"Latitude of 2nd point\",4,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433,\n" + " ID[\"EPSG\",9122]]],\n" " PARAMETER[\"False easting\",5,\n" " LENGTHUNIT[\"metre\",1],\n" " ID[\"EPSG\",8806]],\n" @@ -2453,8 +2483,8 @@ TEST(operation, imw_polyconic_export) { WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()), "PROJECTION[\"International_Map_of_the_World_Polyconic\"],\n" "PARAMETER[\"central_meridian\",1],\n" - "PARAMETER[\"standard_parallel_1\",3],\n" - "PARAMETER[\"standard_parallel_2\",4],\n" + "PARAMETER[\"Latitude_Of_1st_Point\",3],\n" + "PARAMETER[\"Latitude_Of_2nd_Point\",4],\n" "PARAMETER[\"false_easting\",5],\n" "PARAMETER[\"false_northing\",6]"); } @@ -2467,7 +2497,8 @@ TEST(operation, krovak_north_oriented_export) { Angle(78.5), Scale(0.9999), Length(5), Length(6)); EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=krovak +lat_0=49.5 +lon_0=42.5 +k=0.9999 +x_0=5 +y_0=6"); + "+proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.2881397222222 " + "+k=0.9999 +x_0=5 +y_0=6"); EXPECT_EQ( conv->exportToWKT(WKTFormatter::create().get()), @@ -2517,7 +2548,8 @@ TEST(operation, krovak_export) { Angle(78.5), Scale(0.9999), Length(5), Length(6)); EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=krovak +axis=swu +lat_0=49.5 +lon_0=42.5 +k=0.9999 +x_0=5 " + "+proj=krovak +axis=swu +lat_0=49.5 +lon_0=42.5 " + "+alpha=30.2881397222222 +k=0.9999 +x_0=5 " "+y_0=6"); EXPECT_EQ( @@ -2706,6 +2738,36 @@ TEST(operation, wkt1_import_mercator_variant_A) { // --------------------------------------------------------------------------- +TEST(operation, wkt1_import_mercator_variant_A_that_is_variant_B) { + // Adresses https://trac.osgeo.org/gdal/ticket/3026 + auto wkt = "PROJCS[\"test\",\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS 1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563]],\n" + " PRIMEM[\"Greenwich\",0],\n" + " UNIT[\"degree\",0.0174532925199433]],\n" + " PROJECTION[\"Mercator_1SP\"],\n" + " PARAMETER[\"latitude_of_origin\",-1],\n" + " PARAMETER[\"central_meridian\",2],\n" + " PARAMETER[\"scale_factor\",1],\n" + " PARAMETER[\"false_easting\",3],\n" + " PARAMETER[\"false_northing\",4],\n" + " UNIT[\"metre\",1]]"; + auto obj = WKTParser().createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + + auto conversion = crs->derivingConversion(); + auto convRef = Conversion::createMercatorVariantB( + PropertyMap().set(IdentifiedObject::NAME_KEY, "unnamed"), Angle(-1), + Angle(2), Length(3), Length(4)); + + EXPECT_TRUE(conversion->isEquivalentTo(convRef.get(), + IComparable::Criterion::EQUIVALENT)); +} + +// --------------------------------------------------------------------------- + TEST(operation, mercator_variant_B_export) { auto conv = Conversion::createMercatorVariantB( PropertyMap(), Angle(1), Angle(2), Length(3), Length(4)); @@ -2742,6 +2804,70 @@ TEST(operation, mercator_variant_B_export) { // --------------------------------------------------------------------------- +TEST(operation, odd_mercator_1sp_with_non_null_latitude) { + auto obj = WKTParser().createFromWKT( + "PROJCS[\"unnamed\",\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4326\"]],\n" + " PROJECTION[\"Mercator_1SP\"],\n" + " PARAMETER[\"latitude_of_origin\",30],\n" + " PARAMETER[\"central_meridian\",0],\n" + " PARAMETER[\"scale_factor\",0.99],\n" + " PARAMETER[\"false_easting\",0],\n" + " PARAMETER[\"false_northing\",0],\n" + " UNIT[\"metre\",1],\n" + " AXIS[\"Easting\",EAST],\n" + " AXIS[\"Northing\",NORTH]]"); + + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + + EXPECT_THROW(crs->exportToPROJString(PROJStringFormatter::create().get()), + FormattingException); +} + +// --------------------------------------------------------------------------- + +TEST(operation, odd_mercator_2sp_with_latitude_of_origin) { + auto obj = WKTParser().createFromWKT( + "PROJCS[\"unnamed\",\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4326\"]],\n" + " PROJECTION[\"Mercator_2SP\"],\n" + " PARAMETER[\"standard_parallel_1\",30],\n" + " PARAMETER[\"latitude_of_origin\",40],\n" + " PARAMETER[\"central_meridian\",0],\n" + " PARAMETER[\"false_easting\",0],\n" + " PARAMETER[\"false_northing\",0],\n" + " UNIT[\"metre\",1],\n" + " AXIS[\"Easting\",EAST],\n" + " AXIS[\"Northing\",NORTH]]"); + + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + + EXPECT_THROW(crs->exportToPROJString(PROJStringFormatter::create().get()), + FormattingException); +} + +// --------------------------------------------------------------------------- + TEST(operation, webmerc_export) { auto conv = Conversion::createPopularVisualisationPseudoMercator( PropertyMap(), Angle(0), Angle(2), Length(3), Length(4)); @@ -2790,8 +2916,6 @@ TEST(operation, webmerc_export) { " AUTHORITY[\"EPSG\",\"8901\"]],\n" " UNIT[\"degree\",0.0174532925199433,\n" " AUTHORITY[\"EPSG\",\"9122\"]],\n" - " AXIS[\"Latitude\",NORTH],\n" - " AXIS[\"Longitude\",EAST],\n" " AUTHORITY[\"EPSG\",\"4326\"]],\n" " PROJECTION[\"Mercator_1SP\"],\n" " PARAMETER[\"central_meridian\",2],\n" @@ -2994,9 +3118,7 @@ TEST(operation, webmerc_import_from_WKT2_EPSG_3785_deprecated) { " SPHEROID[\"Popular Visualisation Sphere\",6378137,0],\n" " TOWGS84[0,0,0,0,0,0,0]],\n" " PRIMEM[\"Greenwich\",0],\n" - " UNIT[\"degree\",0.0174532925199433],\n" - " AXIS[\"Latitude\",NORTH],\n" - " AXIS[\"Longitude\",EAST]],\n" + " UNIT[\"degree\",0.0174532925199433]],\n" " PROJECTION[\"Mercator_1SP\"],\n" " PARAMETER[\"central_meridian\",0],\n" " PARAMETER[\"scale_factor\",1],\n" @@ -3730,7 +3852,7 @@ TEST(operation, conversion_inverse) { " ID[\"EPSG\",8807]]]"); EXPECT_EQ(inv->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=pipeline +step +inv +proj=tmerc +lat_0=1 +lon_0=2 +k_0=3 " + "+proj=pipeline +step +inv +proj=tmerc +lat_0=1 +lon_0=2 +k=3 " "+x_0=4 +y_0=5"); EXPECT_TRUE(inv->isEquivalentTo(inv.get())); @@ -5463,7 +5585,7 @@ TEST(operation, compoundCRS_to_compoundCRS_with_vertical_transform) { compound2); ASSERT_TRUE(op != nullptr); EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=pipeline +step +inv +proj=tmerc +lat_0=1 +lon_0=2 +k_0=3 " + "+proj=pipeline +step +inv +proj=tmerc +lat_0=1 +lon_0=2 +k=3 " "+x_0=4 +y_0=5 +ellps=WGS84 +step " "+proj=vgridshift +grids=bla.gtx +multiplier=0.001 +step " "+proj=utm +zone=32 " @@ -5473,7 +5595,7 @@ TEST(operation, compoundCRS_to_compoundCRS_with_vertical_transform) { formatter->setUseETMercForTMerc(true); EXPECT_EQ(op->exportToPROJString(formatter.get()), "+proj=pipeline +step +inv +proj=etmerc +lat_0=1 +lon_0=2 " - "+k_0=3 +x_0=4 +y_0=5 +ellps=WGS84 +step " + "+k=3 +x_0=4 +y_0=5 +ellps=WGS84 +step " "+proj=vgridshift +grids=bla.gtx +multiplier=0.001 +step " "+proj=utm +zone=32 " "+ellps=WGS84"); @@ -5485,7 +5607,7 @@ TEST(operation, compoundCRS_to_compoundCRS_with_vertical_transform) { "+proj=pipeline +step +inv +proj=utm +zone=32 +ellps=WGS84 " "+step +inv +proj=vgridshift +grids=bla.gtx " "+multiplier=0.001 +step +proj=etmerc +lat_0=1 +lon_0=2 " - "+k_0=3 +x_0=4 +y_0=5 +ellps=WGS84"); + "+k=3 +x_0=4 +y_0=5 +ellps=WGS84"); } auto opInverse = CoordinateOperationFactory::create()->createOperation(