From 2f38bbbc0aab502b51041402682c81727ee8c088 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 9 Dec 2024 00:41:33 +0100 Subject: [PATCH] createOperations(): fine tune for 'GDA94 + AHD height' to 'GDA2020 to AVWS height' --- .../operation/coordinateoperationfactory.cpp | 40 +++++++++++-------- test/unit/test_operationfactory.cpp | 38 ++++++++++++++++++ 2 files changed, 61 insertions(+), 17 deletions(-) diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index 078a0196d3..34418018d2 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -6452,17 +6452,23 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( if (srcGeog == nullptr || dstGeog == nullptr) { return; } - - // Use PointMotionOperations if appropriate and available + const bool srcGeogIsSameAsDstGeog = srcGeog->_isEquivalentTo( + dstGeog.get(), util::IComparable::Criterion::EQUIVALENT); const auto &authFactory = context.context->getAuthorityFactory(); auto dbContext = authFactory ? authFactory->databaseContext().as_nullable() : nullptr; + const auto intermGeogSrc = srcGeog->promoteTo3D(std::string(), dbContext); + const auto intermGeogDst = + srcGeogIsSameAsDstGeog ? intermGeogSrc + : dstGeog->promoteTo3D(std::string(), dbContext); + const auto opsGeogSrcToGeogDst = createOperations( + intermGeogSrc, sourceEpoch, intermGeogDst, targetEpoch, context); + // Use PointMotionOperations if appropriate and available if (authFactory && sourceEpoch.has_value() && targetEpoch.has_value() && !sourceEpoch->coordinateEpoch()._isEquivalentTo( targetEpoch->coordinateEpoch()) && - srcGeog->_isEquivalentTo(dstGeog.get(), - util::IComparable::Criterion::EQUIVALENT)) { + srcGeogIsSameAsDstGeog) { const auto pmoSrc = authFactory->getPointMotionOperationsFor( NN_NO_CHECK(srcGeog), true); if (!pmoSrc.empty()) { @@ -6588,9 +6594,20 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( // If we didn't find a non-ballpark transformation between // the 2 vertical CRS, then try through intermediate geographic // CRS + // Do that although when the geographic CRS of the source and + // target CRS are not the same, but only if they have a 3D + // known version, and there is a non-ballpark transformation + // between them. + // This helps for "GDA94 + AHD height" to "GDA2020 + AVWS + // height" going through GDA94 3D and GDA2020 3D bTryThroughIntermediateGeogCRS = (verticalTransforms.size() == 1 && - verticalTransforms.front()->hasBallparkTransformation()); + verticalTransforms.front()->hasBallparkTransformation()) || + (!srcGeogIsSameAsDstGeog && + !intermGeogSrc->identifiers().empty() && + !intermGeogDst->identifiers().empty() && + !opsGeogSrcToGeogDst.empty() && + !opsGeogSrcToGeogDst.front()->hasBallparkTransformation()); } } } @@ -6648,14 +6665,6 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( // WGS 84 + EGM96 --> ETRS89 + Belfast height where // there is a geoid model for EGM96 referenced to WGS 84 // and a geoid model for Belfast height referenced to ETRS89 - const auto intermGeogSrc = - srcGeog->promoteTo3D(std::string(), dbContext); - const bool intermGeogSrcIsSameAsIntermGeogDst = - srcGeog->_isEquivalentTo(dstGeog.get()); - const auto intermGeogDst = - intermGeogSrcIsSameAsIntermGeogDst - ? intermGeogSrc - : dstGeog->promoteTo3D(std::string(), dbContext); const auto opsSrcToGeog = createOperations( sourceCRS, sourceEpoch, intermGeogSrc, sourceEpoch, context); const auto opsGeogToTarget = createOperations( @@ -6695,9 +6704,6 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( double bestAccuracy = -1; size_t bestStepCount = 0; if (hasNonTrivialSrcTransf && hasNonTrivialTargetTransf) { - const auto opsGeogSrcToGeogDst = - createOperations(intermGeogSrc, sourceEpoch, intermGeogDst, - targetEpoch, context); // In first pass, exclude (horizontal) ballpark operations, but // accept them in second pass. for (int pass = 0; pass <= 1 && res.empty(); pass++) { @@ -6726,7 +6732,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( try { res.emplace_back( ConcatenatedOperation::createComputeMetadata( - intermGeogSrcIsSameAsIntermGeogDst + srcGeogIsSameAsDstGeog ? std::vector< CoordinateOperationNNPtr>{op1, op3} diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index ad38690999..755c749017 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -5992,6 +5992,44 @@ TEST(operation, compoundCRS_to_compoundCRS_context_helmert) { // --------------------------------------------------------------------------- +TEST(operation, compoundCRS_to_compoundCRS_GDA94_AHD_to_GDA2020_AVWS) { + auto dbContext = DatabaseContext::create(); + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + // GDA94 + AHD height + authFactory->createCoordinateReferenceSystem("9464"), + // GDA2020 + AVWS height + authFactory->createCoordinateReferenceSystem("9462"), ctxt); + ASSERT_GE(list.size(), 1U); + EXPECT_EQ(list[0]->nameStr(), + "Inverse of 'GDA94 to GDA94 + AHD height (1)' + " + "GDA94 to GDA2020 (1) + " + "GDA2020 to AVWS height (2)"); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=vgridshift +grids=au_ga_AUSGeoid09_V1.01.tif " + "+multiplier=1 " + "+step +proj=cart +ellps=GRS80 " + "+step +proj=helmert +x=0.06155 +y=-0.01087 +z=-0.04019 " + "+rx=-0.0394924 +ry=-0.0327221 +rz=-0.0328979 +s=-0.009994 " + "+convention=coordinate_frame " + "+step +inv +proj=cart +ellps=GRS80 " + "+step +inv +proj=vgridshift +grids=au_ga_AGQG_20201120.tif " + "+multiplier=1 " + "+step +proj=unitconvert +xy_in=rad +xy_out=deg " + "+step +proj=axisswap +order=2,1"); +} + +// --------------------------------------------------------------------------- + TEST( operation, compoundCRS_to_compoundCRS_concatenated_operation_with_two_vert_transformation) {