Merge pull request #2222 from rouault/fix_2217

Fixes related to CompoundCRS and BoundCRS

5 files changed, 515 insertions, 149 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 0f5e4c4d..3ffd208f 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -2043,15 +2043,19 @@ SingleOperation::gridsNeeded(const io::DatabaseContextPtr &databaseContext,
         if (opParamvalue) {
             const auto &value = opParamvalue->parameterValue();
             if (value->type() == ParameterValue::Type::FILENAME) {
-                GridDescription desc;
-                desc.shortName = value->valueFile();
-                if (databaseContext) {
-                    databaseContext->lookForGridInfo(
-                        desc.shortName, considerKnownGridsAsAvailable,
-                        desc.fullName, desc.packageName, desc.url,
-                        desc.directDownload, desc.openLicense, desc.available);
+                const auto gridNames = split(value->valueFile(), ",");
+                for (const auto &gridName : gridNames) {
+                    GridDescription desc;
+                    desc.shortName = gridName;
+                    if (databaseContext) {
+                        databaseContext->lookForGridInfo(
+                            desc.shortName, considerKnownGridsAsAvailable,
+                            desc.fullName, desc.packageName, desc.url,
+                            desc.directDownload, desc.openLicense,
+                            desc.available);
+                    }
+                    res.insert(desc);
                 }
-                res.insert(desc);
             }
         }
     }
@@ -6527,7 +6531,7 @@ struct Transformation::Private {
 
     TransformationPtr forwardOperation_{};
 
-    TransformationNNPtr registerInv(util::BaseObjectNNPtr thisIn,
+    TransformationNNPtr registerInv(const Transformation *thisIn,
                                     TransformationNNPtr invTransform);
 };
 //! @endcond
@@ -7996,12 +8000,11 @@ createApproximateInverseIfPossible(const Transformation *op) {
 
 //! @cond Doxygen_Suppress
 TransformationNNPtr
-Transformation::Private::registerInv(util::BaseObjectNNPtr thisIn,
+Transformation::Private::registerInv(const Transformation *thisIn,
                                      TransformationNNPtr invTransform) {
-    invTransform->d->forwardOperation_ =
-        util::nn_dynamic_pointer_cast<Transformation>(thisIn);
+    invTransform->d->forwardOperation_ = thisIn->shallowClone().as_nullable();
     invTransform->setHasBallparkTransformation(
-        invTransform->d->forwardOperation_->hasBallparkTransformation());
+        thisIn->hasBallparkTransformation());
     return invTransform;
 }
 //! @endcond
@@ -8039,12 +8042,24 @@ TransformationNNPtr Transformation::inverseAsTransformation() const {
             parameterValueNumericAsSI(EPSG_CODE_PARAMETER_Y_AXIS_TRANSLATION);
         double z =
             parameterValueNumericAsSI(EPSG_CODE_PARAMETER_Z_AXIS_TRANSLATION);
+        auto properties = createPropertiesForInverse(this, false, false);
         return d->registerInv(
-            shared_from_this(),
-            createGeocentricTranslations(
-                createPropertiesForInverse(this, false, false), l_targetCRS,
-                l_sourceCRS, negate(x), negate(y), negate(z),
-                coordinateOperationAccuracies()));
+            this, create(properties, l_targetCRS, l_sourceCRS, nullptr,
+                         createMethodMapNameEPSGCode(
+                             useOperationMethodEPSGCodeIfPresent(
+                                 properties, methodEPSGCode)),
+                         VectorOfParameters{
+                             createOpParamNameEPSGCode(
+                                 EPSG_CODE_PARAMETER_X_AXIS_TRANSLATION),
+                             createOpParamNameEPSGCode(
+                                 EPSG_CODE_PARAMETER_Y_AXIS_TRANSLATION),
+                             createOpParamNameEPSGCode(
+                                 EPSG_CODE_PARAMETER_Z_AXIS_TRANSLATION),
+                         },
+                         createParams(common::Length(negate(x)),
+                                      common::Length(negate(y)),
+                                      common::Length(negate(z))),
+                         coordinateOperationAccuracies()));
     }
 
     if (methodEPSGCode == EPSG_CODE_METHOD_MOLODENSKY ||
@@ -8062,14 +8077,14 @@ TransformationNNPtr Transformation::inverseAsTransformation() const {
 
         if (methodEPSGCode == EPSG_CODE_METHOD_ABRIDGED_MOLODENSKY) {
             return d->registerInv(
-                shared_from_this(),
+                this,
                 createAbridgedMolodensky(
                     createPropertiesForInverse(this, false, false), l_targetCRS,
                     l_sourceCRS, negate(x), negate(y), negate(z), negate(da),
                     negate(df), coordinateOperationAccuracies()));
         } else {
             return d->registerInv(
-                shared_from_this(),
+                this,
                 createMolodensky(createPropertiesForInverse(this, false, false),
                                  l_targetCRS, l_sourceCRS, negate(x), negate(y),
                                  negate(z), negate(da), negate(df),
@@ -8082,10 +8097,9 @@ TransformationNNPtr Transformation::inverseAsTransformation() const {
             parameterValueMeasure(EPSG_CODE_PARAMETER_LONGITUDE_OFFSET);
         const common::Angle newOffset(negate(offset.value()), offset.unit());
         return d->registerInv(
-            shared_from_this(),
-            createLongitudeRotation(
-                createPropertiesForInverse(this, false, false), l_targetCRS,
-                l_sourceCRS, newOffset));
+            this, createLongitudeRotation(
+                      createPropertiesForInverse(this, false, false),
+                      l_targetCRS, l_sourceCRS, newOffset));
     }
 
     if (methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC2D_OFFSETS) {
@@ -8100,11 +8114,10 @@ TransformationNNPtr Transformation::inverseAsTransformation() const {
                                           offsetLong.unit());
 
         return d->registerInv(
-            shared_from_this(),
-            createGeographic2DOffsets(
-                createPropertiesForInverse(this, false, false), l_targetCRS,
-                l_sourceCRS, newOffsetLat, newOffsetLong,
-                coordinateOperationAccuracies()));
+            this, createGeographic2DOffsets(
+                      createPropertiesForInverse(this, false, false),
+                      l_targetCRS, l_sourceCRS, newOffsetLat, newOffsetLong,
+                      coordinateOperationAccuracies()));
     }
 
     if (methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC3D_OFFSETS) {
@@ -8124,11 +8137,10 @@ TransformationNNPtr Transformation::inverseAsTransformation() const {
                                              offsetHeight.unit());
 
         return d->registerInv(
-            shared_from_this(),
-            createGeographic3DOffsets(
-                createPropertiesForInverse(this, false, false), l_targetCRS,
-                l_sourceCRS, newOffsetLat, newOffsetLong, newOffsetHeight,
-                coordinateOperationAccuracies()));
+            this, createGeographic3DOffsets(
+                      createPropertiesForInverse(this, false, false),
+                      l_targetCRS, l_sourceCRS, newOffsetLat, newOffsetLong,
+                      newOffsetHeight, coordinateOperationAccuracies()));
     }
 
     if (methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC2D_WITH_HEIGHT_OFFSETS) {
@@ -8148,11 +8160,10 @@ TransformationNNPtr Transformation::inverseAsTransformation() const {
                                              offsetHeight.unit());
 
         return d->registerInv(
-            shared_from_this(),
-            createGeographic2DWithHeightOffsets(
-                createPropertiesForInverse(this, false, false), l_targetCRS,
-                l_sourceCRS, newOffsetLat, newOffsetLong, newOffsetHeight,
-                coordinateOperationAccuracies()));
+            this, createGeographic2DWithHeightOffsets(
+                      createPropertiesForInverse(this, false, false),
+                      l_targetCRS, l_sourceCRS, newOffsetLat, newOffsetLong,
+                      newOffsetHeight, coordinateOperationAccuracies()));
     }
 
     if (methodEPSGCode == EPSG_CODE_METHOD_VERTICAL_OFFSET) {
@@ -8163,7 +8174,7 @@ TransformationNNPtr Transformation::inverseAsTransformation() const {
                                              offsetHeight.unit());
 
         return d->registerInv(
-            shared_from_this(),
+            this,
             createVerticalOffset(createPropertiesForInverse(this, false, false),
                                  l_targetCRS, l_sourceCRS, newOffsetHeight,
                                  coordinateOperationAccuracies()));
@@ -8173,18 +8184,17 @@ TransformationNNPtr Transformation::inverseAsTransformation() const {
         const double convFactor = parameterValueNumericAsSI(
             EPSG_CODE_PARAMETER_UNIT_CONVERSION_SCALAR);
         return d->registerInv(
-            shared_from_this(),
-            createChangeVerticalUnit(
-                createPropertiesForInverse(this, false, false), l_targetCRS,
-                l_sourceCRS, common::Scale(1.0 / convFactor),
-                coordinateOperationAccuracies()));
+            this, createChangeVerticalUnit(
+                      createPropertiesForInverse(this, false, false),
+                      l_targetCRS, l_sourceCRS, common::Scale(1.0 / convFactor),
+                      coordinateOperationAccuracies()));
     }
 
 #ifdef notdef
     // We don't need that currently, but we might...
     if (methodEPSGCode == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL) {
         return d->registerInv(
-            shared_from_this(),
+            this,
             createHeightDepthReversal(
                 createPropertiesForInverse(this, false, false), l_targetCRS,
                 l_sourceCRS, coordinateOperationAccuracies()));
@@ -8989,20 +8999,41 @@ static void ThrowExceptionNotGeodeticGeographic(const char *trfrm_name) {
 
 // ---------------------------------------------------------------------------
 
-static void setupPROJGeodeticSourceCRS(io::PROJStringFormatter *formatter,
-                                       const crs::CRSNNPtr &crs, bool addPushV3,
-                                       const char *trfrm_name) {
-    auto sourceCRSGeog = dynamic_cast<const crs::GeographicCRS *>(crs.get());
-    if (!sourceCRSGeog) {
-        auto compoundCRS = dynamic_cast<const crs::CompoundCRS *>(crs.get());
+// If crs is a geographic CRS, or a compound CRS of a geographic CRS,
+// or a compoundCRS of a bound CRS of a geographic CRS, return that
+// geographic CRS
+static crs::GeographicCRSPtr
+extractGeographicCRSIfGeographicCRSOrEquivalent(const crs::CRSNNPtr &crs) {
+    auto geogCRS = util::nn_dynamic_pointer_cast<crs::GeographicCRS>(crs);
+    if (!geogCRS) {
+        auto compoundCRS = util::nn_dynamic_pointer_cast<crs::CompoundCRS>(crs);
         if (compoundCRS) {
             const auto &components = compoundCRS->componentReferenceSystems();
             if (!components.empty()) {
-                sourceCRSGeog = dynamic_cast<const crs::GeographicCRS *>(
-                    components[0].get());
+                geogCRS = util::nn_dynamic_pointer_cast<crs::GeographicCRS>(
+                    components[0]);
+                if (!geogCRS) {
+                    auto boundCRS =
+                        util::nn_dynamic_pointer_cast<crs::BoundCRS>(
+                            components[0]);
+                    if (boundCRS) {
+                        geogCRS =
+                            util::nn_dynamic_pointer_cast<crs::GeographicCRS>(
+                                boundCRS->baseCRS());
+                    }
+                }
             }
         }
     }
+    return geogCRS;
+}
+
+// ---------------------------------------------------------------------------
+
+static void setupPROJGeodeticSourceCRS(io::PROJStringFormatter *formatter,
+                                       const crs::CRSNNPtr &crs, bool addPushV3,
+                                       const char *trfrm_name) {
+    auto sourceCRSGeog = extractGeographicCRSIfGeographicCRSOrEquivalent(crs);
     if (sourceCRSGeog) {
         formatter->startInversion();
         sourceCRSGeog->_exportToPROJString(formatter);
@@ -9030,17 +9061,7 @@ static void setupPROJGeodeticSourceCRS(io::PROJStringFormatter *formatter,
 static void setupPROJGeodeticTargetCRS(io::PROJStringFormatter *formatter,
                                        const crs::CRSNNPtr &crs, bool addPopV3,
                                        const char *trfrm_name) {
-    auto targetCRSGeog = dynamic_cast<const crs::GeographicCRS *>(crs.get());
-    if (!targetCRSGeog) {
-        auto compoundCRS = dynamic_cast<const crs::CompoundCRS *>(crs.get());
-        if (compoundCRS) {
-            const auto &components = compoundCRS->componentReferenceSystems();
-            if (!components.empty()) {
-                targetCRSGeog = dynamic_cast<const crs::GeographicCRS *>(
-                    components[0].get());
-            }
-        }
-    }
+    auto targetCRSGeog = extractGeographicCRSIfGeographicCRSOrEquivalent(crs);
     if (targetCRSGeog) {
         formatter->addStep("cart");
         formatter->setCurrentStepInverted(true);
@@ -9540,14 +9561,14 @@ void Transformation::_exportToPROJString(
                                                          : CTABLE2Filename;
     if (!hGridShiftFilename.empty()) {
         auto sourceCRSGeog =
-            dynamic_cast<const crs::GeographicCRS *>(sourceCRS().get());
+            extractGeographicCRSIfGeographicCRSOrEquivalent(sourceCRS());
         if (!sourceCRSGeog) {
             throw io::FormattingException(
                 concat("Can apply ", methodName, " only to GeographicCRS"));
         }
 
         auto targetCRSGeog =
-            dynamic_cast<const crs::GeographicCRS *>(targetCRS().get());
+            extractGeographicCRSIfGeographicCRSOrEquivalent(targetCRS());
         if (!targetCRSGeog) {
             throw io::FormattingException(
                 concat("Can apply ", methodName, " only to GeographicCRS"));
@@ -13191,7 +13212,7 @@ CoordinateOperationFactory::Private::createOperations(
         return applyInverse(createOperations(targetCRS, sourceCRS, context));
     }
 
-    // boundCRS to boundCRS using the same geographic hubCRS
+    // boundCRS to boundCRS
     if (boundSrc && boundDst) {
         createOperationsBoundToBound(sourceCRS, targetCRS, context, boundSrc,
                                      boundDst, res);
@@ -14327,6 +14348,7 @@ void CoordinateOperationFactory::Private::createOperationsBoundToBound(
 
     ENTER_FUNCTION();
 
+    // BoundCRS to BoundCRS of horizontal CRS using the same (geographic) hub
     const auto &hubSrc = boundSrc->hubCRS();
     auto hubSrcGeog = dynamic_cast<const crs::GeographicCRS *>(hubSrc.get());
     const auto &hubDst = boundDst->hubCRS();
@@ -14375,6 +14397,28 @@ void CoordinateOperationFactory::Private::createOperationsBoundToBound(
         }
     }
 
+    // BoundCRS to BoundCRS of vertical CRS using the same vertical datum
+    // ==> ignore the bound transformation
+    auto baseOfBoundSrcAsVertCRS =
+        dynamic_cast<crs::VerticalCRS *>(boundSrc->baseCRS().get());
+    auto baseOfBoundDstAsVertCRS =
+        dynamic_cast<crs::VerticalCRS *>(boundDst->baseCRS().get());
+    if (baseOfBoundSrcAsVertCRS && baseOfBoundDstAsVertCRS) {
+        const auto datumSrc = baseOfBoundSrcAsVertCRS->datum();
+        const auto datumDst = baseOfBoundDstAsVertCRS->datum();
+        if (datumSrc && datumDst &&
+            datumSrc->nameStr() == datumDst->nameStr() &&
+            (datumSrc->nameStr() != "unknown" ||
+             boundSrc->transformation()->_isEquivalentTo(
+                 boundDst->transformation().get(),
+                 util::IComparable::Criterion::EQUIVALENT))) {
+            res = createOperations(boundSrc->baseCRS(), boundDst->baseCRS(),
+                                   context);
+            return;
+        }
+    }
+
+    // BoundCRS to BoundCRS of vertical CRS
     auto vertCRSOfBaseOfBoundSrc = boundSrc->baseCRS()->extractVerticalCRS();
     auto vertCRSOfBaseOfBoundDst = boundDst->baseCRS()->extractVerticalCRS();
     if (hubSrcGeog && hubDstGeog &&
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index c23bd29b..ec342a8f 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -4620,18 +4620,9 @@ BoundCRS::create(const CRSNNPtr &baseCRSIn, const CRSNNPtr &hubCRSIn,
 BoundCRSNNPtr
 BoundCRS::createFromTOWGS84(const CRSNNPtr &baseCRSIn,
                             const std::vector<double> &TOWGS84Parameters) {
-
-    auto geodCRS = baseCRSIn->extractGeodeticCRS();
-    auto targetCRS =
-        geodCRS.get() == nullptr ||
-                dynamic_cast<const crs::GeographicCRS *>(geodCRS.get())
-            ? util::nn_static_pointer_cast<crs::CRS>(
-                  crs::GeographicCRS::EPSG_4326)
-            : util::nn_static_pointer_cast<crs::CRS>(
-                  crs::GeodeticCRS::EPSG_4978);
-    return create(
-        baseCRSIn, targetCRS,
-        operation::Transformation::createTOWGS84(baseCRSIn, TOWGS84Parameters));
+    auto transf =
+        operation::Transformation::createTOWGS84(baseCRSIn, TOWGS84Parameters);
+    return create(baseCRSIn, transf->targetCRS(), transf);
 }
 
 // ---------------------------------------------------------------------------
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 97f741be..9d9fa4a3 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -3987,6 +3987,79 @@ WKTParser::Private::buildParametricDatum(const WKTNodeNNPtr &node) {
 
 // ---------------------------------------------------------------------------
 
+static CRSNNPtr
+createBoundCRSSourceTransformationCRS(const crs::CRSPtr &sourceCRS,
+                                      const crs::CRSPtr &targetCRS) {
+    CRSPtr sourceTransformationCRS;
+    if (dynamic_cast<GeographicCRS *>(targetCRS.get())) {
+        GeographicCRSPtr sourceGeographicCRS =
+            sourceCRS->extractGeographicCRS();
+        sourceTransformationCRS = sourceGeographicCRS;
+        if (sourceGeographicCRS) {
+            if (sourceGeographicCRS->datum() != nullptr &&
+                sourceGeographicCRS->primeMeridian()
+                        ->longitude()
+                        .getSIValue() != 0.0) {
+                sourceTransformationCRS =
+                    GeographicCRS::create(
+                        util::PropertyMap().set(
+                            common::IdentifiedObject::NAME_KEY,
+                            sourceGeographicCRS->nameStr() +
+                                " (with Greenwich prime meridian)"),
+                        datum::GeodeticReferenceFrame::create(
+                            util::PropertyMap().set(
+                                common::IdentifiedObject::NAME_KEY,
+                                sourceGeographicCRS->datum()->nameStr() +
+                                    " (with Greenwich prime meridian)"),
+                            sourceGeographicCRS->datum()->ellipsoid(),
+                            util::optional<std::string>(),
+                            datum::PrimeMeridian::GREENWICH),
+                        sourceGeographicCRS->coordinateSystem())
+                        .as_nullable();
+            }
+        } else {
+            auto vertSourceCRS =
+                std::dynamic_pointer_cast<VerticalCRS>(sourceCRS);
+            if (!vertSourceCRS) {
+                throw ParsingException(
+                    "Cannot find GeographicCRS or VerticalCRS in sourceCRS");
+            }
+            const auto &axis = vertSourceCRS->coordinateSystem()->axisList()[0];
+            if (axis->unit() == common::UnitOfMeasure::METRE &&
+                &(axis->direction()) == &AxisDirection::UP) {
+                sourceTransformationCRS = sourceCRS;
+            } else {
+                std::string sourceTransformationCRSName(
+                    vertSourceCRS->nameStr());
+                if (ends_with(sourceTransformationCRSName, " (ftUS)")) {
+                    sourceTransformationCRSName.resize(
+                        sourceTransformationCRSName.size() - strlen(" (ftUS)"));
+                }
+                if (ends_with(sourceTransformationCRSName, " depth")) {
+                    sourceTransformationCRSName.resize(
+                        sourceTransformationCRSName.size() - strlen(" depth"));
+                }
+                if (!ends_with(sourceTransformationCRSName, " height")) {
+                    sourceTransformationCRSName += " height";
+                }
+                sourceTransformationCRS =
+                    VerticalCRS::create(
+                        PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                          sourceTransformationCRSName),
+                        vertSourceCRS->datum(), vertSourceCRS->datumEnsemble(),
+                        VerticalCS::createGravityRelatedHeight(
+                            common::UnitOfMeasure::METRE))
+                        .as_nullable();
+            }
+        }
+    } else {
+        sourceTransformationCRS = sourceCRS;
+    }
+    return NN_NO_CHECK(sourceTransformationCRS);
+}
+
+// ---------------------------------------------------------------------------
+
 CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) {
     const auto *nodeP = node->GP();
     auto &datumNode =
@@ -4111,16 +4184,18 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) {
                     gridName != "g2012a_conus.gtx,g2012a_alaska.gtx,"
                                 "g2012a_guam.gtx,g2012a_hawaii.gtx,"
                                 "g2012a_puertorico.gtx,g2012a_samoa.gtx") {
-                    std::string transformationName(crs->nameStr());
-                    if (!ends_with(transformationName, " height")) {
-                        transformationName += " height";
-                    }
-                    transformationName += " to WGS84 ellipsoidal height";
+                    auto sourceTransformationCRS =
+                        createBoundCRSSourceTransformationCRS(
+                            crs.as_nullable(),
+                            GeographicCRS::EPSG_4979.as_nullable());
                     auto transformation = Transformation::
                         createGravityRelatedHeightToGeographic3D(
-                            PropertyMap().set(IdentifiedObject::NAME_KEY,
-                                              transformationName),
-                            crs, GeographicCRS::EPSG_4979, nullptr, gridName,
+                            PropertyMap().set(
+                                IdentifiedObject::NAME_KEY,
+                                sourceTransformationCRS->nameStr() +
+                                    " to WGS84 ellipsoidal height"),
+                            sourceTransformationCRS, GeographicCRS::EPSG_4979,
+                            nullptr, gridName,
                             std::vector<PositionalAccuracyNNPtr>());
                     return nn_static_pointer_cast<CRS>(BoundCRS::create(
                         crs, GeographicCRS::EPSG_4979, transformation));
@@ -4190,52 +4265,6 @@ CRSNNPtr WKTParser::Private::buildCompoundCRS(const WKTNodeNNPtr &node) {
 
 // ---------------------------------------------------------------------------
 
-static CRSNNPtr
-createBoundCRSSourceTransformationCRS(const crs::CRSPtr &sourceCRS,
-                                      const crs::CRSPtr &targetCRS) {
-    CRSPtr sourceTransformationCRS;
-    if (dynamic_cast<GeographicCRS *>(targetCRS.get())) {
-        GeographicCRSPtr sourceGeographicCRS =
-            sourceCRS->extractGeographicCRS();
-        sourceTransformationCRS = sourceGeographicCRS;
-        if (sourceGeographicCRS) {
-            if (sourceGeographicCRS->datum() != nullptr &&
-                sourceGeographicCRS->primeMeridian()
-                        ->longitude()
-                        .getSIValue() != 0.0) {
-                sourceTransformationCRS =
-                    GeographicCRS::create(
-                        util::PropertyMap().set(
-                            common::IdentifiedObject::NAME_KEY,
-                            sourceGeographicCRS->nameStr() +
-                                " (with Greenwich prime meridian)"),
-                        datum::GeodeticReferenceFrame::create(
-                            util::PropertyMap().set(
-                                common::IdentifiedObject::NAME_KEY,
-                                sourceGeographicCRS->datum()->nameStr() +
-                                    " (with Greenwich prime meridian)"),
-                            sourceGeographicCRS->datum()->ellipsoid(),
-                            util::optional<std::string>(),
-                            datum::PrimeMeridian::GREENWICH),
-                        sourceGeographicCRS->coordinateSystem())
-                        .as_nullable();
-            }
-        } else {
-            sourceTransformationCRS =
-                std::dynamic_pointer_cast<VerticalCRS>(sourceCRS);
-            if (!sourceTransformationCRS) {
-                throw ParsingException(
-                    "Cannot find GeographicCRS or VerticalCRS in sourceCRS");
-            }
-        }
-    } else {
-        sourceTransformationCRS = sourceCRS;
-    }
-    return NN_NO_CHECK(sourceTransformationCRS);
-}
-
-// ---------------------------------------------------------------------------
-
 BoundCRSNNPtr WKTParser::Private::buildBoundCRS(const WKTNodeNNPtr &node) {
     const auto *nodeP = node->GP();
     auto &abridgedNode =
@@ -4515,12 +4544,27 @@ CRSPtr WKTParser::Private::buildCRS(const WKTNodeNNPtr &node) {
     const auto *nodeP = node->GP();
     const std::string &name(nodeP->value());
 
+    const auto applyHorizontalBoundCRSParams = [&](const CRSNNPtr &crs) {
+        if (!toWGS84Parameters_.empty()) {
+            auto ret = BoundCRS::createFromTOWGS84(crs, toWGS84Parameters_);
+            toWGS84Parameters_.clear();
+            return util::nn_static_pointer_cast<CRS>(ret);
+        } else if (!datumPROJ4Grids_.empty()) {
+            auto ret = BoundCRS::createFromNadgrids(crs, datumPROJ4Grids_);
+            datumPROJ4Grids_.clear();
+            return util::nn_static_pointer_cast<CRS>(ret);
+        }
+        return crs;
+    };
+
     if (isGeodeticCRS(name)) {
         if (!isNull(nodeP->lookForChild(WKTConstants::BASEGEOGCRS,
                                         WKTConstants::BASEGEODCRS))) {
-            return buildDerivedGeodeticCRS(node);
+            return util::nn_static_pointer_cast<CRS>(
+                applyHorizontalBoundCRSParams(buildDerivedGeodeticCRS(node)));
         } else {
-            return util::nn_static_pointer_cast<CRS>(buildGeodeticCRS(node));
+            return util::nn_static_pointer_cast<CRS>(
+                applyHorizontalBoundCRSParams(buildGeodeticCRS(node)));
         }
     }
 
@@ -4545,12 +4589,14 @@ CRSPtr WKTParser::Private::buildCRS(const WKTNodeNNPtr &node) {
                     PROJStringParser().createFromPROJString(projString);
                 auto crs = nn_dynamic_pointer_cast<CRS>(projObj);
                 if (crs) {
-                    return crs;
+                    return util::nn_static_pointer_cast<CRS>(
+                        applyHorizontalBoundCRSParams(NN_NO_CHECK(crs)));
                 }
             } catch (const io::ParsingException &) {
             }
         }
-        return util::nn_static_pointer_cast<CRS>(buildProjectedCRS(node));
+        return util::nn_static_pointer_cast<CRS>(
+            applyHorizontalBoundCRSParams(buildProjectedCRS(node)));
     }
 
     if (ci_equal(name, WKTConstants::VERT_CS) ||
@@ -4623,16 +4669,6 @@ BaseObjectNNPtr WKTParser::Private::build(const WKTNodeNNPtr &node) {
 
     auto crs = buildCRS(node);
     if (crs) {
-        if (!toWGS84Parameters_.empty()) {
-            return util::nn_static_pointer_cast<BaseObject>(
-                BoundCRS::createFromTOWGS84(NN_NO_CHECK(crs),
-                                            toWGS84Parameters_));
-        }
-        if (!datumPROJ4Grids_.empty()) {
-            return util::nn_static_pointer_cast<BaseObject>(
-                BoundCRS::createFromNadgrids(NN_NO_CHECK(crs),
-                                             datumPROJ4Grids_));
-        }
         return util::nn_static_pointer_cast<BaseObject>(NN_NO_CHECK(crs));
     }
 
@@ -7594,7 +7630,10 @@ std::set<std::string> PROJStringFormatter::getUsedGridNames() const {
     for (const auto &step : d->steps_) {
         for (const auto &param : step.paramValues) {
             if (param.keyEquals("grids")) {
-                res.insert(param.value);
+                const auto gridNames = split(param.value, ",");
+                for (const auto &gridName : gridNames) {
+                    res.insert(gridName);
+                }
             }
         }
     }
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index dff804e8..7ad9cf4a 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -3442,6 +3442,32 @@ TEST(wkt_parse, WKT1_VERT_DATUM_EXTENSION) {
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, WKT1_VERT_DATUM_EXTENSION_units_ftUS) {
+    auto wkt = "VERT_CS[\"NAVD88 height (ftUS)\","
+               "     VERT_DATUM[\"North American Vertical Datum 1988\",2005,"
+               "         EXTENSION[\"PROJ4_GRIDS\",\"foo.gtx\"],"
+               "         AUTHORITY[\"EPSG\",\"5103\"]],"
+               "     UNIT[\"US survey foot\",0.304800609601219,"
+               "         AUTHORITY[\"EPSG\",\"9003\"]],"
+               "     AXIS[\"Gravity-related height\",UP],"
+               "     AUTHORITY[\"EPSG\",\"6360\"]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<BoundCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    EXPECT_EQ(crs->transformation()->nameStr(),
+              "NAVD88 height to WGS84 ellipsoidal height"); // no (ftUS)
+    auto sourceTransformationCRS = crs->transformation()->sourceCRS();
+    auto sourceTransformationVertCRS =
+        nn_dynamic_pointer_cast<VerticalCRS>(sourceTransformationCRS);
+    EXPECT_EQ(
+        sourceTransformationVertCRS->coordinateSystem()->axisList()[0]->unit(),
+        UnitOfMeasure::METRE);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(wkt_parse, WKT1_DATUM_EXTENSION) {
     auto wkt =
         "PROJCS[\"unnamed\",\n"
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 9b339525..3d002a17 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -7136,6 +7136,61 @@ TEST(operation,
 
 // ---------------------------------------------------------------------------
 
+TEST(operation,
+     compoundCRS_with_boundProjCRS_with_ftus_and_boundVerticalCRS_to_geogCRS) {
+
+    auto wkt =
+        "COMPD_CS[\"NAD_1983_StatePlane_Alabama_West_FIPS_0102_Feet + "
+        "NAVD88 height - Geoid12B (US Feet)\",\n"
+        "    PROJCS[\"NAD_1983_StatePlane_Alabama_West_FIPS_0102_Feet\",\n"
+        "        GEOGCS[\"NAD83\",\n"
+        "            DATUM[\"North_American_Datum_1983\",\n"
+        "                SPHEROID[\"GRS 1980\",6378137,298.257222101,\n"
+        "                    AUTHORITY[\"EPSG\",\"7019\"]],\n"
+        "                TOWGS84[0,0,0,0,0,0,0],\n"
+        "                AUTHORITY[\"EPSG\",\"6269\"]],\n"
+        "            PRIMEM[\"Greenwich\",0],\n"
+        "            UNIT[\"Degree\",0.0174532925199433]],\n"
+        "        PROJECTION[\"Transverse_Mercator\"],\n"
+        "        PARAMETER[\"latitude_of_origin\",30],\n"
+        "        PARAMETER[\"central_meridian\",-87.5],\n"
+        "        PARAMETER[\"scale_factor\",0.999933333333333],\n"
+        "        PARAMETER[\"false_easting\",1968500],\n"
+        "        PARAMETER[\"false_northing\",0],\n"
+        "        UNIT[\"US survey foot\",0.304800609601219,\n"
+        "            AUTHORITY[\"EPSG\",\"9003\"]],\n"
+        "        AXIS[\"Easting\",EAST],\n"
+        "        AXIS[\"Northing\",NORTH],\n"
+        "        AUTHORITY[\"ESRI\",\"102630\"]],\n"
+        "    VERT_CS[\"NAVD88 height (ftUS)\",\n"
+        "        VERT_DATUM[\"North American Vertical Datum 1988\",2005,\n"
+        "            EXTENSION[\"PROJ4_GRIDS\",\"foo.gtx\"],\n"
+        "            AUTHORITY[\"EPSG\",\"5103\"]],\n"
+        "        UNIT[\"US survey foot\",0.304800609601219,\n"
+        "            AUTHORITY[\"EPSG\",\"9003\"]],\n"
+        "        AXIS[\"Gravity-related height\",UP],\n"
+        "        AUTHORITY[\"EPSG\",\"6360\"]]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<CompoundCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        NN_NO_CHECK(crs), GeographicCRS::EPSG_4979);
+    ASSERT_TRUE(op != nullptr);
+
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=unitconvert +xy_in=us-ft +xy_out=m "
+              "+step +inv +proj=tmerc +lat_0=30 +lon_0=-87.5 "
+              "+k=0.999933333333333 +x_0=600000 +y_0=0 +ellps=GRS80 "
+              "+step +proj=unitconvert +z_in=us-ft +z_out=m "
+              "+step +proj=vgridshift +grids=foo.gtx +multiplier=1 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, geocent_to_compoundCRS) {
     auto objSrc = PROJStringParser().createFromPROJString(
         "+proj=geocent +datum=WGS84 +units=m +type=crs");
@@ -7296,6 +7351,217 @@ TEST(operation, compoundCRS_to_compoundCRS_with_bound_crs_in_horiz_and_vert) {
 
 // ---------------------------------------------------------------------------
 
+TEST(
+    operation,
+    compoundCRS_to_compoundCRS_with_bound_crs_in_horiz_and_vert_same_geoidgrids) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=longlat +ellps=GRS67 +nadgrids=@foo.gsb +geoidgrids=@foo.gtx "
+        "+type=crs");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +ellps=GRS80 +nadgrids=@bar.gsb +geoidgrids=@foo.gtx "
+        "+type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst));
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +proj=hgridshift +grids=@foo.gsb "
+              "+step +inv +proj=hgridshift +grids=@bar.gsb "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(
+    operation,
+    compoundCRS_to_compoundCRS_with_bound_crs_in_horiz_and_vert_same_geoidgrids_different_vunits) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=longlat +ellps=GRS67 +nadgrids=@foo.gsb +geoidgrids=@foo.gtx "
+        "+type=crs");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +ellps=GRS80 +nadgrids=@bar.gsb +geoidgrids=@foo.gtx "
+        "+vunits=us-ft +type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst));
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +proj=hgridshift +grids=@foo.gsb "
+              "+step +proj=unitconvert +z_in=m +z_out=us-ft "
+              "+step +inv +proj=hgridshift +grids=@bar.gsb "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(
+    operation,
+    compoundCRS_to_compoundCRS_with_bound_crs_in_horiz_and_vert_same_nadgrids_same_geoidgrids) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=longlat +ellps=GRS67 +nadgrids=@foo.gsb +geoidgrids=@foo.gtx "
+        "+type=crs");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +ellps=GRS80 +nadgrids=@foo.gsb +geoidgrids=@foo.gtx "
+        "+type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst));
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=noop");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(
+    operation,
+    compoundCRS_to_compoundCRS_with_bound_crs_in_horiz_and_vert_same_towgs84_same_geoidgrids) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=longlat +ellps=GRS67 +towgs84=0,0,0 +geoidgrids=@foo.gtx "
+        "+type=crs");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +ellps=GRS80 +towgs84=0,0,0 +geoidgrids=@foo.gtx "
+        "+type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst));
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +proj=push +v_3 "
+              "+step +proj=cart +ellps=GRS67 "
+              "+step +inv +proj=cart +ellps=GRS80 "
+              "+step +proj=pop +v_3 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(
+    operation,
+    compoundCRS_to_compoundCRS_with_bound_crs_in_horiz_and_vert_WKT1_same_geoidgrids_context) {
+    auto objSrc = WKTParser().createFromWKT(
+        "COMPD_CS[\"NAD83 / Alabama West + NAVD88 height - Geoid12B "
+        "(Meters)\",\n"
+        "    PROJCS[\"NAD83 / Alabama West\",\n"
+        "        GEOGCS[\"NAD83\",\n"
+        "            DATUM[\"North_American_Datum_1983\",\n"
+        "                SPHEROID[\"GRS 1980\",6378137,298.257222101,\n"
+        "                    AUTHORITY[\"EPSG\",\"7019\"]],\n"
+        "                TOWGS84[0,0,0,0,0,0,0],\n"
+        "                AUTHORITY[\"EPSG\",\"6269\"]],\n"
+        "            PRIMEM[\"Greenwich\",0,\n"
+        "                AUTHORITY[\"EPSG\",\"8901\"]],\n"
+        "            UNIT[\"degree\",0.0174532925199433,\n"
+        "                AUTHORITY[\"EPSG\",\"9122\"]],\n"
+        "            AUTHORITY[\"EPSG\",\"4269\"]],\n"
+        "        PROJECTION[\"Transverse_Mercator\"],\n"
+        "        PARAMETER[\"latitude_of_origin\",30],\n"
+        "        PARAMETER[\"central_meridian\",-87.5],\n"
+        "        PARAMETER[\"scale_factor\",0.999933333],\n"
+        "        PARAMETER[\"false_easting\",600000],\n"
+        "        PARAMETER[\"false_northing\",0],\n"
+        "        UNIT[\"metre\",1,\n"
+        "            AUTHORITY[\"EPSG\",\"9001\"]],\n"
+        "        AXIS[\"X\",EAST],\n"
+        "        AXIS[\"Y\",NORTH],\n"
+        "        AUTHORITY[\"EPSG\",\"26930\"]],\n"
+        "    VERT_CS[\"NAVD88 height\",\n"
+        "        VERT_DATUM[\"North American Vertical Datum 1988\",2005,\n"
+        "            "
+        "EXTENSION[\"PROJ4_GRIDS\",\"g2012a_alaska.gtx,g2012a_hawaii.gtx,"
+        "g2012a_conus.gtx\"],\n"
+        "            AUTHORITY[\"EPSG\",\"5103\"]],\n"
+        "        UNIT[\"metre\",1,\n"
+        "            AUTHORITY[\"EPSG\",\"9001\"]],\n"
+        "        AXIS[\"Gravity-related height\",UP],\n"
+        "        AUTHORITY[\"EPSG\",\"5703\"]]]");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+    auto objDst = WKTParser().createFromWKT(
+        "COMPD_CS[\"NAD_1983_StatePlane_Alabama_West_FIPS_0102_Feet + NAVD88 "
+        "height - Geoid12B (US Feet)\",\n"
+        "    PROJCS[\"NAD_1983_StatePlane_Alabama_West_FIPS_0102_Feet\",\n"
+        "        GEOGCS[\"NAD83\",\n"
+        "            DATUM[\"North_American_Datum_1983\",\n"
+        "                SPHEROID[\"GRS 1980\",6378137,298.257222101,\n"
+        "                    AUTHORITY[\"EPSG\",\"7019\"]],\n"
+        "                TOWGS84[0,0,0,0,0,0,0],\n"
+        "                AUTHORITY[\"EPSG\",\"6269\"]],\n"
+        "            PRIMEM[\"Greenwich\",0],\n"
+        "            UNIT[\"Degree\",0.0174532925199433]],\n"
+        "        PROJECTION[\"Transverse_Mercator\"],\n"
+        "        PARAMETER[\"latitude_of_origin\",30],\n"
+        "        PARAMETER[\"central_meridian\",-87.5],\n"
+        "        PARAMETER[\"scale_factor\",0.999933333333333],\n"
+        "        PARAMETER[\"false_easting\",1968500],\n"
+        "        PARAMETER[\"false_northing\",0],\n"
+        "        UNIT[\"US survey foot\",0.304800609601219,\n"
+        "            AUTHORITY[\"EPSG\",\"9003\"]],\n"
+        "        AXIS[\"Easting\",EAST],\n"
+        "        AXIS[\"Northing\",NORTH],\n"
+        "        AUTHORITY[\"ESRI\",\"102630\"]],\n"
+        "    VERT_CS[\"NAVD88 height (ftUS)\",\n"
+        "        VERT_DATUM[\"North American Vertical Datum 1988\",2005,\n"
+        "            "
+        "EXTENSION[\"PROJ4_GRIDS\",\"g2012a_alaska.gtx,g2012a_hawaii.gtx,"
+        "g2012a_conus.gtx\"],\n"
+        "            AUTHORITY[\"EPSG\",\"5103\"]],\n"
+        "        UNIT[\"US survey foot\",0.304800609601219,\n"
+        "            AUTHORITY[\"EPSG\",\"9003\"]],\n"
+        "        AXIS[\"Gravity-related height\",UP],\n"
+        "        AUTHORITY[\"EPSG\",\"6360\"]]]");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+
+    auto dbContext = DatabaseContext::create();
+    auto authFactory = AuthorityFactory::create(dbContext, "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(
+        NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst), ctxt);
+    ASSERT_GE(list.size(), 1U);
+    EXPECT_EQ(list[0]->nameStr(),
+              "Inverse of unnamed + "
+              "Transformation from NAD83 to WGS84 + "
+              "NAVD88 height to NAVD88 height (ftUS) + "
+              "Inverse of Transformation from NAD83 to WGS84 + "
+              "unnamed");
+    auto grids = list[0]->gridsNeeded(dbContext, false);
+    EXPECT_TRUE(grids.empty());
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +inv +proj=tmerc +lat_0=30 +lon_0=-87.5 +k=0.999933333 "
+              "+x_0=600000 +y_0=0 +ellps=GRS80 "
+              "+step +proj=unitconvert +z_in=m +z_out=us-ft "
+              "+step +proj=tmerc +lat_0=30 +lon_0=-87.5 +k=0.999933333333333 "
+              "+x_0=600000 +y_0=0 +ellps=GRS80 "
+              "+step +proj=unitconvert +xy_in=m +xy_out=us-ft");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, compoundCRS_to_compoundCRS_context) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), "EPSG");