From 74d98e24fc3607bd1fd5882411392c35607a00f5 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 18 Aug 2023 17:59:17 +0200 Subject: [PATCH] createOperations(): use point motion operations for change of epoch between CoordinateMetadata --- include/proj/coordinateoperation.hpp | 13 + .../operation/coordinateoperationfactory.cpp | 177 +++++++++++--- src/iso19111/operation/singleoperation.cpp | 223 +++++++++++++++++- src/iso19111/operation/transformation.cpp | 2 +- test/unit/test_operation.cpp | 64 +++++ test/unit/test_operationfactory.cpp | 187 +++++++++++++++ 6 files changed, 625 insertions(+), 41 deletions(-) diff --git a/include/proj/coordinateoperation.hpp b/include/proj/coordinateoperation.hpp index 923ded9011..c1e43e9d85 100644 --- a/include/proj/coordinateoperation.hpp +++ b/include/proj/coordinateoperation.hpp @@ -1766,6 +1766,8 @@ class PROJ_GCC_DLL PointMotionOperation : public SingleOperation { PROJ_DLL ~PointMotionOperation() override; //! @endcond + PROJ_DLL const crs::CRSNNPtr &sourceCRS() PROJ_PURE_DECL; + PROJ_DLL CoordinateOperationNNPtr inverse() const override; PROJ_DLL static PointMotionOperationNNPtr @@ -1792,6 +1794,10 @@ class PROJ_GCC_DLL PointMotionOperation : public SingleOperation { PROJ_INTERNAL PointMotionOperationNNPtr shallowClone() const; + PROJ_INTERNAL PointMotionOperationNNPtr + cloneWithEpochs(const common::DataEpoch &sourceEpoch, + const common::DataEpoch &targetEpoch) const; + PROJ_INTERNAL void _exportToPROJString(io::PROJStringFormatter *formatter) const override; // throw(FormattingException) @@ -1814,6 +1820,8 @@ class PROJ_GCC_DLL PointMotionOperation : public SingleOperation { PROJ_INTERNAL CoordinateOperationNNPtr _shallowClone() const override; private: + PROJ_OPAQUE_PRIVATE_DATA + PointMotionOperation &operator=(const PointMotionOperation &) = delete; }; @@ -2097,6 +2105,11 @@ class PROJ_GCC_DLL CoordinateOperationFactory { const coordinates::CoordinateMetadataNNPtr &targetCoordinateMetadata, const CoordinateOperationContextNNPtr &context) const; + PROJ_DLL std::vector createOperations( + const coordinates::CoordinateMetadataNNPtr &sourceCoordinateMetadata, + const coordinates::CoordinateMetadataNNPtr &targetCoordinateMetadata, + const CoordinateOperationContextNNPtr &context) const; + PROJ_DLL static CoordinateOperationFactoryNNPtr create(); protected: diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index 823ba14114..56b7fda5cc 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -547,6 +547,7 @@ struct CoordinateOperationFactory::Private { bool inCreateOperationsWithDatumPivotAntiRecursion = false; bool inCreateOperationsGeogToVertWithAlternativeGeog = false; bool inCreateOperationsGeogToVertWithIntermediateVert = false; + bool ignoreCoordinateEpoch = false; bool skipHorizontalTransformation = false; int nRecLevelCreateOperations = 0; std::map, @@ -624,7 +625,7 @@ struct CoordinateOperationFactory::Private { const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, Private::Context &context, const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst, - std::vector &res); + std::vector &res, bool forceBallpark); static void createOperationsFromSphericalPlanetocentric( const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, @@ -705,7 +706,7 @@ struct CoordinateOperationFactory::Private { std::vector &res, const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, Private::Context &context, const crs::GeographicCRS *geogSrc, - const crs::GeographicCRS *geogDst); + const crs::GeographicCRS *geogDst, bool forceBallpark); static void createOperationsWithDatumPivot( std::vector &res, @@ -2085,17 +2086,16 @@ CoordinateOperationFactory::Private::findsOpsInRegistryWithIntermediate( // --------------------------------------------------------------------------- //! @cond Doxygen_Suppress -static TransformationNNPtr -createBallparkGeographicOffset(const crs::CRSNNPtr &sourceCRS, - const crs::CRSNNPtr &targetCRS, - const io::DatabaseContextPtr &dbContext) { +static TransformationNNPtr createBallparkGeographicOffset( + const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, + const io::DatabaseContextPtr &dbContext, bool forceBallpark) { const crs::GeographicCRS *geogSrc = dynamic_cast(sourceCRS.get()); const crs::GeographicCRS *geogDst = dynamic_cast(targetCRS.get()); const bool isSameDatum = - geogSrc && geogDst && + !forceBallpark && geogSrc && geogDst && geogSrc->datumNonNull(dbContext)->_isEquivalentTo( geogDst->datumNonNull(dbContext).get(), util::IComparable::Criterion::EQUIVALENT, dbContext); @@ -2406,9 +2406,12 @@ createGeodToGeodPROJBased(const crs::CRSNNPtr &geodSrc, util::nn_dynamic_pointer_cast(geodSrc), util::nn_dynamic_pointer_cast(geodDst)); - auto properties = util::PropertyMap().set( - common::IdentifiedObject::NAME_KEY, - buildTransfName(geodSrc->nameStr(), geodDst->nameStr())); + auto properties = + util::PropertyMap() + .set(common::IdentifiedObject::NAME_KEY, + buildTransfName(geodSrc->nameStr(), geodDst->nameStr())) + .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY, + metadata::Extent::WORLD); return createPROJBased(properties, exportable, geodSrc, geodDst, nullptr, {}, false); } @@ -2635,7 +2638,8 @@ std::vector CoordinateOperationFactory::Private::createOperationsGeogToGeog( std::vector &res, const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, Private::Context &context, - const crs::GeographicCRS *geogSrc, const crs::GeographicCRS *geogDst) { + const crs::GeographicCRS *geogSrc, const crs::GeographicCRS *geogDst, + bool forceBallpark) { assert(sourceCRS.get() == geogSrc); assert(targetCRS.get() == geogDst); @@ -2669,9 +2673,11 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog( const auto dbContext = authFactory ? authFactory->databaseContext().as_nullable() : nullptr; - const bool sameDatum = geogSrc->datumNonNull(dbContext)->_isEquivalentTo( - geogDst->datumNonNull(dbContext).get(), - util::IComparable::Criterion::EQUIVALENT, dbContext); + const bool sameDatum = + !forceBallpark && + geogSrc->datumNonNull(dbContext)->_isEquivalentTo( + geogDst->datumNonNull(dbContext).get(), + util::IComparable::Criterion::EQUIVALENT, dbContext); // Do the CRS differ by their axis order ? bool axisReversal2D = false; @@ -2765,8 +2771,8 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog( metadata::Extent::WORLD), datum, dstCS)); - steps.emplace_back( - createBallparkGeographicOffset(sourceCRS, interm_crs, dbContext)); + steps.emplace_back(createBallparkGeographicOffset( + sourceCRS, interm_crs, dbContext, forceBallpark)); steps.emplace_back(Transformation::createLongitudeRotation( util::PropertyMap() @@ -2801,10 +2807,10 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog( metadata::Extent::WORLD), sourceCRS, interm_crs, offset_pm)); steps.emplace_back(createBallparkGeographicOffset( - interm_crs, targetCRS, dbContext)); + interm_crs, targetCRS, dbContext, forceBallpark)); } else { steps.emplace_back(createBallparkGeographicOffset( - sourceCRS, targetCRS, dbContext)); + sourceCRS, targetCRS, dbContext, forceBallpark)); } } @@ -3304,7 +3310,7 @@ CoordinateOperationFactory::Private::createOperations( // Special case if both CRS are geodetic if (geodSrc && geodDst && !derivedSrc && !derivedDst) { createOperationsGeodToGeod(sourceCRS, targetCRS, context, geodSrc, - geodDst, res); + geodDst, res, /*forceBallpark=*/false); return res; } @@ -3515,6 +3521,86 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase( return true; } + // Use PointMotionOperations if appropriate and available + if (geodSrc && geodDst && !context.ignoreCoordinateEpoch && + context.context->getSourceCoordinateEpoch().has_value() && + context.context->getTargetCoordinateEpoch().has_value() && + !context.context->getSourceCoordinateEpoch() + ->coordinateEpoch() + ._isEquivalentTo(context.context->getTargetCoordinateEpoch() + ->coordinateEpoch())) { + struct CoordinateEpochIgnorer { + Context &context; + + explicit CoordinateEpochIgnorer(Context &contextIn) + : context(contextIn) { + assert(!context.ignoreCoordinateEpoch); + context.ignoreCoordinateEpoch = true; + } + + ~CoordinateEpochIgnorer() { context.ignoreCoordinateEpoch = false; } + }; + CoordinateEpochIgnorer guard(context); + + const auto pmoSrc = + context.context->getAuthorityFactory()->getPointMotionOperationsFor( + NN_NO_CHECK( + util::nn_dynamic_pointer_cast(sourceCRS)), + true); + if (!pmoSrc.empty()) { + const auto pmoDst = + context.context->getAuthorityFactory() + ->getPointMotionOperationsFor( + NN_NO_CHECK( + util::nn_dynamic_pointer_cast( + targetCRS)), + true); + if (pmoDst.size() == pmoSrc.size()) { + bool ok = true; + for (size_t i = 0; i < pmoSrc.size(); ++i) { + if (pmoSrc[i]->_isEquivalentTo(pmoDst[i].get())) { + auto pmo = pmoSrc[i]->cloneWithEpochs( + *(context.context->getSourceCoordinateEpoch()), + *(context.context->getTargetCoordinateEpoch())); + std::vector ops; + if (!pmo->sourceCRS()->_isEquivalentTo( + sourceCRS.get(), + util::IComparable::Criterion::EQUIVALENT)) { + auto tmp = createOperations( + sourceCRS, pmo->sourceCRS(), context); + assert(!tmp.empty()); + ops.emplace_back(tmp.front()); + } + ops.emplace_back(pmo); + // pmo->sourceCRS() == pmo->targetCRS() by definition + if (!pmo->sourceCRS()->_isEquivalentTo( + targetCRS.get(), + util::IComparable::Criterion::EQUIVALENT)) { + auto tmp = createOperations(pmo->sourceCRS(), + targetCRS, context); + assert(!tmp.empty()); + ops.emplace_back(tmp.front()); + } + res.emplace_back( + ConcatenatedOperation::createComputeMetadata( + ops, disallowEmptyIntersection)); + } else { + ok = false; + break; + } + } + if (ok) { + std::vector resTmp; + createOperationsGeodToGeod(sourceCRS, targetCRS, context, + geodSrc, geodDst, resTmp, + /*forceBallpark=*/true); + res.insert(res.end(), resTmp.begin(), resTmp.end()); + return true; + } + } + } + } + bool resFindDirectNonEmptyBeforeFiltering = false; res = findOpsInRegistryDirect(sourceCRS, targetCRS, context, resFindDirectNonEmptyBeforeFiltering); @@ -3738,7 +3824,7 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid( opsUnitConvert, tmpCRSWithSrcZ, NN_NO_CHECK(op->sourceCRS()), context, dynamic_cast(tmpCRSWithSrcZ.get()), - opSourceCRSGeog); + opSourceCRSGeog, /*forceBallpark=*/false); assert(opsUnitConvert.size() == 1); opPtr = opsUnitConvert.front().as_nullable(); } @@ -4028,7 +4114,7 @@ std::vector CoordinateOperationFactory::Private:: NN_NO_CHECK(op->sourceCRS()), context, dynamic_cast( tmpCRSWithSrcZ.get()), - tmpCRS); + tmpCRS, /*forceBallpark=*/false); assert(opsUnitConvert.size() == 1); auto concat = ConcatenatedOperation::createComputeMetadata( {opsUnitConvert.front(), op}, disallowEmptyIntersection); @@ -4126,8 +4212,8 @@ void CoordinateOperationFactory::Private:: void CoordinateOperationFactory::Private::createOperationsGeodToGeod( const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, Private::Context &context, const crs::GeodeticCRS *geodSrc, - const crs::GeodeticCRS *geodDst, - std::vector &res) { + const crs::GeodeticCRS *geodDst, std::vector &res, + bool forceBallpark) { ENTER_FUNCTION(); @@ -4162,7 +4248,7 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod( if (geogSrc && geogDst) { createOperationsGeogToGeog(res, sourceCRS, targetCRS, context, geogSrc, - geogDst); + geogDst, forceBallpark); return; } @@ -4187,8 +4273,14 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod( // Same datum ? if (IsSameDatum()) { - res.emplace_back( - Conversion::createGeographicGeocentric(sourceCRS, targetCRS)); + if (forceBallpark) { + auto op = createGeodToGeodPROJBased(sourceCRS, targetCRS); + op->setHasBallparkTransformation(true); + res.emplace_back(op); + } else { + res.emplace_back(Conversion::createGeographicGeocentric( + sourceCRS, targetCRS)); + } } else if (isSrcGeocentric && geogDst) { #if 0 // The below logic was used between PROJ >= 6.0 and < 9.2 @@ -4243,7 +4335,7 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod( // Apply previous case in reverse way std::vector resTmp; createOperationsGeodToGeod(targetCRS, sourceCRS, context, geodDst, - geodSrc, resTmp); + geodSrc, resTmp, forceBallpark); resTmp = applyInverse(resTmp); res.insert(res.end(), resTmp.begin(), resTmp.end()); } @@ -4252,9 +4344,10 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod( } if (isSrcGeocentric && isTargetGeocentric) { - if (sourceCRS->_isEquivalentTo( - targetCRS.get(), util::IComparable::Criterion::EQUIVALENT) || - IsSameDatum()) { + if (!forceBallpark && + (sourceCRS->_isEquivalentTo( + targetCRS.get(), util::IComparable::Criterion::EQUIVALENT) || + IsSameDatum())) { std::string name(NULL_GEOCENTRIC_TRANSLATION); name += " from "; name += sourceCRS->nameStr(); @@ -6415,6 +6508,30 @@ CoordinateOperationFactory::createOperations( // --------------------------------------------------------------------------- +/** \brief Find a list of CoordinateOperation from a source coordinate metadata + * to a target coordinate metadata. + * @param sourceCoordinateMetadata source CoordinateMetadata. + * @param targetCoordinateMetadata target CoordinateMetadata. + * @param context Search context. + * @return a list + * @since 9.4 + */ +std::vector +CoordinateOperationFactory::createOperations( + const coordinates::CoordinateMetadataNNPtr &sourceCoordinateMetadata, + const coordinates::CoordinateMetadataNNPtr &targetCoordinateMetadata, + const CoordinateOperationContextNNPtr &context) const { + auto newContext = context->clone(); + newContext->setSourceCoordinateEpoch( + sourceCoordinateMetadata->coordinateEpoch()); + newContext->setTargetCoordinateEpoch( + targetCoordinateMetadata->coordinateEpoch()); + return createOperations(sourceCoordinateMetadata->crs(), + targetCoordinateMetadata->crs(), newContext); +} + +// --------------------------------------------------------------------------- + /** \brief Instantiate a CoordinateOperationFactory. */ CoordinateOperationFactoryNNPtr CoordinateOperationFactory::create() { diff --git a/src/iso19111/operation/singleoperation.cpp b/src/iso19111/operation/singleoperation.cpp index 46263bdf2e..4af85370d7 100644 --- a/src/iso19111/operation/singleoperation.cpp +++ b/src/iso19111/operation/singleoperation.cpp @@ -4119,6 +4119,40 @@ bool InverseCoordinateOperation::_isEquivalentTo( // --------------------------------------------------------------------------- +//! @cond Doxygen_Suppress +struct PointMotionOperation::Private { + std::shared_ptr> sourceCoordinateEpoch_{ + std::make_shared>()}; + std::shared_ptr> targetCoordinateEpoch_{ + std::make_shared>()}; + + Private() = default; + Private(const Private &) = default; + + void setEpochs(const common::DataEpoch &sourceEpoch, + const common::DataEpoch &targetEpoch) { + sourceCoordinateEpoch_ = + std::make_shared>(sourceEpoch); + targetCoordinateEpoch_ = + std::make_shared>(targetEpoch); + } +}; + +// --------------------------------------------------------------------------- + +// Avoid rounding issues due to year -> second (SI unit) -> year roundtrips +static double getRoundedEpochInDecimalYear(double year) { + // Try to see if the value is close to xxxx.yyy decimal year. + if (std::fabs(1000 * year - std::round(1000 * year)) <= 1e-3) { + year = std::round(1000 * year) / 1000.0; + } + return year; +} + +//! @endcond + +// --------------------------------------------------------------------------- + //! @cond Doxygen_Suppress PointMotionOperation::~PointMotionOperation() = default; //! @endcond @@ -4150,6 +4184,27 @@ PointMotionOperationNNPtr PointMotionOperation::create( crsIn, methodIn, values, accuracies); pmo->assignSelf(pmo); pmo->setProperties(properties); + + const std::string l_name = pmo->nameStr(); + auto pos = l_name.find(" from epoch "); + if (pos != std::string::npos) { + pos += strlen(" from epoch "); + const auto pos2 = l_name.find(" to epoch ", pos); + if (pos2 != std::string::npos) { + const double sourceYear = std::stod(l_name.substr(pos, pos2 - pos)); + const double targetYear = + std::stod(l_name.substr(pos2 + strlen(" to epoch "))); + pmo->d->sourceCoordinateEpoch_ = + std::make_shared>( + common::DataEpoch(common::Measure( + sourceYear, common::UnitOfMeasure::YEAR))); + pmo->d->targetCoordinateEpoch_ = + std::make_shared>( + common::DataEpoch(common::Measure( + targetYear, common::UnitOfMeasure::YEAR))); + } + } + return pmo; } @@ -4202,7 +4257,7 @@ PointMotionOperation::PointMotionOperation( const crs::CRSNNPtr &crsIn, const OperationMethodNNPtr &methodIn, const std::vector &values, const std::vector &accuracies) - : SingleOperation(methodIn) { + : SingleOperation(methodIn), d(internal::make_unique()) { setParameterValues(values); setCRSs(crsIn, crsIn, nullptr); setAccuracies(accuracies); @@ -4211,13 +4266,41 @@ PointMotionOperation::PointMotionOperation( // --------------------------------------------------------------------------- PointMotionOperation::PointMotionOperation(const PointMotionOperation &other) - : CoordinateOperation(other), SingleOperation(other) {} + : CoordinateOperation(other), SingleOperation(other), + d(internal::make_unique(*(other.d))) {} // --------------------------------------------------------------------------- CoordinateOperationNNPtr PointMotionOperation::inverse() const { - return NN_NO_CHECK(std::dynamic_pointer_cast( - shared_from_this().as_nullable())); + auto inverse = shallowClone(); + if (d->sourceCoordinateEpoch_->has_value()) { + // Switch source and target epochs + inverse->d->sourceCoordinateEpoch_ = d->targetCoordinateEpoch_; + inverse->d->targetCoordinateEpoch_ = d->sourceCoordinateEpoch_; + + auto l_name = inverse->nameStr(); + auto pos = l_name.find(" from epoch "); + if (pos != std::string::npos) + l_name.resize(pos); + + const double sourceYear = getRoundedEpochInDecimalYear( + (*(inverse->d->sourceCoordinateEpoch_)) + ->coordinateEpoch() + .convertToUnit(common::UnitOfMeasure::YEAR)); + const double targetYear = getRoundedEpochInDecimalYear( + (*(inverse->d->targetCoordinateEpoch_)) + ->coordinateEpoch() + .convertToUnit(common::UnitOfMeasure::YEAR)); + + l_name += " from epoch "; + l_name += toString(sourceYear); + l_name += " to epoch "; + l_name += toString(targetYear); + util::PropertyMap newProperties; + newProperties.set(IdentifiedObject::NAME_KEY, l_name); + inverse->setProperties(newProperties); + } + return inverse; } // --------------------------------------------------------------------------- @@ -4257,12 +4340,11 @@ PointMotionOperation::substitutePROJAlternativeGridNames( return self; } - auto l_sourceCRS = NN_NO_CHECK(sourceCRS()); auto parameters = std::vector{createOpParamNameEPSGCode( EPSG_CODE_PARAMETER_POINT_MOTION_VELOCITY_GRID_FILE)}; return PointMotionOperation::create( - createSimilarPropertiesOperation(self), l_sourceCRS, + createSimilarPropertiesOperation(self), sourceCRS(), createSimilarPropertiesMethod(method()), parameters, {ParameterValue::createFilename(projFilename)}, coordinateOperationAccuracies()); @@ -4273,6 +4355,16 @@ PointMotionOperation::substitutePROJAlternativeGridNames( // --------------------------------------------------------------------------- +/** \brief Return the source crs::CRS of the operation. + * + * @return the source CRS. + */ +const crs::CRSNNPtr &PointMotionOperation::sourceCRS() PROJ_PURE_DEFN { + return CoordinateOperation::getPrivate()->strongRef_->sourceCRS_; +} + +// --------------------------------------------------------------------------- + //! @cond Doxygen_Suppress PointMotionOperationNNPtr PointMotionOperation::shallowClone() const { @@ -4289,6 +4381,38 @@ CoordinateOperationNNPtr PointMotionOperation::_shallowClone() const { // --------------------------------------------------------------------------- +PointMotionOperationNNPtr PointMotionOperation::cloneWithEpochs( + const common::DataEpoch &sourceEpoch, + const common::DataEpoch &targetEpoch) const { + auto pmo = + PointMotionOperation::nn_make_shared(*this); + + pmo->assignSelf(pmo); + pmo->setCRSs(this, false); + + pmo->d->setEpochs(sourceEpoch, targetEpoch); + + const double sourceYear = getRoundedEpochInDecimalYear( + sourceEpoch.coordinateEpoch().convertToUnit( + common::UnitOfMeasure::YEAR)); + const double targetYear = getRoundedEpochInDecimalYear( + targetEpoch.coordinateEpoch().convertToUnit( + common::UnitOfMeasure::YEAR)); + + auto l_name = nameStr(); + l_name += " from epoch "; + l_name += toString(sourceYear); + l_name += " to epoch "; + l_name += toString(targetYear); + util::PropertyMap newProperties; + newProperties.set(IdentifiedObject::NAME_KEY, l_name); + pmo->setProperties(newProperties); + + return pmo; +} + +// --------------------------------------------------------------------------- + void PointMotionOperation::_exportToWKT(io::WKTFormatter *formatter) const { if (formatter->version() != io::WKTFormatter::Version::WKT2 || !formatter->use2019Keywords()) { @@ -4309,7 +4433,6 @@ void PointMotionOperation::_exportToWKT(io::WKTFormatter *formatter) const { } auto l_sourceCRS = sourceCRS(); - assert(l_sourceCRS); const bool canExportCRSId = !(formatter->idOnTopLevelOnly() && formatter->topLevelHasId()); @@ -4355,10 +4478,90 @@ void PointMotionOperation::_exportToWKT(io::WKTFormatter *formatter) const { // --------------------------------------------------------------------------- void PointMotionOperation::_exportToPROJString( - io::PROJStringFormatter * /*formatter*/) const // throw(FormattingException) + io::PROJStringFormatter *formatter) const // throw(FormattingException) { - throw io::FormattingException( - "CoordinateOperationNNPtr::_exportToPROJString() unimplemented"); + if (formatter->convention() == + io::PROJStringFormatter::Convention::PROJ_4) { + throw io::FormattingException( + "PointMotionOperation cannot be exported as a PROJ.4 string"); + } + + const int methodEPSGCode = method()->getEPSGCode(); + if (methodEPSGCode == + EPSG_CODE_METHOD_POINT_MOTION_BY_GRID_CANADA_NTV2_VEL) { + if (!d->sourceCoordinateEpoch_->has_value()) { + throw io::FormattingException( + "CoordinateOperationNNPtr::_exportToPROJString() unimplemented " + "when source coordinate epoch is missing"); + } + if (!d->targetCoordinateEpoch_->has_value()) { + throw io::FormattingException( + "CoordinateOperationNNPtr::_exportToPROJString() unimplemented " + "when target coordinate epoch is missing"); + } + + auto l_sourceCRS = dynamic_cast(sourceCRS().get()); + + if (!l_sourceCRS->isGeocentric()) { + formatter->startInversion(); + l_sourceCRS->_exportToPROJString(formatter); + formatter->stopInversion(); + + formatter->addStep("cart"); + l_sourceCRS->ellipsoid()->_exportToPROJString(formatter); + } else { + formatter->startInversion(); + l_sourceCRS->addGeocentricUnitConversionIntoPROJString(formatter); + formatter->stopInversion(); + } + + const double sourceYear = getRoundedEpochInDecimalYear( + (*(d->sourceCoordinateEpoch_)) + ->coordinateEpoch() + .convertToUnit(common::UnitOfMeasure::YEAR)); + const double targetYear = getRoundedEpochInDecimalYear( + (*(d->targetCoordinateEpoch_)) + ->coordinateEpoch() + .convertToUnit(common::UnitOfMeasure::YEAR)); + + formatter->addStep("set"); + formatter->addParam("v_4", sourceYear); + formatter->addParam("omit_fwd"); + + formatter->addStep("deformation"); + formatter->addParam("dt", targetYear - sourceYear); + const auto &fileParameter = + parameterValue(EPSG_NAME_PARAMETER_POINT_MOTION_VELOCITY_GRID_FILE, + EPSG_CODE_PARAMETER_POINT_MOTION_VELOCITY_GRID_FILE); + if (fileParameter && + fileParameter->type() == ParameterValue::Type::FILENAME) { + formatter->addParam("grids", fileParameter->valueFile()); + } else { + throw io::FormattingException( + "CoordinateOperationNNPtr::_exportToPROJString(): missing " + "velocity grid file parameter"); + } + + formatter->addStep("set"); + formatter->addParam("v_4", targetYear); + formatter->addParam("omit_inv"); + + if (!l_sourceCRS->isGeocentric()) { + formatter->startInversion(); + formatter->addStep("cart"); + l_sourceCRS->ellipsoid()->_exportToPROJString(formatter); + formatter->stopInversion(); + + l_sourceCRS->_exportToPROJString(formatter); + } else { + l_sourceCRS->addGeocentricUnitConversionIntoPROJString(formatter); + } + + } else { + throw io::FormattingException( + "CoordinateOperationNNPtr::_exportToPROJString() unimplemented for " + "this method"); + } } // --------------------------------------------------------------------------- diff --git a/src/iso19111/operation/transformation.cpp b/src/iso19111/operation/transformation.cpp index d660e269d8..22569c5a5f 100644 --- a/src/iso19111/operation/transformation.cpp +++ b/src/iso19111/operation/transformation.cpp @@ -1851,7 +1851,7 @@ void Transformation::_exportToPROJString( return; } - throw io::FormattingException("Unimplemented"); + throw io::FormattingException("Unimplemented " + nameStr()); } } // namespace operation diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 98fbd2c490..16b2b7eadf 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -5777,3 +5777,67 @@ TEST(operation, export_of_boundCRS_with_proj_string_method) { "+step +proj=axisswap +order=2,1 " "+step +proj=unitconvert +xy_in=rad +xy_out=deg"); } + +// --------------------------------------------------------------------------- + +TEST(operation, PointMotionOperation_with_epochs) { + auto wkt = + "POINTMOTIONOPERATION[\"Canada velocity grid v7 from epoch 2010 to " + "epoch 2002\",\n" + " SOURCECRS[\n" + " GEOGCRS[\"NAD83(CSRS)v7\",\n" + " DATUM[\"North American Datum of 1983 (CSRS) version 7\",\n" + " ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n" + " LENGTHUNIT[\"metre\",1]]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " CS[ellipsoidal,3],\n" + " AXIS[\"geodetic latitude (Lat)\",north,\n" + " ORDER[1],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " AXIS[\"geodetic longitude (Lon)\",east,\n" + " ORDER[2],\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " AXIS[\"ellipsoidal height (h)\",up,\n" + " ORDER[3],\n" + " LENGTHUNIT[\"metre\",1]],\n" + " ID[\"EPSG\",8254]]],\n" + " METHOD[\"Point motion by grid (Canada NTv2_Vel)\",\n" + " ID[\"EPSG\",1070]],\n" + " PARAMETERFILE[\"Point motion velocity grid " + "file\",\"ca_nrc_NAD83v70VG.tif\"],\n" + " OPERATIONACCURACY[0.01],\n" + " ID[\"DERIVED_FROM(EPSG)\",9483],\n" + " REMARK[\"File initially published with name cvg70.cvb, later " + "renamed to NAD83v70VG.gvb with no change of content.\"]]"; + + auto obj = WKTParser().createFromWKT(wkt); + auto pmo = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(pmo != nullptr); + EXPECT_EQ(pmo->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=set +v_4=2010 +omit_fwd " + "+step +proj=deformation +dt=-8 +grids=ca_nrc_NAD83v70VG.tif " + "+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"); + + EXPECT_EQ(pmo->inverse()->nameStr(), + "Canada velocity grid v7 from epoch 2002 to epoch 2010"); + EXPECT_EQ( + pmo->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=GRS80 " + "+step +proj=set +v_4=2002 +omit_fwd " + "+step +proj=deformation +dt=8 +grids=ca_nrc_NAD83v70VG.tif " + "+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"); +} diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index fa8102ff2b..82a58a2c6f 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -37,6 +37,7 @@ #include "proj/common.hpp" #include "proj/coordinateoperation.hpp" +#include "proj/coordinates.hpp" #include "proj/coordinatesystem.hpp" #include "proj/crs.hpp" #include "proj/datum.hpp" @@ -52,6 +53,7 @@ #include using namespace osgeo::proj::common; +using namespace osgeo::proj::coordinates; using namespace osgeo::proj::crs; using namespace osgeo::proj::cs; using namespace osgeo::proj::datum; @@ -9307,3 +9309,188 @@ TEST(operation, "+step +proj=affine +xoff=20 " "+step +proj=axisswap +order=2,1"); } + +// --------------------------------------------------------------------------- + +TEST(operation, createOperation_point_motion_operation_geog2D) { + auto dbContext = DatabaseContext::create(); + auto factory = AuthorityFactory::create(dbContext, "EPSG"); + // "NAD83(CSRS)v7" + auto crs = factory->createCoordinateReferenceSystem("8255"); + auto crs_2002 = CoordinateMetadata::create(crs, 2002.0, dbContext); + auto crs_2010 = CoordinateMetadata::create(crs, 2010.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_2002, crs_2010, ctxt); + ASSERT_EQ(list.size(), 2U); + EXPECT_FALSE(list[0]->hasBallparkTransformation()); + 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=set +v_4=2002 +omit_fwd " + "+step +proj=deformation +dt=8 +grids=ca_nrc_NAD83v70VG.tif " + "+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"); + EXPECT_TRUE(list[1]->hasBallparkTransformation()); + EXPECT_EQ(list[1]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=noop"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, createOperation_point_motion_operation_geog3D) { + auto dbContext = DatabaseContext::create(); + auto factory = AuthorityFactory::create(dbContext, "EPSG"); + // "NAD83(CSRS)v7" + auto crs = factory->createCoordinateReferenceSystem("8254"); + auto crs_2002 = CoordinateMetadata::create(crs, 2002.0, dbContext); + auto crs_2010 = CoordinateMetadata::create(crs, 2010.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_2002, crs_2010, ctxt); + ASSERT_EQ(list.size(), 2U); + EXPECT_FALSE(list[0]->hasBallparkTransformation()); + 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=set +v_4=2002 +omit_fwd " + "+step +proj=deformation +dt=8 +grids=ca_nrc_NAD83v70VG.tif " + "+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"); + EXPECT_TRUE(list[1]->hasBallparkTransformation()); + EXPECT_EQ(list[1]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=noop"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, createOperation_point_motion_operation_geocentric) { + auto dbContext = DatabaseContext::create(); + auto factory = AuthorityFactory::create(dbContext, "EPSG"); + // "NAD83(CSRS)v7" + auto crs = factory->createCoordinateReferenceSystem("8253"); + auto crs_2002 = CoordinateMetadata::create(crs, 2002.0, dbContext); + auto crs_2010 = CoordinateMetadata::create(crs, 2010.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_2002, crs_2010, ctxt); + ASSERT_EQ(list.size(), 2U); + EXPECT_FALSE(list[0]->hasBallparkTransformation()); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=set +v_4=2002 +omit_fwd " + "+step +proj=deformation +dt=8 +grids=ca_nrc_NAD83v70VG.tif " + "+step +proj=set +v_4=2010 +omit_inv"); + EXPECT_TRUE(list[1]->hasBallparkTransformation()); + EXPECT_EQ(list[1]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=noop"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, createOperation_point_motion_operation_geocentric_to_geog3D) { + auto dbContext = DatabaseContext::create(); + auto factory = AuthorityFactory::create(dbContext, "EPSG"); + // "NAD83(CSRS)v7" + auto crs_geocentric = factory->createCoordinateReferenceSystem("8253"); + auto crs_2002 = + CoordinateMetadata::create(crs_geocentric, 2002.0, dbContext); + auto crs_geog3d = factory->createCoordinateReferenceSystem("8254"); + auto crs_2010 = CoordinateMetadata::create(crs_geog3d, 2010.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_2002, crs_2010, ctxt); + ASSERT_EQ(list.size(), 2U); + EXPECT_FALSE(list[0]->hasBallparkTransformation()); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=set +v_4=2002 +omit_fwd " + "+step +proj=deformation +dt=8 +grids=ca_nrc_NAD83v70VG.tif " + "+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"); + EXPECT_TRUE(list[1]->hasBallparkTransformation()); + EXPECT_EQ(list[1]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+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, + createOperation_point_motion_operation_NAD83_CSRS_v7_TO_ITRF2014) { + auto dbContext = DatabaseContext::create(); + auto factory = AuthorityFactory::create(dbContext, "EPSG"); + // "NAD83(CSRS)v7" + auto sourceCRS = factory->createCoordinateReferenceSystem("8254"); + auto crs_2002 = CoordinateMetadata::create(sourceCRS, 2002.0, dbContext); + // ITRF2014 + auto targetCRS = factory->createCoordinateReferenceSystem("7912"); + auto crs_2005 = CoordinateMetadata::create(targetCRS, 2005.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_2002, crs_2005, ctxt); + ASSERT_EQ(list.size(), 2U); + EXPECT_FALSE(list[0]->hasBallparkTransformation()); + 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=set +v_4=2002 +omit_fwd " + "+step +proj=deformation +dt=3 +grids=ca_nrc_NAD83v70VG.tif " + "+step +proj=set +v_4=2005 +omit_inv " + "+step +inv +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 +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"); + EXPECT_TRUE(list[1]->hasBallparkTransformation()); + EXPECT_EQ(list[1]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=noop"); +}