Skip to content

Commit

Permalink
Merge pull request #3943 from OSGeo/backport-3940-to-9.3
Browse files Browse the repository at this point in the history 
[Backport 9.3] createOperations(): fix bad PROJ pipeline when converting between Geog3D with non-metre height to CompoundCRS
  • Loading branch information
@rouault
rouault committed Nov 4, 2023
2 parents aaad5eb + 3b0008a commit a29c310
Show file tree
Hide file tree
Showing 7 changed files with 571 additions and 63 deletions.
4 changes: 4 additions & 0 deletions include/proj/coordinateoperation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -726,6 +726,10 @@ class PROJ_GCC_DLL SingleOperation : virtual public CoordinateOperation {
createOperationParameterValueFromInterpolationCRS(int methodEPSGCode,
int crsEPSGCode);

PROJ_INTERNAL static void
exportToPROJStringChangeVerticalUnit(io::PROJStringFormatter *formatter,
double convFactor);

private:
PROJ_OPAQUE_PRIVATE_DATA
SingleOperation &operator=(const SingleOperation &other) = delete;
Expand Down
208 changes: 206 additions & 2 deletions src/iso19111/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8695,8 +8695,8 @@ const std::string &PROJStringFormatter::toString() const {
}

// +step +proj=unitconvert +xy_in=X1 +z_in=Z1 +xy_out=X2 +z_out=Z2
// +step +proj=unitconvert +z_in=Z2 +z_out=Z3
// ==> step +proj=unitconvert +xy_in=X1 +z_in=Z1 +xy_out=X2
// +step +proj=unitconvert +z_in=Z2 +z_out=Z3
// ==> +step +proj=unitconvert +xy_in=X1 +z_in=Z1 +xy_out=X2
// +z_out=Z3
if (prevStep.name == "unitconvert" &&
curStep.name == "unitconvert" && !prevStep.inverted &&
Expand Down Expand Up @@ -8727,6 +8727,157 @@ const std::string &PROJStringFormatter::toString() const {
continue;
}

// +step +proj=unitconvert +z_in=Z1 +z_out=Z2
// +step +proj=unitconvert +xy_in=X1 +z_in=Z2 +xy_out=X2 +z_out=Z3
// ==> +step +proj=unitconvert +xy_in=X1 +z_in=Z1 +xy_out=X2
// +z_out=Z3
if (prevStep.name == "unitconvert" &&
curStep.name == "unitconvert" && !prevStep.inverted &&
!curStep.inverted && prevStep.paramValues.size() == 2 &&
curStep.paramValues.size() == 4 &&
prevStep.paramValues[0].keyEquals("z_in") &&
prevStep.paramValues[1].keyEquals("z_out") &&
curStep.paramValues[0].keyEquals("xy_in") &&
curStep.paramValues[1].keyEquals("z_in") &&
curStep.paramValues[2].keyEquals("xy_out") &&
curStep.paramValues[3].keyEquals("z_out") &&
prevStep.paramValues[1].value == curStep.paramValues[1].value) {
auto xy_in = curStep.paramValues[0].value;
auto z_in = prevStep.paramValues[0].value;
auto xy_out = curStep.paramValues[2].value;
auto z_out = curStep.paramValues[3].value;

iterCur->paramValues.clear();
iterCur->paramValues.emplace_back(
Step::KeyValue("xy_in", xy_in));
iterCur->paramValues.emplace_back(Step::KeyValue("z_in", z_in));
iterCur->paramValues.emplace_back(
Step::KeyValue("xy_out", xy_out));
iterCur->paramValues.emplace_back(
Step::KeyValue("z_out", z_out));

deletePrevIter();
continue;
}

// +step +proj=unitconvert +xy_in=X1 +z_in=Z1 +xy_out=X2 +z_out=Z2
// +step +proj=unitconvert +xy_in=X2 +xy_out=X3
// ==> +step +proj=unitconvert +xy_in=X1 +z_in=Z1 +xy_out=X3
// +z_out=Z2
if (prevStep.name == "unitconvert" &&
curStep.name == "unitconvert" && !prevStep.inverted &&
!curStep.inverted && prevStep.paramValues.size() == 4 &&
curStep.paramValues.size() == 2 &&
prevStep.paramValues[0].keyEquals("xy_in") &&
prevStep.paramValues[1].keyEquals("z_in") &&
prevStep.paramValues[2].keyEquals("xy_out") &&
prevStep.paramValues[3].keyEquals("z_out") &&
curStep.paramValues[0].keyEquals("xy_in") &&
curStep.paramValues[1].keyEquals("xy_out") &&
prevStep.paramValues[2].value == curStep.paramValues[0].value) {
auto xy_in = prevStep.paramValues[0].value;
auto z_in = prevStep.paramValues[1].value;
auto xy_out = curStep.paramValues[1].value;
auto z_out = prevStep.paramValues[3].value;

iterCur->paramValues.clear();
iterCur->paramValues.emplace_back(
Step::KeyValue("xy_in", xy_in));
iterCur->paramValues.emplace_back(Step::KeyValue("z_in", z_in));
iterCur->paramValues.emplace_back(
Step::KeyValue("xy_out", xy_out));
iterCur->paramValues.emplace_back(
Step::KeyValue("z_out", z_out));

deletePrevIter();
continue;
}

// clang-format off
// A bit odd. Used to simplify geog3d_feet -> EPSG:6318+6360
// of https://github.com/OSGeo/PROJ/issues/3938
// where we get originally
// +step +proj=unitconvert +xy_in=deg +z_in=ft +xy_out=rad +z_out=us-ft
// +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m
// and want it simplified as:
// +step +proj=unitconvert +xy_in=deg +z_in=ft +xy_out=deg +z_out=us-ft
//
// More generally:
// +step +proj=unitconvert +xy_in=X1 +z_in=Z1 +xy_out=X2 +z_out=Z2
// +step +proj=unitconvert +xy_in=X2 +z_in=Z3 +xy_out=X3 +z_out=Z3
// ==> +step +proj=unitconvert +xy_in=X1 +z_in=Z1 +xy_out=X3 +z_out=Z2
// clang-format on
if (prevStep.name == "unitconvert" &&
curStep.name == "unitconvert" && !prevStep.inverted &&
!curStep.inverted && prevStep.paramValues.size() == 4 &&
curStep.paramValues.size() == 4 &&
prevStep.paramValues[0].keyEquals("xy_in") &&
prevStep.paramValues[1].keyEquals("z_in") &&
prevStep.paramValues[2].keyEquals("xy_out") &&
prevStep.paramValues[3].keyEquals("z_out") &&
curStep.paramValues[0].keyEquals("xy_in") &&
curStep.paramValues[1].keyEquals("z_in") &&
curStep.paramValues[2].keyEquals("xy_out") &&
curStep.paramValues[3].keyEquals("z_out") &&
prevStep.paramValues[2].value == curStep.paramValues[0].value &&
curStep.paramValues[1].value == curStep.paramValues[3].value) {
auto xy_in = prevStep.paramValues[0].value;
auto z_in = prevStep.paramValues[1].value;
auto xy_out = curStep.paramValues[2].value;
auto z_out = prevStep.paramValues[3].value;

iterCur->paramValues.clear();
iterCur->paramValues.emplace_back(
Step::KeyValue("xy_in", xy_in));
iterCur->paramValues.emplace_back(Step::KeyValue("z_in", z_in));
iterCur->paramValues.emplace_back(
Step::KeyValue("xy_out", xy_out));
iterCur->paramValues.emplace_back(
Step::KeyValue("z_out", z_out));

deletePrevIter();
continue;
}

// clang-format off
// Variant of above
// +step +proj=unitconvert +xy_in=X1 +z_in=Z1 +xy_out=X2 +z_out=Z1
// +step +proj=unitconvert +xy_in=X2 +z_in=Z2 +xy_out=X3 +z_out=Z3
// ==> +step +proj=unitconvert +xy_in=X1 +z_in=Z2 +xy_out=X3 +z_out=Z3
// clang-format on
if (prevStep.name == "unitconvert" &&
curStep.name == "unitconvert" && !prevStep.inverted &&
!curStep.inverted && prevStep.paramValues.size() == 4 &&
curStep.paramValues.size() == 4 &&
prevStep.paramValues[0].keyEquals("xy_in") &&
prevStep.paramValues[1].keyEquals("z_in") &&
prevStep.paramValues[2].keyEquals("xy_out") &&
prevStep.paramValues[3].keyEquals("z_out") &&
curStep.paramValues[0].keyEquals("xy_in") &&
curStep.paramValues[1].keyEquals("z_in") &&
curStep.paramValues[2].keyEquals("xy_out") &&
curStep.paramValues[3].keyEquals("z_out") &&
prevStep.paramValues[1].value ==
prevStep.paramValues[3].value &&
curStep.paramValues[0].value == prevStep.paramValues[2].value) {
auto xy_in = prevStep.paramValues[0].value;
auto z_in = curStep.paramValues[1].value;
auto xy_out = curStep.paramValues[2].value;
auto z_out = curStep.paramValues[3].value;

iterCur->paramValues.clear();
iterCur->paramValues.emplace_back(
Step::KeyValue("xy_in", xy_in));
iterCur->paramValues.emplace_back(Step::KeyValue("z_in", z_in));
iterCur->paramValues.emplace_back(
Step::KeyValue("xy_out", xy_out));
iterCur->paramValues.emplace_back(
Step::KeyValue("z_out", z_out));

deletePrevIter();
continue;
}

// unitconvert (1), axisswap order=2,1, unitconvert(2) ==>
// axisswap order=2,1, unitconvert (1), unitconvert(2) which
// will get further optimized by previous case
Expand Down Expand Up @@ -9275,6 +9426,59 @@ const std::string &PROJStringFormatter::toString() const {
}
}

{
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();

const auto deletePrevAndCurIter = [&steps, &iterPrev, &iterCur]() {
iterCur = steps.erase(iterPrev, std::next(iterCur));
if (iterCur != steps.begin())
iterCur = std::prev(iterCur);
if (iterCur == steps.begin() && iterCur != steps.end())
++iterCur;
};

// axisswap order=2,1 followed by itself is a no-op
if (curStep.name == "axisswap" && prevStep.name == "axisswap" &&
curStepParamCount == 1 && prevStepParamCount == 1 &&
curStep.paramValues[0].equals("order", "2,1") &&
prevStep.paramValues[0].equals("order", "2,1")) {
deletePrevAndCurIter();
continue;
}

// 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;
}
}

if (steps.size() > 1 ||
(steps.size() == 1 &&
(steps.front().inverted || steps.front().hasKey("omit_inv") ||
Expand Down
24 changes: 1 addition & 23 deletions src/iso19111/operation/conversion.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4144,29 +4144,7 @@ void Conversion::_exportToPROJString(
"requires an input and output vertical CRS");
}
}
auto uom = common::UnitOfMeasure(std::string(), convFactor,
common::UnitOfMeasure::Type::LINEAR)
.exportToPROJString();
auto reverse_uom =
convFactor == 0.0
? std::string()
: common::UnitOfMeasure(std::string(), 1.0 / convFactor,
common::UnitOfMeasure::Type::LINEAR)
.exportToPROJString();
if (uom == "m") {
// do nothing
} else if (!uom.empty()) {
formatter->addStep("unitconvert");
formatter->addParam("z_in", uom);
formatter->addParam("z_out", "m");
} else if (!reverse_uom.empty()) {
formatter->addStep("unitconvert");
formatter->addParam("z_in", "m");
formatter->addParam("z_out", reverse_uom);
} else {
formatter->addStep("affine");
formatter->addParam("s33", convFactor);
}
exportToPROJStringChangeVerticalUnit(formatter, convFactor);
bConversionDone = true;
bEllipsoidParametersDone = true;
} else if (methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC_TOPOCENTRIC) {
Expand Down
4 changes: 0 additions & 4 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2250,16 +2250,12 @@ struct MyPROJStringExportableHorizVertical final
// cppcheck-suppress functionStatic
_exportToPROJString(io::PROJStringFormatter *formatter) const override {

formatter->pushOmitZUnitConversion();

horizTransform->_exportToPROJString(formatter);

formatter->startInversion();
geogDst->addAngularUnitConvertAndAxisSwap(formatter);
formatter->stopInversion();

formatter->popOmitZUnitConversion();

formatter->pushOmitHorizontalConversionInVertTransformation();
verticalTransform->_exportToPROJString(formatter);
formatter->popOmitHorizontalConversionInVertTransformation();
Expand Down
69 changes: 45 additions & 24 deletions src/iso19111/operation/singleoperation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3001,6 +3001,50 @@ static void setupPROJGeodeticTargetCRS(io::PROJStringFormatter *formatter,

// ---------------------------------------------------------------------------

/* static */
void SingleOperation::exportToPROJStringChangeVerticalUnit(
io::PROJStringFormatter *formatter, double convFactor) {

const auto uom = common::UnitOfMeasure(std::string(), convFactor,
common::UnitOfMeasure::Type::LINEAR)
.exportToPROJString();
const auto reverse_uom =
convFactor == 0.0
? std::string()
: common::UnitOfMeasure(std::string(), 1.0 / convFactor,
common::UnitOfMeasure::Type::LINEAR)
.exportToPROJString();
if (uom == "m") {
// do nothing
} else if (!uom.empty()) {
formatter->addStep("unitconvert");
formatter->addParam("z_in", uom);
formatter->addParam("z_out", "m");
} else if (!reverse_uom.empty()) {
formatter->addStep("unitconvert");
formatter->addParam("z_in", "m");
formatter->addParam("z_out", reverse_uom);
} else if (fabs(convFactor -
common::UnitOfMeasure::FOOT.conversionToSI() /
common::UnitOfMeasure::US_FOOT.conversionToSI()) <
1e-10) {
formatter->addStep("unitconvert");
formatter->addParam("z_in", "ft");
formatter->addParam("z_out", "us-ft");
} else if (fabs(convFactor -
common::UnitOfMeasure::US_FOOT.conversionToSI() /
common::UnitOfMeasure::FOOT.conversionToSI()) < 1e-10) {
formatter->addStep("unitconvert");
formatter->addParam("z_in", "us-ft");
formatter->addParam("z_out", "ft");
} else {
formatter->addStep("affine");
formatter->addParam("s33", convFactor);
}
}

// ---------------------------------------------------------------------------

bool SingleOperation::exportToPROJStringGeneric(
io::PROJStringFormatter *formatter) const {
const int methodEPSGCode = method()->getEPSGCode();
Expand Down Expand Up @@ -3124,30 +3168,7 @@ bool SingleOperation::exportToPROJStringGeneric(
if (methodEPSGCode == EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT) {
const double convFactor = parameterValueNumericAsSI(
EPSG_CODE_PARAMETER_UNIT_CONVERSION_SCALAR);
const auto uom =
common::UnitOfMeasure(std::string(), convFactor,
common::UnitOfMeasure::Type::LINEAR)
.exportToPROJString();
const auto reverse_uom =
convFactor == 0.0
? std::string()
: common::UnitOfMeasure(std::string(), 1.0 / convFactor,
common::UnitOfMeasure::Type::LINEAR)
.exportToPROJString();
if (uom == "m") {
// do nothing
} else if (!uom.empty()) {
formatter->addStep("unitconvert");
formatter->addParam("z_in", uom);
formatter->addParam("z_out", "m");
} else if (!reverse_uom.empty()) {
formatter->addStep("unitconvert");
formatter->addParam("z_in", "m");
formatter->addParam("z_out", reverse_uom);
} else {
formatter->addStep("affine");
formatter->addParam("s33", convFactor);
}
exportToPROJStringChangeVerticalUnit(formatter, convFactor);
return true;
}

Expand Down
Loading

0 comments on commit a29c310

Please sign in to comment.