Adjust createBoundCRSToWGS84IfPossible() and operation filtering (for POSGAR 2007 to WGS84 issues) (fixes #2356)

- We make createBoundCRSToWGS84IfPossible() more restrictive. If there are more
  than one Helmert transformation from the CRS to WGS 84 covering the area of
  use of the CRS, we do not create a BoundCRS / +towgs84
- In createOperations() filtering, we are less aggressive in discarding operations
  that have the same area of use but worse accuracy. We do it only if they involve
  more transformation steps.

We now get:
```
$ projinfo EPSG:5340 -o PROJ
PROJ.4 string:
+proj=longlat +ellps=GRS80 +no_defs +type=crs

$ projinfo -s EPSG:5340 -t EPSG:4326 --spatial-test intersects --summary
Candidate operations found: 2
EPSG:9264, POSGAR 2007 to WGS 84 (2), 0.5 m, Argentina
EPSG:5351, POSGAR 2007 to WGS 84 (1), 1.0 m, Argentina
```

7 files changed, 77 insertions, 151 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 10b1baf4..0e61d2ee 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -11407,6 +11407,23 @@ static size_t getStepCount(const CoordinateOperationNNPtr &op) {
 
 // ---------------------------------------------------------------------------
 
+// Return number of steps that are transformations (and not conversions)
+static size_t getTransformationStepCount(const CoordinateOperationNNPtr &op) {
+    auto concat = dynamic_cast<const ConcatenatedOperation *>(op.get());
+    size_t stepCount = 1;
+    if (concat) {
+        stepCount = 0;
+        for (const auto &subOp : concat->operations()) {
+            if (dynamic_cast<const Conversion *>(subOp.get()) == nullptr) {
+                stepCount++;
+            }
+        }
+    }
+    return stepCount;
+}
+
+// ---------------------------------------------------------------------------
+
 static bool isNullTransformation(const std::string &name) {
     if (name.find(" + ") != std::string::npos)
         return false;
@@ -11442,8 +11459,7 @@ struct FilterResults {
         // results
         // ...
         removeSyntheticNullTransforms();
-        if (context->getDiscardSuperseded())
-            removeUninterestingOps();
+        removeUninterestingOps();
         removeDuplicateOps();
         removeSyntheticNullTransforms();
         return *this;
@@ -11758,50 +11774,23 @@ struct FilterResults {
     void removeUninterestingOps() {
 
         // Eliminate operations that bring nothing, ie for a given area of use,
-        // do not keep operations that have greater accuracy. Actually we must
-        // be a bit more subtle than that, and take into account grid
-        // availability
+        // do not keep operations that have similar or worse accuracy, but
+        // involve more (non conversion) steps
         std::vector<CoordinateOperationNNPtr> resTemp;
         metadata::ExtentPtr lastExtent;
         double lastAccuracy = -1;
-        bool lastHasGrids = false;
-        bool lastGridsAvailable = true;
-        std::set<std::set<std::string>> setOfSetOfGrids;
         size_t lastStepCount = 0;
         CoordinateOperationPtr lastOp;
 
         bool first = true;
-        const auto gridAvailabilityUse = context->getGridAvailabilityUse();
         for (const auto &op : res) {
             const auto curAccuracy = getAccuracy(op);
             bool dummy = false;
             const auto curExtent = getExtent(op, true, dummy);
-            bool curHasGrids = false;
-            bool curGridsAvailable = true;
-            std::set<std::string> curSetOfGrids;
-
-            const auto curStepCount = getStepCount(op);
-
-            if (context->getAuthorityFactory()) {
-                const auto gridsNeeded = op->gridsNeeded(
-                    context->getAuthorityFactory()->databaseContext(),
-                    gridAvailabilityUse ==
-                        CoordinateOperationContext::GridAvailabilityUse::
-                            KNOWN_AVAILABLE);
-                for (const auto &gridDesc : gridsNeeded) {
-                    curHasGrids = true;
-                    curSetOfGrids.insert(gridDesc.shortName);
-                    if (!gridDesc.available) {
-                        curGridsAvailable = false;
-                    }
-                }
-            }
+            const auto curStepCount = getTransformationStepCount(op);
 
             if (first) {
                 resTemp.emplace_back(op);
-
-                lastHasGrids = curHasGrids;
-                lastGridsAvailable = curGridsAvailable;
                 first = false;
             } else {
                 if (lastOp->_isEquivalentTo(op.get())) {
@@ -11812,66 +11801,19 @@ struct FilterResults {
                      (curExtent && lastExtent &&
                       curExtent->contains(NN_NO_CHECK(lastExtent)) &&
                       lastExtent->contains(NN_NO_CHECK(curExtent))));
-                if (((curAccuracy > lastAccuracy && lastAccuracy >= 0) ||
+                if (((curAccuracy >= lastAccuracy && lastAccuracy >= 0) ||
                      (curAccuracy < 0 && lastAccuracy >= 0)) &&
-                    sameExtent) {
-                    // If that set of grids has always been used for that
-                    // extent,
-                    // no need to add them again
-                    if (setOfSetOfGrids.find(curSetOfGrids) !=
-                        setOfSetOfGrids.end()) {
-                        continue;
-                    }
-
-                    const bool sameNameOrEmptyName =
-                        ((!curExtent && !lastExtent) ||
-                         (curExtent && lastExtent &&
-                          !curExtent->description()->empty() &&
-                          *(curExtent->description()) ==
-                              *(lastExtent->description())));
-
-                    // If we have already found a operation without grids for
-                    // that extent, no need to add any lower accuracy operation
-                    if (!lastHasGrids && sameNameOrEmptyName) {
-                        continue;
-                    }
-                    // If we had only operations involving grids, but one
-                    // past operation had available grids, no need to add
-                    // the new one.
-                    if (curHasGrids && curGridsAvailable &&
-                        lastGridsAvailable) {
-                        continue;
-                    }
-                } else if (curAccuracy == lastAccuracy && sameExtent) {
-                    if (curStepCount > lastStepCount) {
-                        continue;
-                    }
+                    sameExtent && curStepCount > lastStepCount) {
+                    continue;
                 }
 
                 resTemp.emplace_back(op);
-
-                if (sameExtent) {
-                    if (!curHasGrids) {
-                        lastHasGrids = false;
-                    }
-                    if (curGridsAvailable) {
-                        lastGridsAvailable = true;
-                    }
-                } else {
-                    setOfSetOfGrids.clear();
-
-                    lastHasGrids = curHasGrids;
-                    lastGridsAvailable = curGridsAvailable;
-                }
             }
 
             lastOp = op.as_nullable();
             lastStepCount = curStepCount;
             lastExtent = curExtent;
             lastAccuracy = curAccuracy;
-            if (!curSetOfGrids.empty()) {
-                setOfSetOfGrids.insert(curSetOfGrids);
-            }
         }
         res = std::move(resTemp);
     }
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 88420c8a..e96b3cc9 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -378,7 +378,8 @@ VerticalCRSPtr CRS::extractVerticalCRS() const {
  * a +towgs84 parameter or a WKT1:GDAL string with a TOWGS node.
  *
  * This method will fetch the GeographicCRS of this CRS and find a
- * transformation to EPSG:4326 using the domain of the validity of the main CRS.
+ * transformation to EPSG:4326 using the domain of the validity of the main CRS,
+ * and there's only one Helmert transformation.
  *
  * @return a CRS.
  */
@@ -456,6 +457,7 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
             auto list =
                 operation::CoordinateOperationFactory::create()
                     ->createOperations(NN_NO_CHECK(geodCRS), hubCRS, ctxt);
+            CRSPtr candidateBoundCRS;
             for (const auto &op : list) {
                 auto transf =
                     util::nn_dynamic_pointer_cast<operation::Transformation>(
@@ -466,8 +468,13 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
                     } catch (const std::exception &) {
                         continue;
                     }
-                    return util::nn_static_pointer_cast<CRS>(BoundCRS::create(
-                        thisAsCRS, hubCRS, NN_NO_CHECK(transf)));
+                    if (candidateBoundCRS) {
+                        candidateBoundCRS = nullptr;
+                        break;
+                    }
+                    candidateBoundCRS =
+                        BoundCRS::create(thisAsCRS, hubCRS, NN_NO_CHECK(transf))
+                            .as_nullable();
                 } else {
                     auto concatenated =
                         dynamic_cast<const operation::ConcatenatedOperation *>(
@@ -499,15 +506,23 @@ CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
                                     } catch (const std::exception &) {
                                         continue;
                                     }
-                                    return util::nn_static_pointer_cast<CRS>(
+                                    if (candidateBoundCRS) {
+                                        candidateBoundCRS = nullptr;
+                                        break;
+                                    }
+                                    candidateBoundCRS =
                                         BoundCRS::create(thisAsCRS, hubCRS,
-                                                         NN_NO_CHECK(transf)));
+                                                         NN_NO_CHECK(transf))
+                                            .as_nullable();
                                 }
                             }
                         }
                     }
                 }
             }
+            if (candidateBoundCRS) {
+                return NN_NO_CHECK(candidateBoundCRS);
+            }
         } catch (const std::exception &) {
         }
     }
diff --git a/test/cli/proj_outIGNF.dist b/test/cli/proj_outIGNF.dist
index c12b883b..35d053fa 100644
--- a/test/cli/proj_outIGNF.dist
+++ b/test/cli/proj_outIGNF.dist
@@ -8,7 +8,7 @@
  311552.5340 1906457.4840  0.0000	358799.172	6342652.486 0.000
  960488.4138 1910172.8812  0.0000	1007068.686	6340907.237 0.000
  600000.0000 1699510.8340  0.0000	645204.279	6133556.746 0.000
-1203792.5981 626873.17210  0.0000	*	* inf
+1203792.5981 626873.17210  0.0000	1238837.253	5057451.037 0.000
 +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX
  600000.0000 2600545.4523  0.0000	179040.148	5610495.275 0.000
  135638.3592 2418760.4094  0.0000	-303729.363	5410118.356 0.000
@@ -17,7 +17,7 @@
  311552.5340 1906457.4840  0.0000	-96825.465	4909184.136 0.000
  960488.4138 1910172.8812  0.0000	523880.019	4909191.141 0.000
  600000.0000 1699510.8340  0.0000	179047.633	4708817.007 0.000
-1203792.5981 626873.17210  0.0000	*	* inf
+1203792.5981 626873.17210  0.0000	658259.467	3623786.764 0.000
 +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX
 2d20'11.4239243" 50d23'59.7718445" 0.0	179040.151	5610495.281 0.000
 -3d57'49.4051448" 48d35'59.7121716" 0.0	-303729.365	5410118.352 0.000
diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist
index 8ece2e62..5e1f1dd5 100644
--- a/test/cli/testprojinfo_out.dist
+++ b/test/cli/testprojinfo_out.dist
@@ -1044,8 +1044,9 @@ Testing -s +proj=longlat +datum=WGS84 +geoidgrids=@foo.gtx +type=crs -t EPSG:432
 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=@foo.gtx +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 
 Testing -s "GDA94" -t "WGS 84 (G1762)" --spatial-test intersects --summary. Should include transformations through ITRF2008 and GDA2020
-Candidate operations found: 6
+Candidate operations found: 7
 unknown id, GDA94 to GDA2020 (1) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.21 m, Australia - GDA
+unknown id, GDA94 to WGS 84 (1) + WGS 84 to WGS 84 (G1762), 5 m, Australia - GDA
 unknown id, Conversion from GDA94 (geog2D) to GDA94 (geocentric) + Inverse of ITRF2008 to GDA94 (1) + Inverse of WGS 84 (G1762) to ITRF2008 (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.04 m, Australia - onshore and EEZ
 unknown id, GDA94 to GDA2020 (2) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Australia - onshore
 unknown id, GDA94 to GDA2020 (3) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Australia - onshore
@@ -1053,16 +1054,19 @@ unknown id, GDA94 to GDA2020 (5) + Conversion from GDA2020 (geog2D) to GDA2020 (
 unknown id, GDA94 to GDA2020 (4) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Christmas Island - onshore
 
 Testing -s "AGD66" -t "WGS 84 (G1762)" --spatial-test intersects --summary. Should include a transformation through GDA2020
-Candidate operations found: 11
+Candidate operations found: 14
 unknown id, AGD66 to WGS 84 (18) + WGS 84 to WGS 84 (G1762), 5 m, Australia - offshore
 unknown id, AGD66 to WGS 84 (16) + WGS 84 to WGS 84 (G1762), 7 m, Australia - onshore
+unknown id, AGD66 to WGS 84 (20) + WGS 84 to WGS 84 (G1762), 11 m, Australia - onshore
 unknown id, AGD66 to WGS 84 (15) + WGS 84 to WGS 84 (G1762), 3 m, Australia - Northern Territory
 unknown id, AGD66 to WGS 84 (13) + WGS 84 to WGS 84 (G1762), 3 m, Australia - New South Wales and Victoria
 unknown id, AGD66 to WGS 84 (21) + WGS 84 to WGS 84 (G1762), 7 m, Papua New Guinea - mainland onshore
 unknown id, AGD66 to WGS 84 (14) + WGS 84 to WGS 84 (G1762), 3 m, Australia - Tasmania
 unknown id, AGD66 to WGS 84 (19) + WGS 84 to WGS 84 (G1762), 4 m, Papua New Guinea - PFTB
+unknown id, AGD66 to WGS 84 (22) + WGS 84 to WGS 84 (G1762), 6 m, Papua New Guinea - PFTB
 unknown id, AGD66 to WGS 84 (23) + WGS 84 to WGS 84 (G1762), 6 m, Papua New Guinea - North Fly
 unknown id, AGD66 to GDA2020 (1) + Conversion from GDA2020 (geog2D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), 0.25 m, Australia - Australian Capital Territory
+unknown id, AGD66 to WGS 84 (12) + WGS 84 to WGS 84 (G1762), 3 m, Australia - Australian Capital Territory
 unknown id, AGD66 to WGS 84 (17) + WGS 84 to WGS 84 (G1762), 3 m, Australia - onshore
 unknown id, Ballpark geographic offset from AGD66 to WGS 84 (G1762), unknown accuracy, World, has ballpark transformation
 
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index b9a6cdd5..564833db 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -540,7 +540,7 @@ TEST_F(CApi, proj_as_proj_string_approx_tmerc_option_yes) {
 // ---------------------------------------------------------------------------
 
 TEST_F(CApi, proj_crs_create_bound_crs_to_WGS84) {
-    auto crs = proj_create_from_database(m_ctxt, "EPSG", "3844",
+    auto crs = proj_create_from_database(m_ctxt, "EPSG", "4807",
                                          PJ_CATEGORY_CRS, false, nullptr);
     ObjectKeeper keeper(crs);
     ASSERT_NE(crs, nullptr);
@@ -552,10 +552,8 @@ TEST_F(CApi, proj_crs_create_bound_crs_to_WGS84) {
     auto proj_4 = proj_as_proj_string(m_ctxt, res, PJ_PROJ_4, nullptr);
     ASSERT_NE(proj_4, nullptr);
     EXPECT_EQ(std::string(proj_4),
-              "+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 "
-              "+y_0=500000 +ellps=krass "
-              "+towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69 "
-              "+units=m +no_defs +type=crs");
+              "+proj=longlat +ellps=clrk80ign +pm=paris "
+              "+towgs84=-168,-60,320,0,0,0,0 +no_defs +type=crs");
 
     auto base_crs = proj_get_source_crs(m_ctxt, res);
     ObjectKeeper keeper_base_crs(base_crs);
@@ -572,8 +570,7 @@ TEST_F(CApi, proj_crs_create_bound_crs_to_WGS84) {
     std::vector<double> values(7, 0);
     EXPECT_TRUE(proj_coordoperation_get_towgs84_values(m_ctxt, transf,
                                                        values.data(), 7, true));
-    auto expected = std::vector<double>{2.329, -147.042, -92.08, -0.309,
-                                        0.325, 0.497,    5.69};
+    auto expected = std::vector<double>{-168, -60, 320, 0, 0, 0, 0};
     EXPECT_EQ(values, expected);
 
     auto res2 = proj_crs_create_bound_crs(m_ctxt, base_crs, hub_crs, transf);
@@ -1485,7 +1482,7 @@ TEST_F(CApi, proj_create_operations_discard_superseded) {
     ASSERT_NE(res, nullptr);
     ObjListKeeper keeper_res(res);
 
-    EXPECT_EQ(proj_list_get_count(res), 2);
+    EXPECT_EQ(proj_list_get_count(res), 4);
 }
 
 // ---------------------------------------------------------------------------
@@ -1594,7 +1591,7 @@ TEST_F(CApi, proj_create_operations_with_pivot) {
         auto res = proj_create_operations(m_ctxt, source_crs, target_crs, ctxt);
         ASSERT_NE(res, nullptr);
         ObjListKeeper keeper_res(res);
-        EXPECT_EQ(proj_list_get_count(res), 7);
+        EXPECT_EQ(proj_list_get_count(res), 8);
         auto op = proj_list_get(m_ctxt, res, 0);
         ASSERT_NE(op, nullptr);
         ObjectKeeper keeper_op(op);
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index caee0646..0e340560 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -5422,40 +5422,18 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) {
     {
         // Pulkovo 42 Romania
         auto crs_3844 = factory->createCoordinateReferenceSystem("3844");
-        auto bound = crs_3844->createBoundCRSToWGS84IfPossible(
-            dbContext, CoordinateOperationContext::IntermediateCRSUse::NEVER);
-        EXPECT_NE(bound, crs_3844);
-        EXPECT_EQ(bound->createBoundCRSToWGS84IfPossible(
+        EXPECT_EQ(crs_3844->createBoundCRSToWGS84IfPossible(
                       dbContext,
                       CoordinateOperationContext::IntermediateCRSUse::NEVER),
-                  bound);
-        auto boundCRS = nn_dynamic_pointer_cast<BoundCRS>(bound);
-        ASSERT_TRUE(boundCRS != nullptr);
-        EXPECT_EQ(
-            boundCRS->exportToPROJString(PROJStringFormatter::create().get()),
-            "+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 "
-            "+y_0=500000 +ellps=krass "
-            "+towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69 "
-            "+units=m +no_defs +type=crs");
+                  crs_3844);
     }
     {
         // Pulkovo 42 Poland
         auto crs_2171 = factory->createCoordinateReferenceSystem("2171");
-        auto bound = crs_2171->createBoundCRSToWGS84IfPossible(
-            dbContext, CoordinateOperationContext::IntermediateCRSUse::NEVER);
-        EXPECT_NE(bound, crs_2171);
-        EXPECT_EQ(bound->createBoundCRSToWGS84IfPossible(
+        EXPECT_EQ(crs_2171->createBoundCRSToWGS84IfPossible(
                       dbContext,
                       CoordinateOperationContext::IntermediateCRSUse::NEVER),
-                  bound);
-        auto boundCRS = nn_dynamic_pointer_cast<BoundCRS>(bound);
-        ASSERT_TRUE(boundCRS != nullptr);
-        EXPECT_EQ(
-            boundCRS->exportToPROJString(PROJStringFormatter::create().get()),
-            "+proj=sterea +lat_0=50.625 +lon_0=21.0833333333333 "
-            "+k=0.9998 +x_0=4637000 +y_0=5647000 +ellps=krass "
-            "+towgs84=33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84 "
-            "+units=m +no_defs +type=crs");
+                  crs_2171);
     }
     {
         // NTF (Paris)
@@ -5547,6 +5525,15 @@ TEST(crs, crs_createBoundCRSToWGS84IfPossible) {
             ::testing::UnitTest::GetInstance()->elapsed_time();
         EXPECT_LE(time_after - time_before, 500);
     }
+    {
+        // POSGAR 2007: it has 2 helmert shifts to WGS84 (#2356). Don't take
+        // an arbitrary one
+        auto crs_5340 = factory->createCoordinateReferenceSystem("5340");
+        EXPECT_EQ(crs_5340->createBoundCRSToWGS84IfPossible(
+                      dbContext,
+                      CoordinateOperationContext::IntermediateCRSUse::NEVER),
+                  crs_5340);
+    }
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 0238638d..eca5d59b 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -4599,7 +4599,7 @@ TEST(operation, geogCRS_to_geogCRS_context_NAD27_to_WGS84) {
         authFactory->createCoordinateReferenceSystem("4267"), // NAD27
         authFactory->createCoordinateReferenceSystem("4326"), // WGS84
         ctxt);
-    ASSERT_EQ(list.size(), 78U);
+    ASSERT_EQ(list.size(), 79U);
     EXPECT_EQ(list[0]->nameStr(),
               "NAD27 to WGS 84 (33)"); // 1.0 m, Canada - NAD27
     EXPECT_EQ(list[1]->nameStr(),
@@ -5050,7 +5050,7 @@ TEST(operation, geogCRS_to_geogCRS_context_concatenated_operation) {
         authFactory->createCoordinateReferenceSystem("4807"), // NTF(Paris)
         authFactory->createCoordinateReferenceSystem("4171"), // RGF93
         ctxt);
-    ASSERT_EQ(list.size(), 5U);
+    ASSERT_EQ(list.size(), 4U);
 
     EXPECT_EQ(list[0]->nameStr(), "NTF (Paris) to RGF93 (1)");
     EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
@@ -5158,7 +5158,7 @@ TEST(operation, geogCRS_to_geogCRS_CH1903_to_CH1903plus_context) {
         authFactory->createCoordinateReferenceSystem("4149"), // CH1903
         authFactory->createCoordinateReferenceSystem("4150"), // CH1903+
         ctxt);
-    ASSERT_TRUE(list.size() == 2U || list.size() == 3U);
+    ASSERT_TRUE(list.size() == 1U);
 
     EXPECT_EQ(list[0]->nameStr(), "CH1903 to CH1903+ (1)");
     EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
@@ -5167,25 +5167,6 @@ TEST(operation, geogCRS_to_geogCRS_CH1903_to_CH1903plus_context) {
               "+step +proj=hgridshift +grids=ch_swisstopo_CHENyx06a.tif "
               "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
               "+step +proj=axisswap +order=2,1");
-
-    if (list.size() == 2U) {
-        // Grids not there
-        EXPECT_EQ(list[1]->nameStr(),
-                  "CH1903 to ETRS89 (1) + Inverse of CH1903+ to ETRS89 (1)");
-        EXPECT_EQ(
-            list[1]->exportToPROJString(PROJStringFormatter::create().get()),
-            "+proj=noop");
-    } else {
-        // Grids available
-        EXPECT_EQ(list[1]->nameStr(),
-                  "CH1903 to ETRS89 (2) + Inverse of CH1903+ to ETRS89 (1)");
-
-        EXPECT_EQ(list[2]->nameStr(),
-                  "CH1903 to ETRS89 (1) + Inverse of CH1903+ to ETRS89 (1)");
-        EXPECT_EQ(
-            list[2]->exportToPROJString(PROJStringFormatter::create().get()),
-            "+proj=noop");
-    }
 }
 
 // ---------------------------------------------------------------------------
@@ -6069,7 +6050,7 @@ TEST(operation, projCRS_to_projCRS_context_compatible_area) {
         authFactory->createCoordinateReferenceSystem(
             "2171"), // Pulkovo 42 Poland I
         ctxt);
-    ASSERT_EQ(list.size(), 1U);
+    ASSERT_EQ(list.size(), 2U);
     EXPECT_EQ(list[0]->nameStr(),
               "Inverse of UTM zone 34N + Inverse of Pulkovo 1942(58) to WGS 84 "
               "(1) + Poland zone I");
@@ -6088,7 +6069,7 @@ TEST(operation, projCRS_to_projCRS_context_compatible_area_bis) {
             "3844"), // Pulkovo 42 Stereo 70 (Romania)
         authFactory->createCoordinateReferenceSystem("32634"), // UTM 34
         ctxt);
-    ASSERT_EQ(list.size(), 1U);
+    ASSERT_EQ(list.size(), 3U);
     EXPECT_EQ(list[0]->nameStr(), "Inverse of Stereo 70 + "
                                   "Pulkovo 1942(58) to WGS 84 "
                                   "(19) + UTM zone 34N");
@@ -6107,7 +6088,7 @@ TEST(operation, projCRS_to_projCRS_context_one_incompatible_area) {
         authFactory->createCoordinateReferenceSystem(
             "2171"), // Pulkovo 42 Poland I
         ctxt);
-    ASSERT_EQ(list.size(), 1U);
+    ASSERT_EQ(list.size(), 2U);
     EXPECT_EQ(list[0]->nameStr(),
               "Inverse of UTM zone 31N + Inverse of Pulkovo 1942(58) to WGS 84 "
               "(1) + Poland zone I");