From 4097baf8e370922817924fa29e338d175d54f78d Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 24 May 2024 04:11:48 +0200 Subject: [PATCH] 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 2de95d9eda..f20c5f363b 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 b9dd783a01..43e957371c 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -2105,15 +2105,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); @@ -7146,15 +7148,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( @@ -10745,7 +10750,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)\","