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
```

3 files changed, 25 insertions, 60 deletions

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");