WKT/PROJJSON importer: improve heuristics to guess the celestial body…

… name by using the ellipsoid name (fixes OSGeo#4148)

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;

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

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<const Ellipsoid *>(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;
 }
 
 // ---------------------------------------------------------------------------

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];
 }
 
 // ---------------------------------------------------------------------------

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_) {

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<GeodeticReferenceFrame>(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<GeodeticReferenceFrame>(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)\","