Add logic to handle interpolationCRS in PROJBasedOperation and ConcatenatedOperation

3 files changed, 81 insertions, 33 deletions

diff --git a/include/proj/internal/coordinateoperation_internal.hpp b/include/proj/internal/coordinateoperation_internal.hpp
index 69f30933..322d71b7 100644
--- a/include/proj/internal/coordinateoperation_internal.hpp
+++ b/include/proj/internal/coordinateoperation_internal.hpp
@@ -275,6 +275,7 @@ class PROJBasedOperation : public SingleOperation {
     create(const util::PropertyMap &properties,
            const io::IPROJStringExportableNNPtr &projExportable, bool inverse,
            const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+           const crs::CRSPtr &interpolationCRS,
            const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies,
            bool hasRoughTransformation);
 
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 3577b6ff..82ec4f30 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -1753,11 +1753,12 @@ static SingleOperationNNPtr createPROJBased(
     const util::PropertyMap &properties,
     const io::IPROJStringExportableNNPtr &projExportable,
     const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    const crs::CRSPtr &interpolationCRS,
     const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies,
     bool hasBallparkTransformation) {
     return util::nn_static_pointer_cast<SingleOperation>(
         PROJBasedOperation::create(properties, projExportable, false, sourceCRS,
-                                   targetCRS, accuracies,
+                                   targetCRS, interpolationCRS, accuracies,
                                    hasBallparkTransformation));
 }
 //! @endcond
@@ -10209,9 +10210,31 @@ ConcatenatedOperationNNPtr ConcatenatedOperation::create(
             "ConcatenatedOperation must have at least 2 operations");
     }
     crs::CRSPtr lastTargetCRS;
+
+    crs::CRSPtr interpolationCRS;
+    bool interpolationCRSValid = true;
     for (size_t i = 0; i < operationsIn.size(); i++) {
         auto l_sourceCRS = operationsIn[i]->sourceCRS();
         auto l_targetCRS = operationsIn[i]->targetCRS();
+
+        if (interpolationCRSValid) {
+            auto subOpInterpCRS = operationsIn[i]->interpolationCRS();
+            if (interpolationCRS == nullptr)
+                interpolationCRS = subOpInterpCRS;
+            else if ((subOpInterpCRS == nullptr &&
+                      interpolationCRS != nullptr) ||
+                     (subOpInterpCRS != nullptr &&
+                      interpolationCRS == nullptr) ||
+                     (subOpInterpCRS != nullptr &&
+                      interpolationCRS != nullptr &&
+                      !(subOpInterpCRS->isEquivalentTo(
+                          interpolationCRS.get(),
+                          util::IComparable::Criterion::EQUIVALENT)))) {
+                interpolationCRS = nullptr;
+                interpolationCRSValid = false;
+            }
+        }
+
         if (l_sourceCRS == nullptr || l_targetCRS == nullptr) {
             throw InvalidOperation("At least one of the operation lacks a "
                                    "source and/or target CRS");
@@ -10249,7 +10272,8 @@ ConcatenatedOperationNNPtr ConcatenatedOperation::create(
     op->assignSelf(op);
     op->setProperties(properties);
     op->setCRSs(NN_NO_CHECK(operationsIn[0]->sourceCRS()),
-                NN_NO_CHECK(operationsIn.back()->targetCRS()), nullptr);
+                NN_NO_CHECK(operationsIn.back()->targetCRS()),
+                interpolationCRS);
     op->setAccuracies(accuracies);
 #ifdef DEBUG_CONCATENATED_OPERATION
     {
@@ -10361,6 +10385,15 @@ void ConcatenatedOperation::fixStepsDirection(
             } else if (i >= 1) {
                 l_sourceCRS = operationsInOut[i - 1]->targetCRS();
                 if (l_sourceCRS) {
+                    derivedCRS = dynamic_cast<const crs::DerivedCRS *>(
+                        l_sourceCRS.get());
+                    if (conv->isEquivalentTo(
+                            derivedCRS->derivingConversion().get(),
+                            util::IComparable::Criterion::EQUIVALENT)) {
+                        op->setCRSs(concatOpTargetCRS, NN_NO_CHECK(l_sourceCRS),
+                                    nullptr);
+                        op = op->inverse();
+                    }
                     op->setCRSs(NN_NO_CHECK(l_sourceCRS), concatOpTargetCRS,
                                 nullptr);
                 }
@@ -12581,7 +12614,8 @@ createGeodToGeodPROJBased(const crs::CRSNNPtr &geodSrc,
     auto properties = util::PropertyMap().set(
         common::IdentifiedObject::NAME_KEY,
         buildTransfName(geodSrc->nameStr(), geodDst->nameStr()));
-    return createPROJBased(properties, exportable, geodSrc, geodDst, {}, false);
+    return createPROJBased(properties, exportable, geodSrc, geodDst, nullptr,
+                           {}, false);
 }
 
 // ---------------------------------------------------------------------------
@@ -12612,7 +12646,7 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased(
                            NN_NO_CHECK(extent));
         }
         return createPROJBased(
-            properties, exportable, sourceCRS, targetCRS,
+            properties, exportable, sourceCRS, targetCRS, nullptr,
             verticalTransform->coordinateOperationAccuracies(),
             verticalTransform->hasBallparkTransformation());
     } else {
@@ -12643,7 +12677,7 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased(
         }
 
         return createPROJBased(
-            properties, exportable, sourceCRS, targetCRS, accuracies,
+            properties, exportable, sourceCRS, targetCRS, nullptr, accuracies,
             horizTransform->hasBallparkTransformation() ||
                 verticalTransform->hasBallparkTransformation());
     }
@@ -12703,7 +12737,7 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased(
     }
 
     return createPROJBased(properties, exportable, sourceCRS, targetCRS,
-                           accuracies, hasBallparkTransformation);
+                           nullptr, accuracies, hasBallparkTransformation);
 }
 
 //! @endcond
@@ -15303,19 +15337,14 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
 
     for (const auto &verticalTransform : verticalTransforms) {
         auto interpolationGeogCRS = NN_NO_CHECK(srcGeog);
-        auto transformationVerticalTransform =
-            dynamic_cast<const Transformation *>(verticalTransform.get());
-        if (transformationVerticalTransform) {
-            auto interpTransformCRS =
-                transformationVerticalTransform->interpolationCRS();
-            if (interpTransformCRS) {
-                auto nn_interpTransformCRS = NN_NO_CHECK(interpTransformCRS);
-                if (dynamic_cast<const crs::GeographicCRS *>(
-                        nn_interpTransformCRS.get())) {
-                    interpolationGeogCRS = NN_NO_CHECK(
-                        util::nn_dynamic_pointer_cast<crs::GeographicCRS>(
-                            nn_interpTransformCRS));
-                }
+        auto interpTransformCRS = verticalTransform->interpolationCRS();
+        if (interpTransformCRS) {
+            auto nn_interpTransformCRS = NN_NO_CHECK(interpTransformCRS);
+            if (dynamic_cast<const crs::GeographicCRS *>(
+                    nn_interpTransformCRS.get())) {
+                interpolationGeogCRS = NN_NO_CHECK(
+                    util::nn_dynamic_pointer_cast<crs::GeographicCRS>(
+                        nn_interpTransformCRS));
             }
         } else {
             auto compSrc0BoundCrs =
@@ -15756,6 +15785,7 @@ PROJBasedOperationNNPtr PROJBasedOperation::create(
     const util::PropertyMap &properties,
     const io::IPROJStringExportableNNPtr &projExportable, bool inverse,
     const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+    const crs::CRSPtr &interpolationCRS,
     const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies,
     bool hasBallparkTransformation) {
 
@@ -15777,7 +15807,7 @@ PROJBasedOperationNNPtr PROJBasedOperation::create(
     auto op = PROJBasedOperation::nn_make_shared<PROJBasedOperation>(method);
     op->assignSelf(op);
     op->projString_ = projString;
-    op->setCRSs(sourceCRS, targetCRS, nullptr);
+    op->setCRSs(sourceCRS, targetCRS, interpolationCRS);
     op->setProperties(
         addDefaultNameIfNeeded(properties, "PROJ-based coordinate operation"));
     op->setAccuracies(accuracies);
@@ -15791,13 +15821,14 @@ PROJBasedOperationNNPtr PROJBasedOperation::create(
 
 CoordinateOperationNNPtr PROJBasedOperation::inverse() const {
 
-    if (projStringExportable_) {
+    if (projStringExportable_ && sourceCRS() && targetCRS()) {
         return util::nn_static_pointer_cast<CoordinateOperation>(
             PROJBasedOperation::create(
                 createPropertiesForInverse(this, false, false),
                 NN_NO_CHECK(projStringExportable_), !inverse_,
                 NN_NO_CHECK(targetCRS()), NN_NO_CHECK(sourceCRS()),
-                coordinateOperationAccuracies(), hasBallparkTransformation()));
+                interpolationCRS(), coordinateOperationAccuracies(),
+                hasBallparkTransformation()));
     }
 
     auto formatter = io::PROJStringFormatter::create();
@@ -15813,6 +15844,10 @@ CoordinateOperationNNPtr PROJBasedOperation::inverse() const {
     auto op = PROJBasedOperation::create(
         createPropertiesForInverse(this, false, false), formatter->toString(),
         targetCRS(), sourceCRS(), coordinateOperationAccuracies());
+    if (sourceCRS() && targetCRS()) {
+        op->setCRSs(NN_NO_CHECK(targetCRS()), NN_NO_CHECK(sourceCRS()),
+                    interpolationCRS());
+    }
     op->setHasBallparkTransformation(hasBallparkTransformation());
     return util::nn_static_pointer_cast<CoordinateOperation>(op);
 }
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index 8d9742d5..08b627b2 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -3348,7 +3348,8 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation(
                "method_auth_name, method_code, method_name, "
                "source_crs_auth_name, source_crs_code, target_crs_auth_name, "
                "target_crs_code, "
-               "accuracy";
+               "interpolation_crs_auth_name, interpolation_crs_code, "
+               "operation_version, accuracy, deprecated";
         constexpr int N_MAX_PARAMS = 7;
         for (int i = 1; i <= N_MAX_PARAMS; ++i) {
             buffer << ", param" << i << "_auth_name";
@@ -3358,7 +3359,7 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation(
             buffer << ", param" << i << "_uom_auth_name";
             buffer << ", param" << i << "_uom_code";
         }
-        buffer << ", operation_version, deprecated FROM other_transformation "
+        buffer << " FROM other_transformation "
                   "WHERE auth_name = ? AND code = ?";
 
         auto res = d->runWithCodeParam(buffer.str(), code);
@@ -3379,7 +3380,12 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation(
             const auto &source_crs_code = row[idx++];
             const auto &target_crs_auth_name = row[idx++];
             const auto &target_crs_code = row[idx++];
+            const auto &interpolation_crs_auth_name = row[idx++];
+            const auto &interpolation_crs_code = row[idx++];
+            const auto &operation_version = row[idx++];
             const auto &accuracy = row[idx++];
+            const auto &deprecated_str = row[idx++];
+            const bool deprecated = deprecated_str == "1";
 
             const size_t base_param_idx = idx;
             std::vector<operation::OperationParameterNNPtr> parameters;
@@ -3410,10 +3416,7 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation(
                     common::Measure(normalized_value, uom)));
             }
             idx = base_param_idx + 6 * N_MAX_PARAMS;
-
-            const auto &operation_version = row[idx++];
-            const auto &deprecated_str = row[idx++];
-            const bool deprecated = deprecated_str == "1";
+            (void)idx;
             assert(idx == row.size());
 
             auto sourceCRS =
@@ -3422,6 +3425,13 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation(
             auto targetCRS =
                 d->createFactory(target_crs_auth_name)
                     ->createCoordinateReferenceSystem(target_crs_code);
+            auto interpolationCRS =
+                interpolation_crs_auth_name.empty()
+                    ? nullptr
+                    : d->createFactory(interpolation_crs_auth_name)
+                          ->createCoordinateReferenceSystem(
+                              interpolation_crs_code)
+                          .as_nullable();
 
             auto props = d->createPropertiesSearchUsages(
                 type, code, name, deprecated, description);
@@ -3438,8 +3448,10 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation(
 
             if (method_auth_name == "PROJ") {
                 if (method_code == "PROJString") {
-                    return operation::SingleOperation::createPROJBased(
+                    auto op = operation::SingleOperation::createPROJBased(
                         props, method_name, sourceCRS, targetCRS, accuracies);
+                    op->setCRSs(sourceCRS, targetCRS, interpolationCRS);
+                    return op;
                 } else if (method_code == "WKT") {
                     auto op = util::nn_dynamic_pointer_cast<
                         operation::CoordinateOperation>(
@@ -3448,7 +3460,7 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation(
                         throw FactoryException("WKT string does not express a "
                                                "coordinate operation");
                     }
-                    op->setCRSs(sourceCRS, targetCRS, nullptr);
+                    op->setCRSs(sourceCRS, targetCRS, interpolationCRS);
                     return NN_NO_CHECK(op);
                 }
             }
@@ -3466,13 +3478,13 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation(
                     method_code_int == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL) {
                     auto op = operation::Conversion::create(props, propsMethod,
                                                             parameters, values);
-                    op->setCRSs(sourceCRS, targetCRS, nullptr);
+                    op->setCRSs(sourceCRS, targetCRS, interpolationCRS);
                     return op;
                 }
             }
             return operation::Transformation::create(
-                props, sourceCRS, targetCRS, nullptr, propsMethod, parameters,
-                values, accuracies);
+                props, sourceCRS, targetCRS, interpolationCRS, propsMethod,
+                parameters, values, accuracies);
 
         } catch (const std::exception &ex) {
             throw buildFactoryException("transformation", code, ex);