Merge pull request #1454 from rouault/fix_ETRS89_3D_to_proj_string_with_geoidgrids_nadgrids

createOperations(): fix case of ETRS89 3D to proj string with nadgrids and geoidgrids

1 files changed, 53 insertions, 4 deletions

diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index a38e9df2..5ac1609d 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -5903,6 +5903,35 @@ TEST(operation, boundCRS_with_basecrs_with_extent_to_geogCRS) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, ETRS89_3D_to_proj_string_with_geoidgrids_nadgrids) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    // ETRS89 3D
+    auto src = authFactory->createCoordinateReferenceSystem("4937");
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 "
+        "+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel "
+        "+nadgrids=rdtrans2008.gsb +geoidgrids=naptrans2008.gtx +units=m "
+        "+type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        src, NN_NO_CHECK(dst), ctxt);
+    ASSERT_EQ(list.size(), 1U);
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +inv +proj=vgridshift +grids=naptrans2008.gtx "
+              "+multiplier=1 "
+              "+step +inv +proj=hgridshift +grids=rdtrans2008.gsb "
+              "+step +proj=sterea +lat_0=52.1561605555556 "
+              "+lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 "
+              "+y_0=463000 +ellps=bessel");
+}
+
+// ---------------------------------------------------------------------------
+
 static VerticalCRSNNPtr createVerticalCRS() {
     PropertyMap propertiesVDatum;
     propertiesVDatum.set(Identifier::CODESPACE_KEY, "EPSG")
@@ -5941,8 +5970,8 @@ static BoundCRSNNPtr createBoundVerticalCRS() {
     auto vertCRS = createVerticalCRS();
     auto transformation =
         Transformation::createGravityRelatedHeightToGeographic3D(
-            PropertyMap(), vertCRS, GeographicCRS::EPSG_4979, "egm08_25.gtx",
-            std::vector<PositionalAccuracyNNPtr>());
+            PropertyMap(), vertCRS, GeographicCRS::EPSG_4979, nullptr,
+            "egm08_25.gtx", std::vector<PositionalAccuracyNNPtr>());
     return BoundCRS::create(vertCRS, GeographicCRS::EPSG_4979, transformation);
 }
 
@@ -6727,6 +6756,22 @@ TEST(operation, compoundCRS_from_WKT2_no_id_to_geogCRS_3D_context) {
     auto list =
         CoordinateOperationFactory::create()->createOperations(src, dst, ctxt);
     ASSERT_GE(list.size(), 1U);
+
+    {
+        // Important here is vgridshift before hgridshift
+        auto op_proj =
+            list[0]->exportToPROJString(PROJStringFormatter::create().get());
+        EXPECT_EQ(
+            op_proj,
+            "+proj=pipeline +step +inv +proj=sterea +lat_0=52.1561605555556 "
+            "+lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 "
+            "+ellps=bessel "
+            "+step +proj=vgridshift +grids=naptrans2008.gtx +multiplier=1 "
+            "+step +proj=hgridshift +grids=rdtrans2008.gsb "
+            "+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
+            "+step +proj=axisswap +order=2,1");
+    }
+
     auto wkt2 =
         "COMPOUNDCRS[\"unknown\",\n"
         "  PROJCRS[\"unknown\",\n"
@@ -6759,8 +6804,12 @@ TEST(operation, compoundCRS_from_WKT2_no_id_to_geogCRS_3D_context) {
     for (size_t i = 0; i < list.size(); i++) {
         const auto &op = list[i];
         const auto &op2 = list2[i];
-        EXPECT_TRUE(
-            op->isEquivalentTo(op2.get(), IComparable::Criterion::EQUIVALENT));
+        auto op_proj =
+            op->exportToPROJString(PROJStringFormatter::create().get());
+        auto op2_proj =
+            op2->exportToPROJString(PROJStringFormatter::create().get());
+        EXPECT_EQ(op_proj, op2_proj) << "op=" << op->nameStr()
+                                     << " op2=" << op2->nameStr();
     }
 }