diff --git a/src/iso19111/operation/concatenatedoperation.cpp b/src/iso19111/operation/concatenatedoperation.cpp index ea089ab75c..7f6622b886 100644 --- a/src/iso19111/operation/concatenatedoperation.cpp +++ b/src/iso19111/operation/concatenatedoperation.cpp @@ -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; } @@ -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); + } } // --------------------------------------------------------------------------- diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index 342886885c..7808be957b 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -3075,21 +3075,35 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot( bool resNonEmptyBeforeFiltering; // Deal with potential epoch change - std::vector opsEpochChange; + std::vector opsEpochChangeSrc; + std::vector 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( 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( + candidateDstGeod)), + true); + if (!pmoDst.empty()) { + opsEpochChangeDst = createOperations( + candidateDstGeod, sourceEpoch, candidateDstGeod, + targetEpoch, context); + } } } @@ -3107,9 +3121,9 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot( assert(!opsThird.empty()); const CoordinateOperationNNPtr &opThird(opsThird[0]); - for (size_t iEpochChange = 0; - iEpochChange < std::max(1, opsEpochChange.size()); - ++iEpochChange) { + const auto nIters = std::max( + 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 @@ -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) { diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 6c77e25504..144a01d652 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -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 " @@ -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"); @@ -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"); @@ -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 " @@ -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 " @@ -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()); } @@ -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 " @@ -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"); } // --------------------------------------------------------------------------- @@ -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 " @@ -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"); } // ---------------------------------------------------------------------------