Skip to content

Commit

Permalink
Merge pull request #3199 from OSGeo/backport-3193-to-9.0
Browse files Browse the repository at this point in the history 
[Backport 9.0] createOperations(): fix/improve result of 'BD72 + Ostend height ' to 'WGS84+EGM96 height' (fixes #3191)
  • Loading branch information
@rouault
rouault committed May 13, 2022
2 parents 5f7e401 + 3c5ef75 commit c58e68d
Show file tree
Hide file tree
Showing 4 changed files with 337 additions and 46 deletions.
72 changes: 36 additions & 36 deletions src/iso19111/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8228,6 +8228,42 @@ const std::string &PROJStringFormatter::toString() const {
}
}

// detect a step and its inverse
if (curStep.inverted != prevStep.inverted &&
curStep.name == prevStep.name &&
curStepParamCount == prevStepParamCount) {
bool allSame = true;
for (size_t j = 0; j < curStepParamCount; j++) {
if (curStep.paramValues[j] != prevStep.paramValues[j]) {
allSame = false;
break;
}
}
if (allSame) {
deletePrevAndCurIter();
continue;
}
}

++iterCur;
}
}

{
auto iterCur = steps.begin();
if (iterCur != steps.end()) {
++iterCur;
}
while (iterCur != steps.end()) {

assert(iterCur != steps.begin());
auto iterPrev = std::prev(iterCur);
auto &prevStep = *iterPrev;
auto &curStep = *iterCur;

const auto curStepParamCount = curStep.paramValues.size();
const auto prevStepParamCount = prevStep.paramValues.size();

// +step +proj=hgridshift +grids=grid_A
// +step +proj=vgridshift [...] <== curStep
// +step +inv +proj=hgridshift +grids=grid_A
Expand Down Expand Up @@ -8266,42 +8302,6 @@ const std::string &PROJStringFormatter::toString() const {
}
}

// detect a step and its inverse
if (curStep.inverted != prevStep.inverted &&
curStep.name == prevStep.name &&
curStepParamCount == prevStepParamCount) {
bool allSame = true;
for (size_t j = 0; j < curStepParamCount; j++) {
if (curStep.paramValues[j] != prevStep.paramValues[j]) {
allSame = false;
break;
}
}
if (allSame) {
deletePrevAndCurIter();
continue;
}
}

++iterCur;
}
}

{
auto iterCur = steps.begin();
if (iterCur != steps.end()) {
++iterCur;
}
while (iterCur != steps.end()) {

assert(iterCur != steps.begin());
auto iterPrev = std::prev(iterCur);
auto &prevStep = *iterPrev;
auto &curStep = *iterCur;

const auto curStepParamCount = curStep.paramValues.size();
const auto prevStepParamCount = prevStep.paramValues.size();

// +step +proj=unitconvert +xy_in=rad +xy_out=deg
// +step +proj=axisswap +order=2,1
// +step +proj=push +v_1 +v_2
Expand Down
23 changes: 20 additions & 3 deletions src/iso19111/operation/concatenatedoperation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -203,13 +203,30 @@ ConcatenatedOperationNNPtr ConcatenatedOperation::create(
}
lastTargetCRS = l_targetCRS;
}

// When chaining VerticalCRS -> GeographicCRS -> VerticalCRS, use
// GeographicCRS as the interpolationCRS
const auto l_sourceCRS = NN_NO_CHECK(operationsIn[0]->sourceCRS());
const auto l_targetCRS = NN_NO_CHECK(operationsIn.back()->targetCRS());
if (operationsIn.size() == 2 && interpolationCRS == nullptr &&
dynamic_cast<const crs::VerticalCRS *>(l_sourceCRS.get()) != nullptr &&
dynamic_cast<const crs::VerticalCRS *>(l_targetCRS.get()) != nullptr) {
const auto geog1 = dynamic_cast<crs::GeographicCRS *>(
operationsIn[0]->targetCRS().get());
const auto geog2 = dynamic_cast<crs::GeographicCRS *>(
operationsIn[1]->sourceCRS().get());
if (geog1 != nullptr && geog2 != nullptr &&
geog1->_isEquivalentTo(geog2,
util::IComparable::Criterion::EQUIVALENT)) {
interpolationCRS = operationsIn[0]->targetCRS();
}
}

auto op = ConcatenatedOperation::nn_make_shared<ConcatenatedOperation>(
operationsIn);
op->assignSelf(op);
op->setProperties(properties);
op->setCRSs(NN_NO_CHECK(operationsIn[0]->sourceCRS()),
NN_NO_CHECK(operationsIn.back()->targetCRS()),
interpolationCRS);
op->setCRSs(l_sourceCRS, l_targetCRS, interpolationCRS);
op->setAccuracies(accuracies);
#ifdef DEBUG_CONCATENATED_OPERATION
{
Expand Down
108 changes: 101 additions & 7 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5454,14 +5454,17 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
return false;
};

const bool hasNonTrivialSrcTransf =
!opsSrcToGeog.empty() &&
(!opsSrcToGeog.front()->hasBallparkTransformation() ||
hasKnownGrid(opsSrcToGeog.front()));
const auto hasNonTrivialTransf =
[&hasKnownGrid](const std::vector<CoordinateOperationNNPtr> &ops) {
return !ops.empty() &&
(!ops.front()->hasBallparkTransformation() ||
hasKnownGrid(ops.front()));
};

const bool hasNonTrivialSrcTransf = hasNonTrivialTransf(opsSrcToGeog);
const bool hasNonTrivialTargetTransf =
!opsGeogToTarget.empty() &&
(!opsGeogToTarget.front()->hasBallparkTransformation() ||
hasKnownGrid(opsGeogToTarget.front()));
hasNonTrivialTransf(opsGeogToTarget);
double bestAccuracy = -1;
if (hasNonTrivialSrcTransf && hasNonTrivialTargetTransf) {
const auto opsGeogSrcToGeogDst =
createOperations(intermGeogSrc, intermGeogDst, context);
Expand Down Expand Up @@ -5502,6 +5505,12 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
op2,
op3},
disallowEmptyIntersection));
const double accuracy = getAccuracy(res.back());
if (accuracy >= 0 &&
(bestAccuracy < 0 ||
accuracy < bestAccuracy)) {
bestAccuracy = accuracy;
}
} catch (const std::exception &) {
}
}
Expand All @@ -5510,6 +5519,91 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
}
}

const auto createOpsInTwoSteps =
[&res,
bestAccuracy](const std::vector<CoordinateOperationNNPtr> &ops1,
const std::vector<CoordinateOperationNNPtr> &ops2) {
std::vector<CoordinateOperationNNPtr> res2;
double bestAccuracy2 = -1;

// In first pass, exclude (horizontal) ballpark operations, but
// accept them in second pass.
for (int pass = 0; pass <= 1 && res2.empty(); pass++) {
for (const auto &op1 : ops1) {
if (pass == 0 && op1->hasBallparkTransformation()) {
// std::cerr << "excluded " << op1->nameStr() <<
// std::endl;
continue;
}
if (op1->nameStr().find(
BALLPARK_VERTICAL_TRANSFORMATION) !=
std::string::npos) {
continue;
}
for (const auto &op2 : ops2) {
if (pass == 0 && op2->hasBallparkTransformation()) {
// std::cerr << "excluded " << op2->nameStr() <<
// std::endl;
continue;
}
if (op2->nameStr().find(
BALLPARK_VERTICAL_TRANSFORMATION) !=
std::string::npos) {
continue;
}
try {
res2.emplace_back(
ConcatenatedOperation::
createComputeMetadata(
{op1, op2},
disallowEmptyIntersection));
const double accuracy =
getAccuracy(res2.back());
if (accuracy >= 0 &&
(bestAccuracy2 < 0 ||
accuracy < bestAccuracy2)) {
bestAccuracy2 = accuracy;
}
} catch (const std::exception &) {
}
}
}
}

// Keep the results of this new attempt, if there are better
// than the previous ones
if (bestAccuracy2 >= 0 &&
(bestAccuracy < 0 || bestAccuracy2 < bestAccuracy)) {
res = std::move(res2);
}
};

// If the promoted-to-3D source geographic CRS is not a known object,
// transformations from it to another 3D one may not be relevant,
// so try doing source -> geogDst 3D -> dest, if geogDst 3D is a known
// object
if (!srcGeog->identifiers().empty() &&
intermGeogSrc->identifiers().empty() &&
!intermGeogDst->identifiers().empty() &&
hasNonTrivialTargetTransf) {
const auto opsSrcToIntermGeog =
createOperations(sourceCRS, intermGeogDst, context);
if (hasNonTrivialTransf(opsSrcToIntermGeog)) {
createOpsInTwoSteps(opsSrcToIntermGeog, opsGeogToTarget);
}
}
// Symetrical situation with the promoted-to-3D target geographic CRS
else if (!dstGeog->identifiers().empty() &&
intermGeogDst->identifiers().empty() &&
!intermGeogSrc->identifiers().empty() &&
hasNonTrivialSrcTransf) {
const auto opsIntermGeogToDst =
createOperations(intermGeogSrc, targetCRS, context);
if (hasNonTrivialTransf(opsIntermGeogToDst)) {
createOpsInTwoSteps(opsSrcToGeog, opsIntermGeogToDst);
}
}

if (!res.empty()) {
return;
}
Expand Down
Loading

0 comments on commit c58e68d

Please sign in to comment.