Fix issue with WGS 84 + EGM2008 -> Amersfoort + NAP height (master only)

Fixes https://lists.osgeo.org/pipermail/proj/2026-February/011959.html

Output of projinfo EPSG:4326+3855 EPSG:7415 --spatial-test intersects
---------------------------------------------------------------------

9.7.1
+++++

Operation No. 1:

```
unknown id, Inverse of WGS 84 to EGM2008 height (1) +
            Inverse of ETRS89 to WGS 84 (1) +
            ETRS89 to NAP height (2) +
            Inverse of Amersfoort to ETRS89 (9) +
            RD New,
            1.114 m, Netherlands - onshore, including Waddenzee, Dutch Wadden Islands and 12-mile offshore coastal zone., at least one grid missing

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1
  +step +inv +proj=vgridshift +grids=nl_nsgi_nlgeo2018.tif +multiplier=1
  +step +inv +proj=hgridshift +grids=nl_nsgi_rdtrans2018.tif
  +step +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079
        +x_0=155000 +y_0=463000 +ellps=bessel
```

master before this change
+++++++++++++++++++++++++

Operation No. 1:

```
unknown id, Inverse of WGS 84 to EGM2008 height (1) +
            Inverse of ETRS89-NLD [AGRS2010] to WGS 84 (1) +
            ETRS89-NLD [AGRS2010] to NAP height (2) +
            Inverse of Amersfoort to ETRS89-NLD [AGRS2010] (9) +
            RD New,
            1.123 m, Netherlands - onshore and offshore., at least one grid missing

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1
  +step +inv +proj=vgridshift +grids=nl_nsgi_nlgeo2018.tif +multiplier=1
  +step +inv +proj=hgridshift +grids=nl_nsgi_rdtrans2018.tif
  +step +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079
        +x_0=155000 +y_0=463000 +ellps=bessel

    OPERATIONACCURACY[1.123],
    USAGE[
        SCOPE["unknown"],
        AREA["Netherlands - onshore and offshore."],
        BBOX[50.75,2.53,55.77,7.22]]]
```

Operation No. 3:
```
unknown id, Inverse of WGS 84 to EGM2008 height (1) +
            Inverse of Amersfoort to WGS 84 (4) +
            ETRS89-NLD [AGRS2010] to NAP height (2) using Amersfoort to ETRS89-NLD [AGRS2010] (9) +
            RD New,
            1.123 m, Netherlands - onshore, including Waddenzee, Dutch Wadden Islands and 12-mile offshore coastal zone., at least one grid missing

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1
  +step +proj=cart +ellps=WGS84
  +step +inv +proj=helmert +x=565.4171 +y=50.3319 +z=465.5524
        +rx=0.398957388243134 +ry=-0.343987817378283 +rz=1.87740163998045
        +s=4.0725 +convention=coordinate_frame
  +step +inv +proj=cart +ellps=bessel
  +step +proj=push +v_1 +v_2
  +step +proj=hgridshift +grids=nl_nsgi_rdtrans2018.tif +omit_inv
  +step +inv +proj=vgridshift +grids=nl_nsgi_nlgeo2018.tif +multiplier=1
  +step +inv +proj=hgridshift +grids=nl_nsgi_rdtrans2018.tif +omit_fwd
  +step +proj=pop +v_1 +v_2
  +step +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079
        +x_0=155000 +y_0=463000 +ellps=bessel
    OPERATIONACCURACY[1.123],
    USAGE[
        SCOPE["unknown"],
        AREA["Netherlands - onshore, including Waddenzee, Dutch Wadden Islands and 12-mile offshore coastal zone."],
        BBOX[50.75,3.2,53.7,7.22]]]
```

Analysis
--------

Due to the ETRS89 datum ensemble changes, PROJ synthetized a new operation (operation 3)
that has the same accuracy than the "correct" (operation 1), and a smaller area of
use, thus it was selected. However that choice was not super relevant compared to 1,
because it relied on using Amersfoort CRS in 3D (with ellipsoidal height), through
the Helmert transformation "Amersfoort to WGS 84 (4)"). Given that there's no entry
in the database for a 3D Amersfoort CRS, this possibility should be discarded when
there is an alternative that uses only known objects.

Author: rouault

src/iso19111/operation/coordinateoperationfactory.cpp

@@ -7033,15 +7033,15 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
                     }
                 }
 
-                // Keep the results of this new attempt, if there are better
+                // Use the results of this new attempt if they are better
                 // than the previous ones
                 if (bestAccuracy2 >= 0 &&
                     (bestAccuracy < 0 || (bestAccuracy2 < bestAccuracy ||
                                           (bestAccuracy2 == bestAccuracy &&
                                            bestStepCount2 < bestStepCount)))) {
                     bestAccuracy = bestAccuracy2;
                     bestStepCount = bestStepCount2;
-                    res.insert(res.end(), res2.begin(), res2.end());
+                    res = std::move(res2);
                 }
             };
 

test/unit/test_operationfactory.cpp

@@ -6585,6 +6585,71 @@ TEST(operation, compoundCRS_to_compoundCRS_issue_3328) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, compoundCRS_to_compoundCRS_WGS84_EGM2008_to_RD_new_NAP_height) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        // WGS 84 + EGM2008 height
+        authFactory->createCoordinateReferenceSystem("9518"),
+        // Amersfoort / RD New + NAP height
+        authFactory->createCoordinateReferenceSystem("7415"), ctxt);
+
+    ASSERT_EQ(list.size(), 4U);
+
+    EXPECT_FALSE(list[0]->hasBallparkTransformation());
+    EXPECT_EQ(list[0]->nameStr(),
+              "Inverse of WGS 84 to EGM2008 height (1) + "
+              "Inverse of ETRS89-NLD [AGRS2010] to WGS 84 (1) + "
+              "ETRS89-NLD [AGRS2010] to NAP height (2) + "
+              "Inverse of Amersfoort to ETRS89-NLD [AGRS2010] (9) + RD New");
+    EXPECT_EQ(
+        list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+        "+proj=pipeline "
+        "+step +proj=axisswap +order=2,1 "
+        "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+        "+step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1 "
+        "+step +inv +proj=vgridshift +grids=nl_nsgi_nlgeo2018.tif "
+        "+multiplier=1 "
+        "+step +inv +proj=hgridshift +grids=nl_nsgi_rdtrans2018.tif "
+        "+step +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 "
+        "+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel");
+    ASSERT_EQ(list[0]->coordinateOperationAccuracies().size(), 1U);
+    EXPECT_EQ(list[0]->coordinateOperationAccuracies()[0]->value(), "1.123");
+
+    EXPECT_FALSE(list[1]->hasBallparkTransformation());
+    EXPECT_EQ(list[1]->nameStr(),
+              "Inverse of WGS 84 to EGM2008 height (1) + "
+              "Inverse of ETRS89-NLD [AGRS2010] to WGS 84 (1) + "
+              "ETRS89-NLD [AGRS2010] to NAP height (2) + "
+              "Inverse of Amersfoort to ETRS89-NLD [AGRS2010] (8) + RD New");
+    ASSERT_EQ(list[1]->coordinateOperationAccuracies().size(), 1U);
+    EXPECT_EQ(list[1]->coordinateOperationAccuracies()[0]->value(), "1.373");
+
+    // Using not available "WGS 84 to EGM2008 height (2)" with 1' EGM2008 grid
+    EXPECT_FALSE(list[2]->hasBallparkTransformation());
+    EXPECT_EQ(list[2]->nameStr(),
+              "Inverse of WGS 84 to EGM2008 height (2) + "
+              "Inverse of ETRS89-NLD [AGRS2010] to WGS 84 (1) + "
+              "ETRS89-NLD [AGRS2010] to NAP height (2) + "
+              "Inverse of Amersfoort to ETRS89-NLD [AGRS2010] (9) + RD New");
+
+    // Using not available "WGS 84 to EGM2008 height (2)" with 1' EGM2008 grid
+    EXPECT_FALSE(list[3]->hasBallparkTransformation());
+    EXPECT_EQ(list[3]->nameStr(),
+              "Inverse of WGS 84 to EGM2008 height (2) + "
+              "Inverse of ETRS89-NLD [AGRS2010] to WGS 84 (1) + "
+              "ETRS89-NLD [AGRS2010] to NAP height (2) + "
+              "Inverse of Amersfoort to ETRS89-NLD [AGRS2010] (8) + RD New");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(
     operation,
     compoundCRS_to_compoundCRS_concatenated_operation_with_two_vert_transformation_and_ballpark_geog) {