From 0100ad08e12a1464865f171a353dc49be43ae766 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 4 Sep 2024 20:57:04 +0200 Subject: [PATCH 1/2] AuthorityFactory::createBetweenGeodeticCRSWithDatumBasedIntermediates(): avoid using intermediate CRS of ancient era when dealing with modern source/target CRS --- src/iso19111/factory.cpp | 57 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 52 insertions(+), 5 deletions(-) diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index a0332d3ae9..bb44c97f26 100644 --- a/src/iso19111/factory.cpp +++ b/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(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 steps; if (!sourceCRS->isEquivalentTo( From 529904c26313978632026a6c25e25896bbc66ab3 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 4 Sep 2024 20:57:40 +0200 Subject: [PATCH 2/2] createOperation(): tune so that ITRF2000->ETRS89 does not return only NKG grid based operations but also time-dependent Helmert --- .../operation/coordinateoperationfactory.cpp | 89 ++++++++++++++++++- test/cli/test_projinfo.yaml | 3 +- test/unit/test_operationfactory.cpp | 43 ++++++++- 3 files changed, 128 insertions(+), 7 deletions(-) diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index c6a3dbcbb8..5b98fd3bc2 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/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 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(); } diff --git a/test/cli/test_projinfo.yaml b/test/cli/test_projinfo.yaml index 58f1f1c1dd..af540c95bb 100644 --- a/test/cli/test_projinfo.yaml +++ b/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 diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 83beebaa3c..6bdf952de0 100644 --- a/test/unit/test_operationfactory.cpp +++ b/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); +}