Skip to content

Commit 9911bad

Browse files
rouaultrouault
authored
Merge pull request #4624 from rouault/fix_4618
createOperationsCompoundToGeog(): improvement to make "PNG94 / PNGMG94 zone 54 + Kumul 34 height" to "WGS 84 (G2139)" perform vertical transformation
2 parents b0ace64 + da59216 commit 9911bad

File tree

2 files changed

+266
-92
lines changed

2 files changed

+266
-92
lines changed

‎src/iso19111/operation/coordinateoperationfactory.cpp‎

Lines changed: 185 additions & 92 deletions
Original file line numberDiff line numberDiff line change
@@ -6090,9 +6090,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
60906090
}
60916091

60926092
// Only do a vertical transformation if the target CRS is 3D.
6093-
const auto dstSingle = dynamic_cast<crs::SingleCRS *>(targetCRS.get());
6094-
if (dstSingle &&
6095-
dstSingle->coordinateSystem()->axisList().size() == 2) {
6093+
if (geogDst->coordinateSystem()->axisList().size() <= 2) {
60966094
auto tmp = createOperations(componentsSrc[0], sourceEpoch,
60976095
targetCRS, targetEpoch, context);
60986096
for (const auto &op : tmp) {
@@ -6109,6 +6107,42 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
61096107
horizTransforms = createOperations(componentsSrc[0], sourceEpoch,
61106108
targetCRS, targetEpoch, context);
61116109
}
6110+
6111+
const auto hasAtLeastOneOpWithNonTrivialAndWithAllGrids =
6112+
[&dbContext](
6113+
const std::vector<CoordinateOperationNNPtr> &verticalTransforms,
6114+
CoordinateOperationContext::GridAvailabilityUse
6115+
gridAvailabilityUse,
6116+
bool ignoreMissingGrids,
6117+
bool *foundRegisteredTransform = nullptr) {
6118+
if (foundRegisteredTransform)
6119+
*foundRegisteredTransform = false;
6120+
for (const auto &op : verticalTransforms) {
6121+
if (hasIdentifiers(op) && dbContext) {
6122+
bool missingGrid = false;
6123+
if (!ignoreMissingGrids) {
6124+
const auto gridsNeeded = op->gridsNeeded(
6125+
dbContext,
6126+
gridAvailabilityUse ==
6127+
CoordinateOperationContext::
6128+
GridAvailabilityUse::KNOWN_AVAILABLE);
6129+
for (const auto &gridDesc : gridsNeeded) {
6130+
if (!gridDesc.available) {
6131+
missingGrid = true;
6132+
break;
6133+
}
6134+
}
6135+
}
6136+
if (foundRegisteredTransform)
6137+
*foundRegisteredTransform = true;
6138+
if (!missingGrid) {
6139+
return true;
6140+
}
6141+
}
6142+
}
6143+
return false;
6144+
};
6145+
61126146
std::vector<CoordinateOperationNNPtr> verticalTransforms;
61136147

61146148
if (componentsSrc.size() >= 2 &&
@@ -6127,7 +6161,6 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
61276161
context.skipHorizontalTransformation = false;
61286162
}
61296163
};
6130-
SetSkipHorizontalTransform setSkipHorizontalTransform(context);
61316164

61326165
struct SetGeogCRSOfVertCRS {
61336166
Context &context;
@@ -6148,39 +6181,23 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
61486181
};
61496182
SetGeogCRSOfVertCRS setGeogCRSOfVertCRS(context, srcGeogCRS);
61506183

6151-
verticalTransforms = createOperations(
6152-
componentsSrc[1], util::optional<common::DataEpoch>(),
6153-
targetCRS->promoteTo3D(std::string(), dbContext),
6154-
util::optional<common::DataEpoch>(), context);
6155-
bool foundRegisteredTransformWithAllGridsAvailable = false;
6184+
{
6185+
SetSkipHorizontalTransform setSkipHorizontalTransform(context);
6186+
6187+
verticalTransforms = createOperations(
6188+
componentsSrc[1], util::optional<common::DataEpoch>(),
6189+
targetCRS, util::optional<common::DataEpoch>(), context);
6190+
}
61566191
const auto gridAvailabilityUse =
61576192
context.context->getGridAvailabilityUse();
61586193
const bool ignoreMissingGrids =
61596194
gridAvailabilityUse ==
61606195
CoordinateOperationContext::GridAvailabilityUse::
61616196
IGNORE_GRID_AVAILABILITY;
6162-
for (const auto &op : verticalTransforms) {
6163-
if (hasIdentifiers(op) && dbContext) {
6164-
bool missingGrid = false;
6165-
if (!ignoreMissingGrids) {
6166-
const auto gridsNeeded = op->gridsNeeded(
6167-
dbContext,
6168-
gridAvailabilityUse ==
6169-
CoordinateOperationContext::
6170-
GridAvailabilityUse::KNOWN_AVAILABLE);
6171-
for (const auto &gridDesc : gridsNeeded) {
6172-
if (!gridDesc.available) {
6173-
missingGrid = true;
6174-
break;
6175-
}
6176-
}
6177-
}
6178-
if (!missingGrid) {
6179-
foundRegisteredTransformWithAllGridsAvailable = true;
6180-
break;
6181-
}
6182-
}
6183-
}
6197+
bool foundRegisteredTransformWithAllGridsAvailable =
6198+
hasAtLeastOneOpWithNonTrivialAndWithAllGrids(
6199+
verticalTransforms, gridAvailabilityUse,
6200+
ignoreMissingGrids);
61846201
if (srcGeogCRS &&
61856202
!srcGeogCRS->_isEquivalentTo(
61866203
geogDst, util::IComparable::Criterion::EQUIVALENT) &&
@@ -6190,40 +6207,22 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
61906207
->demoteTo2D(std::string(), dbContext)
61916208
->promoteTo3D(
61926209
std::string(), dbContext,
6193-
geogDst->coordinateSystem()->axisList().size() == 3
6194-
? geogDst->coordinateSystem()->axisList()[2]
6195-
: cs::VerticalCS::createGravityRelatedHeight(
6196-
common::UnitOfMeasure::METRE)
6197-
->axisList()[0]);
6198-
auto verticalTransformsTmp = createOperations(
6199-
componentsSrc[1], util::optional<common::DataEpoch>(),
6200-
geogCRSTmp, util::optional<common::DataEpoch>(), context);
6201-
bool foundRegisteredTransform = false;
6202-
bool foundRegisteredTransformWithAllGridsAvailable2 = false;
6203-
for (const auto &op : verticalTransformsTmp) {
6204-
if (hasIdentifiers(op) && dbContext) {
6205-
bool missingGrid = false;
6206-
if (!ignoreMissingGrids) {
6207-
const auto gridsNeeded = op->gridsNeeded(
6208-
dbContext,
6209-
gridAvailabilityUse ==
6210-
CoordinateOperationContext::
6211-
GridAvailabilityUse::KNOWN_AVAILABLE);
6212-
for (const auto &gridDesc : gridsNeeded) {
6213-
if (!gridDesc.available) {
6214-
missingGrid = true;
6215-
break;
6216-
}
6217-
}
6218-
}
6219-
foundRegisteredTransform = true;
6220-
if (!missingGrid) {
6221-
foundRegisteredTransformWithAllGridsAvailable2 =
6222-
true;
6223-
break;
6224-
}
6225-
}
6210+
geogDst->coordinateSystem()->axisList()[2]);
6211+
std::vector<CoordinateOperationNNPtr> verticalTransformsTmp;
6212+
{
6213+
SetSkipHorizontalTransform setSkipHorizontalTransform(
6214+
context);
6215+
6216+
verticalTransformsTmp = createOperations(
6217+
componentsSrc[1], util::optional<common::DataEpoch>(),
6218+
geogCRSTmp, util::optional<common::DataEpoch>(),
6219+
context);
62266220
}
6221+
bool foundRegisteredTransform = false;
6222+
const bool foundRegisteredTransformWithAllGridsAvailable2 =
6223+
hasAtLeastOneOpWithNonTrivialAndWithAllGrids(
6224+
verticalTransformsTmp, gridAvailabilityUse,
6225+
ignoreMissingGrids, &foundRegisteredTransform);
62276226
if (foundRegisteredTransformWithAllGridsAvailable2 &&
62286227
!foundRegisteredTransformWithAllGridsAvailable) {
62296228
verticalTransforms = std::move(verticalTransformsTmp);
@@ -6233,6 +6232,94 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
62336232
verticalTransformsTmp.end());
62346233
}
62356234
}
6235+
6236+
// Case for EPSG:5550+7651 ("PNG94 / PNGMG94 zone 54 + Kumul 34
6237+
// height") to EPSG:9754 ("WGS 84 (G2139)")
6238+
// (https://github.com/OSGeo/PROJ/issues/4618)
6239+
// where there is a vertical transformation between Kumul 34 height
6240+
// and WGS84 going through EGM96 and the transformation between
6241+
// PNG94 and WGS 84 (G2139) goes through WGS84 We can then get a
6242+
// good result by doing EPSG:5550+7651 -> EPSG:4979 and EPSG:4979 ->
6243+
// EPSG:9754
6244+
const auto hasNonBallparkLambda =
6245+
[](const CoordinateOperationNNPtr &op) {
6246+
return !op->hasBallparkTransformation();
6247+
};
6248+
const bool onlyVerticalBallparkTransformation =
6249+
std::find_if(verticalTransforms.begin(),
6250+
verticalTransforms.end(),
6251+
hasNonBallparkLambda) == verticalTransforms.end();
6252+
if (onlyVerticalBallparkTransformation) {
6253+
std::map<std::string, crs::CRSNNPtr> candidateIntermGeogCRS;
6254+
for (const auto &op : horizTransforms) {
6255+
if (!op->hasBallparkTransformation()) {
6256+
auto concatenatedOp =
6257+
dynamic_cast<const ConcatenatedOperation *>(
6258+
op.get());
6259+
if (concatenatedOp) {
6260+
for (const auto &subOp :
6261+
concatenatedOp->operations()) {
6262+
auto subOpTargetCRS = subOp->targetCRS();
6263+
if (subOpTargetCRS) {
6264+
auto tmpCRS = subOpTargetCRS->promoteTo3D(
6265+
std::string(), dbContext);
6266+
if (!tmpCRS->identifiers().empty() &&
6267+
!tmpCRS->_isEquivalentTo(
6268+
targetCRS.get(),
6269+
util::IComparable::Criterion::
6270+
EQUIVALENT) &&
6271+
candidateIntermGeogCRS.find(
6272+
tmpCRS->nameStr()) ==
6273+
candidateIntermGeogCRS.end()) {
6274+
candidateIntermGeogCRS.insert(
6275+
std::make_pair(tmpCRS->nameStr(),
6276+
std::move(tmpCRS)));
6277+
}
6278+
}
6279+
}
6280+
}
6281+
}
6282+
}
6283+
6284+
for (const auto &[name, geogCRSTmp] : candidateIntermGeogCRS) {
6285+
(void)name;
6286+
const auto verticalTransformsTmp = createOperations(
6287+
componentsSrc[1], util::optional<common::DataEpoch>(),
6288+
geogCRSTmp, util::optional<common::DataEpoch>(),
6289+
context);
6290+
const bool onlyVerticalBallparkTransformationTmp =
6291+
std::find_if(verticalTransformsTmp.begin(),
6292+
verticalTransformsTmp.end(),
6293+
hasNonBallparkLambda) ==
6294+
verticalTransformsTmp.end();
6295+
if (!onlyVerticalBallparkTransformationTmp) {
6296+
const auto ops1 = createOperations(
6297+
sourceCRS, util::optional<common::DataEpoch>(),
6298+
geogCRSTmp, util::optional<common::DataEpoch>(),
6299+
context);
6300+
const auto ops2 = createOperations(
6301+
geogCRSTmp, util::optional<common::DataEpoch>(),
6302+
targetCRS, util::optional<common::DataEpoch>(),
6303+
context);
6304+
for (const auto &op1 : ops1) {
6305+
if (!op1->hasBallparkTransformation()) {
6306+
for (const auto &op2 : ops2) {
6307+
if (!op2->hasBallparkTransformation()) {
6308+
try {
6309+
res.emplace_back(
6310+
ConcatenatedOperation::
6311+
createComputeMetadata(
6312+
{op1, op2},
6313+
disallowEmptyIntersection));
6314+
} catch (const std::exception &) {
6315+
}
6316+
}
6317+
}
6318+
}
6319+
}
6320+
}
6321+
}
6322+
}
62366323
}
62376324

62386325
if (horizTransforms.empty() || verticalTransforms.empty()) {
@@ -6850,9 +6937,9 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
68506937
}
68516938

68526939
const auto createOpsInTwoSteps =
6853-
[&res, bestAccuracy,
6854-
bestStepCount](const std::vector<CoordinateOperationNNPtr> &ops1,
6855-
const std::vector<CoordinateOperationNNPtr> &ops2) {
6940+
[&res, &bestAccuracy, &bestStepCount](
6941+
const std::vector<CoordinateOperationNNPtr> &ops1,
6942+
const std::vector<CoordinateOperationNNPtr> &ops2) {
68566943
std::vector<CoordinateOperationNNPtr> res2;
68576944
double bestAccuracy2 = -1;
68586945
size_t bestStepCount2 = 0;
@@ -6912,35 +6999,41 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
69126999
(bestAccuracy < 0 || (bestAccuracy2 < bestAccuracy ||
69137000
(bestAccuracy2 == bestAccuracy &&
69147001
bestStepCount2 < bestStepCount)))) {
6915-
res = std::move(res2);
7002+
bestAccuracy = bestAccuracy2;
7003+
bestStepCount = bestStepCount2;
7004+
res.insert(res.end(), res2.begin(), res2.end());
69167005
}
69177006
};
69187007

6919-
// If the promoted-to-3D source geographic CRS is not a known object,
6920-
// transformations from it to another 3D one may not be relevant,
6921-
// so try doing source -> geogDst 3D -> dest, if geogDst 3D is a known
6922-
// object
6923-
if (!srcGeog->identifiers().empty() &&
6924-
intermGeogSrc->identifiers().empty() &&
6925-
!intermGeogDst->identifiers().empty() &&
6926-
hasNonTrivialTargetTransf) {
6927-
const auto opsSrcToIntermGeog = createOperations(
6928-
sourceCRS, util::optional<common::DataEpoch>(), intermGeogDst,
6929-
util::optional<common::DataEpoch>(), context);
6930-
if (hasNonTrivialTransf(opsSrcToIntermGeog)) {
6931-
createOpsInTwoSteps(opsSrcToIntermGeog, opsGeogToTarget);
7008+
if (res.empty() || intermGeogSrc->identifiers().empty() ||
7009+
intermGeogDst->identifiers().empty()) {
7010+
// If the promoted-to-3D source geographic CRS is not a known
7011+
// object, transformations from it to another 3D one may not be
7012+
// relevant, so try doing source -> geogDst 3D -> dest, if geogDst
7013+
// 3D is a known object
7014+
if (!srcGeog->identifiers().empty() &&
7015+
!intermGeogDst->identifiers().empty() &&
7016+
hasNonTrivialTargetTransf) {
7017+
const auto opsSrcToIntermGeog = createOperations(
7018+
sourceCRS, util::optional<common::DataEpoch>(),
7019+
intermGeogDst, util::optional<common::DataEpoch>(),
7020+
context);
7021+
if (hasNonTrivialTransf(opsSrcToIntermGeog)) {
7022+
createOpsInTwoSteps(opsSrcToIntermGeog, opsGeogToTarget);
7023+
}
69327024
}
6933-
}
6934-
// Symmetrical situation with the promoted-to-3D target geographic CRS
6935-
else if (!dstGeog->identifiers().empty() &&
6936-
intermGeogDst->identifiers().empty() &&
6937-
!intermGeogSrc->identifiers().empty() &&
6938-
hasNonTrivialSrcTransf) {
6939-
const auto opsIntermGeogToDst = createOperations(
6940-
intermGeogSrc, util::optional<common::DataEpoch>(), targetCRS,
6941-
util::optional<common::DataEpoch>(), context);
6942-
if (hasNonTrivialTransf(opsIntermGeogToDst)) {
6943-
createOpsInTwoSteps(opsSrcToGeog, opsIntermGeogToDst);
7025+
7026+
// Symmetrical situation with the promoted-to-3D target geographic
7027+
// CRS
7028+
if (!dstGeog->identifiers().empty() &&
7029+
!intermGeogSrc->identifiers().empty() &&
7030+
hasNonTrivialSrcTransf) {
7031+
const auto opsIntermGeogToDst = createOperations(
7032+
intermGeogSrc, util::optional<common::DataEpoch>(),
7033+
targetCRS, util::optional<common::DataEpoch>(), context);
7034+
if (hasNonTrivialTransf(opsIntermGeogToDst)) {
7035+
createOpsInTwoSteps(opsSrcToGeog, opsIntermGeogToDst);
7036+
}
69447037
}
69457038
}
69467039

0 commit comments

Comments
 (0)