@@ -626,6 +626,15 @@ struct CoordinateOperationFactory::Private {
626626 const crs::BoundCRS *boundSrc, const crs::BoundCRS *boundDst,
627627 std::vector<CoordinateOperationNNPtr> &res);
628628
629+ static std::vector<CoordinateOperationNNPtr>
630+ createOperationsEnsembleCRSToOtherGeodCRS (
631+ const crs::CRSNNPtr &sourceCRS,
632+ const util::optional<common::DataEpoch> &sourceEpoch,
633+ const crs::CRSNNPtr &targetCRS,
634+ const util::optional<common::DataEpoch> &targetEpoch,
635+ Private::Context &context, const crs::GeodeticCRS *geodSrc,
636+ const crs::GeodeticCRS *geodDst);
637+
629638 static bool createOperationsFromDatabase (
630639 const crs::CRSNNPtr &sourceCRS,
631640 const util::optional<common::DataEpoch> &sourceEpoch,
@@ -4031,6 +4040,116 @@ useOnlyReplacedOperations(const std::vector<CoordinateOperationNNPtr> &ops) {
40314040
40324041// ---------------------------------------------------------------------------
40334042
4043+ // Below helps for example when doing ETRS89 to Pulkovo 1942(58) which has
4044+ // direct but imprecise operations. By trying to use members of the ETRS89
4045+ // datum ensemble, we can get precise operations using ETRS89-ROU [ETRF2000]
4046+ // for example.
4047+ std::vector<CoordinateOperationNNPtr>
4048+ CoordinateOperationFactory::Private::createOperationsEnsembleCRSToOtherGeodCRS (
4049+ const crs::CRSNNPtr &sourceCRS,
4050+ const util::optional<common::DataEpoch> &sourceEpoch,
4051+ const crs::CRSNNPtr &targetCRS,
4052+ const util::optional<common::DataEpoch> &targetEpoch,
4053+ Private::Context &context, const crs::GeodeticCRS *geodSrc,
4054+ const crs::GeodeticCRS *geodDst) {
4055+ const auto getNatureOf = [](const crs::GeodeticCRS &crs) -> std::string {
4056+ if (dynamic_cast <const crs::GeographicCRS *>(&crs) != nullptr ) {
4057+ if (crs.coordinateSystem ()->axisList ().size () == 2 )
4058+ return " geographic2D" ;
4059+ else
4060+ return " geographic3D" ;
4061+ } else {
4062+ return " geodetic" ;
4063+ }
4064+ };
4065+
4066+ const auto &authFactory = context.context ->getAuthorityFactory ();
4067+ const auto srcEnsemble = geodSrc->datumEnsemble ();
4068+ assert (srcEnsemble);
4069+ const auto &dstDomains = geodDst->domains ();
4070+ const auto &areaOfInterest = context.context ->getAreaOfInterest ();
4071+ const bool is3D = geodSrc->coordinateSystem ()->axisList ().size () == 3 &&
4072+ geodDst->coordinateSystem ()->axisList ().size () == 3 ;
4073+ const char *intermediateType = is3D ? " geographic 3D" : " geographic 2D" ;
4074+
4075+ std::vector<CoordinateOperationNNPtr> newOps;
4076+ const auto &dstExtent = dstDomains[0 ]->domainOfValidity ();
4077+ if (dstExtent) {
4078+ for (const auto &srcDatumCandidate : srcEnsemble->datums ()) {
4079+ const auto &srcDatumDomains = srcDatumCandidate->domains ();
4080+ if (srcDatumDomains.empty ())
4081+ continue ;
4082+
4083+ const auto &srcDatumExtent = srcDatumDomains[0 ]->domainOfValidity ();
4084+ if (!srcDatumExtent ||
4085+ !srcDatumExtent->intersects (NN_NO_CHECK (dstExtent))) {
4086+ continue ;
4087+ }
4088+
4089+ if (areaOfInterest &&
4090+ !srcDatumExtent->intersects (NN_NO_CHECK (areaOfInterest))) {
4091+ continue ;
4092+ }
4093+
4094+ const auto srcDatumGeogCRSList =
4095+ authFactory->createGeodeticCRSFromDatum (
4096+ NN_NO_CHECK (
4097+ std::static_pointer_cast<datum::GeodeticReferenceFrame>(
4098+ srcDatumCandidate.as_nullable ())),
4099+ std::string (), intermediateType);
4100+ for (const auto &srcDatumGeogCRS : srcDatumGeogCRSList) {
4101+ auto ops = createOperations (srcDatumGeogCRS, sourceEpoch,
4102+ targetCRS, targetEpoch, context);
4103+ if (getNatureOf (*(srcDatumGeogCRS.get ())) ==
4104+ getNatureOf (*geodSrc) &&
4105+ srcDatumGeogCRS->coordinateSystem ()->_isEquivalentTo (
4106+ geodSrc->coordinateSystem ().get (),
4107+ util::IComparable::Criterion::EQUIVALENT )) {
4108+ for (const auto &op : ops) {
4109+ if (!op->hasBallparkTransformation ()) {
4110+ auto newOp = op->shallowClone ();
4111+ auto interpolationCRS = newOp->interpolationCRS ();
4112+ if (interpolationCRS &&
4113+ interpolationCRS->_isEquivalentTo (
4114+ srcDatumGeogCRS.get (),
4115+ util::IComparable::Criterion::EQUIVALENT )) {
4116+ newOp->setInterpolationCRS (sourceCRS);
4117+ }
4118+ setCRSs (newOp.get (), sourceCRS, targetCRS);
4119+ newOps.push_back (newOp);
4120+ }
4121+ }
4122+ } else {
4123+ // Expected to get only one as this is a null
4124+ // transformation between an ensemble and one
4125+ // of its member
4126+ const auto opsFirst =
4127+ createOperations (sourceCRS, sourceEpoch,
4128+ srcDatumGeogCRS, targetEpoch, context);
4129+ if (!opsFirst.empty () &&
4130+ !opsFirst.front ()->hasBallparkTransformation ()) {
4131+ const auto &opFirst = opsFirst.front ();
4132+ for (const auto &op : ops) {
4133+ if (!op->hasBallparkTransformation ()) {
4134+ newOps.push_back (
4135+ ConcatenatedOperation::
4136+ createComputeMetadata (
4137+ {opFirst, op},
4138+ context
4139+ .disallowEmptyIntersection ()));
4140+ }
4141+ }
4142+ }
4143+ }
4144+ }
4145+ }
4146+ }
4147+
4148+ return newOps;
4149+ }
4150+
4151+ // ---------------------------------------------------------------------------
4152+
40344153bool CoordinateOperationFactory::Private::createOperationsFromDatabase (
40354154 const crs::CRSNNPtr &sourceCRS,
40364155 const util::optional<common::DataEpoch> &sourceEpoch,
@@ -4360,181 +4479,24 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
43604479 };
43614480 AntiRecursionGuard guard (context);
43624481
4363- const auto getNatureOf =
4364- [](const crs::GeodeticCRS &crs) -> std::string {
4365- if (dynamic_cast <const crs::GeographicCRS *>(&crs) != nullptr ) {
4366- if (crs.coordinateSystem ()->axisList ().size () == 2 )
4367- return " geographic2D" ;
4368- else
4369- return " geographic3D" ;
4370- } else {
4371- return " geodetic" ;
4372- }
4373- };
4374-
4375- const auto &authFactory = context.context ->getAuthorityFactory ();
43764482 const auto srcEnsemble = geodSrc->datumEnsemble ();
43774483 const auto &srcDomains = geodSrc->domains ();
43784484 const auto dstEnsemble = geodDst->datumEnsemble ();
43794485 const auto &dstDomains = geodDst->domains ();
4380- const auto &areaOfInterest = context.context ->getAreaOfInterest ();
4381- const bool is3D = geodSrc->coordinateSystem ()->axisList ().size () == 3 &&
4382- geodDst->coordinateSystem ()->axisList ().size () == 3 ;
4383- const char *intermediateType = is3D ? " geographic 3D" : " geographic 2D" ;
43844486 if (srcEnsemble && !dstDomains.empty ()) {
4385- const auto &dstExtent = dstDomains[0 ]->domainOfValidity ();
4386- if (dstExtent) {
4387- for (const auto &srcDatumCandidate : srcEnsemble->datums ()) {
4388- const auto &srcDatumDomains = srcDatumCandidate->domains ();
4389- if (srcDatumDomains.empty ())
4390- continue ;
4391-
4392- const auto &srcDatumExtent =
4393- srcDatumDomains[0 ]->domainOfValidity ();
4394- if (!srcDatumExtent ||
4395- !srcDatumExtent->intersects (NN_NO_CHECK (dstExtent))) {
4396- continue ;
4397- }
4398-
4399- if (areaOfInterest && !srcDatumExtent->intersects (
4400- NN_NO_CHECK (areaOfInterest))) {
4401- continue ;
4402- }
4403-
4404- const auto srcDatumGeogCRSList =
4405- authFactory->createGeodeticCRSFromDatum (
4406- NN_NO_CHECK (std::static_pointer_cast<
4407- datum::GeodeticReferenceFrame>(
4408- srcDatumCandidate.as_nullable ())),
4409- std::string (), intermediateType);
4410- for (const auto &srcDatumGeogCRS : srcDatumGeogCRSList) {
4411- auto ops =
4412- createOperations (srcDatumGeogCRS, sourceEpoch,
4413- targetCRS, targetEpoch, context);
4414- if (getNatureOf (*(srcDatumGeogCRS.get ())) ==
4415- getNatureOf (*geodSrc) &&
4416- srcDatumGeogCRS->coordinateSystem ()
4417- ->_isEquivalentTo (
4418- geodSrc->coordinateSystem ().get (),
4419- util::IComparable::Criterion::EQUIVALENT )) {
4420- for (const auto &op : ops) {
4421- if (!op->hasBallparkTransformation ()) {
4422- auto newOp = op->shallowClone ();
4423- auto interpolationCRS =
4424- newOp->interpolationCRS ();
4425- if (interpolationCRS &&
4426- interpolationCRS->_isEquivalentTo (
4427- srcDatumGeogCRS.get (),
4428- util::IComparable::Criterion::
4429- EQUIVALENT )) {
4430- newOp->setInterpolationCRS (sourceCRS);
4431- }
4432- setCRSs (newOp.get (), sourceCRS, targetCRS);
4433- res.push_back (newOp);
4434- }
4435- }
4436- } else {
4437- // Expected to get only one as this is a null
4438- // transformation between an ensemble and one
4439- // of its member
4440- const auto opsFirst = createOperations (
4441- sourceCRS, sourceEpoch, srcDatumGeogCRS,
4442- targetEpoch, context);
4443- if (!opsFirst.empty () &&
4444- !opsFirst.front ()
4445- ->hasBallparkTransformation ()) {
4446- const auto &opFirst = opsFirst.front ();
4447- for (const auto &op : ops) {
4448- if (!op->hasBallparkTransformation ()) {
4449- res.push_back (
4450- ConcatenatedOperation::createComputeMetadata (
4451- {opFirst, op},
4452- context
4453- .disallowEmptyIntersection ()));
4454- }
4455- }
4456- }
4457- }
4458- }
4459- }
4460- }
4487+ std::vector<CoordinateOperationNNPtr> newOps =
4488+ createOperationsEnsembleCRSToOtherGeodCRS (
4489+ sourceCRS, sourceEpoch, targetCRS, targetEpoch, context,
4490+ geodSrc, geodDst);
4491+ res.insert (res.end (), std::make_move_iterator (newOps.begin ()),
4492+ std::make_move_iterator (newOps.end ()));
44614493 } else if (dstEnsemble && !srcDomains.empty ()) {
4462- // This is the symmetric case of the previous one
4463-
4464- const auto &srcExtent = srcDomains[0 ]->domainOfValidity ();
4465- if (srcExtent) {
4466- for (const auto &dstDatumCandidate : dstEnsemble->datums ()) {
4467- const auto &dstDatumDomains = dstDatumCandidate->domains ();
4468- if (dstDatumDomains.empty ())
4469- continue ;
4470-
4471- const auto &dstDatumExtent =
4472- dstDatumDomains[0 ]->domainOfValidity ();
4473- if (!dstDatumExtent ||
4474- !dstDatumExtent->intersects (NN_NO_CHECK (srcExtent))) {
4475- continue ;
4476- }
4477-
4478- if (areaOfInterest && !dstDatumExtent->intersects (
4479- NN_NO_CHECK (areaOfInterest))) {
4480- continue ;
4481- }
4482-
4483- const auto dstDatumGeogCRSList =
4484- authFactory->createGeodeticCRSFromDatum (
4485- NN_NO_CHECK (std::static_pointer_cast<
4486- datum::GeodeticReferenceFrame>(
4487- dstDatumCandidate.as_nullable ())),
4488- std::string (), intermediateType);
4489- for (const auto &dstDatumGeogCRS : dstDatumGeogCRSList) {
4490- const auto ops = createOperations (
4491- sourceCRS, sourceEpoch, dstDatumGeogCRS,
4492- targetEpoch, context);
4493- if (getNatureOf (*(dstDatumGeogCRS.get ())) ==
4494- getNatureOf (*geodDst) &&
4495- dstDatumGeogCRS->coordinateSystem ()
4496- ->_isEquivalentTo (
4497- geodDst->coordinateSystem ().get (),
4498- util::IComparable::Criterion::EQUIVALENT )) {
4499- for (const auto &op : ops) {
4500- if (!op->hasBallparkTransformation ()) {
4501- auto newOp = op->shallowClone ();
4502- auto interpolationCRS =
4503- newOp->interpolationCRS ();
4504- if (interpolationCRS &&
4505- interpolationCRS->_isEquivalentTo (
4506- dstDatumGeogCRS.get (),
4507- util::IComparable::Criterion::
4508- EQUIVALENT )) {
4509- newOp->setInterpolationCRS (targetCRS);
4510- }
4511- setCRSs (newOp.get (), sourceCRS, targetCRS);
4512- res.push_back (newOp);
4513- }
4514- }
4515- } else {
4516- // Expected to get only one as this is a null
4517- // transformation between an ensemble and one
4518- // of its member
4519- const auto opsLast = createOperations (
4520- dstDatumGeogCRS, targetEpoch, targetCRS,
4521- targetEpoch, context);
4522- if (!opsLast.empty () &&
4523- !opsLast.front ()->hasBallparkTransformation ()) {
4524- const auto &opLast = opsLast.front ();
4525- for (const auto &op : ops) {
4526- if (!op->hasBallparkTransformation ()) {
4527- res.push_back (
4528- ConcatenatedOperation::createComputeMetadata (
4529- {op, opLast},
4530- context
4531- .disallowEmptyIntersection ()));
4532- }
4533- }
4534- }
4535- }
4536- }
4537- }
4494+ std::vector<CoordinateOperationNNPtr> newOps =
4495+ createOperationsEnsembleCRSToOtherGeodCRS (
4496+ targetCRS, targetEpoch, sourceCRS, sourceEpoch, context,
4497+ geodDst, geodSrc);
4498+ for (const auto &newOp : newOps) {
4499+ res.push_back (newOp->inverse ());
45384500 }
45394501 }
45404502
0 commit comments