ProjectedCRS::identify(): do not return CRS whose ellipsoid is totally different from the input one (#4635)

Author: rouault

src/iso19111/crs.cpp

@@ -5039,6 +5039,8 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
             }
             hasNonMatchingId = true;
         } else if (!insignificantName) {
+            const double semiMajorAxis =
+                ellipsoid->semiMajorAxis().getSIValue();
             for (int ipass = 0; ipass < 2; ipass++) {
                 const bool approximateMatch = ipass == 1;
                 auto objects = authorityFactory->createObjectsFromNameEx(
@@ -5053,8 +5055,14 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                         thisName.c_str(), pairObjName.second.c_str());
                     foundEquivalentName |= eqName;
 
-                    if (addCRS(crsNN, eqName, false).second == 100) {
-                        return res;
+                    if (std::fabs(semiMajorAxis - crsNN->baseCRS()
+                                                      ->ellipsoid()
+                                                      ->semiMajorAxis()
+                                                      .getSIValue()) <=
+                        1e-4 * semiMajorAxis) {
+                        if (addCRS(crsNN, eqName, false).second == 100) {
+                            return res;
+                        }
                     }
                 }
                 if (!res.empty()) {

src/iso19111/factory.cpp

@@ -10008,13 +10008,27 @@ AuthorityFactory::createProjectedCRSFromExisting(
         }
     }
 
-    std::string sql("SELECT projected_crs.auth_name, projected_crs.code, "
-                    "projected_crs.name FROM projected_crs "
-                    "JOIN conversion_table conv ON "
-                    "projected_crs.conversion_auth_name = conv.auth_name AND "
-                    "projected_crs.conversion_code = conv.code WHERE "
-                    "projected_crs.deprecated = 0 AND ");
+    std::string sql(
+        "SELECT projected_crs.auth_name, projected_crs.code, "
+        "projected_crs.name FROM projected_crs "
+        "JOIN conversion_table conv ON "
+        "projected_crs.conversion_auth_name = conv.auth_name AND "
+        "projected_crs.conversion_code = conv.code "
+        "JOIN geodetic_crs gcrs ON "
+        "gcrs.auth_name = projected_crs.geodetic_crs_auth_name AND "
+        "gcrs.code = projected_crs.geodetic_crs_code "
+        "JOIN geodetic_datum datum ON "
+        "datum.auth_name = gcrs.datum_auth_name AND "
+        "datum.code = gcrs.datum_code "
+        "JOIN ellipsoid ellps ON "
+        "ellps.auth_name = datum.ellipsoid_auth_name AND "
+        "ellps.code = datum.ellipsoid_code "
+        "WHERE "
+        "abs(ellps.semi_major_axis - ?) <= 1e-4 * ellps.semi_major_axis AND "
+        "projected_crs.deprecated = 0 AND ");
     ListOfParams params;
+    params.emplace_back(
+        toString(crs->baseCRS()->ellipsoid()->semiMajorAxis().value()));
     if (!candidatesGeodCRS.empty()) {
         sql += buildSqlLookForAuthNameCode(candidatesGeodCRS, params,
                                            "projected_crs.geodetic_crs_");

test/unit/test_crs.cpp

@@ -3765,6 +3765,51 @@ TEST(crs, projectedCRS_identify_db) {
         EXPECT_EQ(*(res.front().first->identifiers()[0]->codeSpace()), "EPSG");
         EXPECT_EQ(res.front().second, 100);
     }
+    {
+        // The ellipsoid definition uses a unusual value for the inverse
+        // flattening Check that we only identify CRS whose ellipsoid shape is
+        // close to the input one (in particular no IAU_2015...)
+        auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(
+            "PROJCRS[\"Transverse Mercator\",\n"
+            "    BASEGEOGCRS[\"DE_DHDN (whole country, 2001) to ETRS89\",\n"
+            "        DATUM[\"DE_DHDN (whole country, 2001) to ETRS89\",\n"
+            "            ELLIPSOID[\"Bessel\",6377397.155,299.15281310608,\n"
+            "                LENGTHUNIT[\"metre\",1,\n"
+            "                    ID[\"EPSG\",9001]]]],\n"
+            "        PRIMEM[\"Greenwich\",0,\n"
+            "            ANGLEUNIT[\"degree\",0.0174532925199433,\n"
+            "                ID[\"EPSG\",9122]]]],\n"
+            "    CONVERSION[\"Transverse Mercator\",\n"
+            "        METHOD[\"Transverse Mercator\",\n"
+            "            ID[\"EPSG\",9807]],\n"
+            "        PARAMETER[\"Latitude of natural origin\",0,\n"
+            "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+            "            ID[\"EPSG\",8801]],\n"
+            "        PARAMETER[\"Longitude of natural origin\",12,\n"
+            "            ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+            "            ID[\"EPSG\",8802]],\n"
+            "        PARAMETER[\"Scale factor at natural origin\",1,\n"
+            "            SCALEUNIT[\"unity\",1],\n"
+            "            ID[\"EPSG\",8805]],\n"
+            "        PARAMETER[\"False easting\",4500000,\n"
+            "            LENGTHUNIT[\"metre\",1],\n"
+            "            ID[\"EPSG\",8806]],\n"
+            "        PARAMETER[\"False northing\",0,\n"
+            "            LENGTHUNIT[\"metre\",1],\n"
+            "            ID[\"EPSG\",8807]]],\n"
+            "    CS[Cartesian,2],\n"
+            "        AXIS[\"easting\",east,\n"
+            "            ORDER[1],\n"
+            "            LENGTHUNIT[\"meters\",1]],\n"
+            "        AXIS[\"northing\",north,\n"
+            "            ORDER[2],\n"
+            "            LENGTHUNIT[\"meters\",1]]]");
+        auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        auto factoryAll = AuthorityFactory::create(dbContext, std::string());
+        auto res = crs->identify(factoryAll);
+        EXPECT_EQ(res.size(), 13U);
+    }
 }
 
 // ---------------------------------------------------------------------------