From d37884e00f27907464675558b2da529729299dcf Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 18 Mar 2022 23:01:50 +0100 Subject: createOperations(): fix transformation involving CompoundCRS, ToWGS84 and PROJ4_GRIDS Fix issue reported in https://lists.osgeo.org/pipermail/gdal-dev/2022-March/055587.html --- .../operation/coordinateoperationfactory.cpp | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) (limited to 'src') 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(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 -- cgit v1.2.3