createOperations(): fix inconsistent chaining exception when transforming from BoundCRS of projected CRS based on NTF Paris to BoundCRS of geog CRS NTF Paris. Fixes https://github.com/OSGeo/gdal/issues/3273

1 files changed, 35 insertions, 0 deletions

diff --git a/src/iso19111/operation/concatenatedoperation.cpp b/src/iso19111/operation/concatenatedoperation.cpp
index 1c65a24b..ce4b015a 100644
--- 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];
     }