Merge pull request #4244 from rouault/ETRF2000_ETRS89

createOperation(): tune so that ITRF2000->ETRS89 does not return only NKG grid based operations but also time-dependent Helmert

src/iso19111/factory.cpp

@@ -636,14 +636,11 @@ SQLiteHandleCache::getHandle(const std::string &path, PJ_CONTEXT *ctx) {
                 // to acquire the mutex in invalidateHandles().
                 SQLiteHandleCache::get().sMutex_.lock();
             },
-            []() {
-                SQLiteHandleCache::get().sMutex_.unlock();
-            },
+            []() { SQLiteHandleCache::get().sMutex_.unlock(); },
             []() {
                 SQLiteHandleCache::get().sMutex_.unlock();
                 SQLiteHandleCache::get().invalidateHandles();
-            }
-        );
+            });
     }
 #endif
 
@@ -8142,6 +8139,25 @@ AuthorityFactory::createBetweenGeodeticCRSWithDatumBasedIntermediates(
         }
     }
 
+    std::string sourceDatumPubDate;
+    const auto sourceDatum = sourceGeodCRS->datumNonNull(d->context());
+    if (sourceDatum->publicationDate().has_value()) {
+        sourceDatumPubDate = sourceDatum->publicationDate()->toString();
+    }
+
+    std::string targetDatumPubDate;
+    const auto targetDatum = targetGeodCRS->datumNonNull(d->context());
+    if (targetDatum->publicationDate().has_value()) {
+        targetDatumPubDate = targetDatum->publicationDate()->toString();
+    }
+
+    const std::string mostAncientDatumPubDate =
+        (!targetDatumPubDate.empty() &&
+         (sourceDatumPubDate.empty() ||
+          targetDatumPubDate < sourceDatumPubDate))
+            ? targetDatumPubDate
+            : sourceDatumPubDate;
+
     auto opFactory = operation::CoordinateOperationFactory::create();
     for (const auto &pair : candidates) {
         const auto &trfmSource = pair.first;
@@ -8169,6 +8185,37 @@ AuthorityFactory::createBetweenGeodeticCRSWithDatumBasedIntermediates(
             continue;
         }
 
+        // Skip operations using a datum that is older than the source or
+        // target datum (e.g to avoid ED50 to WGS84 to go through NTF)
+        if (!mostAncientDatumPubDate.empty()) {
+            const auto isOlderCRS = [this, &mostAncientDatumPubDate](
+                                        const crs::CRSPtr &crs) {
+                const auto geogCRS =
+                    dynamic_cast<const crs::GeodeticCRS *>(crs.get());
+                if (geogCRS) {
+                    const auto datum = geogCRS->datumNonNull(d->context());
+                    // Hum, theoretically we'd want to check
+                    // datum->publicationDate()->toString() <
+                    // mostAncientDatumPubDate but that would exclude doing
+                    // IG05/12 Intermediate CRS to ITRF2014 through ITRF2008,
+                    // since IG05/12 Intermediate CRS has been published later
+                    // than ITRF2008. So use a cut of date for ancient vs
+                    // "modern" era.
+                    constexpr const char *CUT_OFF_DATE = "1900-01-01";
+                    if (datum->publicationDate().has_value() &&
+                        datum->publicationDate()->toString() < CUT_OFF_DATE &&
+                        mostAncientDatumPubDate > CUT_OFF_DATE) {
+                        return true;
+                    }
+                }
+                return false;
+            };
+
+            if (isOlderCRS(op1Source) || isOlderCRS(op1Target) ||
+                isOlderCRS(op2Source) || isOlderCRS(op2Target))
+                continue;
+        }
+
         std::vector<operation::CoordinateOperationNNPtr> steps;
 
         if (!sourceCRS->isEquivalentTo(

src/iso19111/operation/coordinateoperationfactory.cpp

@@ -3940,16 +3940,99 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
         doFilterAndCheckPerfectOp = !res.empty();
     }
 
-    if (res.empty() && !context.inCreateOperationsWithDatumPivotAntiRecursion &&
+    // 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 CRSs.
+    // But only do that if we have a relatively small number of solutions.
+    // Otherwise we might just spend time appending more dubious solutions.
+    bool tooSmallAreas = false;
+    // 10: arbitrary threshold so that particular case triggers for
+    // EPSG:9989 (ITRF2000) to EPSG:4937 (ETRS89), but not for
+    // "WGS 84" to "NAD83(CSRS)v2", to avoid excessive computation time.
+    if (!res.empty() && res.size() < 10) {
+        const auto &areaOfInterest = context.context->getAreaOfInterest();
+        double targetArea = 0.0;
+        if (areaOfInterest) {
+            targetArea = getPseudoArea(NN_NO_CHECK(areaOfInterest));
+        } else if (context.extent1) {
+            if (context.extent2) {
+                auto intersection =
+                    context.extent1->intersection(NN_NO_CHECK(context.extent2));
+                if (intersection)
+                    targetArea = getPseudoArea(NN_NO_CHECK(intersection));
+            } else {
+                targetArea = getPseudoArea(NN_NO_CHECK(context.extent1));
+            }
+        } else if (context.extent2) {
+            targetArea = getPseudoArea(NN_NO_CHECK(context.extent2));
+        }
+        if (targetArea > 0) {
+            tooSmallAreas = true;
+            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);
+                    }
+                }
+
+                constexpr double HALF_RATIO = 0.5;
+                if (area > HALF_RATIO * targetArea) {
+                    tooSmallAreas = false;
+                    break;
+                }
+            }
+        }
+    }
+
+    if ((res.empty() || tooSmallAreas) &&
+        !context.inCreateOperationsWithDatumPivotAntiRecursion &&
         !resFindDirectNonEmptyBeforeFiltering && geodSrc && geodDst &&
         !sameGeodeticDatum && context.context->getIntermediateCRS().empty() &&
         context.context->getAllowUseIntermediateCRS() !=
             CoordinateOperationContext::IntermediateCRSUse::NEVER) {
         // Currently triggered by "IG05/12 Intermediate CRS" to ITRF2014
+
+        std::set<std::string> oSetNames;
+        if (tooSmallAreas) {
+            for (const auto &op : res) {
+                oSetNames.insert(op->nameStr());
+            }
+        }
         auto resWithIntermediate = findsOpsInRegistryWithIntermediate(
             sourceCRS, targetCRS, context, true);
-        res.insert(res.end(), resWithIntermediate.begin(),
-                   resWithIntermediate.end());
+        if (tooSmallAreas && !res.empty()) {
+            // Only insert operations we didn't already find
+            for (const auto &op : resWithIntermediate) {
+                if (oSetNames.find(op->nameStr()) == oSetNames.end()) {
+                    res.emplace_back(op);
+                }
+            }
+        } else {
+            res.insert(res.end(), resWithIntermediate.begin(),
+                       resWithIntermediate.end());
+        }
         doFilterAndCheckPerfectOp = !res.empty();
     }
 

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_ITRF2000_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(
+        // ITRF2000
+        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);
+}