BoundCRS::identify(): improvements to discard CRS that aren't relevant (fixes #1801)

Fix for
```
projinfo --identify "+proj=utm +zone=48 +a=6377276.345 +b=6356075.41314024 +towgs84=198,881,317,0,0,0,0 +units=m +no_defs +type=crs"
```

to only return BoundCRS of EPSG:3148: 70 %

Previously it also returned EPSG:23948 and EPSG:24048 whose projected CRS-only
parts where likely matches, but those 2 CRSs don't have a +towgs84=198,881,317,0,0,0,0,
so discard them.

3 files changed, 61 insertions, 30 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 0469cae1..2c3e38ac 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -10189,7 +10189,9 @@ bool ConcatenatedOperation::_isEquivalentTo(
     const util::IComparable *other,
     util::IComparable::Criterion criterion) const {
     auto otherCO = dynamic_cast<const ConcatenatedOperation *>(other);
-    if (otherCO == nullptr || !ObjectUsage::_isEquivalentTo(other, criterion)) {
+    if (otherCO == nullptr ||
+        (criterion == util::IComparable::Criterion::STRICT &&
+         !ObjectUsage::_isEquivalentTo(other, criterion))) {
         return false;
     }
     const auto &steps = operations();
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 46aeee20..f231d072 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -4553,10 +4553,32 @@ std::list<std::pair<CRSNNPtr, int>>
 BoundCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
     typedef std::pair<CRSNNPtr, int> Pair;
     std::list<Pair> res;
+    std::list<Pair> resMatchOfTransfToWGS84;
     if (authorityFactory &&
         d->hubCRS_->_isEquivalentTo(GeographicCRS::EPSG_4326.get(),
                                     util::IComparable::Criterion::EQUIVALENT)) {
         auto resTemp = d->baseCRS_->identify(authorityFactory);
+
+        std::string refTransfPROJString;
+        bool refTransfPROJStringValid = false;
+        auto refTransf = d->transformation_->normalizeForVisualization();
+        try {
+            refTransfPROJString = refTransf->exportToPROJString(
+                io::PROJStringFormatter::create().get());
+            refTransfPROJString = replaceAll(
+                refTransfPROJString,
+                " +rx=0 +ry=0 +rz=0 +s=0 +convention=position_vector", "");
+            refTransfPROJStringValid = true;
+        } catch (const std::exception &) {
+        }
+        bool refIsNullTransform = false;
+        if (isTOWGS84Compatible()) {
+            auto params = transformation()->getTOWGS84Parameters();
+            if (params == std::vector<double>{0, 0, 0, 0, 0, 0, 0}) {
+                refIsNullTransform = true;
+            }
+        }
+
         for (const auto &pair : resTemp) {
             const auto &candidateBaseCRS = pair.first;
             auto projCRS =
@@ -4567,54 +4589,46 @@ BoundCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
             if (geodCRS) {
                 auto context = operation::CoordinateOperationContext::create(
                     authorityFactory, nullptr, 0.0);
+                context->setSpatialCriterion(
+                    operation::CoordinateOperationContext::SpatialCriterion::
+                        PARTIAL_INTERSECTION);
                 auto ops =
                     operation::CoordinateOperationFactory::create()
                         ->createOperations(NN_NO_CHECK(geodCRS),
                                            GeographicCRS::EPSG_4326, context);
-                std::string refTransfPROJString;
-                bool refTransfPROJStringValid = false;
-                try {
-                    refTransfPROJString =
-                        d->transformation_->exportToPROJString(
-                            io::PROJStringFormatter::create().get());
-                    refTransfPROJStringValid = true;
-                    if (refTransfPROJString == "+proj=axisswap +order=2,1") {
-                        refTransfPROJString.clear();
-                    }
-                } catch (const std::exception &) {
-                }
+
                 bool foundOp = false;
                 for (const auto &op : ops) {
+                    auto opNormalized = op->normalizeForVisualization();
                     std::string opTransfPROJString;
                     bool opTransfPROJStringValid = false;
                     if (op->nameStr().find("Ballpark geographic") == 0) {
-                        if (isTOWGS84Compatible()) {
-                            auto params =
-                                transformation()->getTOWGS84Parameters();
-                            if (params ==
-                                std::vector<double>{0, 0, 0, 0, 0, 0, 0}) {
-                                res.emplace_back(create(candidateBaseCRS,
-                                                        d->hubCRS_,
-                                                        transformation()),
-                                                 pair.second);
-                                foundOp = true;
-                                break;
-                            }
+                        if (refIsNullTransform) {
+                            res.emplace_back(create(candidateBaseCRS,
+                                                    d->hubCRS_,
+                                                    transformation()),
+                                             pair.second);
+                            foundOp = true;
+                            break;
                         }
                         continue;
                     }
                     try {
-                        opTransfPROJString = op->exportToPROJString(
+                        opTransfPROJString = opNormalized->exportToPROJString(
                             io::PROJStringFormatter::create().get());
                         opTransfPROJStringValid = true;
+                        opTransfPROJString = replaceAll(
+                            opTransfPROJString, " +rx=0 +ry=0 +rz=0 +s=0 "
+                                                "+convention=position_vector",
+                            "");
                     } catch (const std::exception &) {
                     }
                     if ((refTransfPROJStringValid && opTransfPROJStringValid &&
                          refTransfPROJString == opTransfPROJString) ||
-                        op->_isEquivalentTo(
-                            d->transformation_.get(),
+                        opNormalized->_isEquivalentTo(
+                            refTransf.get(),
                             util::IComparable::Criterion::EQUIVALENT)) {
-                        res.emplace_back(
+                        resMatchOfTransfToWGS84.emplace_back(
                             create(candidateBaseCRS, d->hubCRS_,
                                    NN_NO_CHECK(util::nn_dynamic_pointer_cast<
                                                operation::Transformation>(op))),
@@ -4631,7 +4645,7 @@ BoundCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
             }
         }
     }
-    return res;
+    return !resMatchOfTransfToWGS84.empty() ? resMatchOfTransfToWGS84 : res;
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index 33d67e0a..df1c257d 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -3825,6 +3825,21 @@ TEST(crs, boundCRS_identify_db) {
             WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
         EXPECT_TRUE(wkt.find("32122") != std::string::npos) << wkt;
     }
+
+    {
+        // Identify from a PROJ string with +towgs84
+        auto obj = PROJStringParser().createFromPROJString(
+            "+proj=utm +zone=48 +a=6377276.345 +b=6356075.41314024 "
+            "+towgs84=198,881,317,0,0,0,0 +units=m +no_defs +type=crs");
+        auto crs = nn_dynamic_pointer_cast<BoundCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        auto res = crs->identify(factoryEPSG);
+        ASSERT_EQ(res.size(), 1U);
+        auto boundCRS = dynamic_cast<const BoundCRS *>(res.front().first.get());
+        ASSERT_TRUE(boundCRS != nullptr);
+        EXPECT_EQ(boundCRS->baseCRS()->getEPSGCode(), 3148);
+        EXPECT_EQ(res.front().second, 70);
+    }
 }
 
 // ---------------------------------------------------------------------------