aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/iso19111/operation/coordinateoperation_private.hpp1
-rw-r--r--src/iso19111/operation/coordinateoperationfactory.cpp58
-rw-r--r--src/iso19111/operation/transformation.cpp1
-rw-r--r--test/unit/test_operationfactory.cpp78
4 files changed, 122 insertions, 16 deletions
diff --git a/src/iso19111/operation/coordinateoperation_private.hpp b/src/iso19111/operation/coordinateoperation_private.hpp
index 42bedd4b..2dc70414 100644
--- a/src/iso19111/operation/coordinateoperation_private.hpp
+++ b/src/iso19111/operation/coordinateoperation_private.hpp
@@ -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 {
diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp
index 44d9b570..4832c7b8 100644
--- a/src/iso19111/operation/coordinateoperationfactory.cpp
+++ b/src/iso19111/operation/coordinateoperationfactory.cpp
@@ -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
@@ -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()
@@ -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) {
@@ -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;
@@ -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
@@ -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());
@@ -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());
diff --git a/src/iso19111/operation/transformation.cpp b/src/iso19111/operation/transformation.cpp
index 984f2756..392d6026 100644
--- a/src/iso19111/operation/transformation.cpp
+++ b/src/iso19111/operation/transformation.cpp
@@ -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 &&
diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
index a91d1906..8cc4fa40 100644
--- a/test/unit/test_operationfactory.cpp
+++ b/test/unit/test_operationfactory.cpp
@@ -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");
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 17:15:40 +0000