From 35a343fb3069d1b4c7474ccf7b32d742db82d3f7 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 24 May 2024 04:11:48 +0200 Subject: [PATCH 1/2] WKT/PROJJSON importer: improve heuristics to guess the celestial body name by using the ellipsoid name (fixes #4148) --- include/proj/datum.hpp | 3 ++- include/proj/internal/datum_internal.hpp | 7 ++++++ src/iso19111/datum.cpp | 27 +++++++++++++++++------- src/iso19111/factory.cpp | 18 ++++++++++++---- src/iso19111/io.cpp | 25 +++++++++++++--------- test/unit/test_io.cpp | 26 +++++++++++++++++++++++ 6 files changed, 83 insertions(+), 23 deletions(-) diff --git a/include/proj/datum.hpp b/include/proj/datum.hpp index 5e81330f67..3fdc4fcdd8 100644 --- a/include/proj/datum.hpp +++ b/include/proj/datum.hpp @@ -338,7 +338,8 @@ class PROJ_GCC_DLL Ellipsoid final : public common::IdentifiedObject, //! @endcond PROJ_INTERNAL static std::string - guessBodyName(const io::DatabaseContextPtr &dbContext, double a); + guessBodyName(const io::DatabaseContextPtr &dbContext, double a, + const std::string &ellpsName = std::string()); PROJ_INTERNAL bool lookForProjWellKnownEllps(std::string &projEllpsName, std::string &ellpsName) const; diff --git a/include/proj/internal/datum_internal.hpp b/include/proj/internal/datum_internal.hpp index 899679f669..8cee02547f 100644 --- a/include/proj/internal/datum_internal.hpp +++ b/include/proj/internal/datum_internal.hpp @@ -35,4 +35,11 @@ #define UNKNOWN_ENGINEERING_DATUM "Unknown engineering datum" +constexpr const char *NON_EARTH_BODY = "Non-Earth body"; + +// Mars (2015) - Sphere uses R=3396190 +// and Mars polar radius (as used by HIRISE JPEG2000) is 3376200m +// which is a 0.59% relative difference. +constexpr double REL_ERROR_FOR_SAME_CELESTIAL_BODY = 0.007; + #endif // DATUM_INTERNAL_HH_INCLUDED diff --git a/src/iso19111/datum.cpp b/src/iso19111/datum.cpp index 7a1fb612fb..a8a246bcc9 100644 --- a/src/iso19111/datum.cpp +++ b/src/iso19111/datum.cpp @@ -1179,24 +1179,35 @@ bool Ellipsoid::_isEquivalentTo(const util::IComparable *other, // --------------------------------------------------------------------------- std::string Ellipsoid::guessBodyName(const io::DatabaseContextPtr &dbContext, - double a) { - // Mars (2015) - Sphere uses R=3396190 - // and Mars polar radius (as used by HIRISE JPEG2000) is 3376200m - // which is a 0.59% relative difference. - constexpr double relError = 0.007; + double a, const std::string &ellpsName) { constexpr double earthMeanRadius = 6375000.0; - if (std::fabs(a - earthMeanRadius) < relError * earthMeanRadius) { + if (std::fabs(a - earthMeanRadius) < + REL_ERROR_FOR_SAME_CELESTIAL_BODY * earthMeanRadius) { return Ellipsoid::EARTH; } if (dbContext) { try { auto factory = io::AuthorityFactory::create(NN_NO_CHECK(dbContext), std::string()); - return factory->identifyBodyFromSemiMajorAxis(a, relError); + if (!ellpsName.empty()) { + auto matches = factory->createObjectsFromName( + ellpsName, {io::AuthorityFactory::ObjectType::ELLIPSOID}, + true, 1); + if (!matches.empty()) { + auto ellps = + static_cast(matches.front().get()); + if (std::fabs(a - ellps->semiMajorAxis().getSIValue()) < + REL_ERROR_FOR_SAME_CELESTIAL_BODY * a) { + return ellps->celestialBody(); + } + } + } + return factory->identifyBodyFromSemiMajorAxis( + a, REL_ERROR_FOR_SAME_CELESTIAL_BODY); } catch (const std::exception &) { } } - return "Non-Earth body"; + return NON_EARTH_BODY; } // --------------------------------------------------------------------------- diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index be211e809a..c8ea85a6d5 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -4524,20 +4524,30 @@ std::string AuthorityFactory::identifyBodyFromSemiMajorAxis(double semi_major_axis, double tolerance) const { auto res = - d->run("SELECT name, (ABS(semi_major_axis - ?) / semi_major_axis ) " - "AS rel_error FROM celestial_body WHERE rel_error <= ?", + d->run("SELECT DISTINCT name, " + "(ABS(semi_major_axis - ?) / semi_major_axis ) AS rel_error " + "FROM celestial_body WHERE rel_error <= ? " + "ORDER BY rel_error, name", {semi_major_axis, tolerance}); if (res.empty()) { throw FactoryException("no match found"); } + constexpr int IDX_NAME = 0; if (res.size() > 1) { + constexpr int IDX_REL_ERROR = 1; + // If the first object has a relative error of 0 and the next one + // a non-zero error, then use the first one. + if (res.front()[IDX_REL_ERROR] == "0" && + (*std::next(res.begin()))[IDX_REL_ERROR] != "0") { + return res.front()[IDX_NAME]; + } for (const auto &row : res) { - if (row[0] != res.front()[0]) { + if (row[IDX_NAME] != res.front()[IDX_NAME]) { throw FactoryException("more than one match found"); } } } - return res.front()[0]; + return res.front()[IDX_NAME]; } // --------------------------------------------------------------------------- diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 67d7340733..a65ce729f0 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -2083,15 +2083,17 @@ EllipsoidNNPtr WKTParser::Private::buildEllipsoid(const WKTNodeNNPtr &node) { Scale invFlattening(invFlatteningChild->GP()->value() == "\"inf\"" ? 0 : asDouble(invFlatteningChild)); - const auto celestialBody( - Ellipsoid::guessBodyName(dbContext_, semiMajorAxis.getSIValue())); + const auto ellpsProperties = buildProperties(node); + std::string ellpsName; + ellpsProperties.getStringValue(IdentifiedObject::NAME_KEY, ellpsName); + const auto celestialBody(Ellipsoid::guessBodyName( + dbContext_, semiMajorAxis.getSIValue(), ellpsName)); if (invFlattening.getSIValue() == 0) { - return Ellipsoid::createSphere(buildProperties(node), semiMajorAxis, + return Ellipsoid::createSphere(ellpsProperties, semiMajorAxis, celestialBody); } else { return Ellipsoid::createFlattenedSphere( - buildProperties(node), semiMajorAxis, invFlattening, - celestialBody); + ellpsProperties, semiMajorAxis, invFlattening, celestialBody); } } catch (const std::exception &e) { throw buildRethrow(__FUNCTION__, e); @@ -7122,15 +7124,18 @@ PrimeMeridianNNPtr JSONParser::buildPrimeMeridian(const json &j) { EllipsoidNNPtr JSONParser::buildEllipsoid(const json &j) { if (j.contains("semi_major_axis")) { auto semiMajorAxis = getLength(j, "semi_major_axis"); - const auto celestialBody( - Ellipsoid::guessBodyName(dbContext_, semiMajorAxis.getSIValue())); + const auto ellpsProperties = buildProperties(j); + std::string ellpsName; + ellpsProperties.getStringValue(IdentifiedObject::NAME_KEY, ellpsName); + const auto celestialBody(Ellipsoid::guessBodyName( + dbContext_, semiMajorAxis.getSIValue(), ellpsName)); if (j.contains("semi_minor_axis")) { - return Ellipsoid::createTwoAxis(buildProperties(j), semiMajorAxis, + return Ellipsoid::createTwoAxis(ellpsProperties, semiMajorAxis, getLength(j, "semi_minor_axis"), celestialBody); } else if (j.contains("inverse_flattening")) { return Ellipsoid::createFlattenedSphere( - buildProperties(j), semiMajorAxis, + ellpsProperties, semiMajorAxis, Scale(getNumber(j, "inverse_flattening")), celestialBody); } else { throw ParsingException( @@ -10717,7 +10722,7 @@ PrimeMeridianNNPtr PROJStringParser::Private::buildPrimeMeridian(Step &step) { std::string PROJStringParser::Private::guessBodyName(double a) { auto ret = Ellipsoid::guessBodyName(dbContext_, a); - if (ret == "Non-Earth body" && dbContext_ == nullptr && ctx_ != nullptr) { + if (ret == NON_EARTH_BODY && dbContext_ == nullptr && ctx_ != nullptr) { dbContext_ = ctx_->get_cpp_context()->getDatabaseContext().as_nullable(); if (dbContext_) { diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 984aef93a4..4ccaff1772 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -266,6 +266,32 @@ TEST(wkt_parse, datum_no_pm_not_earth) { // --------------------------------------------------------------------------- +TEST(wkt_parse, guess_celestial_body_from_ellipsoid_name) { + auto obj = WKTParser() + .attachDatabaseContext(DatabaseContext::create()) + .createFromWKT("DATUM[\"unnamed\",\n" + " ELLIPSOID[\"Ananke\",10000,0,\n" + " LENGTHUNIT[\"metre\",1]]]"); + auto datum = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(datum != nullptr); + EXPECT_EQ(datum->ellipsoid()->celestialBody(), "Ananke"); +} + +// --------------------------------------------------------------------------- + +TEST(wkt_parse, guess_celestial_body_from_ellipsoid_name_false_positive) { + auto obj = WKTParser() + .attachDatabaseContext(DatabaseContext::create()) + .createFromWKT("DATUM[\"unnamed\",\n" + " ELLIPSOID[\"Ananke\",999999,0,\n" + " LENGTHUNIT[\"metre\",1]]]"); + auto datum = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(datum != nullptr); + EXPECT_EQ(datum->ellipsoid()->celestialBody(), "Non-Earth body"); +} + +// --------------------------------------------------------------------------- + TEST(wkt_parse, dynamic_geodetic_reference_frame) { auto obj = WKTParser().createFromWKT( "GEOGCRS[\"WGS 84 (G1762)\"," From 38f218da43ba9c928cdc77adac042b558eced5f5 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 24 May 2024 04:12:23 +0200 Subject: [PATCH 2/2] createOperations(): improve detection of compatible/incompatible celestial bodies (fixes #4148) --- .../operation/coordinateoperationfactory.cpp | 39 ++++--- test/unit/test_operationfactory.cpp | 109 ++++++++++++++++++ 2 files changed, 132 insertions(+), 16 deletions(-) diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index 9d37b60506..da6dd3fcf6 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -38,6 +38,7 @@ #include "proj/metadata.hpp" #include "proj/util.hpp" +#include "proj/internal/datum_internal.hpp" #include "proj/internal/internal.hpp" #include "proj/internal/io_internal.hpp" #include "proj/internal/tracing.hpp" @@ -4442,28 +4443,34 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod( ENTER_FUNCTION(); - if (geodSrc->ellipsoid()->celestialBody() != - geodDst->ellipsoid()->celestialBody()) { + const auto &srcEllps = geodSrc->ellipsoid(); + const auto &dstEllps = geodDst->ellipsoid(); + if (srcEllps->celestialBody() == dstEllps->celestialBody() && + srcEllps->celestialBody() != NON_EARTH_BODY) { + // Same celestial body, with a known name (that is not the generic + // NON_EARTH_BODY used when we can't guess it) ==> compatible + } else if ((srcEllps->celestialBody() != dstEllps->celestialBody() && + srcEllps->celestialBody() != NON_EARTH_BODY && + dstEllps->celestialBody() != NON_EARTH_BODY) || + std::fabs(srcEllps->semiMajorAxis().getSIValue() - + dstEllps->semiMajorAxis().getSIValue()) > + REL_ERROR_FOR_SAME_CELESTIAL_BODY * + dstEllps->semiMajorAxis().getSIValue()) { const char *envVarVal = getenv("PROJ_IGNORE_CELESTIAL_BODY"); - if (envVarVal == nullptr) { + if (envVarVal == nullptr || ci_equal(envVarVal, "NO") || + ci_equal(envVarVal, "FALSE") || ci_equal(envVarVal, "OFF")) { std::string osMsg( "Source and target ellipsoid do not belong to the same " "celestial body ("); - osMsg += geodSrc->ellipsoid()->celestialBody(); + osMsg += srcEllps->celestialBody(); osMsg += " vs "; - osMsg += geodDst->ellipsoid()->celestialBody(); - osMsg += "). You may override this check by setting the " - "PROJ_IGNORE_CELESTIAL_BODY environment variable to YES."; - throw util::UnsupportedOperationException(osMsg.c_str()); - } else if (ci_equal(envVarVal, "NO") || ci_equal(envVarVal, "FALSE") || - ci_equal(envVarVal, "OFF")) { - std::string osMsg( - "Source and target ellipsoid do not belong to the same " - "celestial body ("); - osMsg += geodSrc->ellipsoid()->celestialBody(); - osMsg += " vs "; - osMsg += geodDst->ellipsoid()->celestialBody(); + osMsg += dstEllps->celestialBody(); osMsg += ")."; + if (envVarVal == nullptr) { + osMsg += " You may override this check by setting the " + "PROJ_IGNORE_CELESTIAL_BODY environment variable " + "to YES."; + } throw util::UnsupportedOperationException(osMsg.c_str()); } } diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 976e8e7ddd..d3548a3ca9 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -637,6 +637,115 @@ TEST(operation, geogCRS_to_geogCRS_context_datum_ensemble) { // --------------------------------------------------------------------------- +TEST(operation, geogCRS_to_geogCRS_context_incompatible_celestial_body) { + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), std::string()); + auto authFactoryEPSG = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto authFactoryIAU_2015 = + AuthorityFactory::create(DatabaseContext::create(), "IAU_2015"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + EXPECT_THROW( + CoordinateOperationFactory::create()->createOperations( + authFactoryEPSG->createCoordinateReferenceSystem("4326"), // WGS 84 + authFactoryIAU_2015->createCoordinateReferenceSystem( + "51200"), // Ananke + ctxt), + UnsupportedOperationException); +} + +// --------------------------------------------------------------------------- + +TEST(operation, + geogCRS_to_geogCRS_context_incompatible_celestial_body_but_same_radius) { + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), std::string()); + auto authFactoryESRI = + AuthorityFactory::create(DatabaseContext::create(), "ESRI"); + auto authFactoryIAU_2015 = + AuthorityFactory::create(DatabaseContext::create(), "IAU_2015"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + EXPECT_THROW(CoordinateOperationFactory::create()->createOperations( + authFactoryESRI->createCoordinateReferenceSystem( + "104936"), // GCS_Pan_2000 + authFactoryIAU_2015->createCoordinateReferenceSystem( + "51200"), // Ananke + ctxt), + UnsupportedOperationException); +} + +// --------------------------------------------------------------------------- + +TEST(operation, geogCRS_to_geogCRS_compatible_unknown_celestial_body) { + auto dbContext = DatabaseContext::create(); + auto authFactory = AuthorityFactory::create(dbContext, std::string()); + auto objSrc = + createFromUserInput("+proj=longlat +R=10000 +type=crs", dbContext); + auto srcCRS = nn_dynamic_pointer_cast(objSrc); + auto objDst = + createFromUserInput("+proj=longlat +R=10001 +type=crs", dbContext); + auto dstCRS = nn_dynamic_pointer_cast(objDst); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(srcCRS), NN_NO_CHECK(dstCRS), ctxt); + ASSERT_EQ(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=noop"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, geogCRS_to_geogCRS_incompatible_unknown_celestial_body) { + auto dbContext = DatabaseContext::create(); + auto authFactory = AuthorityFactory::create(dbContext, std::string()); + auto objSrc = + createFromUserInput("+proj=longlat +R=10000 +type=crs", dbContext); + auto srcCRS = nn_dynamic_pointer_cast(objSrc); + auto objDst = + createFromUserInput("+proj=longlat +R=99999 +type=crs", dbContext); + auto dstCRS = nn_dynamic_pointer_cast(objDst); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + EXPECT_THROW(CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(srcCRS), NN_NO_CHECK(dstCRS), ctxt), + UnsupportedOperationException); +} + +// --------------------------------------------------------------------------- + +TEST(operation, + geogCRS_to_geogCRS_compatible_celestial_body_through_semi_major_axis) { + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), std::string()); + auto authFactoryIAU_2015 = + AuthorityFactory::create(DatabaseContext::create(), "IAU_2015"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + + auto dst_wkt = "GEOGCRS[\"unknown\",\n" + " DATUM[\"unknown\",\n" + " ELLIPSOID[\"unknown\",9999,0,\n" // Ananke is 10000 + " LENGTHUNIT[\"metre\",1]]],\n" + " PRIMEM[\"Reference_Meridian\",0,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " CS[ellipsoidal,2],\n" + " AXIS[\"geodetic latitude (Lat)\",north,\n" + " ORDER[1],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " AXIS[\"geodetic longitude (Lon)\",east,\n" + " ORDER[2],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]]]"; + auto dstObj = WKTParser().createFromWKT(dst_wkt); + auto dstCRS = nn_dynamic_pointer_cast(dstObj); + + auto list = CoordinateOperationFactory::create()->createOperations( + authFactoryIAU_2015->createCoordinateReferenceSystem("51200"), // Ananke + NN_NO_CHECK(dstCRS), ctxt); + ASSERT_EQ(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=noop"); +} + +// --------------------------------------------------------------------------- + TEST(operation, geogCRS_to_derived_geogCRS_3D) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");