WKT CONCATENATEDOPERATION parsing: allow CONVERSION steps and reverse operations when neededs (fixes #1197)

3 files changed, 215 insertions, 132 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index c531a8bc..d00be93f 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -8794,6 +8794,21 @@ ConcatenatedOperation::operations() const {
 
 // ---------------------------------------------------------------------------
 
+//! @cond Doxygen_Suppress
+static bool compareStepCRS(const crs::CRS *a, const crs::CRS *b) {
+    const auto &aIds = a->identifiers();
+    const auto &bIds = b->identifiers();
+    if (aIds.size() == 1 && bIds.size() == 1 &&
+        aIds[0]->code() == bIds[0]->code() &&
+        *aIds[0]->codeSpace() == *bIds[0]->codeSpace()) {
+        return true;
+    }
+    return a->_isEquivalentTo(b, util::IComparable::Criterion::EQUIVALENT);
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
 /** \brief Instantiate a ConcatenatedOperation
  *
  * @param properties See \ref general_properties. At minimum the name should
@@ -8823,16 +8838,7 @@ ConcatenatedOperationNNPtr ConcatenatedOperation::create(
                                    "source and/or target CRS");
         }
         if (i >= 1) {
-            const auto &sourceCRSIds = l_sourceCRS->identifiers();
-            const auto &targetCRSIds = lastTargetCRS->identifiers();
-            if (sourceCRSIds.size() == 1 && targetCRSIds.size() == 1 &&
-                sourceCRSIds[0]->code() == targetCRSIds[0]->code() &&
-                *sourceCRSIds[0]->codeSpace() ==
-                    *targetCRSIds[0]->codeSpace()) {
-                // same id --> ok
-            } else if (!l_sourceCRS->_isEquivalentTo(
-                           lastTargetCRS.get(),
-                           util::IComparable::Criterion::EQUIVALENT)) {
+            if (!compareStepCRS(l_sourceCRS.get(), lastTargetCRS.get())) {
                 throw InvalidOperation(
                     "Inconsistent chaining of CRS in operations");
             }
@@ -8852,6 +8858,137 @@ ConcatenatedOperationNNPtr ConcatenatedOperation::create(
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+void ConcatenatedOperation::fixStepsDirection(
+    const crs::CRSNNPtr &concatOpSourceCRS,
+    const crs::CRSNNPtr &concatOpTargetCRS,
+    std::vector<CoordinateOperationNNPtr> &operationsInOut) {
+
+    // Set of heuristics to assign CRS to steps, and possibly reverse them.
+
+    for (size_t i = 0; i < operationsInOut.size(); ++i) {
+        auto &op = operationsInOut[i];
+        auto l_sourceCRS = op->sourceCRS();
+        auto l_targetCRS = op->targetCRS();
+        auto conv = dynamic_cast<const Conversion *>(op.get());
+        if (conv && i == 0 && !l_sourceCRS && !l_targetCRS) {
+            auto derivedCRS =
+                dynamic_cast<const crs::DerivedCRS *>(concatOpSourceCRS.get());
+            if (derivedCRS) {
+                if (i + 1 < operationsInOut.size()) {
+                    // use the sourceCRS of the next operation as our target CRS
+                    l_targetCRS = operationsInOut[i + 1]->sourceCRS();
+                    // except if it looks like the next operation should
+                    // actually be reversed !!!
+                    if (l_targetCRS &&
+                        !compareStepCRS(l_targetCRS.get(),
+                                        derivedCRS->baseCRS().get()) &&
+                        operationsInOut[i + 1]->targetCRS() &&
+                        compareStepCRS(
+                            operationsInOut[i + 1]->targetCRS().get(),
+                            derivedCRS->baseCRS().get())) {
+                        l_targetCRS = operationsInOut[i + 1]->targetCRS();
+                    }
+                }
+                if (!l_targetCRS) {
+                    l_targetCRS = derivedCRS->baseCRS().as_nullable();
+                }
+                auto invConv =
+                    util::nn_dynamic_pointer_cast<InverseConversion>(op);
+                auto nn_targetCRS = NN_NO_CHECK(l_targetCRS);
+                if (invConv) {
+                    invConv->inverse()->setCRSs(nn_targetCRS, concatOpSourceCRS,
+                                                nullptr);
+                    op->setCRSs(concatOpSourceCRS, nn_targetCRS, nullptr);
+                } else {
+                    op->setCRSs(nn_targetCRS, concatOpSourceCRS, nullptr);
+                    op = op->inverse();
+                }
+            } else if (i + 1 < operationsInOut.size()) {
+                l_targetCRS = operationsInOut[i + 1]->sourceCRS();
+                if (l_targetCRS) {
+                    op->setCRSs(concatOpSourceCRS, NN_NO_CHECK(l_targetCRS),
+                                nullptr);
+                }
+            }
+        } else if (conv && i + 1 == operationsInOut.size() && !l_sourceCRS &&
+                   !l_targetCRS) {
+            auto derivedCRS =
+                dynamic_cast<const crs::DerivedCRS *>(concatOpTargetCRS.get());
+            if (derivedCRS) {
+                if (i >= 1) {
+                    // use the sourceCRS of the previous operation as our source
+                    // CRS
+                    l_sourceCRS = operationsInOut[i - 1]->targetCRS();
+                    // except if it looks like the previous operation should
+                    // actually be reversed !!!
+                    if (l_sourceCRS &&
+                        !compareStepCRS(l_sourceCRS.get(),
+                                        derivedCRS->baseCRS().get()) &&
+                        operationsInOut[i - 1]->sourceCRS() &&
+                        compareStepCRS(
+                            operationsInOut[i - 1]->sourceCRS().get(),
+                            derivedCRS->baseCRS().get())) {
+                        l_targetCRS = operationsInOut[i - 1]->sourceCRS();
+                    }
+                }
+                if (!l_sourceCRS) {
+                    l_sourceCRS = derivedCRS->baseCRS().as_nullable();
+                }
+                op->setCRSs(NN_NO_CHECK(l_sourceCRS), concatOpTargetCRS,
+                            nullptr);
+            } else if (i >= 1) {
+                l_sourceCRS = operationsInOut[i - 1]->targetCRS();
+                if (l_sourceCRS) {
+                    op->setCRSs(NN_NO_CHECK(l_sourceCRS), concatOpTargetCRS,
+                                nullptr);
+                }
+            }
+        } else if (conv && i > 0 && i < operationsInOut.size() - 1) {
+            // For an intermediate conversion, use the target CRS of the
+            // previous step and the source CRS of the next step
+            l_sourceCRS = operationsInOut[i - 1]->targetCRS();
+            l_targetCRS = operationsInOut[i + 1]->sourceCRS();
+            if (l_sourceCRS && l_targetCRS) {
+                op->setCRSs(NN_NO_CHECK(l_sourceCRS), NN_NO_CHECK(l_targetCRS),
+                            nullptr);
+            }
+        } else if (!conv && l_sourceCRS && l_targetCRS) {
+            // Transformations might be mentioned in their foward directions,
+            // whereas we should instead use the reverse path.
+            auto prevOpTarget = (i == 0) ? concatOpSourceCRS.as_nullable()
+                                         : operationsInOut[i - 1]->targetCRS();
+            if (!compareStepCRS(l_sourceCRS.get(), prevOpTarget.get()) &&
+                compareStepCRS(l_targetCRS.get(), prevOpTarget.get())) {
+                op = op->inverse();
+            }
+        }
+    }
+
+    if (!operationsInOut.empty()) {
+        auto l_sourceCRS = operationsInOut.front()->sourceCRS();
+        if (l_sourceCRS &&
+            !compareStepCRS(l_sourceCRS.get(), concatOpSourceCRS.get())) {
+            throw InvalidOperation("The source CRS of the first step of "
+                                   "concatenated operation is not the same "
+                                   "as the source CRS of the concantenated "
+                                   "operation itself");
+        }
+
+        auto l_targetCRS = operationsInOut.back()->targetCRS();
+        if (l_targetCRS &&
+            !compareStepCRS(l_targetCRS.get(), concatOpTargetCRS.get())) {
+            throw InvalidOperation("The target CRS of the last step of "
+                                   "concatenated operation is not the same "
+                                   "as the target CRS of the concantenated "
+                                   "operation itself");
+        }
+    }
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 static std::string computeConcatenatedName(
     const std::vector<CoordinateOperationNNPtr> &flattenOps) {
     std::string name;
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index aaf8d109..8454966a 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -3126,117 +3126,12 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation(
                                                     std::string()));
             }
 
-            // In case the operation is a conversion (we hope this is the
-            // forward case!)
-            if (!operations[0]->sourceCRS() || !operations[0]->targetCRS()) {
-                if (!operations[1]->sourceCRS()) {
-                    throw FactoryException(
-                        "chaining of conversion not supported");
-                }
-                operations[0]->setCRSs(
-                    d->createFactory(source_crs_auth_name)
-                        ->createCoordinateReferenceSystem(source_crs_code),
-                    NN_NO_CHECK(operations[1]->sourceCRS()), nullptr);
-            }
-
-            // Some concatenated operations, like 8443, might actually chain
-            // reverse operations rather than forward operations.
-            {
-                const auto &op0SrcId =
-                    operations[0]->sourceCRS()->identifiers()[0];
-                if (op0SrcId->code() != source_crs_code ||
-                    *op0SrcId->codeSpace() != source_crs_auth_name) {
-                    operations[0] = operations[0]->inverse();
-                }
-            }
-
-            {
-                const auto &op0SrcId =
-                    operations[0]->sourceCRS()->identifiers()[0];
-                if (op0SrcId->code() != source_crs_code ||
-                    *op0SrcId->codeSpace() != source_crs_auth_name) {
-                    throw FactoryException(
-                        "Source CRS of first operation in concatenated "
-                        "operation " +
-                        code + " does not match source CRS of "
-                               "concatenated operation");
-                }
-            }
-
-            // In case the operation is a conversion (we hope this is the
-            // forward case!)
-            if (!operations[1]->sourceCRS() || !operations[1]->targetCRS()) {
-                if (step3_auth_name.empty()) {
-                    operations[1]->setCRSs(
-                        NN_NO_CHECK(operations[0]->targetCRS()),
-                        d->createFactory(target_crs_auth_name)
-                            ->createCoordinateReferenceSystem(target_crs_code),
-                        nullptr);
-                } else {
-                    if (!operations[2]->sourceCRS()) {
-                        throw FactoryException(
-                            "chaining of conversion not supported");
-                    }
-                    operations[1]->setCRSs(
-                        NN_NO_CHECK(operations[0]->targetCRS()),
-                        NN_NO_CHECK(operations[2]->sourceCRS()), nullptr);
-                }
-            }
-
-            const auto &op1SrcId = operations[1]->sourceCRS()->identifiers()[0];
-            const auto &op0TargetId =
-                operations[0]->targetCRS()->identifiers()[0];
-            while (true) {
-                if (step3_auth_name.empty()) {
-                    const auto &opLastTargetId =
-                        operations.back()->targetCRS()->identifiers()[0];
-                    if (opLastTargetId->code() == target_crs_code &&
-                        *opLastTargetId->codeSpace() == target_crs_auth_name) {
-                        // in case we have only 2 steps, and
-                        // step2.targetCRS == concatenate.targetCRS do nothing,
-                        // but ConcatenatedOperation::create() will ultimately
-                        // check that step1.targetCRS == step2.sourceCRS
-                        break;
-                    }
-                }
-                if (op1SrcId->code() != op0TargetId->code() ||
-                    *op1SrcId->codeSpace() != *op0TargetId->codeSpace()) {
-                    operations[1] = operations[1]->inverse();
-                }
-                break;
-            }
-
-            if (!step3_auth_name.empty()) {
-
-                const auto &op2Src = operations[2]->sourceCRS();
-                // In case the operation is a conversion (we hope this is the
-                // forward case!)
-                if (!op2Src || !operations[2]->targetCRS()) {
-                    operations[2]->setCRSs(
-                        NN_NO_CHECK(operations[1]->targetCRS()),
-                        d->createFactory(target_crs_auth_name)
-                            ->createCoordinateReferenceSystem(target_crs_code),
-                        nullptr);
-                }
-
-                const auto &op2SrcId = op2Src->identifiers()[0];
-                const auto &op1TargetId =
-                    operations[1]->targetCRS()->identifiers()[0];
-                if (op2SrcId->code() != op1TargetId->code() ||
-                    *op2SrcId->codeSpace() != *op1TargetId->codeSpace()) {
-                    operations[2] = operations[2]->inverse();
-                }
-            }
-
-            const auto &opLastTargetId =
-                operations.back()->targetCRS()->identifiers()[0];
-            if (opLastTargetId->code() != target_crs_code ||
-                *opLastTargetId->codeSpace() != target_crs_auth_name) {
-                throw FactoryException(
-                    "Target CRS of last operation in concatenated operation " +
-                    code +
-                    " doest not match target CRS of concatenated operation");
-            }
+            operation::ConcatenatedOperation::fixStepsDirection(
+                d->createFactory(source_crs_auth_name)
+                    ->createCoordinateReferenceSystem(source_crs_code),
+                d->createFactory(target_crs_auth_name)
+                    ->createCoordinateReferenceSystem(target_crs_code),
+                operations);
 
             auto props =
                 d->createProperties(code, name, deprecated,
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 4ca5a7f0..6a2c3e1a 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -1153,9 +1153,11 @@ struct WKTParser::Private {
 
     BaseObjectNNPtr build(const WKTNodeNNPtr &node);
 
-    IdentifierPtr buildId(const WKTNodeNNPtr &node, bool tolerant = true);
+    IdentifierPtr buildId(const WKTNodeNNPtr &node, bool tolerant,
+                          bool removeInverseOf);
 
-    PropertyMap &buildProperties(const WKTNodeNNPtr &node);
+    PropertyMap &buildProperties(const WKTNodeNNPtr &node,
+                                 bool removeInverseOf = false);
 
     ObjectDomainPtr buildObjectDomain(const WKTNodeNNPtr &node);
 
@@ -1396,11 +1398,16 @@ double WKTParser::Private::asDouble(const WKTNodeNNPtr &node) {
 // ---------------------------------------------------------------------------
 
 IdentifierPtr WKTParser::Private::buildId(const WKTNodeNNPtr &node,
-                                          bool tolerant) {
+                                          bool tolerant, bool removeInverseOf) {
     const auto *nodeP = node->GP();
     const auto &nodeChidren = nodeP->children();
     if (nodeChidren.size() >= 2) {
         auto codeSpace = stripQuotes(nodeChidren[0]);
+        if (removeInverseOf && starts_with(codeSpace, "INVERSE(") &&
+            codeSpace.back() == ')') {
+            codeSpace = codeSpace.substr(strlen("INVERSE("));
+            codeSpace.resize(codeSpace.size() - 1);
+        }
         auto code = stripQuotes(nodeChidren[1]);
         auto &citationNode = nodeP->lookForChild(WKTConstants::CITATION);
         auto &uriNode = nodeP->lookForChild(WKTConstants::URI);
@@ -1444,7 +1451,8 @@ IdentifierPtr WKTParser::Private::buildId(const WKTNodeNNPtr &node,
 
 // ---------------------------------------------------------------------------
 
-PropertyMap &WKTParser::Private::buildProperties(const WKTNodeNNPtr &node) {
+PropertyMap &WKTParser::Private::buildProperties(const WKTNodeNNPtr &node,
+                                                 bool removeInverseOf) {
 
     if (propertyCount_ == MAX_PROPERTY_SIZE) {
         throw ParsingException("MAX_PROPERTY_SIZE reached");
@@ -1460,6 +1468,10 @@ PropertyMap &WKTParser::Private::buildProperties(const WKTNodeNNPtr &node) {
     if (!nodeChildren.empty()) {
         const auto &nodeName(nodeP->value());
         auto name(stripQuotes(nodeChildren[0]));
+        if (removeInverseOf && starts_with(name, "Inverse of ")) {
+            name = name.substr(strlen("Inverse of "));
+        }
+
         if (ends_with(name, " (deprecated)")) {
             name.resize(name.size() - strlen(" (deprecated)"));
             properties->set(common::IdentifiedObject::DEPRECATED_KEY, true);
@@ -1512,7 +1524,7 @@ PropertyMap &WKTParser::Private::buildProperties(const WKTNodeNNPtr &node) {
         const auto &subNodeName(subNode->GP()->value());
         if (ci_equal(subNodeName, WKTConstants::ID) ||
             ci_equal(subNodeName, WKTConstants::AUTHORITY)) {
-            auto id = buildId(subNode);
+            auto id = buildId(subNode, true, removeInverseOf);
             if (id) {
                 identifiers->add(NN_NO_CHECK(id));
             }
@@ -1990,7 +2002,7 @@ GeodeticReferenceFrameNNPtr WKTParser::Private::buildGeodeticReferenceFrame(
                 auto &idNode = nodeP->lookForChild(WKTConstants::AUTHORITY);
                 if (!isNull(idNode)) {
                     try {
-                        auto id = buildId(idNode);
+                        auto id = buildId(idNode, true, false);
                         auto authFactory2 = AuthorityFactory::create(
                             NN_NO_CHECK(dbContext_), *id->codeSpace());
                         auto dbDatum =
@@ -2850,8 +2862,22 @@ WKTParser::Private::buildConversion(const WKTNodeNNPtr &node,
     consumeParameters(node, false, parameters, values, defaultLinearUnit,
                       defaultAngularUnit);
 
-    return Conversion::create(buildProperties(node),
-                              buildProperties(methodNode), parameters, values);
+    auto &convProps = buildProperties(node);
+    auto &methodProps = buildProperties(methodNode);
+    std::string convName;
+    std::string methodName;
+    if (convProps.getStringValue(IdentifiedObject::NAME_KEY, convName) &&
+        methodProps.getStringValue(IdentifiedObject::NAME_KEY, methodName) &&
+        starts_with(convName, "Inverse of ") &&
+        starts_with(methodName, "Inverse of ")) {
+
+        auto &invConvProps = buildProperties(node, true);
+        auto &invMethodProps = buildProperties(methodNode, true);
+        return NN_NO_CHECK(util::nn_dynamic_pointer_cast<Conversion>(
+            Conversion::create(invConvProps, invMethodProps, parameters, values)
+                ->inverse()));
+    }
+    return Conversion::create(convProps, methodProps, parameters, values);
 }
 
 // ---------------------------------------------------------------------------
@@ -2918,8 +2944,28 @@ WKTParser::Private::buildCoordinateOperation(const WKTNodeNNPtr &node) {
 
 ConcatenatedOperationNNPtr
 WKTParser::Private::buildConcatenatedOperation(const WKTNodeNNPtr &node) {
+
+    const auto *nodeP = node->GP();
+    auto &sourceCRSNode = nodeP->lookForChild(WKTConstants::SOURCECRS);
+    if (/*isNull(sourceCRSNode) ||*/ sourceCRSNode->GP()->childrenSize() != 1) {
+        ThrowMissing(WKTConstants::SOURCECRS);
+    }
+    auto sourceCRS = buildCRS(sourceCRSNode->GP()->children()[0]);
+    if (!sourceCRS) {
+        throw ParsingException("Invalid content in SOURCECRS node");
+    }
+
+    auto &targetCRSNode = nodeP->lookForChild(WKTConstants::TARGETCRS);
+    if (/*isNull(targetCRSNode) ||*/ targetCRSNode->GP()->childrenSize() != 1) {
+        ThrowMissing(WKTConstants::TARGETCRS);
+    }
+    auto targetCRS = buildCRS(targetCRSNode->GP()->children()[0]);
+    if (!targetCRS) {
+        throw ParsingException("Invalid content in TARGETCRS node");
+    }
+
     std::vector<CoordinateOperationNNPtr> operations;
-    for (const auto &childNode : node->GP()->children()) {
+    for (const auto &childNode : nodeP->children()) {
         if (ci_equal(childNode->GP()->value(), WKTConstants::STEP)) {
             if (childNode->GP()->childrenSize() != 1) {
                 throw ParsingException("Invalid content in STEP node");
@@ -2932,6 +2978,10 @@ WKTParser::Private::buildConcatenatedOperation(const WKTNodeNNPtr &node) {
             operations.emplace_back(NN_NO_CHECK(op));
         }
     }
+
+    ConcatenatedOperation::fixStepsDirection(
+        NN_NO_CHECK(sourceCRS), NN_NO_CHECK(targetCRS), operations);
+
     try {
         return ConcatenatedOperation::create(
             buildProperties(node), operations,
@@ -4213,7 +4263,7 @@ BaseObjectNNPtr WKTParser::Private::build(const WKTNodeNNPtr &node) {
     if (ci_equal(name, WKTConstants::ID) ||
         ci_equal(name, WKTConstants::AUTHORITY)) {
         return util::nn_static_pointer_cast<BaseObject>(
-            NN_NO_CHECK(buildId(node, false)));
+            NN_NO_CHECK(buildId(node, false, false)));
     }
 
     throw ParsingException(concat("unhandled keyword: ", name));
@@ -5914,7 +5964,8 @@ PROJStringParser::Private::buildDatum(const Step &step,
     PropertyMap grfMap;
 
     // It is arguable that we allow the prime meridian of a datum defined by
-    // its name to be overridden, but this is found at least in a regression test
+    // its name to be overridden, but this is found at least in a regression
+    // test
     // of GDAL. So let's keep the ellipsoid part of the datum in that case and
     // use the specified prime meridian.
     const auto overridePmIfNeeded =