CoordinateOperationFactory::Private::createOperationsWithDatumPivot()…

…: use candidateSrcGeod and candidateDstGeod of same type in priority

src/iso19111/operation/coordinateoperationfactory.cpp

@@ -109,6 +109,17 @@ namespace operation {
 static std::string objectAsStr(const common::IdentifiedObject *obj) {
     std::string ret(obj->nameStr());
     const auto &ids = obj->identifiers();
+    if (const auto *geogCRS = dynamic_cast<const crs::GeographicCRS *>(obj)) {
+        if (geogCRS->coordinateSystem()->axisList().size() == 3U)
+            ret += " (geographic3D)";
+        else
+            ret += " (geographic2D)";
+    }
+    if (const auto *geodCRS = dynamic_cast<const crs::GeodeticCRS *>(obj)) {
+        if (geodCRS->isGeocentric()) {
+            ret += " (geocentric)";
+        }
+    }
     if (!ids.empty()) {
         ret += " (";
         ret += (*ids[0]->codeSpace()) + ":" + ids[0]->code();
@@ -3100,12 +3111,16 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
     // If RGF93GEO is returned before then we go through WGS84 and use
     // instead a Helmert transformation.
     //
-    // Actually, in the general case, we do the lookup in 2 passes with the 2
+    // Actually, in the general case, we do the lookup in 3 passes with the 2
     // above steps in each pass:
     // - one first pass where we only consider direct transformations (no
     //   other intermediate CRS)
     // - a second pass where we allow transformation through another
-    //   intermediate CRS.
+    //   intermediate CRS, but we make sure the candidate geodetic CRS are of
+    //   the same type
+    // - a third where we allow transformation through another
+    //   intermediate CRS, where the candidate geodetic CRS are of different
+    //   type.
     // ... but when transforming between 2 IGNF CRS, we do just one single pass
     // by allowing directly all transformation. There is no strong reason for
     // that particular case, except that otherwise we'd get different results
@@ -3121,22 +3136,45 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
         const auto &ids = crs->identifiers();
         return !ids.empty() && *(ids.front()->codeSpace()) == "IGNF";
     };
-    const int nIters = (isIGNF(sourceCRS) && isIGNF(targetCRS)) ? 1 : 2;
+    const int nIters = (isIGNF(sourceCRS) && isIGNF(targetCRS)) ? 1 : 3;
+
+    const auto getType = [](const crs::GeodeticCRSNNPtr &crs) {
+        if (auto geogCRS =
+                dynamic_cast<const crs::GeographicCRS *>(crs.get())) {
+            if (geogCRS->coordinateSystem()->axisList().size() == 3)
+                return 1;
+            return 0;
+        }
+        return 2;
+    };
+
     for (int iter = 0; iter < nIters; ++iter) {
-        const bool useOnlyDirectRegistryOp = (iter == 0 && nIters == 2);
+        const bool useOnlyDirectRegistryOp = (iter == 0 && nIters == 3);
         for (const auto &candidateSrcGeod : candidatesSrcGeod) {
             if (candidateSrcGeod->nameStr() == sourceCRS->nameStr()) {
+                const auto typeSource =
+                    (iter >= 1) ? getType(candidateSrcGeod) : -1;
                 auto sourceSrcGeodModified(sourceAndTargetAre3D
                                                ? candidateSrcGeod->promoteTo3D(
                                                      std::string(), dbContext)
                                                : candidateSrcGeod);
                 for (const auto &candidateDstGeod : candidatesDstGeod) {
                     if (candidateDstGeod->nameStr() == targetCRS->nameStr()) {
+                        if (iter == 1) {
+                            if (typeSource != getType(candidateDstGeod)) {
+                                continue;
+                            }
+                        } else if (iter == 2) {
+                            if (typeSource == getType(candidateDstGeod)) {
+                                continue;
+                            }
+                        }
 #ifdef TRACE_CREATE_OPERATIONS
-                        ENTER_BLOCK("try " + objectAsStr(sourceCRS.get()) +
-                                    "->" + objectAsStr(candidateSrcGeod.get()) +
-                                    "->" + objectAsStr(candidateDstGeod.get()) +
-                                    "->" + objectAsStr(targetCRS.get()) + ")");
+                        ENTER_BLOCK("iter=" + toString(iter) + ", try " +
+                                    objectAsStr(sourceCRS.get()) + "->" +
+                                    objectAsStr(candidateSrcGeod.get()) + "->" +
+                                    objectAsStr(candidateDstGeod.get()) + "->" +
+                                    objectAsStr(targetCRS.get()) + ")");
 #endif
                         const auto opsFirst = createOperations(
                             sourceCRS, sourceSrcGeodModified, context);

test/unit/test_operationfactory.cpp

@@ -9633,3 +9633,33 @@ TEST(operation, createOperation_Geographic3D_Offset_by_velocity_grid) {
               "+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
               "+step +proj=axisswap +order=2,1");
 }
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, createOperation_test_createOperationsWithDatumPivot_iter_1) {
+    // Test
+    // CoordinateOperationFactory::Private::createOperationsWithDatumPivot()
+    // iter=1, ie getType(candidateSrcGeod) == getType(candidateDstGeod)
+    auto dbContext = DatabaseContext::create();
+    auto factoryEPSG = AuthorityFactory::create(dbContext, "EPSG");
+    // NAD83(CSRS)v2 (2D)
+    auto sourceCRS = factoryEPSG->createCoordinateReferenceSystem("8237");
+    // NAD83(CSRS)v3 (2D)
+    auto targetCRS = factoryEPSG->createCoordinateReferenceSystem("8240");
+    auto ctxt = CoordinateOperationContext::create(
+        AuthorityFactory::create(dbContext, std::string()), nullptr, 0);
+    // Do *NOT* set SpatialCriterion::PARTIAL_INTERSECTION, otherwise we'd
+    // get the direct operations
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        sourceCRS, targetCRS, ctxt);
+    ASSERT_GE(list.size(), 1U);
+    EXPECT_FALSE(list[0]->hasBallparkTransformation());
+    EXPECT_STREQL(list[0]->nameStr().c_str(),
+                  "Inverse of NAD83(CSRS)v8 to NAD83(CSRS)v2 (1) + "
+                  "NAD83(CSRS)v8 to NAD83(CSRS)v3 (1)");
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=noop");
+}