createOperations(): fix Compound to Geog3D/Projected3D CRS with non-metre ellipsoidal height

src/iso19111/operation/coordinateoperationfactory.cpp

@@ -3476,8 +3476,8 @@ std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
 
 std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
     createOperationsGeogToVertWithAlternativeGeog(
-        const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS
-        const crs::CRSNNPtr &targetCRS,      // vertical CRS
+        const crs::CRSNNPtr &sourceCRS, // geographic CRS
+        const crs::CRSNNPtr &targetCRS, // vertical CRS
         Private::Context &context) {
 
     ENTER_FUNCTION();
@@ -3501,10 +3501,39 @@ std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
     // Generally EPSG has operations from GeogCrs to VertCRS
     auto ops = findOpsInRegistryDirectTo(targetCRS, context);
 
+    const auto geogCRS =
+        dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
+    assert(geogCRS);
+    const auto &srcAxisList = geogCRS->coordinateSystem()->axisList();
     for (const auto &op : ops) {
-        const auto tmpCRS = op->sourceCRS();
-        if (tmpCRS && dynamic_cast<const crs::GeographicCRS *>(tmpCRS.get())) {
-            res.emplace_back(op);
+        const auto tmpCRS =
+            dynamic_cast<const crs::GeographicCRS *>(op->sourceCRS().get());
+        if (tmpCRS) {
+            if (srcAxisList.size() == 3 &&
+                srcAxisList[2]->unit().conversionToSI() != 1) {
+
+                const auto &authFactory =
+                    context.context->getAuthorityFactory();
+                const auto dbContext =
+                    authFactory->databaseContext().as_nullable();
+                auto tmpCRSWithSrcZ =
+                    tmpCRS->demoteTo2D(std::string(), dbContext)
+                        ->promoteTo3D(std::string(), dbContext, srcAxisList[2]);
+
+                std::vector<CoordinateOperationNNPtr> opsUnitConvert;
+                createOperationsGeogToGeog(
+                    opsUnitConvert, tmpCRSWithSrcZ,
+                    NN_NO_CHECK(op->sourceCRS()), context,
+                    dynamic_cast<const crs::GeographicCRS *>(
+                        tmpCRSWithSrcZ.get()),
+                    tmpCRS);
+                assert(opsUnitConvert.size() == 1);
+                auto concat = ConcatenatedOperation::createComputeMetadata(
+                    {opsUnitConvert.front(), op}, disallowEmptyIntersection);
+                res.emplace_back(concat);
+            } else {
+                res.emplace_back(op);
+            }
         }
     }
 
@@ -4564,11 +4593,18 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
                 !srcGeogCRS->_isEquivalentTo(
                     geogDst, util::IComparable::Criterion::EQUIVALENT) &&
                 !srcGeogCRS->is2DPartOf3D(NN_NO_CHECK(geogDst), dbContext)) {
-                auto verticalTransformsTmp = createOperations(
-                    componentsSrc[1],
+                auto geogCRSTmp =
                     NN_NO_CHECK(srcGeogCRS)
-                        ->promoteTo3D(std::string(), dbContext),
-                    context);
+                        ->demoteTo2D(std::string(), dbContext)
+                        ->promoteTo3D(
+                            std::string(), dbContext,
+                            geogDst->coordinateSystem()->axisList().size() == 3
+                                ? geogDst->coordinateSystem()->axisList()[2]
+                                : cs::VerticalCS::createGravityRelatedHeight(
+                                      common::UnitOfMeasure::METRE)
+                                      ->axisList()[0]);
+                auto verticalTransformsTmp =
+                    createOperations(componentsSrc[1], geogCRSTmp, context);
                 bool foundRegisteredTransform = false;
                 foundRegisteredTransformWithAllGridsAvailable = false;
                 for (const auto &op : verticalTransformsTmp) {

test/unit/test_operationfactory.cpp

@@ -5867,6 +5867,38 @@ TEST(operation, createOperation_on_crs_with_bound_crs_and_wktext) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, compoundCRS_to_proj_string_with_non_metre_height) {
+    auto objSrc =
+        createFromUserInput("EPSG:6318+5703", DatabaseContext::create(), false);
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +ellps=GRS80 +vunits=us-ft +type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    auto list =
+        CoordinateOperationFactory::create()->createOperations(
+            NN_NO_CHECK(src), NN_NO_CHECK(dst), ctxt);
+    ASSERT_GT(list.size(), 1U);
+    // What is important to check here is the vertical unit conversion
+    EXPECT_EQ(
+        list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+        "+proj=pipeline "
+        "+step +proj=axisswap +order=2,1 "
+        "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+        "+step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 "
+        "+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=us-ft");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, createOperation_ossfuzz_18587) {
     auto objSrc =
         createFromUserInput("EPSG:4326", DatabaseContext::create(), false);