createOperation(): tune so that ETRF2000->ETRS89 does not return only…

… NKG grid based operations but also time-dependent Helmert

src/iso19111/operation/coordinateoperationfactory.cpp

@@ -3922,6 +3922,66 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
         foundInstantiableOp = true;
     }
 
+    // Browse through candidate operations and check if their area of use
+    // is sufficiently large compared to the area of use of the
+    // source/target CRS. If it is not, we might need to extend the lookup
+    // to using intermediate pivots.
+    const auto &areaOfInterest = context.context->getAreaOfInterest();
+    double largestArea = 0;
+    for (const auto &op : res) {
+        bool dummy = false;
+        auto extentOp = getExtent(op, true, dummy);
+        double area = 0.0;
+        if (extentOp) {
+            if (areaOfInterest) {
+                area = getPseudoArea(
+                    extentOp->intersection(NN_NO_CHECK(areaOfInterest)));
+            } else if (context.extent1 && context.extent2) {
+                auto x = extentOp->intersection(NN_NO_CHECK(context.extent1));
+                auto y = extentOp->intersection(NN_NO_CHECK(context.extent2));
+                area =
+                    getPseudoArea(x) + getPseudoArea(y) -
+                    ((x && y) ? getPseudoArea(x->intersection(NN_NO_CHECK(y)))
+                              : 0.0);
+            } else if (context.extent1) {
+                area = getPseudoArea(
+                    extentOp->intersection(NN_NO_CHECK(context.extent1)));
+            } else if (context.extent2) {
+                area = getPseudoArea(
+                    extentOp->intersection(NN_NO_CHECK(context.extent2)));
+            } else {
+                area = getPseudoArea(extentOp);
+            }
+        }
+
+        largestArea = std::max(largestArea, area);
+    }
+
+    constexpr double HALF_RATIO = 0.5;
+    bool tooSmallAreas = false;
+    if (areaOfInterest) {
+        tooSmallAreas =
+            (largestArea <
+             HALF_RATIO * getPseudoArea(NN_NO_CHECK(areaOfInterest)));
+    } else if (context.extent1) {
+        if (context.extent2) {
+            auto intersection =
+                context.extent1->intersection(NN_NO_CHECK(context.extent2));
+            if (intersection)
+                tooSmallAreas =
+                    largestArea <
+                    HALF_RATIO * getPseudoArea(NN_NO_CHECK(intersection));
+        } else {
+            tooSmallAreas =
+                (largestArea <
+                 HALF_RATIO * getPseudoArea(NN_NO_CHECK(context.extent1)));
+        }
+    } else if (context.extent2) {
+        tooSmallAreas =
+            (largestArea <
+             HALF_RATIO * getPseudoArea(NN_NO_CHECK(context.extent2)));
+    }
+
     // NAD27 to NAD83 has tens of results already. No need to look
     // for a pivot
     if (!sameGeodeticDatum &&
@@ -3940,7 +4000,8 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
         doFilterAndCheckPerfectOp = !res.empty();
     }
 
-    if (res.empty() && !context.inCreateOperationsWithDatumPivotAntiRecursion &&
+    if ((res.empty() || tooSmallAreas) &&
+        !context.inCreateOperationsWithDatumPivotAntiRecursion &&
         !resFindDirectNonEmptyBeforeFiltering && geodSrc && geodDst &&
         !sameGeodeticDatum && context.context->getIntermediateCRS().empty() &&
         context.context->getAllowUseIntermediateCRS() !=

test/cli/test_projinfo.yaml

@@ -1376,7 +1376,8 @@ tests:
   comment: "Undocumented option: --normalize-axis-order"
   out: |
     +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=stere +lat_0=90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 +ellps=WGS84
-- args: -s EPSG:4936 -t EPSG:4978 --spatial-test intersects --summary
+- args: -s EPSG:4936 -t EPSG:4978 --spatial-test intersects --summary --bbox 5,54,6,55
+  comment: "WGS 84 to ETRS89 (2) uses a transformation method not supported by PROJ currently (time-specific Helmert), and thus must be sorted last"
   out: |
     Candidate operations found: 2
     unknown id, Ballpark geocentric translation from ETRS89 to WGS 84, unknown accuracy, World, has ballpark transformation

test/unit/test_operationfactory.cpp

@@ -1192,12 +1192,14 @@ TEST(operation, geogCRS_to_geogCRS_context_lonlat_vs_latlon_crs) {
         authFactory->createCoordinateReferenceSystem("10310"), // RGNC15
         ctxt);
     ASSERT_EQ(list.size(), 3U);
-    // Check that we get direct transformation, and not through WGS 84
+    EXPECT_EQ(list[0]->nameStr(),
+              "RGNC91-93 to WGS 84 (1) + Inverse of RGNC15 to WGS 84 (1)");
+    // Check that we get direct transformation, and not only through WGS 84
     // The difficulty here is that the transformation is registered between
     // "RGNC91-93 (lon-lat)" et "RGNC15 (lon-lat)"
-    EXPECT_EQ(list[0]->nameStr(), "axis order change (2D) + RGNC91-93 to "
-                                  "RGNC15 (2) + axis order change (2D)");
     EXPECT_EQ(list[1]->nameStr(), "axis order change (2D) + RGNC91-93 to "
+                                  "RGNC15 (2) + axis order change (2D)");
+    EXPECT_EQ(list[2]->nameStr(), "axis order change (2D) + RGNC91-93 to "
                                   "RGNC15 (1) + axis order change (2D)");
 }
 
@@ -10937,3 +10939,38 @@ TEST(operation, createOperation_epsg_10622_local_orthographic) {
               "+alpha=27.7928209333 +k=0.9999968 +x_0=0 +y_0=0 +ellps=GRS80 "
               "+step +proj=unitconvert +xy_in=m +xy_out=us-ft");
 }
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, createOperation_ETRF2000_to_ETRS89) {
+    auto dbContext = DatabaseContext::create();
+    auto authFactory = AuthorityFactory::create(dbContext, std::string());
+    auto authFactoryEPSG = AuthorityFactory::create(dbContext, "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        // ETRF2000
+        authFactoryEPSG->createCoordinateReferenceSystem("9989"),
+        // ETRS89
+        authFactoryEPSG->createCoordinateReferenceSystem("4937"), ctxt);
+
+    // Check that we find a mix of grid-based (NKG ones, with very narrow area
+    // of use, but easier to lookup) and non-grid based operations (wider area
+    // of use, but require more effort to infer)
+    bool foundNonGridBasedOp = false;
+    bool foundGridBaseOp = false;
+    for (const auto &op : list) {
+        if (!op->hasBallparkTransformation()) {
+            if (op->gridsNeeded(dbContext, false).empty())
+                foundNonGridBasedOp = true;
+            else
+                foundGridBaseOp = true;
+        }
+    }
+    EXPECT_TRUE(foundNonGridBasedOp);
+    EXPECT_TRUE(foundGridBaseOp);
+}