Skip to content

Commit

Permalink
createOperations(): use point motion operations for change of epoch b…
Browse files Browse the repository at this point in the history 
…etween CoordinateMetadata
  • Loading branch information
@rouault
rouault committed Aug 18, 2023
1 parent 59941b6 commit 74d98e2
Show file tree
Hide file tree
Showing 6 changed files with 625 additions and 41 deletions.
13 changes: 13 additions & 0 deletions include/proj/coordinateoperation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)

Expand All @@ -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;
};

Expand Down Expand Up @@ -2097,6 +2105,11 @@ class PROJ_GCC_DLL CoordinateOperationFactory {
const coordinates::CoordinateMetadataNNPtr &targetCoordinateMetadata,
const CoordinateOperationContextNNPtr &context) const;

PROJ_DLL std::vector<CoordinateOperationNNPtr> createOperations(
const coordinates::CoordinateMetadataNNPtr &sourceCoordinateMetadata,
const coordinates::CoordinateMetadataNNPtr &targetCoordinateMetadata,
const CoordinateOperationContextNNPtr &context) const;

PROJ_DLL static CoordinateOperationFactoryNNPtr create();

protected:
Expand Down
177 changes: 147 additions & 30 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::pair<io::AuthorityFactory::ObjectType, std::string>,
Expand Down Expand Up @@ -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<CoordinateOperationNNPtr> &res);
std::vector<CoordinateOperationNNPtr> &res, bool forceBallpark);

static void createOperationsFromSphericalPlanetocentric(
const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
Expand Down Expand Up @@ -705,7 +706,7 @@ struct CoordinateOperationFactory::Private {
std::vector<CoordinateOperationNNPtr> &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<CoordinateOperationNNPtr> &res,
Expand Down Expand Up @@ -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<const crs::GeographicCRS *>(sourceCRS.get());
const crs::GeographicCRS *geogDst =
dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
const bool isSameDatum =
geogSrc && geogDst &&
!forceBallpark && geogSrc && geogDst &&
geogSrc->datumNonNull(dbContext)->_isEquivalentTo(
geogDst->datumNonNull(dbContext).get(),
util::IComparable::Criterion::EQUIVALENT, dbContext);
Expand Down Expand Up @@ -2406,9 +2406,12 @@ createGeodToGeodPROJBased(const crs::CRSNNPtr &geodSrc,
util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(geodSrc),
util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(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);
}
Expand Down Expand Up @@ -2635,7 +2638,8 @@ std::vector<CoordinateOperationNNPtr>
CoordinateOperationFactory::Private::createOperationsGeogToGeog(
std::vector<CoordinateOperationNNPtr> &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);
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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()
Expand Down Expand Up @@ -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));
}
}

Expand Down Expand Up @@ -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;
}

Expand Down Expand Up @@ -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<crs::GeodeticCRS>(sourceCRS)),
true);
if (!pmoSrc.empty()) {
const auto pmoDst =
context.context->getAuthorityFactory()
->getPointMotionOperationsFor(
NN_NO_CHECK(
util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(
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<operation::CoordinateOperationNNPtr> 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<CoordinateOperationNNPtr> 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);
Expand Down Expand Up @@ -3738,7 +3824,7 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid(
opsUnitConvert, tmpCRSWithSrcZ, NN_NO_CHECK(op->sourceCRS()),
context,
dynamic_cast<const crs::GeographicCRS *>(tmpCRSWithSrcZ.get()),
opSourceCRSGeog);
opSourceCRSGeog, /*forceBallpark=*/false);
assert(opsUnitConvert.size() == 1);
opPtr = opsUnitConvert.front().as_nullable();
}
Expand Down Expand Up @@ -4028,7 +4114,7 @@ std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
NN_NO_CHECK(op->sourceCRS()), context,
dynamic_cast<const crs::GeographicCRS *>(
tmpCRSWithSrcZ.get()),
tmpCRS);
tmpCRS, /*forceBallpark=*/false);
assert(opsUnitConvert.size() == 1);
auto concat = ConcatenatedOperation::createComputeMetadata(
{opsUnitConvert.front(), op}, disallowEmptyIntersection);
Expand Down Expand Up @@ -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<CoordinateOperationNNPtr> &res) {
const crs::GeodeticCRS *geodDst, std::vector<CoordinateOperationNNPtr> &res,
bool forceBallpark) {

ENTER_FUNCTION();

Expand Down Expand Up @@ -4162,7 +4248,7 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod(

if (geogSrc && geogDst) {
createOperationsGeogToGeog(res, sourceCRS, targetCRS, context, geogSrc,
geogDst);
geogDst, forceBallpark);
return;
}

Expand All @@ -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
Expand Down Expand Up @@ -4243,7 +4335,7 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod(
// Apply previous case in reverse way
std::vector<CoordinateOperationNNPtr> resTmp;
createOperationsGeodToGeod(targetCRS, sourceCRS, context, geodDst,
geodSrc, resTmp);
geodSrc, resTmp, forceBallpark);
resTmp = applyInverse(resTmp);
res.insert(res.end(), resTmp.begin(), resTmp.end());
}
Expand All @@ -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();
Expand Down Expand Up @@ -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<CoordinateOperationNNPtr>
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() {
Expand Down
Loading

0 comments on commit 74d98e2

Please sign in to comment.