ProjectedCRS::identify(): fix identification of EPSG:3059 from 'LKS92_Latvia_TM' name, and other cleanups/improvements (fixes #2214)

3 files changed, 118 insertions, 72 deletions

diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 6b726f34..6f20cf12 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -2236,8 +2236,8 @@ PJ *proj_get_target_crs(PJ_CONTEXT *ctx, const PJ *obj) {
  * </li>
  * <li>90% means that CRS are equivalent, but the names are not exactly the
  * same.</li>
- * <li>70% means that CRS are equivalent), but the names do not match at
- * all.</li>
+ * <li>70% means that CRS are equivalent, but the names are not equivalent.
+ * </li>
  * <li>25% means that the CRS are not equivalent, but there is some similarity
  * in
  * the names.</li>
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index de02eb35..c23bd29b 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -1703,7 +1703,7 @@ static bool hasCodeCompatibleOfAuthorityFactory(
  * same.
  * <li>70% means that CRS are equivalent (equivalent datum and coordinate
  * system),
- * but the names do not match at all.</li>
+ * but the names are not equivalent.</li>
  * <li>60% means that ellipsoid, prime meridian and coordinate systems are
  * equivalent, but the CRS and datum names do not match.</li>
  * <li>25% means that the CRS are not equivalent, but there is some similarity
@@ -2808,7 +2808,7 @@ bool VerticalCRS::_isEquivalentTo(
  * single result is returned.
  * 90% means that CRS are equivalent, but the names are not exactly the same.
  * 70% means that CRS are equivalent (equivalent datum and coordinate system),
- * but the names do not match at all.
+ * but the names are not equivalent.
  * 25% means that the CRS are not equivalent, but there is some similarity in
  * the names.
  *
@@ -3630,8 +3630,11 @@ void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter,
  * single result is returned.
  * 90% means that CRS are equivalent, but the names are not exactly the same.
  * 70% means that CRS are equivalent (equivalent base CRS, conversion and
- * coordinate system), but the names do not match at all.
- * 50% means that CRS have similarity (equivalent base CRS and conversion),
+ * coordinate system), but the names are not equivalent.
+ * 60% means that CRS have strong similarity (equivalent base datum, conversion
+ * and coordinate system), but the names are not equivalent.
+ * 50% means that CRS have similarity (equivalent base ellipsoid and
+ * conversion),
  * but the coordinate system do not match (e.g. different axis ordering or
  * axis unit).
  * 25% means that the CRS are not equivalent, but there is some similarity in
@@ -3655,6 +3658,11 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
 
     std::list<std::pair<GeodeticCRSNNPtr, int>> baseRes;
     const auto &l_baseCRS(baseCRS());
+    const auto l_datum = l_baseCRS->datum();
+    const bool significantNameForDatum =
+        l_datum != nullptr && !ci_starts_with(l_datum->nameStr(), "unknown") &&
+        l_datum->nameStr() != "unnamed";
+    const auto &ellipsoid = l_baseCRS->ellipsoid();
     auto geogCRS = dynamic_cast<const GeographicCRS *>(l_baseCRS.get());
     if (geogCRS &&
         geogCRS->coordinateSystem()->axisOrder() ==
@@ -3742,6 +3750,59 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
         }
     }
 
+    const bool l_implicitCS = CRS::getPrivate()->implicitCS_;
+    const auto addCRS = [&](const ProjectedCRSNNPtr &crs, const bool eqName) {
+        const auto &l_unit = cs->axisList()[0]->unit();
+        if (_isEquivalentTo(crs.get(), util::IComparable::Criterion::
+                                           EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS,
+                            dbContext) ||
+            (l_implicitCS &&
+             l_unit._isEquivalentTo(
+                 crs->coordinateSystem()->axisList()[0]->unit(),
+                 util::IComparable::Criterion::EQUIVALENT) &&
+             l_baseCRS->_isEquivalentTo(
+                 crs->baseCRS().get(), util::IComparable::Criterion::
+                                           EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS,
+                 dbContext) &&
+             derivingConversionRef()->_isEquivalentTo(
+                 crs->derivingConversionRef().get(),
+                 util::IComparable::Criterion::EQUIVALENT, dbContext))) {
+            if (crs->nameStr() == thisName) {
+                res.clear();
+                res.emplace_back(crs, 100);
+            } else {
+                res.emplace_back(crs, eqName ? 90 : 70);
+            }
+        } else if (ellipsoid->_isEquivalentTo(
+                       crs->baseCRS()->ellipsoid().get(),
+                       util::IComparable::Criterion::EQUIVALENT, dbContext) &&
+                   derivingConversionRef()->_isEquivalentTo(
+                       crs->derivingConversionRef().get(),
+                       util::IComparable::Criterion::EQUIVALENT, dbContext)) {
+            if ((l_implicitCS &&
+                 l_unit._isEquivalentTo(
+                     crs->coordinateSystem()->axisList()[0]->unit(),
+                     util::IComparable::Criterion::EQUIVALENT)) ||
+                cs->_isEquivalentTo(crs->coordinateSystem().get(),
+                                    util::IComparable::Criterion::EQUIVALENT,
+                                    dbContext)) {
+                if (!significantNameForDatum ||
+                    (crs->baseCRS()->datum() &&
+                     l_datum->_isEquivalentTo(
+                         crs->baseCRS()->datum().get(),
+                         util::IComparable::Criterion::EQUIVALENT))) {
+                    res.emplace_back(crs, 70);
+                } else {
+                    res.emplace_back(crs, 60);
+                }
+            } else {
+                res.emplace_back(crs, 50);
+            }
+        } else {
+            res.emplace_back(crs, 25);
+        }
+    };
+
     if (authorityFactory) {
 
         const bool unsignificantName = thisName.empty() ||
@@ -3782,43 +3843,10 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                     const bool eqName = metadata::Identifier::isEquivalentName(
                         thisName.c_str(), crs->nameStr().c_str());
                     foundEquivalentName |= eqName;
-                    if (_isEquivalentTo(
-                            crs.get(), util::IComparable::Criterion::
-                                           EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS,
-                            dbContext)) {
-                        if (crs->nameStr() == thisName) {
-                            res.clear();
-                            res.emplace_back(crsNN, 100);
-                            return res;
-                        }
-                        res.emplace_back(crsNN, eqName ? 90 : 70);
-                    } else if (objects.size() == 1 &&
-                               CRS::getPrivate()->implicitCS_ &&
-                               coordinateSystem()
-                                   ->axisList()[0]
-                                   ->unit()
-                                   ._isEquivalentTo(
-                                       crs->coordinateSystem()
-                                           ->axisList()[0]
-                                           ->unit(),
-                                       util::IComparable::Criterion::
-                                           EQUIVALENT) &&
-                               l_baseCRS->_isEquivalentTo(
-                                   crs->baseCRS().get(),
-                                   util::IComparable::Criterion::
-                                       EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS,
-                                   dbContext) &&
-                               derivingConversionRef()->_isEquivalentTo(
-                                   crs->derivingConversionRef().get(),
-                                   util::IComparable::Criterion::EQUIVALENT,
-                                   dbContext)) {
-                        res.clear();
-                        res.emplace_back(crsNN, crs->nameStr() == thisName
-                                                    ? 100
-                                                    : eqName ? 90 : 70);
+
+                    addCRS(crsNN, eqName);
+                    if (res.back().second == 100) {
                         return res;
-                    } else {
-                        res.emplace_back(crsNN, 25);
                     }
                 }
                 if (!res.empty()) {
@@ -3867,7 +3895,6 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                 shared_from_this().as_nullable()));
             auto candidates =
                 authorityFactory->createProjectedCRSFromExisting(self);
-            const auto &ellipsoid = l_baseCRS->ellipsoid();
             for (const auto &crs : candidates) {
                 const auto &ids = crs->identifiers();
                 assert(!ids.empty());
@@ -3877,30 +3904,7 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                     continue;
                 }
 
-                if (_isEquivalentTo(crs.get(),
-                                    util::IComparable::Criterion::
-                                        EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS,
-                                    dbContext)) {
-                    res.emplace_back(crs, unsignificantName ? 90 : 70);
-                } else if (ellipsoid->_isEquivalentTo(
-                               crs->baseCRS()->ellipsoid().get(),
-                               util::IComparable::Criterion::EQUIVALENT,
-                               dbContext) &&
-                           derivingConversionRef()->_isEquivalentTo(
-                               crs->derivingConversionRef().get(),
-                               util::IComparable::Criterion::EQUIVALENT,
-                               dbContext)) {
-                    if (coordinateSystem()->_isEquivalentTo(
-                            crs->coordinateSystem().get(),
-                            util::IComparable::Criterion::EQUIVALENT,
-                            dbContext)) {
-                        res.emplace_back(crs, 70);
-                    } else {
-                        res.emplace_back(crs, 50);
-                    }
-                } else {
-                    res.emplace_back(crs, 25);
-                }
+                addCRS(crs, unsignificantName);
             }
 
             res.sort(lambdaSort);
@@ -4293,7 +4297,7 @@ bool CompoundCRS::_isEquivalentTo(
  * single result is returned.
  * 90% means that CRS are equivalent, but the names are not exactly the same.
  * 70% means that CRS are equivalent (equivalent horizontal and vertical CRS),
- * but the names do not match at all.
+ * but the names are not equivalent.
  * 25% means that the CRS are not equivalent, but there is some similarity in
  * the names.
  *
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index 270e7f95..c00556a7 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -2303,7 +2303,7 @@ TEST(crs, projectedCRS_identify_db) {
         auto res = crs->identify(factoryEPSG);
         ASSERT_EQ(res.size(), 1U);
         EXPECT_EQ(res.front().first->getEPSGCode(), 3044);
-        EXPECT_EQ(res.front().second, 25);
+        EXPECT_EQ(res.front().second, 50);
     }
     {
         // Identify from a WKT1 string, without explicit axis
@@ -2397,7 +2397,7 @@ TEST(crs, projectedCRS_identify_db) {
         auto res = crs->identify(factoryEPSG);
         EXPECT_EQ(res.size(), 1U);
         EXPECT_EQ(res.front().first->getEPSGCode(), 2327);
-        EXPECT_EQ(res.front().second, 50);
+        EXPECT_EQ(res.front().second, 70);
     }
     {
         // EPSG:6646 as PROJ.4 string, using clrk80 which is pretty generic
@@ -2439,7 +2439,7 @@ TEST(crs, projectedCRS_identify_db) {
         for (const auto &pair : res) {
             if (pair.first->identifiers()[0]->code() == "102039") {
                 found = true;
-                EXPECT_EQ(pair.second, 25);
+                EXPECT_EQ(pair.second, 50);
                 break;
             }
         }
@@ -2574,6 +2574,48 @@ TEST(crs, projectedCRS_identify_db) {
         EXPECT_EQ(res.front().first->getEPSGCode(), 2154);
         EXPECT_EQ(res.front().second, 70);
     }
+    {
+        // Test identification of LKS92_Latvia_TM (#2214)
+        auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(
+            "PROJCS[\"LKS92_Latvia_TM\",GEOGCS[\"GCS_LKS92\","
+            "DATUM[\"D_Latvia_1992\","
+            "SPHEROID[\"GRS_1980\",6378137,298.257222101]],"
+            "PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]],"
+            "PROJECTION[\"Transverse_Mercator\"],"
+            "PARAMETER[\"latitude_of_origin\",0],"
+            "PARAMETER[\"central_meridian\",24],"
+            "PARAMETER[\"scale_factor\",0.9996],"
+            "PARAMETER[\"false_easting\",500000],"
+            "PARAMETER[\"false_northing\",-6000000],"
+            "UNIT[\"Meter\",1]]");
+        auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        auto res = crs->identify(factoryEPSG);
+        ASSERT_EQ(res.size(), 1U);
+        EXPECT_EQ(res.front().first->getEPSGCode(), 3059);
+        EXPECT_EQ(res.front().second, 90);
+    }
+    {
+        // Test identification of CRS where everything but datum names matches
+        auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(
+            "PROJCS[\"WGS_1984_UTM_Zone_31N\",GEOGCS[\"GCS_WGS_1984\","
+            "DATUM[\"wrong_datum_name\","
+            "SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],"
+            "PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],"
+            "PROJECTION[\"Transverse_Mercator\"],"
+            "PARAMETER[\"False_Easting\",500000.0],"
+            "PARAMETER[\"False_Northing\",0.0],"
+            "PARAMETER[\"Central_Meridian\",3.0],"
+            "PARAMETER[\"Scale_Factor\",0.9996],"
+            "PARAMETER[\"Latitude_Of_Origin\",0.0],"
+            "UNIT[\"Meter\",1.0]]");
+        auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        auto res = crs->identify(factoryEPSG);
+        ASSERT_EQ(res.size(), 1U);
+        EXPECT_EQ(res.front().first->getEPSGCode(), 32631);
+        EXPECT_EQ(res.front().second, 60);
+    }
 }
 
 // ---------------------------------------------------------------------------
@@ -4121,7 +4163,7 @@ TEST(crs, boundCRS_identify_db) {
         EXPECT_EQ(boundCRS->baseCRS()->getEPSGCode(), 2106);
         EXPECT_EQ(boundCRS->transformation()->nameStr(),
                   "NZGD2000 to WGS 84 (1)");
-        EXPECT_EQ(res.front().second, 50);
+        EXPECT_EQ(res.front().second, 70);
     }
 
     {