createOperations(): fix transformation involving CompoundCRS, ToWGS84 and PROJ4_GRIDS

Fix issue reported in https://lists.osgeo.org/pipermail/gdal-dev/2022-March/055587.html

2 files changed, 77 insertions, 0 deletions

diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp
index a6a2c986..d5869331 100644
--- a/src/iso19111/operation/coordinateoperationfactory.cpp
+++ b/src/iso19111/operation/coordinateoperationfactory.cpp
@@ -5097,6 +5097,28 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
                         componentsSrc[0], NN_NO_CHECK(interpolationGeogCRS),
                         context);
 
+                    // e.g when doing COMPOUND_CRS[
+                    //      NAD83(CRS)+TOWGS84[0,0,0],
+                    //      CGVD28 height + EXTENSION["PROJ4_GRIDS","HT2_0.gtx"]
+                    // to NAD83(CRS) 3D
+                    const auto boundSrc =
+                        dynamic_cast<crs::BoundCRS *>(componentsSrc[0].get());
+                    if (boundSrc &&
+                        boundSrc->baseCRS()->isEquivalentTo(
+                            targetCRS->demoteTo2D(std::string(), dbContext)
+                                .get(),
+                            util::IComparable::Criterion::EQUIVALENT) &&
+                        boundSrc->hubCRS()->isEquivalentTo(
+                            interpolationGeogCRS
+                                ->demoteTo2D(std::string(), dbContext)
+                                .get(),
+                            util::IComparable::Criterion::EQUIVALENT)) {
+                        // Make sure to use the same horizontal transformation
+                        // (likely a null shift)
+                        interpToTargetOps = applyInverse(srcToInterpOps);
+                        return;
+                    }
+
                     // But do the interpolation CRS to targetCRS in 3D
                     // to have proper ellipsoid height transformation.
                     // We need to force the vertical axis of this 3D'ified
diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
index bb69e347..01d0da9b 100644
--- a/test/unit/test_operationfactory.cpp
+++ b/test/unit/test_operationfactory.cpp
@@ -5969,6 +5969,61 @@ TEST(operation, compoundCRS_of_vertCRS_with_geoid_model_to_geogCRS) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation,
+     compoundCRS_of_horizCRS_with_TOWGS84_vertCRS_with_geoid_model_to_geogCRS) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    auto wkt = "COMPD_CS[\"NAD83(CSRS) + CGVD28 height - HT2_0\",\n"
+               "    GEOGCS[\"NAD83(CSRS)\",\n"
+               "        DATUM[\"NAD83_Canadian_Spatial_Reference_System\",\n"
+               "            SPHEROID[\"GRS 1980\",6378137,298.257222101,\n"
+               "                AUTHORITY[\"EPSG\",\"7019\"]],\n"
+               "            TOWGS84[0,0,0,0,0,0,0],\n"
+               "            AUTHORITY[\"EPSG\",\"6140\"]],\n"
+               "        PRIMEM[\"Greenwich\",0,\n"
+               "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+               "        UNIT[\"degree\",0.0174532925199433,\n"
+               "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
+               "        AUTHORITY[\"EPSG\",\"4617\"]],\n"
+               "    VERT_CS[\"CGVD28 height - HT2_0\",\n"
+               "        VERT_DATUM[\"Canadian Geodetic Vertical Datum of "
+               "1928\",2005,\n"
+               "            EXTENSION[\"PROJ4_GRIDS\",\"HT2_0.gtx\"],\n"
+               "            AUTHORITY[\"EPSG\",\"5114\"]],\n"
+               "        UNIT[\"metre\",1,\n"
+               "            AUTHORITY[\"EPSG\",\"9001\"]],\n"
+               "        AXIS[\"Gravity-related height\",UP],\n"
+               "        AUTHORITY[\"EPSG\",\"5713\"]]]";
+    auto srcObj =
+        createFromUserInput(wkt, authFactory->databaseContext(), false);
+    auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
+    ASSERT_TRUE(src != nullptr);
+    // NAD83(CSRS) 3D
+    auto dst = authFactory->createCoordinateReferenceSystem("4955");
+
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        NN_NO_CHECK(src), dst, ctxt);
+    ASSERT_EQ(list.size(), 1U);
+    auto op_proj =
+        list[0]->exportToPROJString(PROJStringFormatter::create().get());
+    EXPECT_EQ(op_proj, "+proj=pipeline "
+                       "+step +proj=push +v_1 +v_2 "
+                       "+step +proj=axisswap +order=2,1 "
+                       "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+                       "+step +proj=vgridshift +grids=HT2_0.gtx +multiplier=1 "
+                       "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+                       "+step +proj=axisswap +order=2,1 "
+                       "+step +proj=pop +v_1 +v_2");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, compoundCRS_from_WKT2_to_geogCRS_3D_context) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), "EPSG");