Skip to content

Commit 08e5a7e

Browse files
rouaultrouault
committed
Factor code
1 parent 97ae4fe commit 08e5a7e

1 file changed

Lines changed: 162 additions & 199 deletions

File tree

‎src/iso19111/operation/coordinateoperationfactory.cpp‎

Lines changed: 162 additions & 199 deletions
Original file line numberDiff line numberDiff line change
@@ -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,143 @@ 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+
4056+
struct AntiRecursionGuard {
4057+
Context &context;
4058+
4059+
explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) {
4060+
assert(!context.inCreateOperationsGeogToGeogWithAlternativeGeog);
4061+
context.inCreateOperationsGeogToGeogWithAlternativeGeog = true;
4062+
4063+
assert(!context.inCreateOperationsWithIntermediate);
4064+
context.inCreateOperationsWithIntermediate = true;
4065+
4066+
// This one does really help for performance otherwise
4067+
// "projinfo -s EPSG:7789 -t EPSG:4936 --hide-ballpark
4068+
// --summary" would be really slow. hopefully that does not
4069+
// discard candidate operations...
4070+
assert(!context.inCreateOperationsWithDatumPivotAntiRecursion);
4071+
context.inCreateOperationsWithDatumPivotAntiRecursion = true;
4072+
}
4073+
4074+
~AntiRecursionGuard() {
4075+
context.inCreateOperationsGeogToGeogWithAlternativeGeog = false;
4076+
context.inCreateOperationsWithDatumPivotAntiRecursion = false;
4077+
context.inCreateOperationsWithIntermediate = false;
4078+
}
4079+
};
4080+
AntiRecursionGuard guard(context);
4081+
4082+
const auto getNatureOf = [](const crs::GeodeticCRS &crs) -> std::string {
4083+
if (dynamic_cast<const crs::GeographicCRS *>(&crs) != nullptr) {
4084+
if (crs.coordinateSystem()->axisList().size() == 2)
4085+
return "geographic2D";
4086+
else
4087+
return "geographic3D";
4088+
} else {
4089+
return "geodetic";
4090+
}
4091+
};
4092+
4093+
const auto &authFactory = context.context->getAuthorityFactory();
4094+
const auto srcEnsemble = geodSrc->datumEnsemble();
4095+
assert(srcEnsemble);
4096+
const auto &dstDomains = geodDst->domains();
4097+
const auto &areaOfInterest = context.context->getAreaOfInterest();
4098+
const bool is3D = geodSrc->coordinateSystem()->axisList().size() == 3 &&
4099+
geodDst->coordinateSystem()->axisList().size() == 3;
4100+
const char *intermediateType = is3D ? "geographic 3D" : "geographic 2D";
4101+
4102+
std::vector<CoordinateOperationNNPtr> newOps;
4103+
const auto &dstExtent = dstDomains[0]->domainOfValidity();
4104+
if (dstExtent) {
4105+
for (const auto &srcDatumCandidate : srcEnsemble->datums()) {
4106+
const auto &srcDatumDomains = srcDatumCandidate->domains();
4107+
if (srcDatumDomains.empty())
4108+
continue;
4109+
4110+
const auto &srcDatumExtent = srcDatumDomains[0]->domainOfValidity();
4111+
if (!srcDatumExtent ||
4112+
!srcDatumExtent->intersects(NN_NO_CHECK(dstExtent))) {
4113+
continue;
4114+
}
4115+
4116+
if (areaOfInterest &&
4117+
!srcDatumExtent->intersects(NN_NO_CHECK(areaOfInterest))) {
4118+
continue;
4119+
}
4120+
4121+
const auto srcDatumGeogCRSList =
4122+
authFactory->createGeodeticCRSFromDatum(
4123+
NN_NO_CHECK(
4124+
std::static_pointer_cast<datum::GeodeticReferenceFrame>(
4125+
srcDatumCandidate.as_nullable())),
4126+
std::string(), intermediateType);
4127+
for (const auto &srcDatumGeogCRS : srcDatumGeogCRSList) {
4128+
auto ops = createOperations(srcDatumGeogCRS, sourceEpoch,
4129+
targetCRS, targetEpoch, context);
4130+
if (getNatureOf(*(srcDatumGeogCRS.get())) ==
4131+
getNatureOf(*geodSrc) &&
4132+
srcDatumGeogCRS->coordinateSystem()->_isEquivalentTo(
4133+
geodSrc->coordinateSystem().get(),
4134+
util::IComparable::Criterion::EQUIVALENT)) {
4135+
for (const auto &op : ops) {
4136+
if (!op->hasBallparkTransformation()) {
4137+
auto newOp = op->shallowClone();
4138+
auto interpolationCRS = newOp->interpolationCRS();
4139+
if (interpolationCRS &&
4140+
interpolationCRS->_isEquivalentTo(
4141+
srcDatumGeogCRS.get(),
4142+
util::IComparable::Criterion::EQUIVALENT)) {
4143+
newOp->setInterpolationCRS(sourceCRS);
4144+
}
4145+
setCRSs(newOp.get(), sourceCRS, targetCRS);
4146+
newOps.push_back(newOp);
4147+
}
4148+
}
4149+
} else {
4150+
// Expected to get only one as this is a null
4151+
// transformation between an ensemble and one
4152+
// of its member
4153+
const auto opsFirst =
4154+
createOperations(sourceCRS, sourceEpoch,
4155+
srcDatumGeogCRS, targetEpoch, context);
4156+
if (!opsFirst.empty() &&
4157+
!opsFirst.front()->hasBallparkTransformation()) {
4158+
const auto &opFirst = opsFirst.front();
4159+
for (const auto &op : ops) {
4160+
if (!op->hasBallparkTransformation()) {
4161+
newOps.push_back(
4162+
ConcatenatedOperation::
4163+
createComputeMetadata(
4164+
{opFirst, op},
4165+
context
4166+
.disallowEmptyIntersection()));
4167+
}
4168+
}
4169+
}
4170+
}
4171+
}
4172+
}
4173+
}
4174+
4175+
return newOps;
4176+
}
4177+
4178+
// ---------------------------------------------------------------------------
4179+
40344180
bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
40354181
const crs::CRSNNPtr &sourceCRS,
40364182
const util::optional<common::DataEpoch> &sourceEpoch,
@@ -4332,209 +4478,26 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
43324478
!(sourceCRS->nameStr() == "WGS 84" &&
43334479
(targetCRS->nameStr() == "GDA94" ||
43344480
targetCRS->nameStr() == "GDA2020"))) {
4335-
struct AntiRecursionGuard {
4336-
Context &context;
4337-
4338-
explicit AntiRecursionGuard(Context &contextIn)
4339-
: context(contextIn) {
4340-
assert(
4341-
!context.inCreateOperationsGeogToGeogWithAlternativeGeog);
4342-
context.inCreateOperationsGeogToGeogWithAlternativeGeog = true;
43434481

4344-
assert(!context.inCreateOperationsWithIntermediate);
4345-
context.inCreateOperationsWithIntermediate = true;
4346-
4347-
// This one does really help for performance otherwise
4348-
// "projinfo -s EPSG:7789 -t EPSG:4936 --hide-ballpark
4349-
// --summary" would be really slow. hopefully that does not
4350-
// discard candidate operations...
4351-
assert(!context.inCreateOperationsWithDatumPivotAntiRecursion);
4352-
context.inCreateOperationsWithDatumPivotAntiRecursion = true;
4353-
}
4354-
4355-
~AntiRecursionGuard() {
4356-
context.inCreateOperationsGeogToGeogWithAlternativeGeog = false;
4357-
context.inCreateOperationsWithDatumPivotAntiRecursion = false;
4358-
context.inCreateOperationsWithIntermediate = false;
4359-
}
4360-
};
4361-
AntiRecursionGuard guard(context);
4362-
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();
4377-
const auto &srcDomains = geodSrc->domains();
4378-
const auto dstEnsemble = geodDst->datumEnsemble();
43794483
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";
43844484
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-
}
4461-
} 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-
}
4485+
auto newOps = createOperationsEnsembleCRSToOtherGeodCRS(
4486+
sourceCRS, sourceEpoch, targetCRS, targetEpoch, context,
4487+
geodSrc, geodDst);
4488+
res.insert(res.end(), std::make_move_iterator(newOps.begin()),
4489+
std::make_move_iterator(newOps.end()));
4490+
} else {
4491+
// Symmetrical case
4492+
const auto dstEnsemble = geodDst->datumEnsemble();
4493+
const auto &srcDomains = geodSrc->domains();
4494+
if (dstEnsemble && !srcDomains.empty()) {
4495+
auto newOps =
4496+
applyInverse(createOperationsEnsembleCRSToOtherGeodCRS(
4497+
targetCRS, targetEpoch, sourceCRS, sourceEpoch, context,
4498+
geodDst, geodSrc));
4499+
res.insert(res.end(), std::make_move_iterator(newOps.begin()),
4500+
std::make_move_iterator(newOps.end()));
45384501
}
45394502
}
45404503

0 commit comments

Comments
 (0)