Merge pull request #1326 from rouault/fix_1323

createOperation(): fix geocent <--> nadgrids+geoidgrids case (fixes #1323)

3 files changed, 274 insertions, 3 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 575eae39..7ee20758 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -11879,6 +11879,34 @@ CoordinateOperationFactory::Private::createOperations(
             }
         }
 
+        auto vertCRSOfBaseOfBoundSrc =
+            boundSrc->baseCRS()->extractVerticalCRS();
+        auto vertCRSOfBaseOfBoundDst =
+            boundDst->baseCRS()->extractVerticalCRS();
+        if (hubSrcGeog && hubDstGeog &&
+            hubSrcGeog->_isEquivalentTo(
+                hubDstGeog, util::IComparable::Criterion::EQUIVALENT) &&
+            vertCRSOfBaseOfBoundSrc && vertCRSOfBaseOfBoundDst) {
+            auto opsFirst = createOperations(sourceCRS, hubSrc, context);
+            auto opsLast = createOperations(hubSrc, targetCRS, context);
+            if (!opsFirst.empty() && !opsLast.empty()) {
+                for (const auto &opFirst : opsFirst) {
+                    for (const auto &opLast : opsLast) {
+                        try {
+                            res.emplace_back(
+                                ConcatenatedOperation::createComputeMetadata(
+                                    {opFirst, opLast},
+                                    !allowEmptyIntersection));
+                        } catch (const InvalidOperationEmptyIntersection &) {
+                        }
+                    }
+                }
+                if (!res.empty()) {
+                    return res;
+                }
+            }
+        }
+
         return createOperations(boundSrc->baseCRS(), boundDst->baseCRS(),
                                 context);
     }
@@ -11914,11 +11942,40 @@ CoordinateOperationFactory::Private::createOperations(
                 return horizTransforms;
             }
         }
+    } else if (compoundSrc && geodDst) {
+        auto datum = geodDst->datum();
+        if (datum) {
+            auto cs =
+                cs::EllipsoidalCS::createLatitudeLongitudeEllipsoidalHeight(
+                    common::UnitOfMeasure::DEGREE,
+                    common::UnitOfMeasure::METRE);
+            auto intermGeog3DCRS = util::nn_static_pointer_cast<crs::CRS>(
+                crs::GeographicCRS::create(
+                    util::PropertyMap()
+                        .set(common::IdentifiedObject::NAME_KEY,
+                             geodDst->nameStr())
+                        .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
+                             metadata::Extent::WORLD),
+                    NN_NO_CHECK(datum), cs));
+            auto sourceToGeog3DOps =
+                createOperations(sourceCRS, intermGeog3DCRS, context);
+            auto geog3DToTargetOps =
+                createOperations(intermGeog3DCRS, targetCRS, context);
+            if (!geog3DToTargetOps.empty()) {
+                for (const auto &op : sourceToGeog3DOps) {
+                    res.emplace_back(
+                        ConcatenatedOperation::createComputeMetadata(
+                            {op, geog3DToTargetOps.front()},
+                            !allowEmptyIntersection));
+                }
+                return res;
+            }
+        }
     }
 
     // reverse of previous case
     auto compoundDst = dynamic_cast<const crs::CompoundCRS *>(targetCRS.get());
-    if (geogSrc && compoundDst) {
+    if (geodSrc && compoundDst) {
         return applyInverse(createOperations(targetCRS, sourceCRS, context));
     }
 
@@ -11958,6 +12015,21 @@ CoordinateOperationFactory::Private::createOperations(
                                         nn_interpTransformCRS));
                             }
                         }
+                    } else {
+                        auto compSrc0BoundCrs = dynamic_cast<crs::BoundCRS *>(
+                            componentsSrc[0].get());
+                        auto compDst0BoundCrs = dynamic_cast<crs::BoundCRS *>(
+                            componentsDst[0].get());
+                        if (compSrc0BoundCrs && compDst0BoundCrs &&
+                            dynamic_cast<crs::GeographicCRS *>(
+                                compSrc0BoundCrs->hubCRS().get()) &&
+                            compSrc0BoundCrs->hubCRS()->_isEquivalentTo(
+                                compDst0BoundCrs->hubCRS().get())) {
+                            interpolationGeogCRS =
+                                NN_NO_CHECK(util::nn_dynamic_pointer_cast<
+                                            crs::GeographicCRS>(
+                                    compSrc0BoundCrs->hubCRS()));
+                        }
                     }
                     auto opSrcCRSToGeogCRS = createOperations(
                         componentsSrc[0], interpolationGeogCRS, context);
@@ -11982,6 +12054,89 @@ CoordinateOperationFactory::Private::createOperations(
         }
     }
 
+    // '+proj=longlat +ellps=GRS67 +nadgrids=@foo.gsb +type=crs' to
+    // '+proj=longlat +ellps=GRS80 +nadgrids=@bar.gsb +geoidgrids=@bar.gtx
+    // +type=crs'
+    if (boundSrc && compoundDst) {
+        const auto &componentsDst = compoundDst->componentReferenceSystems();
+        if (!componentsDst.empty()) {
+            auto compDst0BoundCrs =
+                dynamic_cast<crs::BoundCRS *>(componentsDst[0].get());
+            if (compDst0BoundCrs) {
+                auto boundSrcHubAsGeogCRS = dynamic_cast<crs::GeographicCRS *>(
+                    boundSrc->hubCRS().get());
+                auto compDst0BoundCrsHubAsGeogCRS =
+                    dynamic_cast<crs::GeographicCRS *>(
+                        compDst0BoundCrs->hubCRS().get());
+                if (boundSrcHubAsGeogCRS && compDst0BoundCrsHubAsGeogCRS) {
+                    const auto &boundSrcHubAsGeogCRSDatum =
+                        boundSrcHubAsGeogCRS->datum();
+                    const auto &compDst0BoundCrsHubAsGeogCRSDatum =
+                        compDst0BoundCrsHubAsGeogCRS->datum();
+                    if (boundSrcHubAsGeogCRSDatum &&
+                        compDst0BoundCrsHubAsGeogCRSDatum &&
+                        boundSrcHubAsGeogCRSDatum->_isEquivalentTo(
+                            compDst0BoundCrsHubAsGeogCRSDatum.get())) {
+                        auto cs = cs::EllipsoidalCS::
+                            createLatitudeLongitudeEllipsoidalHeight(
+                                common::UnitOfMeasure::DEGREE,
+                                common::UnitOfMeasure::METRE);
+                        auto intermGeog3DCRS = util::nn_static_pointer_cast<
+                            crs::CRS>(crs::GeographicCRS::create(
+                            util::PropertyMap()
+                                .set(common::IdentifiedObject::NAME_KEY,
+                                     boundSrcHubAsGeogCRS->nameStr())
+                                .set(
+                                    common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
+                                    metadata::Extent::WORLD),
+                            NN_NO_CHECK(boundSrcHubAsGeogCRSDatum), cs));
+                        auto sourceToGeog3DOps = createOperations(
+                            sourceCRS, intermGeog3DCRS, context);
+                        auto geog3DToTargetOps = createOperations(
+                            intermGeog3DCRS, targetCRS, context);
+                        for (const auto &opSrc : sourceToGeog3DOps) {
+                            for (const auto &opDst : geog3DToTargetOps) {
+                                if (opSrc->targetCRS() && opDst->sourceCRS() &&
+                                    !opSrc->targetCRS()->_isEquivalentTo(
+                                        opDst->sourceCRS().get())) {
+                                    // Shouldn't happen normally, but typically
+                                    // one of them can be 2D and the other 3D
+                                    // due to above createOperations() not
+                                    // exactly setting the expected source and
+                                    // target CRS.
+                                    // So create an adapter operation...
+                                    auto intermOps = createOperations(
+                                        NN_NO_CHECK(opSrc->targetCRS()),
+                                        NN_NO_CHECK(opDst->sourceCRS()),
+                                        context);
+                                    if (!intermOps.empty()) {
+                                        res.emplace_back(
+                                            ConcatenatedOperation::
+                                                createComputeMetadata(
+                                                    {opSrc, intermOps.front(),
+                                                     opDst},
+                                                    !allowEmptyIntersection));
+                                    }
+                                } else {
+                                    res.emplace_back(
+                                        ConcatenatedOperation::
+                                            createComputeMetadata(
+                                                {opSrc, opDst},
+                                                !allowEmptyIntersection));
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    // reverse of previous case
+    if (boundDst && compoundSrc) {
+        return applyInverse(createOperations(targetCRS, sourceCRS, context));
+    }
+
     return res;
 }
 //! @endcond
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 51f93e2f..d5f38fd5 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -5049,7 +5049,7 @@ const std::string &PROJStringFormatter::toString() const {
                     first.paramValues[1].keyEquals("z_in") &&
                     first.paramValues[2].keyEquals("xy_out") &&
                     first.paramValues[3].keyEquals("z_out") &&
-                    second.paramValues[0].keyEquals("xy_in=") &&
+                    second.paramValues[0].keyEquals("xy_in") &&
                     second.paramValues[1].keyEquals("xy_out") &&
                     first.paramValues[0].value == second.paramValues[1].value &&
                     first.paramValues[2].value == second.paramValues[0].value) {
@@ -5074,6 +5074,22 @@ const std::string &PROJStringFormatter::toString() const {
                 break;
             }
 
+            // unitconvert (1), axisswap order=2,1, unitconvert(2)  ==>
+            // axisswap order=2,1, unitconvert (1), unitconvert(2) which
+            // will get further optimized by previous case
+            if (i + 1 < d->steps_.size() && prevStep.name == "unitconvert" &&
+                curStep.name == "axisswap" && curStepParamCount == 1 &&
+                curStep.paramValues[0].equals("order", "2,1")) {
+                auto iterNext = iterCur;
+                ++iterNext;
+                auto &nextStep = *iterNext;
+                if (nextStep.name == "unitconvert") {
+                    std::swap(*iterPrev, *iterCur);
+                    changeDone = true;
+                    break;
+                }
+            }
+
             // axisswap order=2,1 followed by itself is a no-op
             if (curStep.name == "axisswap" && prevStep.name == "axisswap" &&
                 curStepParamCount == 1 && prevStepParamCount == 1 &&
@@ -6707,7 +6723,7 @@ PROJStringParser::Private::buildGeocentricCRS(int iStep, int iUnitConvert) {
 
     auto datum = buildDatum(step, title);
 
-    UnitOfMeasure unit = UnitOfMeasure::METRE;
+    UnitOfMeasure unit = buildUnit(step, "units", "");
     if (iUnitConvert >= 0) {
         auto &stepUnitConvert = steps_[iUnitConvert];
         const std::string *xy_in = &getParamValue(stepUnitConvert, "xy_in");
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index e71ac716..9111b862 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -6026,6 +6026,51 @@ TEST(operation, compoundCRS_with_boundProjCRS_and_boundVerticalCRS_to_geogCRS) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, geocent_to_compoundCRS) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=geocent +datum=WGS84 +units=m +type=crs");
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +ellps=GRS67 +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=pipeline +step +inv +proj=cart +ellps=WGS84 +step +inv "
+              "+proj=vgridshift +grids=@foo.gtx +multiplier=1 +step +inv "
+              "+proj=hgridshift +grids=@foo.gsb +step +proj=unitconvert "
+              "+xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, geocent_to_compoundCRS_context) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    // WGS84 geocentric
+    auto src = authFactory->createCoordinateReferenceSystem("4978");
+    auto objDst = PROJStringParser().createFromPROJString(
+        "+proj=longlat +ellps=GRS67 +nadgrids=@foo.gsb +geoidgrids=@foo.gtx "
+        "+type=crs");
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
+    ASSERT_TRUE(dst != nullptr);
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        src, NN_CHECK_ASSERT(dst), ctxt);
+    ASSERT_EQ(list.size(), 1U);
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline +step +inv +proj=cart +ellps=WGS84 +step +inv "
+              "+proj=vgridshift +grids=@foo.gtx +multiplier=1 +step +inv "
+              "+proj=hgridshift +grids=@foo.gsb +step +proj=unitconvert "
+              "+xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, compoundCRS_to_compoundCRS) {
     auto compound1 = CompoundCRS::create(
         PropertyMap(),
@@ -6116,6 +6161,31 @@ TEST(operation, compoundCRS_to_compoundCRS_with_vertical_transform) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, compoundCRS_to_compoundCRS_with_bound_crs_in_horiz_and_vert) {
+    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=@bar.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 +proj=vgridshift +grids=@foo.gtx +multiplier=1 "
+              "+step +inv +proj=vgridshift +grids=@bar.gtx +multiplier=1 "
+              "+step +inv +proj=hgridshift +grids=@bar.gsb "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, compoundCRS_to_compoundCRS_context) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), "EPSG");
@@ -6336,6 +6406,36 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, boundCRS_to_compoundCRS) {
+    auto objSrc = PROJStringParser().createFromPROJString(
+        "+proj=longlat +ellps=GRS67 +nadgrids=@foo.gsb +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=@bar.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=vgridshift +grids=@bar.gtx +multiplier=1 "
+              "+step +inv +proj=hgridshift +grids=@bar.gsb "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+
+    auto opInverse = CoordinateOperationFactory::create()->createOperation(
+        NN_CHECK_ASSERT(dst), NN_CHECK_ASSERT(src));
+    ASSERT_TRUE(opInverse != nullptr);
+    EXPECT_TRUE(opInverse->inverse()->_isEquivalentTo(op.get()));
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, IGNF_LAMB1_TO_EPSG_4326) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), std::string());