Skip to content

Commit

Permalink
createOperations(): fine tune for 'GDA94 + AHD height' to 'GDA2020 to…
Browse files Browse the repository at this point in the history
… AVWS height'
  • Loading branch information
rouault committed Dec 8, 2024
1 parent a33372b commit 2f38bbb
Show file tree
Hide file tree
Showing 2 changed files with 61 additions and 17 deletions.
40 changes: 23 additions & 17 deletions src/iso19111/operation/coordinateoperationfactory.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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()) {
Expand Down Expand Up @@ -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());
}
}
}
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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++) {
Expand Down Expand Up @@ -6726,7 +6732,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
try {
res.emplace_back(
ConcatenatedOperation::createComputeMetadata(
intermGeogSrcIsSameAsIntermGeogDst
srcGeogIsSameAsDstGeog
? std::vector<
CoordinateOperationNNPtr>{op1,
op3}
Expand Down
38 changes: 38 additions & 0 deletions test/unit/test_operationfactory.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down

0 comments on commit 2f38bbb

Please sign in to comment.