Merge pull request #3124 from OSGeo/backport-3123-to-9.09.0

[Backport 9.0] createOperations(): fix transformation involving CompoundCRS, ToWGS84 and PROJ4_GRIDS

1 files changed, 22 insertions, 0 deletions

diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp
index 20042f22..e8b1a371 100644
--- a/src/iso19111/operation/coordinateoperationfactory.cpp
+++ b/src/iso19111/operation/coordinateoperationfactory.cpp
@@ -5084,6 +5084,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