From 5527b10ed140e20fac8e183317514fd59e4c8b99 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 7 Sep 2021 19:57:36 +0200 Subject: [PATCH 1/4] Support importing/exporting WKT & PROJJSON of 2D axis spherical planetocentric geodetic CRS --- include/proj/coordinatesystem.hpp | 5 ++ include/proj/crs.hpp | 2 + scripts/reference_exported_symbols.txt | 2 + src/iso19111/coordinatesystem.cpp | 21 +++++++ src/iso19111/crs.cpp | 27 ++++++++- src/iso19111/io.cpp | 33 ++++++---- test/unit/test_io.cpp | 84 ++++++++++++++++++++++++++ 7 files changed, 161 insertions(+), 13 deletions(-) diff --git a/include/proj/coordinatesystem.hpp b/include/proj/coordinatesystem.hpp index e16501681a..b40b038d5c 100644 --- a/include/proj/coordinatesystem.hpp +++ b/include/proj/coordinatesystem.hpp @@ -314,6 +314,11 @@ class PROJ_GCC_DLL SphericalCS final : public CoordinateSystem { const CoordinateSystemAxisNNPtr &axis2, const CoordinateSystemAxisNNPtr &axis3); + PROJ_DLL static SphericalCSNNPtr + create(const util::PropertyMap &properties, + const CoordinateSystemAxisNNPtr &axis1, + const CoordinateSystemAxisNNPtr &axis2); + protected: PROJ_INTERNAL explicit SphericalCS( const std::vector &axisIn); diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp index dcab094af2..7fde88c876 100644 --- a/include/proj/crs.hpp +++ b/include/proj/crs.hpp @@ -270,6 +270,8 @@ class PROJ_GCC_DLL GeodeticCRS : virtual public SingleCRS, PROJ_DLL bool isGeocentric() PROJ_PURE_DECL; + PROJ_DLL bool isSphericalPlanetocentric() PROJ_PURE_DECL; + PROJ_DLL static GeodeticCRSNNPtr create(const util::PropertyMap &properties, const datum::GeodeticReferenceFrameNNPtr &datum, diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt index d49a59b220..6cb734e270 100644 --- a/scripts/reference_exported_symbols.txt +++ b/scripts/reference_exported_symbols.txt @@ -152,6 +152,7 @@ osgeo::proj::crs::GeodeticCRS::ellipsoid() const osgeo::proj::crs::GeodeticCRS::~GeodeticCRS() osgeo::proj::crs::GeodeticCRS::identify(std::shared_ptr const&) const osgeo::proj::crs::GeodeticCRS::isGeocentric() const +osgeo::proj::crs::GeodeticCRS::isSphericalPlanetocentric() const osgeo::proj::crs::GeodeticCRS::primeMeridian() const osgeo::proj::crs::GeodeticCRS::velocityModel() const osgeo::proj::crs::GeographicCRS::coordinateSystem() const @@ -225,6 +226,7 @@ osgeo::proj::cs::OrdinalCS::create(osgeo::proj::util::PropertyMap const&, std::v osgeo::proj::cs::OrdinalCS::~OrdinalCS() osgeo::proj::cs::ParametricCS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn > const&) osgeo::proj::cs::ParametricCS::~ParametricCS() +osgeo::proj::cs::SphericalCS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn > const&, dropbox::oxygen::nn > const&) osgeo::proj::cs::SphericalCS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn > const&, dropbox::oxygen::nn > const&, dropbox::oxygen::nn > const&) osgeo::proj::cs::SphericalCS::~SphericalCS() osgeo::proj::cs::TemporalCountCS::create(osgeo::proj::util::PropertyMap const&, dropbox::oxygen::nn > const&) diff --git a/src/iso19111/coordinatesystem.cpp b/src/iso19111/coordinatesystem.cpp index 498e303531..f9db540674 100644 --- a/src/iso19111/coordinatesystem.cpp +++ b/src/iso19111/coordinatesystem.cpp @@ -667,6 +667,27 @@ SphericalCSNNPtr SphericalCS::create(const util::PropertyMap &properties, // --------------------------------------------------------------------------- +/** \brief Instantiate a SphericalCS with 2 axis. + * + * This is an extension to ISO19111 to support (planet)-ocentric CS with + * geocentric latitude. + * + * @param properties See \ref general_properties. + * @param axis1 The first axis. + * @param axis2 The second axis. + * @return a new SphericalCS. + */ +SphericalCSNNPtr SphericalCS::create(const util::PropertyMap &properties, + const CoordinateSystemAxisNNPtr &axis1, + const CoordinateSystemAxisNNPtr &axis2) { + std::vector axis{axis1, axis2}; + auto cs(SphericalCS::nn_make_shared(axis)); + cs->setProperties(properties); + return cs; +} + +// --------------------------------------------------------------------------- + //! @cond Doxygen_Suppress EllipsoidalCS::~EllipsoidalCS() = default; //! @endcond diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 7c8fcd8153..6f5110437a 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -1608,7 +1608,7 @@ GeodeticCRS::velocityModel() PROJ_PURE_DEFN { // --------------------------------------------------------------------------- -/** \brief Return whether the CRS is a geocentric one. +/** \brief Return whether the CRS is a Cartesian geocentric one. * * A geocentric CRS is a geodetic CRS that has a Cartesian coordinate system * with three axis, whose direction is respectively @@ -1629,6 +1629,31 @@ bool GeodeticCRS::isGeocentric() PROJ_PURE_DEFN { // --------------------------------------------------------------------------- +/** \brief Return whether the CRS is a Spherical planetocentric one. + * + * A Spherical planetocentric CRS is a geodetic CRS that has a spherical + * (angular) coordinate system with 2 axis, which represent geocentric latitude/ + * longitude or longitude/geocentric latitude. + * + * Such CRS are typically used in use case that apply to non-Earth bodies. + * + * @return true if the CRS is a Spherical planetocentric CRS. + * + * @since 8.2 + */ +bool GeodeticCRS::isSphericalPlanetocentric() PROJ_PURE_DEFN { + const auto &cs = coordinateSystem(); + const auto &axisList = cs->axisList(); + return axisList.size() == 2 && + dynamic_cast(cs.get()) != nullptr && + ((ci_equal(axisList[0]->nameStr(), "planetocentric latitude") && + ci_equal(axisList[1]->nameStr(), "planetocentric longitude")) || + (ci_equal(axisList[0]->nameStr(), "planetocentric longitude") && + ci_equal(axisList[1]->nameStr(), "planetocentric latitude"))); +} + +// --------------------------------------------------------------------------- + /** \brief Instantiate a GeodeticCRS from a datum::GeodeticReferenceFrame and a * cs::SphericalCS. * diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index a1b06ae740..1373c58304 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -2798,7 +2798,11 @@ WKTParser::Private::buildCS(const WKTNodeNNPtr &node, /* maybe null */ return VerticalCS::create(csMap, axisList[0]); } } else if (ci_equal(csType, "spherical")) { - if (axisCount == 3) { + if (axisCount == 2) { + // Extension to ISO19111 to support (planet)-ocentric CS with + // geocentric latitude + return SphericalCS::create(csMap, axisList[0], axisList[1]); + } else if (axisCount == 3) { return SphericalCS::create(csMap, axisList[0], axisList[1], axisList[2]); } @@ -5945,62 +5949,67 @@ CoordinateSystemNNPtr JSONParser::buildCS(const json &j) { axisList.emplace_back(buildAxis(axis)); } const PropertyMap &csMap = emptyPropertyMap; + const auto axisCount = axisList.size(); if (subtype == "ellipsoidal") { - if (axisList.size() == 2) { + if (axisCount == 2) { return EllipsoidalCS::create(csMap, axisList[0], axisList[1]); } - if (axisList.size() == 3) { + if (axisCount == 3) { return EllipsoidalCS::create(csMap, axisList[0], axisList[1], axisList[2]); } throw ParsingException("Expected 2 or 3 axis"); } if (subtype == "Cartesian") { - if (axisList.size() == 2) { + if (axisCount == 2) { return CartesianCS::create(csMap, axisList[0], axisList[1]); } - if (axisList.size() == 3) { + if (axisCount == 3) { return CartesianCS::create(csMap, axisList[0], axisList[1], axisList[2]); } throw ParsingException("Expected 2 or 3 axis"); } if (subtype == "vertical") { - if (axisList.size() == 1) { + if (axisCount == 1) { return VerticalCS::create(csMap, axisList[0]); } throw ParsingException("Expected 1 axis"); } if (subtype == "spherical") { - if (axisList.size() == 3) { + if (axisCount == 2) { + // Extension to ISO19111 to support (planet)-ocentric CS with + // geocentric latitude + return SphericalCS::create(csMap, axisList[0], axisList[1]); + } else if (axisCount == 3) { return SphericalCS::create(csMap, axisList[0], axisList[1], axisList[2]); } - throw ParsingException("Expected 3 axis"); + throw ParsingException("Expected 2 or 3 axis"); } if (subtype == "ordinal") { return OrdinalCS::create(csMap, axisList); } if (subtype == "parametric") { - if (axisList.size() == 1) { + if (axisCount == 1) { return ParametricCS::create(csMap, axisList[0]); } throw ParsingException("Expected 1 axis"); } if (subtype == "TemporalDateTime") { - if (axisList.size() == 1) { + if (axisCount == 1) { return DateTimeTemporalCS::create(csMap, axisList[0]); } throw ParsingException("Expected 1 axis"); } if (subtype == "TemporalCount") { - if (axisList.size() == 1) { + if (axisCount == 1) { return TemporalCountCS::create(csMap, axisList[0]); } throw ParsingException("Expected 1 axis"); } if (subtype == "TemporalMeasure") { - if (axisList.size() == 1) { + if (axisCount == 1) { return TemporalMeasureCS::create(csMap, axisList[0]); } throw ParsingException("Expected 1 axis"); diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 72ac18689f..96e3a84947 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -894,6 +894,37 @@ TEST(wkt_parse, wkt2_EPSG_4979) { // --------------------------------------------------------------------------- +TEST(wkt_parse, wkt2_spherical_planetocentric) { + const auto wkt = + "GEODCRS[\"Mercury (2015) / Ocentric\",\n" + " DATUM[\"Mercury (2015)\",\n" + " ELLIPSOID[\"Mercury (2015)\",2440530,1075.12334801762,\n" + " LENGTHUNIT[\"metre\",1]],\n" + " ANCHOR[\"Hun Kal: 20.0\"]],\n" + " PRIMEM[\"Reference Meridian\",0,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " CS[spherical,2],\n" + " AXIS[\"planetocentric latitude (U)\",north,\n" + " ORDER[1],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " AXIS[\"planetocentric longitude (V)\",east,\n" + " ORDER[2],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " ID[\"IAU\",19902,2015],\n" + " REMARK[\"Source of IAU Coordinate systems: " + "doi://10.1007/s10569-017-9805-5\"]]"; + auto obj = WKTParser().createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_TRUE(crs->isSphericalPlanetocentric()); + EXPECT_EQ( + crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT2_2019).get()), + wkt); +} + +// --------------------------------------------------------------------------- + static void checkGeocentric(GeodeticCRSPtr crs) { // Explicitly check this is NOT a GeographicCRS EXPECT_TRUE(!dynamic_cast(crs.get())); @@ -11476,6 +11507,59 @@ TEST(json_import, geographic_crs) { // --------------------------------------------------------------------------- +TEST(json_import, spherical_planetocentric) { + const auto json = "{\n" + " \"$schema\": \"foo\",\n" + " \"type\": \"GeodeticCRS\",\n" + " \"name\": \"Mercury (2015) / Ocentric\",\n" + " \"datum\": {\n" + " \"type\": \"GeodeticReferenceFrame\",\n" + " \"name\": \"Mercury (2015)\",\n" + " \"anchor\": \"Hun Kal: 20.0\",\n" + " \"ellipsoid\": {\n" + " \"name\": \"Mercury (2015)\",\n" + " \"semi_major_axis\": 2440530,\n" + " \"inverse_flattening\": 1075.12334801762\n" + " },\n" + " \"prime_meridian\": {\n" + " \"name\": \"Reference Meridian\",\n" + " \"longitude\": 0\n" + " }\n" + " },\n" + " \"coordinate_system\": {\n" + " \"subtype\": \"spherical\",\n" + " \"axis\": [\n" + " {\n" + " \"name\": \"Planetocentric latitude\",\n" + " \"abbreviation\": \"U\",\n" + " \"direction\": \"north\",\n" + " \"unit\": \"degree\"\n" + " },\n" + " {\n" + " \"name\": \"Planetocentric longitude\",\n" + " \"abbreviation\": \"V\",\n" + " \"direction\": \"east\",\n" + " \"unit\": \"degree\"\n" + " }\n" + " ]\n" + " },\n" + " \"id\": {\n" + " \"authority\": \"IAU\",\n" + " \"code\": 19902\n" + " },\n" + " \"remarks\": \"Source of IAU Coordinate systems: " + "doi://10.1007/s10569-017-9805-5\"\n" + "}"; + auto obj = createFromUserInput(json, nullptr); + auto gcrs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(gcrs != nullptr); + EXPECT_TRUE(gcrs->isSphericalPlanetocentric()); + EXPECT_EQ(gcrs->exportToJSON(&(JSONFormatter::create()->setSchema("foo"))), + json); +} + +// --------------------------------------------------------------------------- + TEST(json_import, geographic_crs_errors) { EXPECT_THROW( createFromUserInput( From 85733181ee7c2777139f5d1db94f2beabb737e96 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 8 Sep 2021 12:50:37 +0200 Subject: [PATCH 2/4] createOperations(): deal with spherical planetocentric geodetic CRS This also fixes conversion between geocentric latlong and geodetic latlong with cs2cs. This was dealt with in PR 1093, but in the wrong direction (the geocentric latitude must be <= in absolute value to the geodetic one) The issue here was linked to the semantics of the +geoc specifier, which affects the semantics of the input coordinates in the forward direction (+geoc means that the input coordinate is is a geocentric latitude), which defeats the logic of doing A to B by using the inverse path of A and the forward path of B. --- include/proj/crs.hpp | 11 +- src/apps/cs2cs.cpp | 33 ++-- src/iso19111/crs.cpp | 149 +++++++++-------- src/iso19111/io.cpp | 152 +++++++++++++----- src/iso19111/operation/conversion.cpp | 10 +- .../operation/coordinateoperationfactory.cpp | 123 ++++++++++++++ src/iso19111/operation/transformation.cpp | 49 +++--- test/cli/testvarious | 12 +- test/cli/tv_out.dist | 4 +- test/unit/test_io.cpp | 27 +++- test/unit/test_operationfactory.cpp | 126 ++++++++++++++- 11 files changed, 525 insertions(+), 171 deletions(-) diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp index 7fde88c876..5a8e75aeac 100644 --- a/include/proj/crs.hpp +++ b/include/proj/crs.hpp @@ -310,6 +310,9 @@ class PROJ_GCC_DLL GeodeticCRS : virtual public SingleCRS, PROJ_INTERNAL void addGeocentricUnitConversionIntoPROJString( io::PROJStringFormatter *formatter) const; + PROJ_INTERNAL void + addAngularUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter) const; + PROJ_INTERNAL void _exportToWKT(io::WKTFormatter *formatter) const override; // throw(io::FormattingException) @@ -402,12 +405,10 @@ class PROJ_GCC_DLL GeographicCRS : public GeodeticCRS { PROJ_PRIVATE : //! @cond Doxygen_Suppress - PROJ_INTERNAL void - addAngularUnitConvertAndAxisSwap( - io::PROJStringFormatter *formatter) const; - PROJ_INTERNAL void _exportToPROJString(io::PROJStringFormatter *formatter) - const override; // throw(FormattingException) + PROJ_INTERNAL void + _exportToPROJString(io::PROJStringFormatter *formatter) + const override; // throw(FormattingException) PROJ_INTERNAL void _exportToJSON(io::JSONFormatter *formatter) const override; // throw(FormattingException) diff --git a/src/apps/cs2cs.cpp b/src/apps/cs2cs.cpp index 19f924d4a7..03fd4dbac5 100644 --- a/src/apps/cs2cs.cpp +++ b/src/apps/cs2cs.cpp @@ -59,10 +59,10 @@ static PJ *transformation = nullptr; -static bool srcIsGeog = false; +static bool srcIsLongLat = false; static double srcToRadians = 0.0; -static bool destIsGeog = false; +static bool destIsLongLat = false; static double destToRadians = 0.0; static bool destIsLatLong = false; @@ -158,7 +158,7 @@ static void process(FILE *fid) if (data.u != HUGE_VAL) { - if (srcIsGeog) { + if (srcIsLongLat) { /* dmstor gives values to radians. Convert now to the SRS unit */ data.u /= srcToRadians; @@ -179,7 +179,7 @@ static void process(FILE *fid) if (data.u == HUGE_VAL) /* error output */ fputs(oterr, stdout); - else if (destIsGeog && !oform) { /*ascii DMS output */ + else if (destIsLongLat && !oform) { /*ascii DMS output */ // rtodms() expect radians: convert from the output SRS unit data.u *= destToRadians; @@ -206,7 +206,7 @@ static void process(FILE *fid) } } else { /* x-y or decimal degree ascii output */ - if (destIsGeog) { + if (destIsLongLat) { data.v *= destToRadians * RAD_TO_DEG; data.u *= destToRadians * RAD_TO_DEG; } @@ -239,7 +239,7 @@ static void process(FILE *fid) /************************************************************************/ static PJ *instantiate_crs(const std::string &definition, - bool &isGeog, double &toRadians, + bool &isLongLatCS, double &toRadians, bool &isLatFirst) { PJ *crs = proj_create(nullptr, pj_add_type_crs_if_needed(definition).c_str()); @@ -247,7 +247,7 @@ static PJ *instantiate_crs(const std::string &definition, return nullptr; } - isGeog = false; + isLongLatCS = false; toRadians = 0.0; isLatFirst = false; @@ -259,11 +259,11 @@ static PJ *instantiate_crs(const std::string &definition, type = proj_get_type(crs); } if (type == PJ_TYPE_GEOGRAPHIC_2D_CRS || - type == PJ_TYPE_GEOGRAPHIC_3D_CRS) { + type == PJ_TYPE_GEOGRAPHIC_3D_CRS || + type == PJ_TYPE_GEODETIC_CRS) { auto cs = proj_crs_get_coordinate_system(nullptr, crs); assert(cs); - isGeog = true; const char *axisName = ""; proj_cs_get_axis_info(nullptr, cs, 0, &axisName, // name, @@ -277,6 +277,9 @@ static PJ *instantiate_crs(const std::string &definition, isLatFirst = NS_PROJ::internal::ci_find(std::string(axisName), "latitude") != std::string::npos; + isLongLatCS = isLatFirst || + NS_PROJ::internal::ci_find(std::string(axisName), "longitude") != + std::string::npos; proj_destroy(cs); } @@ -736,7 +739,7 @@ int main(int argc, char **argv) { PJ *src = nullptr; if (!fromStr.empty()) { bool ignored; - src = instantiate_crs(fromStr, srcIsGeog, + src = instantiate_crs(fromStr, srcIsLongLat, srcToRadians, ignored); if (!src) { emess(3, "cannot instantiate source coordinate system"); @@ -745,7 +748,7 @@ int main(int argc, char **argv) { PJ *dst = nullptr; if (!toStr.empty()) { - dst = instantiate_crs(toStr, destIsGeog, + dst = instantiate_crs(toStr, destIsLongLat, destToRadians, destIsLatLong); if (!dst) { emess(3, "cannot instantiate target coordinate system"); @@ -760,7 +763,7 @@ int main(int argc, char **argv) { emess(3, "missing target CRS and source CRS is not a projected CRS"); } - destIsGeog = true; + destIsLongLat = true; } else if (fromStr.empty()) { assert(dst); bool ignored; @@ -770,7 +773,7 @@ int main(int argc, char **argv) { emess(3, "missing source CRS and target CRS is not a projected CRS"); } - srcIsGeog = true; + srcIsLongLat = true; } proj_destroy(src); proj_destroy(dst); @@ -835,13 +838,13 @@ int main(int argc, char **argv) { } /* set input formatting control */ - if (!srcIsGeog) + if (!srcIsLongLat) informat = strtod; else { informat = dmstor; } - if (!destIsGeog && !oform) + if (!destIsLongLat && !oform) oform = "%.2f"; /* process input file list */ diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 6f5110437a..be5938785e 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -1995,6 +1995,61 @@ void GeodeticCRS::addGeocentricUnitConversionIntoPROJString( // --------------------------------------------------------------------------- +//! @cond Doxygen_Suppress +void GeodeticCRS::addAngularUnitConvertAndAxisSwap( + io::PROJStringFormatter *formatter) const { + const auto &axisList = coordinateSystem()->axisList(); + + formatter->addStep("unitconvert"); + formatter->addParam("xy_in", "rad"); + if (axisList.size() == 3 && !formatter->omitZUnitConversion()) { + formatter->addParam("z_in", "m"); + } + { + const auto &unitHoriz = axisList[0]->unit(); + const auto projUnit = unitHoriz.exportToPROJString(); + if (projUnit.empty()) { + formatter->addParam("xy_out", unitHoriz.conversionToSI()); + } else { + formatter->addParam("xy_out", projUnit); + } + } + if (axisList.size() == 3 && !formatter->omitZUnitConversion()) { + const auto &unitZ = axisList[2]->unit(); + auto projVUnit = unitZ.exportToPROJString(); + if (projVUnit.empty()) { + formatter->addParam("z_out", unitZ.conversionToSI()); + } else { + formatter->addParam("z_out", projVUnit); + } + } + + const char *order[2] = {nullptr, nullptr}; + const char *one = "1"; + const char *two = "2"; + for (int i = 0; i < 2; i++) { + const auto &dir = axisList[i]->direction(); + if (&dir == &cs::AxisDirection::WEST) { + order[i] = "-1"; + } else if (&dir == &cs::AxisDirection::EAST) { + order[i] = one; + } else if (&dir == &cs::AxisDirection::SOUTH) { + order[i] = "-2"; + } else if (&dir == &cs::AxisDirection::NORTH) { + order[i] = two; + } + } + if (order[0] && order[1] && (order[0] != one || order[1] != two)) { + formatter->addStep("axisswap"); + char orderStr[10]; + sprintf(orderStr, "%.2s,%.2s", order[0], order[1]); + formatter->addParam("order", orderStr); + } +} +//! @endcond + +// --------------------------------------------------------------------------- + //! @cond Doxygen_Suppress void GeodeticCRS::_exportToPROJString( io::PROJStringFormatter *formatter) const // throw(io::FormattingException) @@ -2007,19 +2062,38 @@ void GeodeticCRS::_exportToPROJString( return; } - if (!isGeocentric()) { - io::FormattingException::Throw( - "GeodeticCRS::exportToPROJString() only " - "supports geocentric coordinate systems"); - } + if (isGeocentric()) { + if (!formatter->getCRSExport()) { + formatter->addStep("cart"); + } else { + formatter->addStep("geocent"); + } - if (!formatter->getCRSExport()) { - formatter->addStep("cart"); + addDatumInfoToPROJString(formatter); + addGeocentricUnitConversionIntoPROJString(formatter); + } else if (isSphericalPlanetocentric()) { + if (!formatter->getCRSExport()) { + formatter->addStep("geoc"); + addDatumInfoToPROJString(formatter); + addAngularUnitConvertAndAxisSwap(formatter); + } else { + io::FormattingException::Throw( + "GeodeticCRS::exportToPROJString() not supported on spherical " + "planetocentric coordinate systems"); + // The below code now works as input to PROJ, but I'm not sure we + // want to propagate this, given that we got cs2cs doing conversion + // in the wrong direction in past versions. + /*formatter->addStep("longlat"); + formatter->addParam("geoc"); + + addDatumInfoToPROJString(formatter);*/ + } } else { - formatter->addStep("geocent"); + io::FormattingException::Throw( + "GeodeticCRS::exportToPROJString() only " + "supports geocentric or spherical planetocentric " + "coordinate systems"); } - addDatumInfoToPROJString(formatter); - addGeocentricUnitConversionIntoPROJString(formatter); } //! @endcond @@ -2883,61 +2957,6 @@ GeographicCRS::demoteTo2D(const std::string &newName, // --------------------------------------------------------------------------- -//! @cond Doxygen_Suppress -void GeographicCRS::addAngularUnitConvertAndAxisSwap( - io::PROJStringFormatter *formatter) const { - const auto &axisList = coordinateSystem()->axisList(); - - formatter->addStep("unitconvert"); - formatter->addParam("xy_in", "rad"); - if (axisList.size() == 3 && !formatter->omitZUnitConversion()) { - formatter->addParam("z_in", "m"); - } - { - const auto &unitHoriz = axisList[0]->unit(); - const auto projUnit = unitHoriz.exportToPROJString(); - if (projUnit.empty()) { - formatter->addParam("xy_out", unitHoriz.conversionToSI()); - } else { - formatter->addParam("xy_out", projUnit); - } - } - if (axisList.size() == 3 && !formatter->omitZUnitConversion()) { - const auto &unitZ = axisList[2]->unit(); - auto projVUnit = unitZ.exportToPROJString(); - if (projVUnit.empty()) { - formatter->addParam("z_out", unitZ.conversionToSI()); - } else { - formatter->addParam("z_out", projVUnit); - } - } - - const char *order[2] = {nullptr, nullptr}; - const char *one = "1"; - const char *two = "2"; - for (int i = 0; i < 2; i++) { - const auto &dir = axisList[i]->direction(); - if (&dir == &cs::AxisDirection::WEST) { - order[i] = "-1"; - } else if (&dir == &cs::AxisDirection::EAST) { - order[i] = one; - } else if (&dir == &cs::AxisDirection::SOUTH) { - order[i] = "-2"; - } else if (&dir == &cs::AxisDirection::NORTH) { - order[i] = two; - } - } - if (order[0] && order[1] && (order[0] != one || order[1] != two)) { - formatter->addStep("axisswap"); - char orderStr[10]; - sprintf(orderStr, "%.2s,%.2s", order[0], order[1]); - formatter->addParam("order", orderStr); - } -} -//! @endcond - -// --------------------------------------------------------------------------- - //! @cond Doxygen_Suppress void GeographicCRS::_exportToPROJString( io::PROJStringFormatter *formatter) const // throw(io::FormattingException) diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 1373c58304..5d7e951e2d 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -8639,10 +8639,10 @@ struct PROJStringParser::Private { PrimeMeridianNNPtr buildPrimeMeridian(Step &step); GeodeticReferenceFrameNNPtr buildDatum(Step &step, const std::string &title); - GeographicCRSNNPtr buildGeographicCRS(int iStep, int iUnitConvert, - int iAxisSwap, bool ignorePROJAxis); + GeodeticCRSNNPtr buildGeodeticCRS(int iStep, int iUnitConvert, + int iAxisSwap, bool ignorePROJAxis); GeodeticCRSNNPtr buildGeocentricCRS(int iStep, int iUnitConvert); - CRSNNPtr buildProjectedCRS(int iStep, GeographicCRSNNPtr geogCRS, + CRSNNPtr buildProjectedCRS(int iStep, const GeodeticCRSNNPtr &geogCRS, int iUnitConvert, int iAxisSwap); CRSNNPtr buildBoundOrCompoundCRSIfNeeded(int iStep, CRSNNPtr crs); UnitOfMeasure buildUnit(Step &step, const std::string &unitsParamName, @@ -8656,6 +8656,9 @@ struct PROJStringParser::Private { EllipsoidalCSNNPtr buildEllipsoidalCS(int iStep, int iUnitConvert, int iAxisSwap, bool ignorePROJAxis); + + SphericalCSNNPtr buildSphericalCS(int iStep, int iUnitConvert, + int iAxisSwap, bool ignorePROJAxis); }; // --------------------------------------------------------------------------- @@ -9282,10 +9285,13 @@ PROJStringParser::Private::processAxisSwap(Step &step, assert(iAxisSwap < 0 || ci_equal(steps_[iAxisSwap].name, "axisswap")); const bool isGeographic = unit.type() == UnitOfMeasure::Type::ANGULAR; + const bool isSpherical = isGeographic && hasParamValue(step, "geoc"); const auto &eastName = - isGeographic ? AxisName::Longitude : AxisName::Easting; - const auto &eastAbbev = - isGeographic ? AxisAbbreviation::lon : AxisAbbreviation::E; + isSpherical ? "Planetocentric longitude" + : isGeographic ? AxisName::Longitude : AxisName::Easting; + const auto &eastAbbev = isSpherical ? "V" + : isGeographic ? AxisAbbreviation::lon + : AxisAbbreviation::E; const auto &eastDir = isGeographic ? AxisDirection::EAST : (axisType == AxisType::NORTH_POLE) @@ -9301,9 +9307,11 @@ PROJStringParser::Private::processAxisSwap(Step &step, : nullMeridian); const auto &northName = - isGeographic ? AxisName::Latitude : AxisName::Northing; - const auto &northAbbev = - isGeographic ? AxisAbbreviation::lat : AxisAbbreviation::N; + isSpherical ? "Planetocentric latitude" + : isGeographic ? AxisName::Latitude : AxisName::Northing; + const auto &northAbbev = isSpherical ? "U" + : isGeographic ? AxisAbbreviation::lat + : AxisAbbreviation::N; const auto &northDir = isGeographic ? AxisDirection::NORTH : (axisType == AxisType::NORTH_POLE) @@ -9320,15 +9328,19 @@ PROJStringParser::Private::processAxisSwap(Step &step, .as_nullable() : nullMeridian); - CoordinateSystemAxisNNPtr west = - createAxis(isGeographic ? AxisName::Longitude : AxisName::Westing, - isGeographic ? AxisAbbreviation::lon : std::string(), - AxisDirection::WEST, unit); + CoordinateSystemAxisNNPtr west = createAxis( + isSpherical ? "Planetocentric longitude" + : isGeographic ? AxisName::Longitude : AxisName::Westing, + isSpherical ? "V" + : isGeographic ? AxisAbbreviation::lon : std::string(), + AxisDirection::WEST, unit); - CoordinateSystemAxisNNPtr south = - createAxis(isGeographic ? AxisName::Latitude : AxisName::Southing, - isGeographic ? AxisAbbreviation::lat : std::string(), - AxisDirection::SOUTH, unit); + CoordinateSystemAxisNNPtr south = createAxis( + isSpherical ? "Planetocentric latitude" + : isGeographic ? AxisName::Latitude : AxisName::Southing, + isSpherical ? "U" + : isGeographic ? AxisAbbreviation::lat : std::string(), + AxisDirection::SOUTH, unit); std::vector axis{east, north}; @@ -9428,6 +9440,42 @@ EllipsoidalCSNNPtr PROJStringParser::Private::buildEllipsoidalCS( // --------------------------------------------------------------------------- +SphericalCSNNPtr PROJStringParser::Private::buildSphericalCS( + int iStep, int iUnitConvert, int iAxisSwap, bool ignorePROJAxis) { + auto &step = steps_[iStep]; + assert(iUnitConvert < 0 || + ci_equal(steps_[iUnitConvert].name, "unitconvert")); + + UnitOfMeasure angularUnit = UnitOfMeasure::DEGREE; + if (iUnitConvert >= 0) { + auto &stepUnitConvert = steps_[iUnitConvert]; + const std::string *xy_in = &getParamValue(stepUnitConvert, "xy_in"); + const std::string *xy_out = &getParamValue(stepUnitConvert, "xy_out"); + if (stepUnitConvert.inverted) { + std::swap(xy_in, xy_out); + } + if (iUnitConvert < iStep) { + std::swap(xy_in, xy_out); + } + if (xy_in->empty() || xy_out->empty() || *xy_in != "rad" || + (*xy_out != "rad" && *xy_out != "deg" && *xy_out != "grad")) { + throw ParsingException("unhandled values for xy_in and/or xy_out"); + } + if (*xy_out == "rad") { + angularUnit = UnitOfMeasure::RADIAN; + } else if (*xy_out == "grad") { + angularUnit = UnitOfMeasure::GRAD; + } + } + + std::vector axis = processAxisSwap( + step, angularUnit, iAxisSwap, AxisType::REGULAR, ignorePROJAxis); + + return SphericalCS::create(emptyPropertyMap, axis[0], axis[1]); +} + +// --------------------------------------------------------------------------- + static double getNumericValue(const std::string ¶mValue, bool *pHasError = nullptr) { try { @@ -9447,7 +9495,7 @@ namespace { template inline void ignoreRetVal(T) {} } // namespace -GeographicCRSNNPtr PROJStringParser::Private::buildGeographicCRS( +GeodeticCRSNNPtr PROJStringParser::Private::buildGeodeticCRS( int iStep, int iUnitConvert, int iAxisSwap, bool ignorePROJAxis) { auto &step = steps_[iStep]; @@ -9462,8 +9510,6 @@ GeographicCRSNNPtr PROJStringParser::Private::buildGeographicCRS( auto props = PropertyMap().set(IdentifiedObject::NAME_KEY, title.empty() ? "unknown" : title); - auto cs = - buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, ignorePROJAxis); if (l_isGeographicStep && (hasUnusedParameters(step) || @@ -9472,7 +9518,17 @@ GeographicCRSNNPtr PROJStringParser::Private::buildGeographicCRS( } props.set("IMPLICIT_CS", true); - return GeographicCRS::create(props, datum, cs); + if (!hasParamValue(step, "geoc")) { + auto cs = + buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, ignorePROJAxis); + + return GeographicCRS::create(props, datum, cs); + } else { + auto cs = + buildSphericalCS(iStep, iUnitConvert, iAxisSwap, ignorePROJAxis); + + return GeodeticCRS::create(props, datum, cs); + } } // --------------------------------------------------------------------------- @@ -9623,8 +9679,10 @@ static bool is_in_stringlist(const std::string &str, const char *stringlist) { // --------------------------------------------------------------------------- -CRSNNPtr PROJStringParser::Private::buildProjectedCRS( - int iStep, GeographicCRSNNPtr geogCRS, int iUnitConvert, int iAxisSwap) { +CRSNNPtr +PROJStringParser::Private::buildProjectedCRS(int iStep, + const GeodeticCRSNNPtr &geodCRS, + int iUnitConvert, int iAxisSwap) { auto &step = steps_[iStep]; const auto mappings = getMappingsFromPROJName(step.name); const MethodMapping *mapping = mappings.empty() ? nullptr : mappings[0]; @@ -9650,7 +9708,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( } if (mapping) { - mapping = selectSphericalOrEllipsoidal(mapping, geogCRS); + mapping = selectSphericalOrEllipsoidal(mapping, geodCRS); } assert(isProjectedStep(step.name)); @@ -9660,7 +9718,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( const auto &title = title_; if (!buildPrimeMeridian(step)->longitude()._isEquivalentTo( - geogCRS->primeMeridian()->longitude(), + geodCRS->primeMeridian()->longitude(), util::IComparable::Criterion::EQUIVALENT)) { throw ParsingException("inconsistent pm values between projectedCRS " "and its base geographicalCRS"); @@ -9926,7 +9984,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( } if (k >= 0 && k <= 1) { const double es = - geogCRS->ellipsoid()->squaredEccentricity(); + geodCRS->ellipsoid()->squaredEccentricity(); if (es < 0 || es == 1) { throw ParsingException("Invalid flattening"); } @@ -10037,10 +10095,16 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( {"PROJ ob_tran o_proj=longlat", "PROJ ob_tran o_proj=lonlat", "PROJ ob_tran o_proj=latlon", "PROJ ob_tran o_proj=latlong"}) { if (starts_with(methodName, substr)) { - return DerivedGeographicCRS::create( - PropertyMap().set(IdentifiedObject::NAME_KEY, "unnamed"), - geogCRS, NN_NO_CHECK(conv), - buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, false)); + auto geogCRS = + util::nn_dynamic_pointer_cast(geodCRS); + if (geogCRS) { + return DerivedGeographicCRS::create( + PropertyMap().set(IdentifiedObject::NAME_KEY, + "unnamed"), + NN_NO_CHECK(geogCRS), NN_NO_CHECK(conv), + buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, + false)); + } } } } @@ -10048,11 +10112,11 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( std::vector axis = processAxisSwap(step, unit, iAxisSwap, axisType, false); - auto csGeogCRS = geogCRS->coordinateSystem(); - auto cs = csGeogCRS->axisList().size() == 2 + auto csGeodCRS = geodCRS->coordinateSystem(); + auto cs = csGeodCRS->axisList().size() == 2 ? CartesianCS::create(emptyPropertyMap, axis[0], axis[1]) : CartesianCS::create(emptyPropertyMap, axis[0], axis[1], - csGeogCRS->axisList()[2]); + csGeodCRS->axisList()[2]); auto props = PropertyMap().set(IdentifiedObject::NAME_KEY, title.empty() ? "unknown" : title); @@ -10066,7 +10130,7 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( bWebMercator ? createPseudoMercator( props.set(IdentifiedObject::NAME_KEY, webMercatorName), cs) - : ProjectedCRS::create(props, geogCRS, NN_NO_CHECK(conv), cs); + : ProjectedCRS::create(props, geodCRS, NN_NO_CHECK(conv), cs); return crs; } @@ -10537,8 +10601,8 @@ PROJStringParser::createFromPROJString(const std::string &projString) { // First run is dry run to mark all recognized/unrecognized tokens for (int iter = 0; iter < 2; iter++) { auto obj = d->buildBoundOrCompoundCRSIfNeeded( - 0, d->buildGeographicCRS(iFirstGeogStep, iFirstUnitConvert, - iFirstAxisSwap, false)); + 0, d->buildGeodeticCRS(iFirstGeogStep, iFirstUnitConvert, + iFirstAxisSwap, false)); if (iter == 1) { return nn_static_pointer_cast(obj); } @@ -10555,14 +10619,14 @@ PROJStringParser::createFromPROJString(const std::string &projString) { iProjStep, d->buildProjectedCRS( iProjStep, - d->buildGeographicCRS(iFirstGeogStep, - iFirstUnitConvert < iFirstGeogStep - ? iFirstUnitConvert - : -1, - iFirstAxisSwap < iFirstGeogStep - ? iFirstAxisSwap - : -1, - true), + d->buildGeodeticCRS(iFirstGeogStep, + iFirstUnitConvert < iFirstGeogStep + ? iFirstUnitConvert + : -1, + iFirstAxisSwap < iFirstGeogStep + ? iFirstAxisSwap + : -1, + true), iFirstUnitConvert < iFirstGeogStep ? iSecondUnitConvert : iFirstUnitConvert, iFirstAxisSwap < iFirstGeogStep ? iSecondAxisSwap diff --git a/src/iso19111/operation/conversion.cpp b/src/iso19111/operation/conversion.cpp index 3d09c9fa22..b8d9c0b750 100644 --- a/src/iso19111/operation/conversion.cpp +++ b/src/iso19111/operation/conversion.cpp @@ -3459,11 +3459,15 @@ void Conversion::_exportToPROJString( } } - srcGeogCRS = std::dynamic_pointer_cast(horiz); - if (srcGeogCRS) { + auto srcGeodCRS = dynamic_cast(horiz.get()); + if (srcGeodCRS) { + srcGeogCRS = std::dynamic_pointer_cast(horiz); + } + if (srcGeodCRS && + (srcGeogCRS || srcGeodCRS->isSphericalPlanetocentric())) { formatter->setOmitProjLongLatIfPossible(true); formatter->startInversion(); - srcGeogCRS->_exportToPROJString(formatter); + srcGeodCRS->_exportToPROJString(formatter); formatter->stopInversion(); formatter->setOmitProjLongLatIfPossible(false); } diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index fc64bd2e4d..60365713a3 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -558,6 +558,17 @@ struct CoordinateOperationFactory::Private { const crs::GeodeticCRS *geodDst, std::vector &res); + static void createOperationsFromSphericalPlanetocentric( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::GeodeticCRS *geodSrc, + std::vector &res); + + static void createOperationsFromBoundOfSphericalPlanetocentric( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::BoundCRS *boundSrc, + const crs::GeodeticCRSNNPtr &geodSrcBase, + std::vector &res); + static void createOperationsDerivedTo( const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, Private::Context &context, const crs::DerivedCRS *derivedSrc, @@ -2989,6 +3000,32 @@ CoordinateOperationFactory::Private::createOperations( return res; } + if (geodSrc && geodSrc->isSphericalPlanetocentric()) { + createOperationsFromSphericalPlanetocentric(sourceCRS, targetCRS, + context, geodSrc, res); + return res; + } else if (geodDst && geodDst->isSphericalPlanetocentric()) { + return applyInverse(createOperations(targetCRS, sourceCRS, context)); + } + + if (boundSrc) { + auto geodSrcBase = util::nn_dynamic_pointer_cast( + boundSrc->baseCRS()); + if (geodSrcBase && geodSrcBase->isSphericalPlanetocentric()) { + createOperationsFromBoundOfSphericalPlanetocentric( + sourceCRS, targetCRS, context, boundSrc, + NN_NO_CHECK(geodSrcBase), res); + return res; + } + } else if (boundDst) { + auto geodDstBase = util::nn_dynamic_pointer_cast( + boundDst->baseCRS()); + if (geodDstBase && geodDstBase->isSphericalPlanetocentric()) { + return applyInverse( + createOperations(targetCRS, sourceCRS, context)); + } + } + // If the source is a derived CRS, then chain the inverse of its // deriving conversion, with transforms from its baseCRS to the // targetCRS @@ -3895,6 +3932,92 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod( // --------------------------------------------------------------------------- +void CoordinateOperationFactory::Private:: + createOperationsFromSphericalPlanetocentric( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::GeodeticCRS *geodSrc, + std::vector &res) { + + ENTER_FUNCTION(); + + // Create an intermediate geographic CRS with the same datum as the + // source spherical planetocentric one + std::string interm_crs_name(geodSrc->nameStr()); + interm_crs_name += " (geographic)"; + auto interm_crs = + util::nn_static_pointer_cast(crs::GeographicCRS::create( + addDomains(util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, interm_crs_name), + geodSrc), + geodSrc->datum(), geodSrc->datumEnsemble(), + cs::EllipsoidalCS::createLatitudeLongitude( + common::UnitOfMeasure::DEGREE))); + + auto opFirst = createGeodToGeodPROJBased(sourceCRS, interm_crs); + auto opsSecond = createOperations(interm_crs, targetCRS, context); + for (const auto &opSecond : opsSecond) { + try { + res.emplace_back(ConcatenatedOperation::createComputeMetadata( + {opFirst, opSecond}, disallowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { + } + } +} + +// --------------------------------------------------------------------------- + +void CoordinateOperationFactory::Private:: + createOperationsFromBoundOfSphericalPlanetocentric( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + Private::Context &context, const crs::BoundCRS *boundSrc, + const crs::GeodeticCRSNNPtr &geodSrcBase, + std::vector &res) { + + ENTER_FUNCTION(); + + // Create an intermediate geographic CRS with the same datum as the + // source spherical planetocentric one + std::string interm_crs_name(geodSrcBase->nameStr()); + interm_crs_name += " (geographic)"; + auto intermGeog = + util::nn_static_pointer_cast(crs::GeographicCRS::create( + addDomains(util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, interm_crs_name), + geodSrcBase.get()), + geodSrcBase->datum(), geodSrcBase->datumEnsemble(), + cs::EllipsoidalCS::createLatitudeLongitude( + common::UnitOfMeasure::DEGREE))); + + // Create an intermediate boundCRS wrapping the above intermediate + // geographic CRS + auto transf = boundSrc->transformation()->shallowClone(); + // keep a reference to the target before patching it with itself + // (this is due to our abuse of passing shared_ptr by reference + auto transfTarget = transf->targetCRS(); + setCRSs(transf.get(), intermGeog, transfTarget); + + auto intermBoundCRS = + crs::BoundCRS::create(intermGeog, boundSrc->hubCRS(), transf); + + auto opFirst = createGeodToGeodPROJBased(geodSrcBase, intermGeog); + setCRSs(opFirst.get(), sourceCRS, intermBoundCRS); + auto opsSecond = createOperations(intermBoundCRS, targetCRS, context); + for (const auto &opSecond : opsSecond) { + try { + auto opSecondClone = opSecond->shallowClone(); + // In theory, we should not need that setCRSs() forcing, but due + // how BoundCRS transformations are implemented currently, we + // need it in practice. + setCRSs(opSecondClone.get(), intermBoundCRS, targetCRS); + res.emplace_back(ConcatenatedOperation::createComputeMetadata( + {opFirst, opSecondClone}, disallowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { + } + } +} + +// --------------------------------------------------------------------------- + void CoordinateOperationFactory::Private::createOperationsDerivedTo( const crs::CRSNNPtr & /*sourceCRS*/, const crs::CRSNNPtr &targetCRS, Private::Context &context, const crs::DerivedCRS *derivedSrc, diff --git a/src/iso19111/operation/transformation.cpp b/src/iso19111/operation/transformation.cpp index 273b636f23..a30e12488b 100644 --- a/src/iso19111/operation/transformation.cpp +++ b/src/iso19111/operation/transformation.cpp @@ -474,13 +474,16 @@ static void getTransformationType(const crs::CRSNNPtr &sourceCRSIn, dynamic_cast(sourceCRSIn.get()); auto targetCRSGeog = dynamic_cast(targetCRSIn.get()); - if (!sourceCRSGeog || !targetCRSGeog) { + if (!(sourceCRSGeog || + (sourceCRSGeod && sourceCRSGeod->isSphericalPlanetocentric())) || + !(targetCRSGeog || + (targetCRSGeod && targetCRSGeod->isSphericalPlanetocentric()))) { throw InvalidOperation("Inconsistent CRS type"); } const auto nSrcAxisCount = - sourceCRSGeog->coordinateSystem()->axisList().size(); + sourceCRSGeod->coordinateSystem()->axisList().size(); const auto nTargetAxisCount = - targetCRSGeog->coordinateSystem()->axisList().size(); + targetCRSGeod->coordinateSystem()->axisList().size(); isGeog2D = nSrcAxisCount == 2 && nTargetAxisCount == 2; isGeog3D = !isGeog2D && nSrcAxisCount >= 2 && nTargetAxisCount >= 2; } @@ -1003,36 +1006,36 @@ TransformationNNPtr Transformation::createTOWGS84( "Invalid number of elements in TOWGS84Parameters"); } - crs::CRSPtr transformSourceCRS = sourceCRSIn->extractGeodeticCRS(); - if (!transformSourceCRS) { + auto transformSourceGeodCRS = sourceCRSIn->extractGeodeticCRS(); + if (!transformSourceGeodCRS) { throw InvalidOperation( "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); - + concat("Transformation from ", + transformSourceGeodCRS->nameStr(), " to WGS84")); + + auto targetCRS = dynamic_cast( + transformSourceGeodCRS.get()) || + transformSourceGeodCRS->isSphericalPlanetocentric() + ? util::nn_static_pointer_cast( + crs::GeographicCRS::EPSG_4326) + : util::nn_static_pointer_cast( + crs::GeodeticCRS::EPSG_4978); + + crs::CRSNNPtr transformSourceCRS = NN_NO_CHECK(transformSourceGeodCRS); if (TOWGS84Parameters.size() == 3) { return createGeocentricTranslations( - properties, NN_NO_CHECK(transformSourceCRS), targetCRS, - TOWGS84Parameters[0], TOWGS84Parameters[1], TOWGS84Parameters[2], - {}); + properties, transformSourceCRS, targetCRS, TOWGS84Parameters[0], + TOWGS84Parameters[1], TOWGS84Parameters[2], {}); } - return createPositionVector(properties, NN_NO_CHECK(transformSourceCRS), - targetCRS, TOWGS84Parameters[0], - TOWGS84Parameters[1], TOWGS84Parameters[2], - TOWGS84Parameters[3], TOWGS84Parameters[4], - TOWGS84Parameters[5], TOWGS84Parameters[6], {}); + return createPositionVector( + properties, transformSourceCRS, targetCRS, TOWGS84Parameters[0], + TOWGS84Parameters[1], TOWGS84Parameters[2], TOWGS84Parameters[3], + TOWGS84Parameters[4], TOWGS84Parameters[5], TOWGS84Parameters[6], {}); } // --------------------------------------------------------------------------- diff --git a/test/cli/testvarious b/test/cli/testvarious index 6e9cd43cca..803403696a 100755 --- a/test/cli/testvarious +++ b/test/cli/testvarious @@ -144,10 +144,10 @@ $EXE +proj=geocent +datum=WGS84 \ EOF # echo "#############################################################" >> ${OUT} -echo Test conversion from geocentric latlong to geodetic latlong >> ${OUT} +echo Test conversion from geodetic latlong to geocentric latlong >> ${OUT} # -$EXE +proj=latlong +datum=WGS84 +geoc \ - +to +proj=latlong +datum=WGS84 \ +$EXE +proj=latlong +datum=WGS84 \ + +to +proj=latlong +datum=WGS84 +geoc \ -E >>${OUT} <> ${OUT} -echo Test conversion from geodetic latlong to geocentric latlong >> ${OUT} +echo Test conversion from geocentric latlong to geodetic latlong >> ${OUT} # -$EXE +proj=latlong +datum=WGS84 \ - +to +proj=latlong +datum=WGS84 +geoc \ +$EXE +proj=latlong +datum=WGS84 +geoc \ + +to +proj=latlong +datum=WGS84 \ -E >>${OUT} <(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_TRUE(crs->isSphericalPlanetocentric()); +#if 1 + EXPECT_THROW( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + FormattingException); +#else + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + input); +#endif +} + +// --------------------------------------------------------------------------- + TEST(io, projparse_projected_wktext) { std::string input("+proj=merc +foo +wktext +type=crs"); auto obj = PROJStringParser().createFromPROJString(input); @@ -10209,13 +10232,13 @@ TEST(io, projparse_init) { } { - auto obj = createFromUserInput("+title=mytitle +geoc +init=epsg:4326", + auto obj = createFromUserInput("+title=mytitle +over +init=epsg:4326", dbContext, true); auto crs = nn_dynamic_pointer_cast(obj); ASSERT_TRUE(crs != nullptr); EXPECT_EQ(crs->nameStr(), "mytitle"); EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()), - "+proj=longlat +geoc +datum=WGS84 +no_defs +type=crs"); + "+proj=longlat +over +datum=WGS84 +no_defs +type=crs"); } { diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 56d9f743ed..ec7d870014 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -6092,7 +6092,7 @@ TEST(operation, createOperation_on_crs_with_canonical_bound_crs) { TEST(operation, createOperation_fallback_to_proj4_strings) { auto objDest = PROJStringParser().createFromPROJString( - "+proj=longlat +geoc +datum=WGS84 +type=crs"); + "+proj=longlat +over +datum=WGS84 +type=crs"); auto dest = nn_dynamic_pointer_cast(objDest); ASSERT_TRUE(dest != nullptr); @@ -6102,7 +6102,7 @@ TEST(operation, createOperation_fallback_to_proj4_strings) { EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline +step +proj=axisswap +order=2,1 " "+step +proj=unitconvert +xy_in=deg +xy_out=rad " - "+step +proj=longlat +geoc +datum=WGS84 " + "+step +proj=longlat +over +datum=WGS84 " "+step +proj=unitconvert +xy_in=rad +xy_out=deg"); } @@ -6110,12 +6110,12 @@ TEST(operation, createOperation_fallback_to_proj4_strings) { TEST(operation, createOperation_fallback_to_proj4_strings_bound_of_geog) { auto objSrc = PROJStringParser().createFromPROJString( - "+proj=longlat +geoc +ellps=GRS80 +towgs84=0,0,0 +type=crs"); + "+proj=longlat +over +ellps=GRS80 +towgs84=0,0,0 +type=crs"); auto src = nn_dynamic_pointer_cast(objSrc); ASSERT_TRUE(src != nullptr); auto objDest = PROJStringParser().createFromPROJString( - "+proj=longlat +geoc +ellps=clrk66 +towgs84=0,0,0 +type=crs"); + "+proj=longlat +over +ellps=clrk66 +towgs84=0,0,0 +type=crs"); auto dest = nn_dynamic_pointer_cast(objDest); ASSERT_TRUE(dest != nullptr); @@ -6125,8 +6125,8 @@ TEST(operation, createOperation_fallback_to_proj4_strings_bound_of_geog) { EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline " "+step +proj=unitconvert +xy_in=deg +xy_out=rad " - "+step +inv +proj=longlat +geoc +ellps=GRS80 +towgs84=0,0,0 " - "+step +proj=longlat +geoc +ellps=clrk66 +towgs84=0,0,0 " + "+step +inv +proj=longlat +over +ellps=GRS80 +towgs84=0,0,0 " + "+step +proj=longlat +over +ellps=clrk66 +towgs84=0,0,0 " "+step +proj=unitconvert +xy_in=rad +xy_out=deg"); } @@ -6832,3 +6832,117 @@ TEST(operation, derivedGeographicCRS_with_to_wgs84_to_geographicCRS) { "+proj=pipeline +step +proj=axisswap +order=2,1 " + pipeline); } } + +// --------------------------------------------------------------------------- + +TEST(operation, createOperation_spherical_ocentric_to_geographic) { + auto objSrc = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +datum=WGS84 +type=crs"); + auto src = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(src != nullptr); + + auto op = CoordinateOperationFactory::create()->createOperation( + NN_CHECK_ASSERT(src), GeographicCRS::EPSG_4326); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +inv +proj=geoc +ellps=WGS84 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg " + "+step +proj=axisswap +order=2,1"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, createOperation_geographic_to_spherical_ocentric) { + auto objDest = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +datum=WGS84 +type=crs"); + auto dest = nn_dynamic_pointer_cast(objDest); + ASSERT_TRUE(dest != nullptr); + + auto op = CoordinateOperationFactory::create()->createOperation( + GeographicCRS::EPSG_4326, NN_CHECK_ASSERT(dest)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=geoc +ellps=WGS84 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, createOperation_bound_of_spherical_ocentric_to_same_type) { + auto objSrc = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +ellps=GRS80 +towgs84=1,2,3 +type=crs"); + auto src = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(src != nullptr); + + auto objDest = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +ellps=clrk66 +towgs84=4,5,6 +type=crs"); + auto dest = nn_dynamic_pointer_cast(objDest); + ASSERT_TRUE(dest != nullptr); + + auto op = CoordinateOperationFactory::create()->createOperation( + NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +inv +proj=geoc +ellps=GRS80 " + "+step +proj=push +v_3 " + "+step +proj=cart +ellps=GRS80 " + "+step +proj=helmert +x=-3 +y=-3 +z=-3 " + "+step +inv +proj=cart +ellps=clrk66 " + "+step +proj=pop +v_3 " + "+step +proj=geoc +ellps=clrk66 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, createOperation_spherical_ocentric_to_projected) { + auto objSrc = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +datum=WGS84 +type=crs"); + auto src = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(src != nullptr); + + auto objDest = PROJStringParser().createFromPROJString( + "+proj=utm +zone=11 +datum=WGS84 +type=crs"); + auto dest = nn_dynamic_pointer_cast(objDest); + ASSERT_TRUE(dest != nullptr); + + auto op = CoordinateOperationFactory::create()->createOperation( + NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +inv +proj=geoc +ellps=WGS84 " + "+step +proj=utm +zone=11 +ellps=WGS84"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, + createOperation_spherical_ocentric_to_projected_of_spherical_ocentric) { + auto objSrc = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +datum=WGS84 +type=crs"); + auto src = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(src != nullptr); + + auto objDest = PROJStringParser().createFromPROJString( + "+proj=utm +geoc +zone=11 +datum=WGS84 +type=crs"); + auto dest = nn_dynamic_pointer_cast(objDest); + ASSERT_TRUE(dest != nullptr); + + auto op = CoordinateOperationFactory::create()->createOperation( + NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +inv +proj=geoc +ellps=WGS84 " + "+step +proj=utm +zone=11 +ellps=WGS84"); +} From bc568fcc99257731a939d93cd0caa4725e6803e4 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 8 Sep 2021 14:29:41 +0200 Subject: [PATCH 3/4] createConversion(): avoid nullptr dereference on a method without parameters --- src/iso19111/operation/conversion.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/iso19111/operation/conversion.cpp b/src/iso19111/operation/conversion.cpp index b8d9c0b750..c1e5cb4431 100644 --- a/src/iso19111/operation/conversion.cpp +++ b/src/iso19111/operation/conversion.cpp @@ -264,7 +264,8 @@ createConversion(const util::PropertyMap &properties, const std::vector &values) { std::vector parameters; - for (int i = 0; mapping->params[i] != nullptr; i++) { + for (int i = 0; mapping->params != nullptr && mapping->params[i] != nullptr; + i++) { const auto *param = mapping->params[i]; auto paramProperties = util::PropertyMap().set( common::IdentifiedObject::NAME_KEY, param->wkt2_name); From 078952e7f078e029d66ab6ca1ed594dfecadd1fc Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 8 Sep 2021 14:34:50 +0200 Subject: [PATCH 4/4] createOperations(): use an explicit conversion operation for geodetic <--> geocentric latitude --- include/proj/coordinateoperation.hpp | 4 + src/iso19111/operation/conversion.cpp | 103 ++++++++++++++++++ .../operation/coordinateoperationfactory.cpp | 38 +++++-- src/iso19111/operation/parammappings.cpp | 2 + src/proj_constants.h | 4 + test/unit/test_operationfactory.cpp | 23 ++++ 6 files changed, 165 insertions(+), 9 deletions(-) diff --git a/include/proj/coordinateoperation.hpp b/include/proj/coordinateoperation.hpp index 3d2ab80025..fc1d5b8a71 100644 --- a/include/proj/coordinateoperation.hpp +++ b/include/proj/coordinateoperation.hpp @@ -1385,6 +1385,10 @@ class PROJ_GCC_DLL Conversion : public SingleOperation { createGeographicGeocentric(const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS); + PROJ_INTERNAL static ConversionNNPtr + createGeographicGeocentricLatitude(const crs::CRSNNPtr &sourceCRS, + const crs::CRSNNPtr &targetCRS); + //! @endcond protected: diff --git a/src/iso19111/operation/conversion.cpp b/src/iso19111/operation/conversion.cpp index c1e5cb4431..f7bb8ec889 100644 --- a/src/iso19111/operation/conversion.cpp +++ b/src/iso19111/operation/conversion.cpp @@ -79,6 +79,7 @@ constexpr double UTM_SCALE_FACTOR = 0.9996; constexpr double UTM_FALSE_EASTING = 500000.0; constexpr double UTM_NORTH_FALSE_NORTHING = 0.0; constexpr double UTM_SOUTH_FALSE_NORTHING = 10000000.0; + //! @endcond // --------------------------------------------------------------------------- @@ -2403,6 +2404,28 @@ Conversion::createGeographicGeocentric(const crs::CRSNNPtr &sourceCRS, // --------------------------------------------------------------------------- +/** \brief Instantiate a conversion between a GeographicCRS and a spherical + * planetocentric GeodeticCRS + * + * This method peforms conversion between geodetic latitude and geocentric + * latitude + * + * @return a new Conversion. + */ +ConversionNNPtr +Conversion::createGeographicGeocentricLatitude(const crs::CRSNNPtr &sourceCRS, + const crs::CRSNNPtr &targetCRS) { + auto properties = util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + buildOpName("Conversion", sourceCRS, targetCRS)); + auto conv = create( + properties, PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE, {}); + conv->setCRSs(sourceCRS, targetCRS, nullptr); + return conv; +} + +// --------------------------------------------------------------------------- + InverseConversion::InverseConversion(const ConversionNNPtr &forward) : Conversion( OperationMethod::create(createPropertiesForInverse(forward->method()), @@ -2509,6 +2532,15 @@ CoordinateOperationNNPtr Conversion::inverse() const { return conv; } + if (method()->nameStr() == + PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE) { + auto conv = + create(createPropertiesForInverse(this, false, false), + PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE, {}); + conv->setCRSs(this, true); + return conv; + } + return InverseConversion::create(NN_NO_CHECK( util::nn_dynamic_pointer_cast(shared_from_this()))); } @@ -3447,6 +3479,77 @@ void Conversion::_exportToPROJString( auto l_sourceCRS = sourceCRS(); auto l_targetCRS = targetCRS(); + if (methodName == PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE) { + + const auto extractGeodeticCRSIfGeodeticCRSOrEquivalent = + [](const crs::CRSPtr &crs) { + auto geodCRS = std::dynamic_pointer_cast(crs); + if (!geodCRS) { + auto compoundCRS = + std::dynamic_pointer_cast(crs); + if (compoundCRS) { + const auto &components = + compoundCRS->componentReferenceSystems(); + if (!components.empty()) { + geodCRS = + util::nn_dynamic_pointer_cast( + components[0]); + if (!geodCRS) { + auto boundCRS = util::nn_dynamic_pointer_cast< + crs::BoundCRS>(components[0]); + if (boundCRS) { + geodCRS = util::nn_dynamic_pointer_cast< + crs::GeodeticCRS>(boundCRS->baseCRS()); + } + } + } + } else { + auto boundCRS = + std::dynamic_pointer_cast(crs); + if (boundCRS) { + geodCRS = + util::nn_dynamic_pointer_cast( + boundCRS->baseCRS()); + } + } + } + return geodCRS; + }; + + auto sourceCRSGeod = dynamic_cast( + extractGeodeticCRSIfGeodeticCRSOrEquivalent(l_sourceCRS).get()); + auto targetCRSGeod = dynamic_cast( + extractGeodeticCRSIfGeodeticCRSOrEquivalent(l_targetCRS).get()); + if (sourceCRSGeod && targetCRSGeod) { + auto sourceCRSGeog = + dynamic_cast(sourceCRSGeod); + auto targetCRSGeog = + dynamic_cast(targetCRSGeod); + bool isSrcGeocentricLat = + sourceCRSGeod->isSphericalPlanetocentric(); + bool isSrcGeographic = sourceCRSGeog != nullptr; + bool isTargetGeocentricLat = + targetCRSGeod->isSphericalPlanetocentric(); + bool isTargetGeographic = targetCRSGeog != nullptr; + if ((isSrcGeocentricLat && isTargetGeographic) || + (isSrcGeographic && isTargetGeocentricLat)) { + + formatter->startInversion(); + sourceCRSGeod->_exportToPROJString(formatter); + formatter->stopInversion(); + + targetCRSGeod->_exportToPROJString(formatter); + + return; + } + } + + throw io::FormattingException("Invalid nature of source and/or " + "targetCRS for Geographic latitude / " + "Geocentric latitude" + "conversion"); + } + crs::GeographicCRSPtr srcGeogCRS; if (!formatter->getCRSExport() && l_sourceCRS && applySourceCRSModifiers) { diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index 60365713a3..1b1cae9b19 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -2993,13 +2993,6 @@ CoordinateOperationFactory::Private::createOperations( } } - // Special case if both CRS are geodetic - if (geodSrc && geodDst && !derivedSrc && !derivedDst) { - createOperationsGeodToGeod(sourceCRS, targetCRS, context, geodSrc, - geodDst, res); - return res; - } - if (geodSrc && geodSrc->isSphericalPlanetocentric()) { createOperationsFromSphericalPlanetocentric(sourceCRS, targetCRS, context, geodSrc, res); @@ -3008,6 +3001,13 @@ CoordinateOperationFactory::Private::createOperations( return applyInverse(createOperations(targetCRS, sourceCRS, context)); } + // Special case if both CRS are geodetic + if (geodSrc && geodDst && !derivedSrc && !derivedDst) { + createOperationsGeodToGeod(sourceCRS, targetCRS, context, geodSrc, + geodDst, res); + return res; + } + if (boundSrc) { auto geodSrcBase = util::nn_dynamic_pointer_cast( boundSrc->baseCRS()); @@ -3940,6 +3940,24 @@ void CoordinateOperationFactory::Private:: ENTER_FUNCTION(); + const auto IsSameDatum = [&context, + &geodSrc](const crs::GeodeticCRS *geodDst) { + const auto &authFactory = context.context->getAuthorityFactory(); + const auto dbContext = + authFactory ? authFactory->databaseContext().as_nullable() + : nullptr; + + return geodSrc->datumNonNull(dbContext)->_isEquivalentTo( + geodDst->datumNonNull(dbContext).get(), + util::IComparable::Criterion::EQUIVALENT); + }; + auto geogDst = dynamic_cast(targetCRS.get()); + if (geogDst && IsSameDatum(geogDst)) { + res.emplace_back(Conversion::createGeographicGeocentricLatitude( + sourceCRS, targetCRS)); + return; + } + // Create an intermediate geographic CRS with the same datum as the // source spherical planetocentric one std::string interm_crs_name(geodSrc->nameStr()); @@ -3953,7 +3971,8 @@ void CoordinateOperationFactory::Private:: cs::EllipsoidalCS::createLatitudeLongitude( common::UnitOfMeasure::DEGREE))); - auto opFirst = createGeodToGeodPROJBased(sourceCRS, interm_crs); + auto opFirst = + Conversion::createGeographicGeocentricLatitude(sourceCRS, interm_crs); auto opsSecond = createOperations(interm_crs, targetCRS, context); for (const auto &opSecond : opsSecond) { try { @@ -3999,7 +4018,8 @@ void CoordinateOperationFactory::Private:: auto intermBoundCRS = crs::BoundCRS::create(intermGeog, boundSrc->hubCRS(), transf); - auto opFirst = createGeodToGeodPROJBased(geodSrcBase, intermGeog); + auto opFirst = + Conversion::createGeographicGeocentricLatitude(geodSrcBase, intermGeog); setCRSs(opFirst.get(), sourceCRS, intermBoundCRS); auto opsSecond = createOperations(intermBoundCRS, targetCRS, context); for (const auto &opSecond : opsSecond) { diff --git a/src/iso19111/operation/parammappings.cpp b/src/iso19111/operation/parammappings.cpp index 9acc9e93a2..df5ae7c5c5 100644 --- a/src/iso19111/operation/parammappings.cpp +++ b/src/iso19111/operation/parammappings.cpp @@ -1340,6 +1340,8 @@ static const MethodMapping otherMethodMappings[] = { {EPSG_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC, EPSG_CODE_METHOD_GEOGRAPHIC_GEOCENTRIC, nullptr, nullptr, nullptr, nullptr}, + {PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE, 0, nullptr, nullptr, + nullptr, nullptr}, {EPSG_NAME_METHOD_LONGITUDE_ROTATION, EPSG_CODE_METHOD_LONGITUDE_ROTATION, nullptr, nullptr, nullptr, paramsLongitudeRotation}, {EPSG_NAME_METHOD_AFFINE_PARAMETRIC_TRANSFORMATION, diff --git a/src/proj_constants.h b/src/proj_constants.h index f67a0583d4..5985f1b1ce 100644 --- a/src/proj_constants.h +++ b/src/proj_constants.h @@ -693,5 +693,9 @@ #define EPSG_NAME_METHOD_GEOGRAPHIC_TOPOCENTRIC "Geographic/topocentric conversions" #define EPSG_CODE_METHOD_GEOGRAPHIC_TOPOCENTRIC 9837 +/* ------------------------------------------------------------------------ */ + +#define PROJ_WKT2_NAME_METHOD_GEOGRAPHIC_GEOCENTRIC_LATITUDE \ + "Geographic latitude / Geocentric latitude" #endif /* PROJ_CONSTANTS_INCLUDED */ diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index ec7d870014..b50542c616 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -6873,6 +6873,29 @@ TEST(operation, createOperation_geographic_to_spherical_ocentric) { // --------------------------------------------------------------------------- +TEST(operation, createOperation_spherical_ocentric_to_geocentric) { + auto objSrc = PROJStringParser().createFromPROJString( + "+proj=longlat +geoc +datum=WGS84 +type=crs"); + auto src = nn_dynamic_pointer_cast(objSrc); + ASSERT_TRUE(src != nullptr); + + auto objDest = PROJStringParser().createFromPROJString( + "+proj=geocent +datum=WGS84 +type=crs"); + auto dest = nn_dynamic_pointer_cast(objDest); + ASSERT_TRUE(dest != nullptr); + + auto op = CoordinateOperationFactory::create()->createOperation( + NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +inv +proj=geoc +ellps=WGS84 " + "+step +proj=cart +ellps=WGS84"); +} + +// --------------------------------------------------------------------------- + TEST(operation, createOperation_bound_of_spherical_ocentric_to_same_type) { auto objSrc = PROJStringParser().createFromPROJString( "+proj=longlat +geoc +ellps=GRS80 +towgs84=1,2,3 +type=crs");