Merge pull request #2265 from rouault/fix_createOperations_bound_to_geod

createOperations(): fix transformations between BoundCRS and Geocentric CRS

3 files changed, 66 insertions, 12 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index e33ad15e..482e6e6d 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -9084,6 +9084,12 @@ extractGeographicCRSIfGeographicCRSOrEquivalent(const crs::CRSNNPtr &crs) {
                     }
                 }
             }
+        } else {
+            auto boundCRS = util::nn_dynamic_pointer_cast<crs::BoundCRS>(crs);
+            if (boundCRS) {
+                geogCRS = util::nn_dynamic_pointer_cast<crs::GeographicCRS>(
+                    boundCRS->baseCRS());
+            }
         }
     }
     return geogCRS;
@@ -11096,10 +11102,9 @@ struct CoordinateOperationFactory::Private {
         const crs::GeographicCRS *geogDst,
         std::vector<CoordinateOperationNNPtr> &res);
 
-    static void createOperationsCompoundToGeod(
+    static void createOperationsToGeod(
         const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
-        Private::Context &context, const crs::CompoundCRS *compoundSrc,
-        const crs::GeodeticCRS *geodDst,
+        Private::Context &context, const crs::GeodeticCRS *geodDst,
         std::vector<CoordinateOperationNNPtr> &res);
 
     static void createOperationsCompoundToCompound(
@@ -13295,14 +13300,18 @@ CoordinateOperationFactory::Private::createOperations(
         return applyInverse(createOperations(targetCRS, sourceCRS, context));
     }
 
+    // Order of comparison between the geogDst vs geodDst is impotant
     if (boundSrc && geogDst) {
         createOperationsBoundToGeog(sourceCRS, targetCRS, context, boundSrc,
                                     geogDst, res);
         return res;
+    } else if (boundSrc && geodDst) {
+        createOperationsToGeod(sourceCRS, targetCRS, context, geodDst, res);
+        return res;
     }
 
     // reverse of previous case
-    if (geogSrc && boundDst) {
+    if (geodSrc && boundDst) {
         return applyInverse(createOperations(targetCRS, sourceCRS, context));
     }
 
@@ -13352,8 +13361,7 @@ CoordinateOperationFactory::Private::createOperations(
                                        compoundSrc, geogDst, res);
         return res;
     } else if (compoundSrc && geodDst) {
-        createOperationsCompoundToGeod(sourceCRS, targetCRS, context,
-                                       compoundSrc, geodDst, res);
+        createOperationsToGeod(sourceCRS, targetCRS, context, geodDst, res);
         return res;
     }
 
@@ -15057,10 +15065,9 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
 
 // ---------------------------------------------------------------------------
 
-void CoordinateOperationFactory::Private::createOperationsCompoundToGeod(
+void CoordinateOperationFactory::Private::createOperationsToGeod(
     const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
-    Private::Context &context, const crs::CompoundCRS * /*compoundSrc*/,
-    const crs::GeodeticCRS *geodDst,
+    Private::Context &context, const crs::GeodeticCRS *geodDst,
     std::vector<CoordinateOperationNNPtr> &res) {
     auto datum = geodDst->datum();
     if (datum) {
@@ -15079,8 +15086,11 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeod(
             createOperations(intermGeog3DCRS, targetCRS, context);
         if (!geog3DToTargetOps.empty()) {
             for (const auto &op : sourceToGeog3DOps) {
+                auto newOp = op->shallowClone();
+                setCRSs(newOp.get(), sourceCRS, intermGeog3DCRS);
                 res.emplace_back(ConcatenatedOperation::createComputeMetadata(
-                    {op, geog3DToTargetOps.front()}, !allowEmptyIntersection));
+                    {newOp, geog3DToTargetOps.front()},
+                    !allowEmptyIntersection));
             }
         }
     }
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 2847d88d..a1918fb7 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -2617,8 +2617,7 @@ TEST_F(CApi, proj_context_get_database_metadata) {
     EXPECT_TRUE(proj_context_get_database_metadata(m_ctxt, "IGNF.VERSION") !=
                 nullptr);
 
-    EXPECT_TRUE(proj_context_get_database_metadata(m_ctxt, "FOO") ==
-                nullptr);
+    EXPECT_TRUE(proj_context_get_database_metadata(m_ctxt, "FOO") == nullptr);
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 6ef6e686..9ebb3794 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -6251,6 +6251,51 @@ TEST(operation, boundCRS_of_geogCRS_to_geogCRS) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, boundCRS_of_geogCRS_to_geodCRS) {
+    auto boundCRS = BoundCRS::createFromTOWGS84(
+        GeographicCRS::EPSG_4807, std::vector<double>{1, 2, 3, 4, 5, 6, 7});
+    auto op = CoordinateOperationFactory::create()->createOperation(
+        boundCRS, GeodeticCRS::EPSG_4978);
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline +step "
+              "+proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=grad +xy_out=rad "
+              "+step +inv +proj=longlat +ellps=clrk80ign +pm=paris "
+              "+step +proj=cart +ellps=clrk80ign "
+              "+step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
+              "+convention=position_vector");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, boundCRS_of_geogCRS_to_geodCRS_not_related_to_hub) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto boundCRS = BoundCRS::createFromTOWGS84(
+        GeographicCRS::EPSG_4807, std::vector<double>{1, 2, 3, 4, 5, 6, 7});
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        boundCRS,
+        // ETRS89 geocentric
+        authFactory->createCoordinateReferenceSystem("4936"), ctxt);
+    ASSERT_EQ(list.size(), 1U);
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline +step "
+              "+proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=grad +xy_out=rad "
+              "+step +inv +proj=longlat +ellps=clrk80ign +pm=paris "
+              "+step +proj=push +v_3 "
+              "+step +proj=cart +ellps=clrk80ign "
+              "+step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
+              "+convention=position_vector "
+              "+step +inv +proj=cart +ellps=GRS80 "
+              "+step +proj=pop +v_3 "
+              "+step +proj=cart +ellps=GRS80");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, boundCRS_of_geogCRS_to_geogCRS_with_area) {
     auto boundCRS = BoundCRS::createFromTOWGS84(
         GeographicCRS::EPSG_4267, std::vector<double>{1, 2, 3, 4, 5, 6, 7});