Skip to content

Commit

Permalink
createOperations(): fix 'ITRF2014 to NAD83(CSRS)v7' type of operation…
Browse files Browse the repository at this point in the history 
…s, i.e. when the point motion operation is on the end, and also make sure that the time-dependent Helmert transformation receives coordinates in the right epoch
  • Loading branch information
@rouault
rouault committed Sep 6, 2023
1 parent a82c1de commit 69a769b
Show file tree
Hide file tree
Showing 3 changed files with 137 additions and 27 deletions.
26 changes: 26 additions & 0 deletions src/iso19111/operation/concatenatedoperation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -706,6 +706,8 @@ CoordinateOperationNNPtr ConcatenatedOperation::inverse() const {
create(properties, inversedOperations, coordinateOperationAccuracies());
op->d->computedName_ = d->computedName_;
op->setHasBallparkTransformation(hasBallparkTransformation());
op->setSourceCoordinateEpoch(targetCoordinateEpoch());
op->setTargetCoordinateEpoch(sourceCoordinateEpoch());
return op;
}

Expand Down Expand Up @@ -838,9 +840,33 @@ CoordinateOperationNNPtr ConcatenatedOperation::_shallowClone() const {
void ConcatenatedOperation::_exportToPROJString(
io::PROJStringFormatter *formatter) const // throw(FormattingException)
{
double sourceYear =
sourceCoordinateEpoch().has_value()
? getRoundedEpochInDecimalYear(
sourceCoordinateEpoch()->coordinateEpoch().convertToUnit(
common::UnitOfMeasure::YEAR))
: 0;
double targetYear =
targetCoordinateEpoch().has_value()
? getRoundedEpochInDecimalYear(
targetCoordinateEpoch()->coordinateEpoch().convertToUnit(
common::UnitOfMeasure::YEAR))
: 0;
if (sourceYear && !targetYear)
targetYear = sourceYear;
else if (targetYear && !sourceYear)
sourceYear = targetYear;
if (sourceYear) {
formatter->addStep("set");
formatter->addParam("v_4", sourceYear);
}
for (const auto &operation : operations()) {
operation->_exportToPROJString(formatter);
}
if (targetYear) {
formatter->addStep("set");
formatter->addParam("v_4", targetYear);
}
}

// ---------------------------------------------------------------------------
Expand Down
60 changes: 39 additions & 21 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3075,21 +3075,35 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
bool resNonEmptyBeforeFiltering;

// Deal with potential epoch change
std::vector<CoordinateOperationNNPtr> opsEpochChange;
std::vector<CoordinateOperationNNPtr> opsEpochChangeSrc;
std::vector<CoordinateOperationNNPtr> opsEpochChangeDst;
if (sourceEpoch.has_value() && targetEpoch.has_value() &&
!sourceEpoch->coordinateEpoch()._isEquivalentTo(
targetEpoch->coordinateEpoch())) {
const auto pmo =
const auto pmoSrc =
context.context->getAuthorityFactory()
->getPointMotionOperationsFor(
NN_NO_CHECK(
util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(
candidateSrcGeod)),
true);
if (!pmo.empty()) {
opsEpochChange =
if (!pmoSrc.empty()) {
opsEpochChangeSrc =
createOperations(candidateSrcGeod, sourceEpoch,
candidateSrcGeod, targetEpoch, context);
} else {
const auto pmoDst =
context.context->getAuthorityFactory()
->getPointMotionOperationsFor(
NN_NO_CHECK(
util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(
candidateDstGeod)),
true);
if (!pmoDst.empty()) {
opsEpochChangeDst = createOperations(
candidateDstGeod, sourceEpoch, candidateDstGeod,
targetEpoch, context);
}
}
}

Expand All @@ -3107,9 +3121,9 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
assert(!opsThird.empty());
const CoordinateOperationNNPtr &opThird(opsThird[0]);

for (size_t iEpochChange = 0;
iEpochChange < std::max<size_t>(1, opsEpochChange.size());
++iEpochChange) {
const auto nIters = std::max<size_t>(
1, std::max(opsEpochChangeSrc.size(), opsEpochChangeDst.size()));
for (size_t iEpochChange = 0; iEpochChange < nIters; ++iEpochChange) {
for (auto &opSecond : opsSecond) {
// Check that it is not a transformation synthetized by
// ourselves
Expand Down Expand Up @@ -3185,25 +3199,29 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
std::string(), dbContext));
}
}
if (isNullFirst) {
auto oldTarget(
NN_CHECK_ASSERT(opSecondCloned->targetCRS()));
setCRSs(opSecondCloned.get(), sourceCRS, oldTarget);
} else {
if (!isNullFirst) {
subOps.emplace_back(opFirst);
}
if (!opsEpochChange.empty()) {
subOps.emplace_back(opsEpochChange[iEpochChange]);
if (!opsEpochChangeSrc.empty()) {
subOps.emplace_back(opsEpochChangeSrc[iEpochChange]);
}
if (isNullThird) {
auto oldSource(
NN_CHECK_ASSERT(opSecondCloned->sourceCRS()));
setCRSs(opSecondCloned.get(), oldSource, targetCRS);
subOps.emplace_back(opSecondCloned);
} else {
subOps.emplace_back(opSecondCloned);
subOps.emplace_back(opSecondCloned);
if (!opsEpochChangeDst.empty()) {
subOps.emplace_back(opsEpochChangeDst[iEpochChange]);
}
if (!isNullThird) {
subOps.emplace_back(opThird);
}

subOps[0] = subOps[0]->shallowClone();
if (subOps[0]->targetCRS())
setCRSs(subOps[0].get(), sourceCRS,
NN_NO_CHECK(subOps[0]->targetCRS()));
subOps.back() = subOps.back()->shallowClone();
if (subOps[0]->sourceCRS())
setCRSs(subOps.back().get(),
NN_NO_CHECK(subOps.back()->sourceCRS()), targetCRS);

#ifdef TRACE_CREATE_OPERATIONS
std::string debugStr;
for (const auto &op : subOps) {
Expand Down
78 changes: 72 additions & 6 deletions test/unit/test_operationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9403,6 +9403,7 @@ TEST(operation, createOperation_point_motion_operation_geog2D) {
EXPECT_FALSE(list[0]->hasBallparkTransformation());
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=set +v_4=2002 "
"+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 "
Expand All @@ -9412,7 +9413,8 @@ TEST(operation, createOperation_point_motion_operation_geog2D) {
"+step +proj=set +v_4=2010 +omit_inv "
"+step +inv +proj=cart +ellps=GRS80 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
"+step +proj=axisswap +order=2,1 "
"+step +proj=set +v_4=2010");
EXPECT_TRUE(list[1]->hasBallparkTransformation());
EXPECT_EQ(list[1]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=noop");
Expand Down Expand Up @@ -9477,10 +9479,21 @@ TEST(operation, createOperation_point_motion_operation_geocentric) {
EXPECT_FALSE(list[0]->hasBallparkTransformation());
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=set +v_4=2002 "
"+step +proj=set +v_4=2002 +omit_fwd "
"+step +proj=deformation +dt=8 +grids=ca_nrc_NAD83v70VG.tif "
"+ellps=GRS80 "
"+step +proj=set +v_4=2010 +omit_inv");
"+step +proj=set +v_4=2010 +omit_inv "
"+step +proj=set +v_4=2010");
EXPECT_EQ(list[0]->inverse()->exportToPROJString(
PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=set +v_4=2010 "
"+step +proj=set +v_4=2010 +omit_fwd "
"+step +proj=deformation +dt=-8 +grids=ca_nrc_NAD83v70VG.tif "
"+ellps=GRS80 "
"+step +proj=set +v_4=2002 +omit_inv "
"+step +proj=set +v_4=2002");
EXPECT_TRUE(list[1]->hasBallparkTransformation());
EXPECT_EQ(list[1]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=noop");
Expand Down Expand Up @@ -9510,13 +9523,15 @@ TEST(operation, createOperation_point_motion_operation_geocentric_to_geog3D) {
EXPECT_FALSE(list[0]->hasBallparkTransformation());
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=set +v_4=2002 "
"+step +proj=set +v_4=2002 +omit_fwd "
"+step +proj=deformation +dt=8 +grids=ca_nrc_NAD83v70VG.tif "
"+ellps=GRS80 "
"+step +proj=set +v_4=2010 +omit_inv "
"+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");
"+step +proj=axisswap +order=2,1 "
"+step +proj=set +v_4=2010");
EXPECT_TRUE(list[1]->hasBallparkTransformation());
EXPECT_EQ(list[1]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
Expand Down Expand Up @@ -9551,6 +9566,7 @@ TEST(operation,
EXPECT_EQ(
list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=set +v_4=2002 "
"+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 "
Expand All @@ -9565,7 +9581,53 @@ TEST(operation,
"+ds=-7.201e-05 +t_epoch=2010 +convention=position_vector "
"+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");
"+step +proj=axisswap +order=2,1 "
"+step +proj=set +v_4=2005");
EXPECT_TRUE(list[1]->hasBallparkTransformation());
}

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

TEST(operation,
createOperation_point_motion_operation_ITRF2014_to_NAD83_CSRS_v7) {
auto dbContext = DatabaseContext::create();
auto factory = AuthorityFactory::create(dbContext, "EPSG");
// ITRF2014
auto sourceCRS = factory->createCoordinateReferenceSystem("7912");
auto crs_2005 = CoordinateMetadata::create(sourceCRS, 2005.0, dbContext);
// "NAD83(CSRS)v7"
auto targetCRS = factory->createCoordinateReferenceSystem("8254");
auto crs_2002 = CoordinateMetadata::create(targetCRS, 2002.0, dbContext);
auto ctxt = CoordinateOperationContext::create(
AuthorityFactory::create(dbContext, std::string()), nullptr, 0);
ctxt->setSpatialCriterion(
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
ctxt->setGridAvailabilityUse(
CoordinateOperationContext::GridAvailabilityUse::
IGNORE_GRID_AVAILABILITY);
auto list = CoordinateOperationFactory::create()->createOperations(
crs_2005, crs_2002, ctxt);
ASSERT_GE(list.size(), 2U);
EXPECT_FALSE(list[0]->hasBallparkTransformation());
EXPECT_EQ(
list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=set +v_4=2005 "
"+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=1.0053 +y=-1.90921 +z=-0.54157 +rx=-0.02678138 "
"+ry=0.00042027 +rz=-0.01093206 +s=0.00036891 +dx=0.00079 +dy=-0.0006 "
"+dz=-0.00144 +drx=-6.667e-05 +dry=0.00075744 +drz=5.133e-05 "
"+ds=-7.201e-05 +t_epoch=2010 +convention=position_vector "
"+step +proj=set +v_4=2005 +omit_fwd "
"+step +proj=deformation +dt=-3 +grids=ca_nrc_NAD83v70VG.tif "
"+ellps=GRS80 "
"+step +proj=set +v_4=2002 +omit_inv "
"+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 "
"+step +proj=set +v_4=2002");
EXPECT_TRUE(list[1]->hasBallparkTransformation());
}

Expand All @@ -9592,6 +9654,7 @@ TEST(operation,
EXPECT_FALSE(list[0]->hasBallparkTransformation());
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=set +v_4=1997 "
"+step +inv +proj=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 "
"+x_0=304800 +y_0=0 +ellps=GRS80 "
"+step +proj=vgridshift +grids=ca_nrc_HT2_1997.tif +multiplier=1 "
Expand All @@ -9603,7 +9666,8 @@ TEST(operation,
"+step +inv +proj=cart +ellps=GRS80 "
"+step +inv +proj=vgridshift +grids=ca_nrc_CGG2013an83.tif "
"+multiplier=1 "
"+step +proj=utm +zone=19 +ellps=GRS80");
"+step +proj=utm +zone=19 +ellps=GRS80 "
"+step +proj=set +v_4=2010");
}

// ---------------------------------------------------------------------------
Expand Down Expand Up @@ -9862,6 +9926,7 @@ TEST(
EXPECT_FALSE(list[0]->hasBallparkTransformation());
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=set +v_4=1997 "
"+step +inv +proj=utm +zone=17 +ellps=GRS80 "
"+step +proj=vgridshift +grids=ca_nrc_CGG2013an83.tif "
"+multiplier=1 "
Expand All @@ -9873,7 +9938,8 @@ TEST(
"+step +inv +proj=cart +ellps=GRS80 "
"+step +inv +proj=vgridshift +grids=ca_nrc_CGG2013an83.tif "
"+multiplier=1 "
"+step +proj=utm +zone=17 +ellps=GRS80");
"+step +proj=utm +zone=17 +ellps=GRS80 "
"+step +proj=set +v_4=2002");
}

// ---------------------------------------------------------------------------
Expand Down

0 comments on commit 69a769b

Please sign in to comment.