Skip to content

Commit 00d0f25

Browse files
rouaultrouault
committed
createOperationsEnsembleCRSToOtherGeodCRS(): better deal with geographic 3D / geocentric CRS
1 parent aad9921 commit 00d0f25

2 files changed

Lines changed: 87 additions & 19 deletions

File tree

‎src/iso19111/operation/coordinateoperationfactory.cpp‎

Lines changed: 85 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -4082,26 +4082,31 @@ CoordinateOperationFactory::Private::createOperationsEnsembleCRSToOtherGeodCRS(
40824082
const auto getNatureOf = [](const crs::GeodeticCRS &crs) -> std::string {
40834083
if (dynamic_cast<const crs::GeographicCRS *>(&crs) != nullptr) {
40844084
if (crs.coordinateSystem()->axisList().size() == 2)
4085-
return "geographic2D";
4085+
return CRS_SUBTYPE_GEOG_2D;
40864086
else
4087-
return "geographic3D";
4087+
return CRS_SUBTYPE_GEOG_3D;
40884088
} else {
4089-
return "geodetic";
4089+
return CRS_SUBTYPE_GEOCENTRIC;
40904090
}
40914091
};
40924092

40934093
const auto &authFactory = context.context->getAuthorityFactory();
40944094
const auto srcEnsemble = geodSrc->datumEnsemble();
40954095
assert(srcEnsemble);
40964096
const auto &dstDomains = geodDst->domains();
4097+
const auto &dstExtent = dstDomains[0]->domainOfValidity();
4098+
if (!dstExtent) {
4099+
return {};
4100+
}
4101+
40974102
const auto &areaOfInterest = context.context->getAreaOfInterest();
40984103
const bool is3D = geodSrc->coordinateSystem()->axisList().size() == 3 &&
40994104
geodDst->coordinateSystem()->axisList().size() == 3;
4100-
const char *intermediateType = is3D ? "geographic 3D" : "geographic 2D";
4101-
41024105
std::vector<CoordinateOperationNNPtr> newOps;
4103-
const auto &dstExtent = dstDomains[0]->domainOfValidity();
4104-
if (dstExtent) {
4106+
for (int iter = 0; iter < (is3D ? 2 : 1); ++iter) {
4107+
const char *intermediateType =
4108+
is3D ? (iter == 0 ? CRS_SUBTYPE_GEOG_3D : CRS_SUBTYPE_GEOCENTRIC)
4109+
: CRS_SUBTYPE_GEOG_2D;
41054110
for (const auto &srcDatumCandidate : srcEnsemble->datums()) {
41064111
const auto &srcDatumDomains = srcDatumCandidate->domains();
41074112
if (srcDatumDomains.empty())
@@ -4125,8 +4130,29 @@ CoordinateOperationFactory::Private::createOperationsEnsembleCRSToOtherGeodCRS(
41254130
srcDatumCandidate.as_nullable())),
41264131
std::string(), intermediateType);
41274132
for (const auto &srcDatumGeogCRS : srcDatumGeogCRSList) {
4128-
auto ops = createOperations(srcDatumGeogCRS, sourceEpoch,
4129-
targetCRS, targetEpoch, context);
4133+
4134+
std::vector<CoordinateOperationNNPtr> ops;
4135+
std::vector<CoordinateOperationNNPtr> opsLastConversion;
4136+
if (getNatureOf(*(srcDatumGeogCRS.get())) ==
4137+
CRS_SUBTYPE_GEOCENTRIC &&
4138+
getNatureOf(*geodDst) == CRS_SUBTYPE_GEOG_3D) {
4139+
auto tmpList = authFactory->createGeodeticCRSFromDatum(
4140+
NN_CHECK_ASSERT(geodDst->datum()), std::string(),
4141+
CRS_SUBTYPE_GEOCENTRIC);
4142+
if (!tmpList.empty()) {
4143+
ops = createOperations(srcDatumGeogCRS, sourceEpoch,
4144+
tmpList.front(), targetEpoch,
4145+
context);
4146+
opsLastConversion =
4147+
createOperations(tmpList.front(), targetEpoch,
4148+
targetCRS, targetEpoch, context);
4149+
}
4150+
}
4151+
if (ops.empty()) {
4152+
opsLastConversion.clear();
4153+
ops = createOperations(srcDatumGeogCRS, sourceEpoch,
4154+
targetCRS, targetEpoch, context);
4155+
}
41304156
if (getNatureOf(*(srcDatumGeogCRS.get())) ==
41314157
getNatureOf(*geodSrc) &&
41324158
srcDatumGeogCRS->coordinateSystem()->_isEquivalentTo(
@@ -4142,6 +4168,12 @@ CoordinateOperationFactory::Private::createOperationsEnsembleCRSToOtherGeodCRS(
41424168
util::IComparable::Criterion::EQUIVALENT)) {
41434169
newOp->setInterpolationCRS(sourceCRS);
41444170
}
4171+
if (!opsLastConversion.empty()) {
4172+
newOp = ConcatenatedOperation::
4173+
createComputeMetadata(
4174+
{newOp, opsLastConversion.front()},
4175+
context.disallowEmptyIntersection());
4176+
}
41454177
setCRSs(newOp.get(), sourceCRS, targetCRS);
41464178
newOps.push_back(newOp);
41474179
}
@@ -4150,20 +4182,55 @@ CoordinateOperationFactory::Private::createOperationsEnsembleCRSToOtherGeodCRS(
41504182
// Expected to get only one as this is a null
41514183
// transformation between an ensemble and one
41524184
// of its member
4153-
const auto opsFirst =
4154-
createOperations(sourceCRS, sourceEpoch,
4155-
srcDatumGeogCRS, sourceEpoch, context);
4185+
std::vector<CoordinateOperationNNPtr> opsFirstConversion;
4186+
std::vector<CoordinateOperationNNPtr> opsFirst;
4187+
if (getNatureOf(*(srcDatumGeogCRS.get())) ==
4188+
CRS_SUBTYPE_GEOCENTRIC &&
4189+
getNatureOf(*geodSrc) == CRS_SUBTYPE_GEOG_3D) {
4190+
auto tmpList = authFactory->createGeodeticCRSFromDatum(
4191+
geodSrc->datumNonNull(
4192+
authFactory->databaseContext()),
4193+
std::string(), CRS_SUBTYPE_GEOCENTRIC);
4194+
if (!tmpList.empty()) {
4195+
opsFirstConversion = createOperations(
4196+
sourceCRS, sourceEpoch, tmpList.front(),
4197+
sourceEpoch, context);
4198+
context
4199+
.inCreateOperationsWithDatumPivotAntiRecursion =
4200+
false;
4201+
opsFirst = createOperations(
4202+
tmpList.front(), sourceEpoch, srcDatumGeogCRS,
4203+
sourceEpoch, context);
4204+
context
4205+
.inCreateOperationsWithDatumPivotAntiRecursion =
4206+
true;
4207+
}
4208+
}
4209+
if (opsFirst.empty()) {
4210+
opsFirstConversion.clear();
4211+
opsFirst = createOperations(sourceCRS, sourceEpoch,
4212+
srcDatumGeogCRS,
4213+
sourceEpoch, context);
4214+
}
41564215
if (!opsFirst.empty() &&
41574216
!opsFirst.front()->hasBallparkTransformation()) {
41584217
const auto &opFirst = opsFirst.front();
41594218
for (const auto &op : ops) {
41604219
if (!op->hasBallparkTransformation()) {
4161-
newOps.push_back(
4162-
ConcatenatedOperation::
4163-
createComputeMetadata(
4164-
{opFirst, op},
4165-
context
4166-
.disallowEmptyIntersection()));
4220+
std::vector<CoordinateOperationNNPtr> tmp;
4221+
if (!opsFirstConversion.empty()) {
4222+
tmp.push_back(opsFirstConversion.front());
4223+
}
4224+
tmp.push_back(opFirst);
4225+
tmp.push_back(op);
4226+
if (!opsLastConversion.empty()) {
4227+
tmp.push_back(opsLastConversion.front());
4228+
}
4229+
auto newOp = ConcatenatedOperation::
4230+
createComputeMetadata(
4231+
tmp,
4232+
context.disallowEmptyIntersection());
4233+
newOps.push_back(newOp);
41674234
}
41684235
}
41694236
}

‎test/cli/test_projinfo.yaml‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1718,7 +1718,8 @@ tests:
17181718
- comment: Quick check of NKG transformations
17191719
args: -s EPSG:7789 -t EPSG:4936 --area EPSG:1080 --summary --hide-ballpark
17201720
out: |
1721-
Candidate operations found: 1
1721+
Candidate operations found: 2
1722+
EPSG:10894, ITRF2014 to ETRS89-DNK (1), 0.006 m, Denmark - onshore and offshore., time-dependent operation
17221723
NKG:ITRF2014_TO_DK, ITRF2014 to ETRS89(DK), 0.01 m, Denmark - onshore and offshore., time-dependent operation
17231724
- args: --dump-db-structure
17241725
head: 5

0 commit comments

Comments
 (0)