createOperations(): fix inconsistent chaining exception when transfor…

…ming from BoundCRS of projected CRS based on NTF Paris to BoundCRS of geog CRS NTF Paris. Fixes OSGeo/gdal#3273
OSGeo · Dec 14, 2020 · 2968770 · 2968770
1 parent cf43868
commit 2968770
Show file tree

Hide file tree

Showing 5 changed files with 94 additions and 39 deletions.
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
@@ -5188,7 +5188,8 @@ BoundCRSNNPtr BoundCRS::createFromNadgrids(const CRSNNPtr &baseCRSIn,
                         " (with Greenwich prime meridian)"),
                 sourceGeographicCRS->datumNonNull(nullptr)->ellipsoid(),
                 util::optional<std::string>(), datum::PrimeMeridian::GREENWICH),
-            sourceGeographicCRS->coordinateSystem());
+            cs::EllipsoidalCS::createLatitudeLongitude(
+                common::UnitOfMeasure::DEGREE));
     }
     std::string transformationName = transformationSourceCRS->nameStr();
     transformationName += " to WGS84";

diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
@@ -4239,7 +4239,8 @@ createBoundCRSSourceTransformationCRS(const crs::CRSPtr &sourceCRS,
                             sourceGeographicCRS->datum()->ellipsoid(),
                             util::optional<std::string>(),
                             datum::PrimeMeridian::GREENWICH),
-                        sourceGeographicCRS->coordinateSystem())
+                        cs::EllipsoidalCS::createLatitudeLongitude(
+                            common::UnitOfMeasure::DEGREE))
                         .as_nullable();
             }
         } else {

diff --git a/src/iso19111/operation/concatenatedoperation.cpp b/src/iso19111/operation/concatenatedoperation.cpp
@@ -446,6 +446,41 @@ CoordinateOperationNNPtr ConcatenatedOperation::createComputeMetadata(
             flattenOps.emplace_back(subOp);
         }
     }
+
+    // Remove consecutive inverse operations
+    if (flattenOps.size() > 2) {
+        std::vector<size_t> indices;
+        for (size_t i = 0; i < flattenOps.size(); ++i)
+            indices.push_back(i);
+        while (true) {
+            bool bHasChanged = false;
+            for (size_t i = 0; i + 1 < indices.size(); ++i) {
+                if (flattenOps[indices[i]]->_isEquivalentTo(
+                        flattenOps[indices[i + 1]]->inverse().get(),
+                        util::IComparable::Criterion::EQUIVALENT) &&
+                    flattenOps[indices[i]]->sourceCRS()->_isEquivalentTo(
+                        flattenOps[indices[i + 1]]->targetCRS().get(),
+                        util::IComparable::Criterion::EQUIVALENT)) {
+                    indices.erase(indices.begin() + i, indices.begin() + i + 2);
+                    bHasChanged = true;
+                    break;
+                }
+            }
+            // We bail out if indices.size() == 2, because potentially
+            // the last 2 remaining ones could auto-cancel, and we would have
+            // to have a special case for that (and this happens in practice).
+            if (!bHasChanged || indices.size() <= 2)
+                break;
+        }
+        if (indices.size() < flattenOps.size()) {
+            std::vector<CoordinateOperationNNPtr> flattenOpsNew;
+            for (size_t i = 0; i < indices.size(); ++i) {
+                flattenOpsNew.emplace_back(flattenOps[indices[i]]);
+            }
+            flattenOps = std::move(flattenOpsNew);
+        }
+    }
+
     if (flattenOps.size() == 1) {
         return flattenOps[0];
     }

diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp
@@ -4255,47 +4255,26 @@ void CoordinateOperationFactory::Private::createOperationsBoundToBound(
     auto hubSrcGeog = dynamic_cast<const crs::GeographicCRS *>(hubSrc.get());
     const auto &hubDst = boundDst->hubCRS();
     auto hubDstGeog = dynamic_cast<const crs::GeographicCRS *>(hubDst.get());
-    auto geogCRSOfBaseOfBoundSrc = boundSrc->baseCRS()->extractGeographicCRS();
-    auto geogCRSOfBaseOfBoundDst = boundDst->baseCRS()->extractGeographicCRS();
     if (hubSrcGeog && hubDstGeog &&
         hubSrcGeog->_isEquivalentTo(hubDstGeog,
-                                    util::IComparable::Criterion::EQUIVALENT) &&
-        geogCRSOfBaseOfBoundSrc && geogCRSOfBaseOfBoundDst) {
-        const bool firstIsNoOp = geogCRSOfBaseOfBoundSrc->_isEquivalentTo(
-            boundSrc->baseCRS().get(),
-            util::IComparable::Criterion::EQUIVALENT);
-        const bool lastIsNoOp = geogCRSOfBaseOfBoundDst->_isEquivalentTo(
-            boundDst->baseCRS().get(),
-            util::IComparable::Criterion::EQUIVALENT);
-        auto opsFirst = createOperations(
-            boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
-        auto opsLast = createOperations(NN_NO_CHECK(geogCRSOfBaseOfBoundDst),
-                                        boundDst->baseCRS(), context);
-        if (!opsFirst.empty() && !opsLast.empty()) {
-            const auto &opSecond = boundSrc->transformation();
-            auto opThird = boundDst->transformation()->inverse();
-            for (const auto &opFirst : opsFirst) {
-                for (const auto &opLast : opsLast) {
-                    try {
-                        std::vector<CoordinateOperationNNPtr> ops;
-                        if (!firstIsNoOp) {
-                            ops.push_back(opFirst);
-                        }
-                        ops.push_back(opSecond);
-                        ops.push_back(opThird);
-                        if (!lastIsNoOp) {
-                            ops.push_back(opLast);
-                        }
-                        res.emplace_back(
-                            ConcatenatedOperation::createComputeMetadata(
-                                ops, disallowEmptyIntersection));
-                    } catch (const InvalidOperationEmptyIntersection &) {
-                    }
+                                    util::IComparable::Criterion::EQUIVALENT)) {
+        auto opsFirst = createOperations(sourceCRS, hubSrc, context);
+        auto opsLast = createOperations(hubSrc, targetCRS, context);
+        for (const auto &opFirst : opsFirst) {
+            for (const auto &opLast : opsLast) {
+                try {
+                    std::vector<CoordinateOperationNNPtr> ops;
+                    ops.push_back(opFirst);
+                    ops.push_back(opLast);
+                    res.emplace_back(
+                        ConcatenatedOperation::createComputeMetadata(
+                            ops, disallowEmptyIntersection));
+                } catch (const InvalidOperationEmptyIntersection &) {
                 }
             }
-            if (!res.empty()) {
-                return;
-            }
+        }
+        if (!res.empty()) {
+            return;
         }
     }
 

diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
@@ -2501,6 +2501,45 @@ TEST(operation, boundCRS_of_projCRS_towgs84_to_boundCRS_of_projCRS_nadgrids) {
 
 // ---------------------------------------------------------------------------
 
+static CRSNNPtr buildCRSFromProjStrThroughWKT(const std::string &projStr) {
+    auto crsFromProj = nn_dynamic_pointer_cast<CRS>(
+        PROJStringParser().createFromPROJString(projStr));
+    if (crsFromProj == nullptr) {
+        throw "crsFromProj == nullptr";
+    }
+    auto crsFromWkt = nn_dynamic_pointer_cast<CRS>(
+        WKTParser().createFromWKT(crsFromProj->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT2_2019).get())));
+    if (crsFromWkt == nullptr) {
+        throw "crsFromWkt == nullptr";
+    }
+    return NN_NO_CHECK(crsFromWkt);
+}
+
+TEST(operation,
+     boundCRS_to_boundCRS_with_base_geog_crs_different_from_source_of_transf) {
+
+    auto src = buildCRSFromProjStrThroughWKT(
+        "+proj=lcc +lat_1=49 +lat_0=49 +lon_0=0 +k_0=0.999877499 +x_0=600000 "
+        "+y_0=200000 +ellps=clrk80ign +pm=paris +towgs84=-168,-60,320,0,0,0,0 "
+        "+units=m +no_defs +type=crs");
+    auto dst = buildCRSFromProjStrThroughWKT(
+        "+proj=longlat +ellps=clrk80ign +pm=paris "
+        "+towgs84=-168,-60,320,0,0,0,0 +no_defs +type=crs");
+
+    auto op = CoordinateOperationFactory::create()->createOperation(src, dst);
+    ASSERT_TRUE(op != nullptr);
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +inv +proj=lcc +lat_1=49 +lat_0=49 +lon_0=0 "
+              "+k_0=0.999877499 +x_0=600000 +y_0=200000 +ellps=clrk80ign "
+              "+pm=paris "
+              "+step +proj=longlat +ellps=clrk80ign +pm=paris "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, boundCRS_with_basecrs_with_extent_to_geogCRS) {
 
     auto wkt =