Skip to content

Commit

Permalink
createOperations(): fix incorrect height transformation between 3D pr…
Browse files Browse the repository at this point in the history 
…omoted RGF93 and CH1903+ (fixes #2541)
  • Loading branch information
@rouault
rouault committed Mar 5, 2021
1 parent f278b5b commit c6a241f
Show file tree
Hide file tree
Showing 4 changed files with 122 additions and 16 deletions.
1 change: 0 additions & 1 deletion src/iso19111/operation/coordinateoperation_private.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,6 @@ struct CoordinateOperation::Private {
util::optional<common::DataEpoch> sourceCoordinateEpoch_{};
util::optional<common::DataEpoch> targetCoordinateEpoch_{};
bool hasBallparkTransformation_ = false;
bool use3DHelmert_ = false;

// do not set this for a ProjectedCRS.definingConversion
struct CRSStrongRef {
Expand Down
58 changes: 44 additions & 14 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2614,9 +2614,13 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
bool isNullFirst) {
const auto opsSecond =
createOperations(candidateSrcGeod, candidateDstGeod, context);
const auto opsThird =
createOperations(candidateDstGeod, targetCRS, context);
const auto opsThird = createOperations(
sourceAndTargetAre3D
? candidateDstGeod->promoteTo3D(std::string(), dbContext)
: candidateDstGeod,
targetCRS, context);
assert(!opsThird.empty());
const CoordinateOperationNNPtr &opThird(opsThird[0]);

for (auto &opSecond : opsSecond) {
// Check that it is not a transformation synthetized by
Expand All @@ -2632,8 +2636,7 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
}

std::vector<CoordinateOperationNNPtr> subOps;
const bool isNullThird =
isNullTransformation(opsThird[0]->nameStr());
const bool isNullThird = isNullTransformation(opThird->nameStr());
CoordinateOperationNNPtr opSecondCloned(
(isNullFirst || isNullThird || sourceAndTargetAre3D)
? opSecond->shallowClone()
Expand Down Expand Up @@ -2665,12 +2668,31 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
}
}
if (sourceAndTargetAre3D) {
opSecondCloned->getPrivate()->use3DHelmert_ = true;
auto invCO = dynamic_cast<InverseCoordinateOperation *>(
opSecondCloned.get());
if (invCO) {
auto invCOForward = invCO->forwardOperation().get();
invCOForward->getPrivate()->use3DHelmert_ = true;

// Force Helmert operations to use the 3D domain, even if the
// ones we found in EPSG are advertized for the 2D domain.
auto concat =
dynamic_cast<ConcatenatedOperation *>(opSecondCloned.get());
if (concat) {
std::vector<CoordinateOperationNNPtr> newSteps;
for (const auto &step : concat->operations()) {
auto newStep = step->shallowClone();
setCRSs(newStep.get(),
newStep->sourceCRS()->promoteTo3D(std::string(),
dbContext),
newStep->targetCRS()->promoteTo3D(std::string(),
dbContext));
newSteps.emplace_back(newStep);
}
opSecondCloned =
ConcatenatedOperation::createComputeMetadata(
newSteps, disallowEmptyIntersection);
} else {
setCRSs(opSecondCloned.get(),
opSecondCloned->sourceCRS()->promoteTo3D(
std::string(), dbContext),
opSecondCloned->targetCRS()->promoteTo3D(
std::string(), dbContext));
}
}
if (isNullFirst) {
Expand All @@ -2685,7 +2707,7 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
subOps.emplace_back(opSecondCloned);
} else {
subOps.emplace_back(opSecondCloned);
subOps.emplace_back(opsThird[0]);
subOps.emplace_back(opThird);
}
#ifdef TRACE_CREATE_OPERATIONS
std::string debugStr;
Expand Down Expand Up @@ -2726,6 +2748,10 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(

for (const auto &candidateSrcGeod : candidatesSrcGeod) {
if (candidateSrcGeod->nameStr() == sourceCRS->nameStr()) {
auto sourceSrcGeodModified(
sourceAndTargetAre3D
? candidateSrcGeod->promoteTo3D(std::string(), dbContext)
: candidateSrcGeod);
for (const auto &candidateDstGeod : candidatesDstGeod) {
if (candidateDstGeod->nameStr() == targetCRS->nameStr()) {
#ifdef TRACE_CREATE_OPERATIONS
Expand All @@ -2734,8 +2760,8 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
objectAsStr(candidateDstGeod.get()) + "->" +
objectAsStr(targetCRS.get()) + ")");
#endif
const auto opsFirst =
createOperations(sourceCRS, candidateSrcGeod, context);
const auto opsFirst = createOperations(
sourceCRS, sourceSrcGeodModified, context);
assert(!opsFirst.empty());
const bool isNullFirst =
isNullTransformation(opsFirst[0]->nameStr());
Expand All @@ -2758,8 +2784,12 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
#ifdef TRACE_CREATE_OPERATIONS
ENTER_BLOCK("");
#endif
auto sourceSrcGeodModified(
sourceAndTargetAre3D
? candidateSrcGeod->promoteTo3D(std::string(), dbContext)
: candidateSrcGeod);
const auto opsFirst =
createOperations(sourceCRS, candidateSrcGeod, context);
createOperations(sourceCRS, sourceSrcGeodModified, context);
assert(!opsFirst.empty());
const bool isNullFirst = isNullTransformation(opsFirst[0]->nameStr());

Expand Down
1 change: 0 additions & 1 deletion src/iso19111/operation/transformation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2595,7 +2595,6 @@ void Transformation::_exportToPROJString(
auto targetCRSGeog =
dynamic_cast<const crs::GeographicCRS *>(targetCRS().get());
const bool addPushPopV3 =
!CoordinateOperation::getPrivate()->use3DHelmert_ &&
((sourceCRSGeog &&
sourceCRSGeog->coordinateSystem()->axisList().size() == 2) ||
(targetCRSGeog &&
Expand Down
78 changes: 78 additions & 0 deletions test/unit/test_operationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1453,6 +1453,84 @@ TEST(operation, geocentricCRS_to_geogCRS_different_datum_context) {

// ---------------------------------------------------------------------------

TEST(operation, geogCRS_3D_to_geogCRS_3D_different_datum_context) {
// Test for https://github.com/OSGeo/PROJ/issues/2541
auto dbContext = DatabaseContext::create();
auto authFactory = AuthorityFactory::create(dbContext, "EPSG");
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
ctxt->setSpatialCriterion(
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
auto list = CoordinateOperationFactory::create()->createOperations(
// RGF93 (3D)
authFactory->createCoordinateReferenceSystem("4965"),
// CH1903+ promoted to 3D
authFactory->createCoordinateReferenceSystem("4150")->promoteTo3D(
std::string(), dbContext),
ctxt);
ASSERT_GE(list.size(), 1U);
EXPECT_EQ(list[0]->nameStr(),
"RGF93 to ETRS89 (1) + Inverse of CH1903+ to ETRS89 (1)");
// Check that there is no +push +v_3
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m "
"+step +proj=cart +ellps=GRS80 "
"+step +proj=helmert +x=-674.374 +y=-15.056 +z=-405.346 "
"+step +inv +proj=cart +ellps=bessel "
"+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
"+step +proj=axisswap +order=2,1");
EXPECT_EQ(list[0]->inverse()->exportToPROJString(
PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m "
"+step +proj=cart +ellps=bessel "
"+step +proj=helmert +x=674.374 +y=15.056 +z=405.346 "
"+step +inv +proj=cart +ellps=GRS80 "
"+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
"+step +proj=axisswap +order=2,1");
}

// ---------------------------------------------------------------------------

TEST(operation, geocentric_to_geogCRS_3D_different_datum_context) {
// Test variant of https://github.com/OSGeo/PROJ/issues/2541
auto dbContext = DatabaseContext::create();
auto authFactory = AuthorityFactory::create(dbContext, "EPSG");
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
ctxt->setSpatialCriterion(
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
auto list = CoordinateOperationFactory::create()->createOperations(
// RGF93 (geocentric)
authFactory->createCoordinateReferenceSystem("4964"),
// CH1903+ promoted to 3D
authFactory->createCoordinateReferenceSystem("4150")->promoteTo3D(
std::string(), dbContext),
ctxt);
ASSERT_GE(list.size(), 1U);
EXPECT_EQ(list[0]->nameStr(),
"Conversion from RGF93 (geocentric) to RGF93 (geog3D) + "
"RGF93 to ETRS89 (1) + "
"Inverse of CH1903+ to ETRS89 (1)");
// Check that there is no +push +v_3
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=helmert +x=-674.374 +y=-15.056 +z=-405.346 "
"+step +inv +proj=cart +ellps=bessel "
"+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
"+step +proj=axisswap +order=2,1");
EXPECT_EQ(list[0]->inverse()->exportToPROJString(
PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m "
"+step +proj=cart +ellps=bessel "
"+step +proj=helmert +x=674.374 +y=15.056 +z=405.346");
}

// ---------------------------------------------------------------------------

TEST(operation, esri_projectedCRS_to_geogCRS_with_ITRF_intermediate_context) {
auto dbContext = DatabaseContext::create();
auto authFactoryEPSG = AuthorityFactory::create(dbContext, "EPSG");
Expand Down

0 comments on commit c6a241f

Please sign in to comment.