Merge pull request #1261 from rouault/make_iso19111_use_push_pop

EPSG Helmert and Molodensky-Badekas methods in the Geographic 2D domain: use the push/pop v_3 operator to preserve the Z component

8 files changed, 334 insertions, 786 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index d36d901f..14b1024d 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -5489,7 +5489,7 @@ void Conversion::_exportToPROJString(
         if (isUTM(zone, north)) {
             bConversionDone = true;
             formatter->addStep("utm");
-            if( useApprox ) {
+            if (useApprox) {
                 formatter->addParam("approx");
             }
             formatter->addParam("zone", zone);
@@ -8138,6 +8138,18 @@ void Transformation::_exportToPROJString(
         double z =
             parameterValueNumericAsSI(EPSG_CODE_PARAMETER_Z_AXIS_TRANSLATION);
 
+        if (methodEPSGCode == EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_2D ||
+            methodEPSGCode ==
+                EPSG_CODE_METHOD_TIME_DEPENDENT_COORDINATE_FRAME_GEOGRAPHIC_2D ||
+            methodEPSGCode == EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_2D ||
+            methodEPSGCode ==
+                EPSG_CODE_METHOD_TIME_DEPENDENT_POSITION_VECTOR_GEOGRAPHIC_2D ||
+            methodEPSGCode ==
+                EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D) {
+            formatter->addStep("push");
+            formatter->addParam("v_3");
+        }
+
         setupPROJGeodeticSourceCRS(formatter, sourceCRS(), "Helmert");
 
         formatter->addStep("helmert");
@@ -8204,6 +8216,18 @@ void Transformation::_exportToPROJString(
 
         setupPROJGeodeticTargetCRS(formatter, targetCRS(), "Helmert");
 
+        if (methodEPSGCode == EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_2D ||
+            methodEPSGCode ==
+                EPSG_CODE_METHOD_TIME_DEPENDENT_COORDINATE_FRAME_GEOGRAPHIC_2D ||
+            methodEPSGCode == EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_2D ||
+            methodEPSGCode ==
+                EPSG_CODE_METHOD_TIME_DEPENDENT_POSITION_VECTOR_GEOGRAPHIC_2D ||
+            methodEPSGCode ==
+                EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D) {
+            formatter->addStep("pop");
+            formatter->addParam("v_3");
+        }
+
         return;
     }
 
@@ -8250,6 +8274,14 @@ void Transformation::_exportToPROJString(
         double pz = parameterValueNumericAsSI(
             EPSG_CODE_PARAMETER_ORDINATE_3_EVAL_POINT);
 
+        if (methodEPSGCode ==
+                EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_PV_GEOGRAPHIC_2D ||
+            methodEPSGCode ==
+                EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_CF_GEOGRAPHIC_2D) {
+            formatter->addStep("push");
+            formatter->addParam("v_3");
+        }
+
         setupPROJGeodeticSourceCRS(formatter, sourceCRS(),
                                    "Molodensky-Badekas");
 
@@ -8273,6 +8305,14 @@ void Transformation::_exportToPROJString(
         setupPROJGeodeticTargetCRS(formatter, targetCRS(),
                                    "Molodensky-Badekas");
 
+        if (methodEPSGCode ==
+                EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_PV_GEOGRAPHIC_2D ||
+            methodEPSGCode ==
+                EPSG_CODE_METHOD_MOLODENSKY_BADEKAS_CF_GEOGRAPHIC_2D) {
+            formatter->addStep("pop");
+            formatter->addParam("v_3");
+        }
+
         return;
     }
 
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 078180a5..6edce579 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -4837,6 +4837,12 @@ const std::string &PROJStringFormatter::toString() const {
                    step.paramValues[1].keyEquals("xy_out") &&
                    step.paramValues[0].value == step.paramValues[1].value) {
             iter = d->steps_.erase(iter);
+        } else if (step.name == "push" && step.inverted) {
+            step.name = "pop";
+            step.inverted = false;
+        } else if (step.name == "pop" && step.inverted) {
+            step.name = "push";
+            step.inverted = false;
         } else {
             ++iter;
         }
@@ -4912,6 +4918,30 @@ const std::string &PROJStringFormatter::toString() const {
                 break;
             }
 
+            // push v_x followed by pop v_x is a no-op.
+            if (curStep.name == "pop" && prevStep.name == "push" &&
+                !curStep.inverted && !prevStep.inverted &&
+                curStepParamCount == 1 && prevStepParamCount == 1 &&
+                curStep.paramValues[0].key == prevStep.paramValues[0].key) {
+                ++iterCur;
+                d->steps_.erase(iterPrev, iterCur);
+                changeDone = true;
+                break;
+            }
+
+            // pop v_x followed by push v_x is, almost, a no-op. For our
+            // purposes,
+            // we consider it as a no-op for better pipeline optimizations.
+            if (curStep.name == "push" && prevStep.name == "pop" &&
+                !curStep.inverted && !prevStep.inverted &&
+                curStepParamCount == 1 && prevStepParamCount == 1 &&
+                curStep.paramValues[0].key == prevStep.paramValues[0].key) {
+                ++iterCur;
+                d->steps_.erase(iterPrev, iterCur);
+                changeDone = true;
+                break;
+            }
+
             // unitconvert (xy) followed by its inverse is a no-op
             if (curStep.name == "unitconvert" &&
                 prevStep.name == "unitconvert" && !curStep.inverted &&
@@ -5072,6 +5102,71 @@ const std::string &PROJStringFormatter::toString() const {
                 }
             }
 
+            // axisswap order=2,1, pop/push v_3, axisswap order=2,1 -> can
+            // suppress axisswap
+            if (i + 1 < d->steps_.size() && prevStep.name == "axisswap" &&
+                (curStep.name == "push" || curStep.name == "pop") &&
+                prevStepParamCount == 1 &&
+                prevStep.paramValues[0].equals("order", "2,1") &&
+                curStepParamCount == 1 && curStep.paramValues[0].key == "v_3") {
+                auto iterNext = iterCur;
+                ++iterNext;
+                auto &nextStep = *iterNext;
+                if (nextStep.name == "axisswap" &&
+                    nextStep.paramValues.size() == 1 &&
+                    nextStep.paramValues[0].equals("order", "2,1")) {
+                    d->steps_.erase(iterPrev);
+                    d->steps_.erase(iterNext);
+                    changeDone = true;
+                    break;
+                }
+            }
+
+            // push v_3, axisswap order=2,1, pop v_3 -> can suppress push/pop
+            if (i + 1 < d->steps_.size() && prevStep.name == "push" &&
+                prevStepParamCount == 1 &&
+                prevStep.paramValues[0].key == "v_3" &&
+                curStep.name == "axisswap" && curStepParamCount == 1 &&
+                curStep.paramValues[0].equals("order", "2,1")) {
+                auto iterNext = iterCur;
+                ++iterNext;
+                auto &nextStep = *iterNext;
+                if (nextStep.name == "pop" &&
+                    nextStep.paramValues.size() == 1 &&
+                    nextStep.paramValues[0].key == "v_3") {
+                    d->steps_.erase(iterPrev);
+                    d->steps_.erase(iterNext);
+                    changeDone = true;
+                    break;
+                }
+            }
+
+            // unitconvert xy_in=A xy_out=B, pop/push v_3, unitconvert xy_in=B
+            // xy_out=A -> can suppress unitconvert
+            if (i + 1 < d->steps_.size() && prevStep.name == "unitconvert" &&
+                (curStep.name == "push" || curStep.name == "pop") &&
+                prevStepParamCount == 2 &&
+                prevStep.paramValues[0].key == "xy_in" &&
+                prevStep.paramValues[1].key == "xy_out" &&
+                curStepParamCount == 1 && curStep.paramValues[0].key == "v_3") {
+                auto iterNext = iterCur;
+                ++iterNext;
+                auto &nextStep = *iterNext;
+                if (nextStep.name == "unitconvert" &&
+                    nextStep.paramValues.size() == 2 &&
+                    nextStep.paramValues[0].key == "xy_in" &&
+                    nextStep.paramValues[1].key == "xy_out" &&
+                    nextStep.paramValues[0].value ==
+                        prevStep.paramValues[1].value &&
+                    nextStep.paramValues[1].value ==
+                        prevStep.paramValues[0].value) {
+                    d->steps_.erase(iterPrev);
+                    d->steps_.erase(iterNext);
+                    changeDone = true;
+                    break;
+                }
+            }
+
             // for practical purposes WGS84 and GRS80 ellipsoids are
             // equivalents (cartesian transform between both lead to differences
             // of the order of 1e-14 deg..).
@@ -5844,14 +5939,6 @@ struct PROJStringParser::Private {
     CRSNNPtr buildBoundOrCompoundCRSIfNeeded(int iStep, CRSNNPtr crs);
     UnitOfMeasure buildUnit(Step &step, const std::string &unitsParamName,
                             const std::string &toMeterParamName);
-    CoordinateOperationNNPtr buildHelmertTransformation(
-        int iStep, int iFirstAxisSwap = -1, int iFirstUnitConvert = -1,
-        int iFirstGeogStep = -1, int iSecondGeogStep = -1,
-        int iSecondAxisSwap = -1, int iSecondUnitConvert = -1);
-    CoordinateOperationNNPtr buildMolodenskyTransformation(
-        int iStep, int iFirstAxisSwap = -1, int iFirstUnitConvert = -1,
-        int iFirstGeogStep = -1, int iSecondGeogStep = -1,
-        int iSecondAxisSwap = -1, int iSecondUnitConvert = -1);
 
     enum class AxisType { REGULAR, NORTH_POLE, SOUTH_POLE };
 
@@ -7170,261 +7257,6 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
     return crs;
 }
 
-// ---------------------------------------------------------------------------
-
-static bool isDatumDefiningParam(const std::string &param) {
-    return (param == "datum" || param == "ellps" || param == "a" ||
-            param == "b" || param == "rf" || param == "f" || param == "R");
-}
-
-// ---------------------------------------------------------------------------
-
-CoordinateOperationNNPtr PROJStringParser::Private::buildHelmertTransformation(
-    int iStep, int iFirstAxisSwap, int iFirstUnitConvert, int iFirstGeogStep,
-    int iSecondGeogStep, int iSecondAxisSwap, int iSecondUnitConvert) {
-    auto &step = steps_[iStep];
-    auto datum = buildDatum(step, std::string());
-    auto cs = CartesianCS::createGeocentric(UnitOfMeasure::METRE);
-
-    auto mapWithUnknownName = createMapWithUnknownName();
-
-    auto sourceCRS =
-        iFirstGeogStep >= 0
-            ? util::nn_static_pointer_cast<crs::CRS>(
-                  buildGeographicCRS(iFirstGeogStep, iFirstUnitConvert,
-                                     iFirstAxisSwap, true, false))
-            : util::nn_static_pointer_cast<crs::CRS>(
-                  GeodeticCRS::create(mapWithUnknownName, datum, cs));
-    auto targetCRS =
-        iSecondGeogStep >= 0
-            ? util::nn_static_pointer_cast<crs::CRS>(
-                  buildGeographicCRS(iSecondGeogStep, iSecondUnitConvert,
-                                     iSecondAxisSwap, true, false))
-            : util::nn_static_pointer_cast<crs::CRS>(
-                  GeodeticCRS::create(mapWithUnknownName, datum, cs));
-
-    double x = 0;
-    double y = 0;
-    double z = 0;
-    double rx = 0;
-    double ry = 0;
-    double rz = 0;
-    double s = 0;
-    double dx = 0;
-    double dy = 0;
-    double dz = 0;
-    double drx = 0;
-    double dry = 0;
-    double drz = 0;
-    double ds = 0;
-    double t_epoch = 0;
-    bool rotationTerms = false;
-    bool timeDependent = false;
-    bool conventionFound = false;
-    bool positionVectorConvention = false;
-
-    struct Params {
-        double *pValue;
-        const char *name;
-        bool *pPresent;
-    };
-    const Params knownParams[] = {
-        {&x, "x", nullptr},
-        {&y, "y", nullptr},
-        {&z, "z", nullptr},
-        {&rx, "rx", &rotationTerms},
-        {&ry, "ry", &rotationTerms},
-        {&rz, "rz", &rotationTerms},
-        {&s, "s", &rotationTerms},
-        {&dx, "dx", &timeDependent},
-        {&dy, "dy", &timeDependent},
-        {&dz, "dz", &timeDependent},
-        {&drx, "drx", &timeDependent},
-        {&dry, "dry", &timeDependent},
-        {&drz, "drz", &timeDependent},
-        {&ds, "ds", &timeDependent},
-        {&t_epoch, "t_epoch", &timeDependent},
-        {nullptr, "exact", nullptr},
-    };
-
-    for (const auto &param : step.paramValues) {
-        if (isDatumDefiningParam(param.key)) {
-            continue;
-        }
-        if (param.key == "convention") {
-            if (param.value == "position_vector") {
-                positionVectorConvention = true;
-                conventionFound = true;
-            } else if (param.value == "coordinate_frame") {
-                positionVectorConvention = false;
-                conventionFound = true;
-            } else {
-                throw ParsingException("unsupported convention");
-            }
-        } else {
-            bool found = false;
-            for (auto &&knownParam : knownParams) {
-                if (param.key == knownParam.name) {
-                    found = true;
-                    if (knownParam.pValue)
-                        *(knownParam.pValue) = getNumericValue(param.value);
-                    if (knownParam.pPresent)
-                        *(knownParam.pPresent) = true;
-                    break;
-                }
-            }
-            if (!found) {
-                throw ParsingException("unsupported keyword for Helmert: " +
-                                       param.key);
-            }
-        }
-    }
-
-    rotationTerms |= timeDependent;
-    if (rotationTerms && !conventionFound) {
-        throw ParsingException("missing convention");
-    }
-
-    std::vector<metadata::PositionalAccuracyNNPtr> emptyAccuracies;
-
-    auto transf = ([&]() {
-        if (!rotationTerms) {
-            return Transformation::createGeocentricTranslations(
-                mapWithUnknownName, sourceCRS, targetCRS, x, y, z,
-                emptyAccuracies);
-        } else if (positionVectorConvention) {
-            if (timeDependent) {
-                return Transformation::createTimeDependentPositionVector(
-                    mapWithUnknownName, sourceCRS, targetCRS, x, y, z, rx, ry,
-                    rz, s, dx, dy, dz, drx, dry, drz, ds, t_epoch,
-                    emptyAccuracies);
-            } else {
-                return Transformation::createPositionVector(
-                    mapWithUnknownName, sourceCRS, targetCRS, x, y, z, rx, ry,
-                    rz, s, emptyAccuracies);
-            }
-        } else {
-            if (timeDependent) {
-                return Transformation::
-                    createTimeDependentCoordinateFrameRotation(
-                        mapWithUnknownName, sourceCRS, targetCRS, x, y, z, rx,
-                        ry, rz, s, dx, dy, dz, drx, dry, drz, ds, t_epoch,
-                        emptyAccuracies);
-            } else {
-                return Transformation::createCoordinateFrameRotation(
-                    mapWithUnknownName, sourceCRS, targetCRS, x, y, z, rx, ry,
-                    rz, s, emptyAccuracies);
-            }
-        }
-    })();
-
-    if (step.inverted) {
-        return util::nn_static_pointer_cast<CoordinateOperation>(
-            transf->inverse());
-    } else {
-        return util::nn_static_pointer_cast<CoordinateOperation>(transf);
-    }
-}
-
-// ---------------------------------------------------------------------------
-
-CoordinateOperationNNPtr
-PROJStringParser::Private::buildMolodenskyTransformation(
-    int iStep, int iFirstAxisSwap, int iFirstUnitConvert, int iFirstGeogStep,
-    int iSecondGeogStep, int iSecondAxisSwap, int iSecondUnitConvert) {
-    auto &step = steps_[iStep];
-
-    double dx = 0;
-    double dy = 0;
-    double dz = 0;
-    double da = 0;
-    double df = 0;
-
-    struct Params {
-        double *pValue;
-        const char *name;
-    };
-    const Params knownParams[] = {
-        {&dx, "dx"}, {&dy, "dy"}, {&dz, "dz"}, {&da, "da"}, {&df, "df"},
-    };
-    bool abridged = false;
-
-    for (const auto &param : step.paramValues) {
-        if (isDatumDefiningParam(param.key)) {
-            continue;
-        } else if (param.key == "abridged") {
-            abridged = true;
-        } else {
-            bool found = false;
-            for (auto &&knownParam : knownParams) {
-                if (param.key == knownParam.name) {
-                    found = true;
-                    if (knownParam.pValue)
-                        *(knownParam.pValue) = getNumericValue(param.value);
-                    break;
-                }
-            }
-            if (!found) {
-                throw ParsingException("unsupported keyword for Molodensky: " +
-                                       param.key);
-            }
-        }
-    }
-
-    auto datum = buildDatum(step, std::string());
-    auto sourceCRS = iFirstGeogStep >= 0
-                         ? buildGeographicCRS(iFirstGeogStep, iFirstUnitConvert,
-                                              iFirstAxisSwap, true, false)
-                         : buildGeographicCRS(iStep, -1, -1, true, false);
-
-    const auto &ellps = sourceCRS->ellipsoid();
-    const double a = ellps->semiMajorAxis().getSIValue();
-    const double rf = ellps->computedInverseFlattening();
-    const double target_a = a + da;
-    const double target_rf = 1.0 / (1.0 / rf + df);
-
-    auto mapWithUnknownName = createMapWithUnknownName();
-
-    auto target_ellipsoid =
-        Ellipsoid::createFlattenedSphere(mapWithUnknownName, Length(target_a),
-                                         Scale(target_rf))
-            ->identify();
-    auto target_datum = GeodeticReferenceFrame::create(
-        mapWithUnknownName, target_ellipsoid, util::optional<std::string>(),
-        PrimeMeridian::GREENWICH);
-
-    auto targetCRS = util::nn_static_pointer_cast<crs::CRS>(
-        iSecondGeogStep >= 0
-            ? buildGeographicCRS(iSecondGeogStep, iSecondUnitConvert,
-                                 iSecondAxisSwap, true, false)
-            : GeographicCRS::create(mapWithUnknownName, target_datum,
-                                    EllipsoidalCS::createLongitudeLatitude(
-                                        UnitOfMeasure::DEGREE)));
-
-    auto sourceCRS_as_CRS = util::nn_static_pointer_cast<crs::CRS>(sourceCRS);
-
-    std::vector<metadata::PositionalAccuracyNNPtr> emptyAccuracies;
-
-    auto transf = ([&]() {
-        if (abridged) {
-            return Transformation::createAbridgedMolodensky(
-                mapWithUnknownName, sourceCRS_as_CRS, targetCRS, dx, dy, dz, da,
-                df, emptyAccuracies);
-        } else {
-            return Transformation::createMolodensky(
-                mapWithUnknownName, sourceCRS_as_CRS, targetCRS, dx, dy, dz, da,
-                df, emptyAccuracies);
-        }
-    })();
-
-    if (step.inverted) {
-        return util::nn_static_pointer_cast<CoordinateOperation>(
-            transf->inverse());
-    } else {
-        return util::nn_static_pointer_cast<CoordinateOperation>(transf);
-    }
-}
-
 //! @endcond
 
 // ---------------------------------------------------------------------------
@@ -7655,10 +7487,6 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
     int iSecondUnitConvert = -1;
     int iFirstAxisSwap = -1;
     int iSecondAxisSwap = -1;
-    int iHelmert = -1;
-    int iFirstCart = -1;
-    int iSecondCart = -1;
-    int iMolodensky = -1;
     bool unexpectedStructure = d->steps_.empty();
     for (int i = 0; i < static_cast<int>(d->steps_.size()); i++) {
         const auto &stepName = d->steps_[i].name;
@@ -7689,27 +7517,6 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
                 unexpectedStructure = true;
                 break;
             }
-        } else if (stepName == "helmert") {
-            if (iHelmert >= 0) {
-                unexpectedStructure = true;
-                break;
-            }
-            iHelmert = i;
-        } else if (stepName == "cart") {
-            if (iFirstCart < 0) {
-                iFirstCart = i;
-            } else if (iSecondCart < 0) {
-                iSecondCart = i;
-            } else {
-                unexpectedStructure = true;
-                break;
-            }
-        } else if (stepName == "molodensky") {
-            if (iMolodensky >= 0) {
-                unexpectedStructure = true;
-                break;
-            }
-            iMolodensky = i;
         } else if (isProjectedStep(stepName)) {
             if (iProjStep >= 0) {
                 unexpectedStructure = true;
@@ -7722,19 +7529,13 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
         }
     }
 
-    if (iHelmert < 0 && iMolodensky < 0 && !d->steps_.empty()) {
+    if (!d->steps_.empty()) {
         // CRS candidate
         if ((d->steps_.size() == 1 &&
              d->getParamValue(d->steps_[0], "type") != "crs") ||
             (d->steps_.size() > 1 && d->getGlobalParamValue("type") != "crs")) {
             unexpectedStructure = true;
         }
-    } else if (iHelmert >= 0 &&
-               (d->hasParamValue(d->steps_[iHelmert], "theta") ||
-                d->hasParamValue(d->steps_[iHelmert], "exact") ||
-                d->hasParamValue(d->steps_[iHelmert], "transpose") ||
-                d->hasParamValue(d->steps_[iHelmert], "towgs84"))) {
-        unexpectedStructure = true;
     }
 
     struct Logger {
@@ -7877,8 +7678,7 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
 
     if (!unexpectedStructure) {
         if (iFirstGeogStep == 0 && !d->steps_[iFirstGeogStep].inverted &&
-            iSecondGeogStep < 0 && iProjStep < 0 && iHelmert < 0 &&
-            iFirstCart < 0 && iMolodensky < 0 &&
+            iSecondGeogStep < 0 && iProjStep < 0 &&
             (iFirstUnitConvert < 0 || iSecondUnitConvert < 0) &&
             (iFirstAxisSwap < 0 || iSecondAxisSwap < 0)) {
             const bool ignoreVUnits = false;
@@ -7895,8 +7695,7 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
         }
         if (iProjStep >= 0 && !d->steps_[iProjStep].inverted &&
             (iFirstGeogStep < 0 || iFirstGeogStep + 1 == iProjStep) &&
-            iMolodensky < 0 && iSecondGeogStep < 0 && iFirstCart < 0 &&
-            iHelmert < 0) {
+            iSecondGeogStep < 0) {
             if (iFirstGeogStep < 0)
                 iFirstGeogStep = iProjStep;
             const bool ignoreVUnits = true;
@@ -7922,47 +7721,6 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
                 }
             }
         }
-
-        if (iProjStep < 0 && iHelmert > 0 && iMolodensky < 0 &&
-            (iFirstGeogStep < 0 || iFirstGeogStep == iFirstCart - 1 ||
-             (iFirstGeogStep == iSecondCart + 1 && iSecondGeogStep < 0)) &&
-            iFirstCart == iHelmert - 1 && iSecondCart == iHelmert + 1 &&
-            (iSecondGeogStep < 0 || iSecondGeogStep == iSecondCart + 1) &&
-            !d->steps_[iFirstCart].inverted &&
-            d->steps_[iSecondCart].inverted && iFirstAxisSwap < iHelmert &&
-            iFirstUnitConvert < iHelmert &&
-            (iSecondAxisSwap < 0 || iSecondAxisSwap > iHelmert) &&
-            (iSecondUnitConvert < 0 || iSecondUnitConvert > iHelmert)) {
-            return d->buildHelmertTransformation(
-                iHelmert, iFirstAxisSwap, iFirstUnitConvert,
-                iFirstGeogStep >= 0 && iFirstGeogStep == iFirstCart - 1
-                    ? iFirstGeogStep
-                    : iFirstCart,
-                iFirstGeogStep == iSecondCart + 1
-                    ? iFirstGeogStep
-                    : iSecondGeogStep == iSecondCart + 1 ? iSecondGeogStep
-                                                         : iSecondCart,
-                iSecondAxisSwap, iSecondUnitConvert);
-        }
-
-        if (iProjStep < 0 && iHelmert < 0 && iMolodensky > 0 &&
-            (iFirstGeogStep < 0 || iFirstGeogStep == iMolodensky - 1 ||
-             (iFirstGeogStep == iMolodensky + 1 && iSecondGeogStep < 0)) &&
-            (iSecondGeogStep < 0 || iSecondGeogStep == iMolodensky + 1) &&
-            iFirstAxisSwap < iMolodensky && iFirstUnitConvert < iMolodensky &&
-            (iSecondAxisSwap < 0 || iSecondAxisSwap > iMolodensky) &&
-            (iSecondUnitConvert < 0 || iSecondUnitConvert > iMolodensky)) {
-            return d->buildMolodenskyTransformation(
-                iMolodensky, iFirstAxisSwap, iFirstUnitConvert,
-                iFirstGeogStep >= 0 && iFirstGeogStep == iMolodensky - 1
-                    ? iFirstGeogStep
-                    : iMolodensky,
-                iFirstGeogStep == iMolodensky + 1
-                    ? iFirstGeogStep
-                    : iSecondGeogStep == iMolodensky + 1 ? iSecondGeogStep
-                                                         : iMolodensky,
-                iSecondAxisSwap, iSecondUnitConvert);
-        }
     }
 
     auto props = PropertyMap();
diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist
index 72e95634..257e9400 100644
--- a/test/cli/tv_out.dist
+++ b/test/cli/tv_out.dist
@@ -7,8 +7,8 @@ Test NAD27 to raw ellipse
 79d00'00.000"W 35d00'00.000"N 0.0	79dW	35dN 0.000
 ##############################################################
 Between two 3parameter approximations on same ellipsoid
-0d00'00.000"W 0d00'00.000"N 0.0	0dE	0dN 4.000
-79d00'00.000"W 45d00'00.000"N 0.0	78d59'59.821"W	44d59'59.983"N 0.540
+0d00'00.000"W 0d00'00.000"N 0.0	0dE	0dN 0.000
+79d00'00.000"W 45d00'00.000"N 0.0	78d59'59.821"W	44d59'59.983"N 0.000
 ##############################################################
 3param to raw ellipsoid on same ellipsoid
 0d00'00.000"W 0d00'00.000"N 0.0	0dE	0dN 0.000
@@ -352,7 +352,7 @@ Test inverse handling
 10 20	-1384841.19	7581707.88 0.00
 ##############################################################
 Test MGI datum gives expected results (#207)
-16.33 48.20	595710.3731286	5357598.4645652 -44.4951085
+16.33 48.20	595710.3731286	5357598.4645652 0.0000000
 ##############################################################
 Test omerc sensitivity with locations 90d from origin(#114)
 56.958381652832 72.8798	-9985.16336453	-227.67701050 0.00000000
@@ -387,10 +387,10 @@ Test bug 244 (2)
 987122.418330284 -14429896.539530909	-140.10000000000	-87.00000000000 0.00000000000
 ##############################################################
 Test bug 245 (use +datum=carthage)
-10 34	592302.9819461	3762148.7340609 -30.3110170
+10 34	592302.9819461	3762148.7340609 0.0000000
 ##############################################################
 Test bug 245 (use expansion of +datum=carthage)
-10 34	592302.9819461	3762148.7340609 -30.3110170
+10 34	592302.9819461	3762148.7340609 0.0000000
 ##############################################################
 Test SCH forward projection
 0.0 0.0	-1977112.0305592	5551475.1418378 6595.7256583
diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp
index b77fecf3..a3cd8eab 100644
--- a/test/unit/gie_self_tests.cpp
+++ b/test/unit/gie_self_tests.cpp
@@ -682,17 +682,17 @@ TEST(gie, proj_create_crs_to_crs_PULKOVO42_ETRS89) {
     EXPECT_NEAR(c.xy.x, 44.999701238, 1e-9);
     EXPECT_NEAR(c.xy.y, 24.998474948, 1e-9);
     EXPECT_EQ(std::string(proj_pj_info(P).definition),
-              "proj=pipeline step proj=axisswap order=2,1 step "
-              "proj=unitconvert xy_in=deg xy_out=rad step proj=cart "
+              "proj=pipeline step proj=push v_3 step proj=axisswap order=2,1 "
+              "step proj=unitconvert xy_in=deg xy_out=rad step proj=cart "
               "ellps=krass step proj=helmert x=2.3287 y=-147.0425 z=-92.0802 "
               "rx=0.3092483 ry=-0.32482185 rz=-0.49729934 s=5.68906266 "
               "convention=coordinate_frame step inv proj=cart ellps=GRS80 step "
               "proj=unitconvert xy_in=rad xy_out=deg step proj=axisswap "
-              "order=2,1");
+              "order=2,1 step proj=pop v_3");
 
     c = proj_trans(P, PJ_INV, c);
-    EXPECT_NEAR(c.xy.x, 45, 1e-9);
-    EXPECT_NEAR(c.xy.y, 25, 1e-9);
+    EXPECT_NEAR(c.xy.x, 45, 1e-8);
+    EXPECT_NEAR(c.xy.y, 25, 1e-8);
 
     c.xyz.x = 45; // Lat
     c.xyz.y = 25; // Long
@@ -711,12 +711,12 @@ TEST(gie, proj_create_crs_to_crs_PULKOVO42_ETRS89) {
     EXPECT_NEAR(c.xy.x, 51.999714150, 1e-9);
     EXPECT_NEAR(c.xy.y, 19.998187811, 1e-9);
     EXPECT_EQ(std::string(proj_pj_info(P).definition),
-              "proj=pipeline step proj=axisswap order=2,1 step "
-              "proj=unitconvert xy_in=deg xy_out=rad step proj=cart "
+              "proj=pipeline step proj=push v_3 step proj=axisswap order=2,1 "
+              "step proj=unitconvert xy_in=deg xy_out=rad step proj=cart "
               "ellps=krass step proj=helmert x=33.4 y=-146.6 z=-76.3 rx=-0.359 "
               "ry=-0.053 rz=0.844 s=-0.84 convention=position_vector step inv "
               "proj=cart ellps=GRS80 step proj=unitconvert xy_in=rad "
-              "xy_out=deg step proj=axisswap order=2,1");
+              "xy_out=deg step proj=axisswap order=2,1 step proj=pop v_3");
 
     proj_destroy(P);
 }
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index d535e412..273a04e6 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -838,13 +838,14 @@ TEST_F(CApi, proj_create_from_database) {
         ASSERT_NE(info.definition, nullptr);
         EXPECT_EQ(
             info.definition,
-            std::string("proj=pipeline step proj=axisswap order=2,1 step "
-                        "proj=unitconvert xy_in=deg xy_out=rad step proj=cart "
-                        "ellps=bessel step proj=helmert x=601.705 y=84.263 "
-                        "z=485.227 rx=-4.7354 ry=-1.3145 rz=-5.393 s=-2.3887 "
-                        "convention=coordinate_frame step inv proj=cart "
-                        "ellps=GRS80 step proj=unitconvert xy_in=rad "
-                        "xy_out=deg step proj=axisswap order=2,1"));
+            std::string("proj=pipeline step proj=push v_3 step proj=axisswap "
+                        "order=2,1 step proj=unitconvert xy_in=deg xy_out=rad "
+                        "step proj=cart ellps=bessel step proj=helmert "
+                        "x=601.705 y=84.263 z=485.227 rx=-4.7354 ry=-1.3145 "
+                        "rz=-5.393 s=-2.3887 convention=coordinate_frame step "
+                        "inv proj=cart ellps=GRS80 step proj=unitconvert "
+                        "xy_in=rad xy_out=deg step proj=axisswap order=2,1 "
+                        "step proj=pop v_3"));
         EXPECT_EQ(info.accuracy, 1);
     }
 }
diff --git a/test/unit/test_factory.cpp b/test/unit/test_factory.cpp
index ce019079..64512e56 100644
--- a/test/unit/test_factory.cpp
+++ b/test/unit/test_factory.cpp
@@ -658,13 +658,13 @@ TEST(factory, AuthorityFactory_createCoordinateOperation_helmert_3) {
                  NoSuchAuthorityCodeException);
     auto op = factory->createCoordinateOperation("1113", false);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +inv "
-              "+proj=longlat +a=6378249.145 +rf=293.4663077 +step +proj=cart "
-              "+a=6378249.145 +rf=293.4663077 +step +proj=helmert +x=-143 "
-              "+y=-90 +z=-294 +step +inv +proj=cart +ellps=WGS84 +step "
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
+              "+inv +proj=longlat +a=6378249.145 +rf=293.4663077 +step "
+              "+proj=cart +a=6378249.145 +rf=293.4663077 +step +proj=helmert "
+              "+x=-143 +y=-90 +z=-294 +step +inv +proj=cart +ellps=WGS84 +step "
               "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
-              "+order=2,1");
+              "+order=2,1 +step +proj=pop +v_3");
 }
 
 // ---------------------------------------------------------------------------
@@ -673,13 +673,13 @@ TEST(factory, AuthorityFactory_createCoordinateOperation_helmert_7_CF) {
     auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
     auto op = factory->createCoordinateOperation("7676", false);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-              "+ellps=bessel +step +proj=helmert +x=577.88891 +y=165.22205 "
-              "+z=391.18289 +rx=-4.9145 +ry=0.94729 +rz=13.05098 +s=7.78664 "
-              "+convention=coordinate_frame +step +inv +proj=cart +ellps=WGS84 "
-              "+step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
-              "+proj=axisswap +order=2,1");
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
+              "+proj=cart +ellps=bessel +step +proj=helmert +x=577.88891 "
+              "+y=165.22205 +z=391.18289 +rx=-4.9145 +ry=0.94729 +rz=13.05098 "
+              "+s=7.78664 +convention=coordinate_frame +step +inv +proj=cart "
+              "+ellps=WGS84 +step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1 +step +proj=pop +v_3");
 }
 
 // ---------------------------------------------------------------------------
@@ -850,14 +850,15 @@ TEST(factory,
     EXPECT_TRUE(so->validateParameters().empty());
 
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-              "+ellps=bessel +step +proj=molobadekas +x=593.032 +y=26 "
-              "+z=478.741 +rx=0.409394387439237 +ry=-0.359705195614311 "
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
+              "+proj=cart +ellps=bessel +step +proj=molobadekas +x=593.032 "
+              "+y=26 +z=478.741 +rx=0.409394387439237 +ry=-0.359705195614311 "
               "+rz=1.86849100035057 +s=4.0772 +px=3903453.148 +py=368135.313 "
               "+pz=5012970.306 +convention=coordinate_frame +step +inv "
               "+proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad "
-              "+xy_out=deg +step +proj=axisswap +order=2,1");
+              "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop "
+              "+v_3");
 }
 
 // ---------------------------------------------------------------------------
@@ -2065,11 +2066,12 @@ TEST_F(FactoryWithTmpDatabase, AuthorityFactory_wkt_based_transformation) {
     ASSERT_EQ(res.size(), 1);
     EXPECT_EQ(res[0]->nameStr(), "My WKT string based op");
     EXPECT_EQ(res[0]->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-              "+ellps=WGS84 +step +proj=helmert +x=1 +y=2 +z=3 +step +inv "
-              "+proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
-              "+xy_out=deg +step +proj=axisswap +order=2,1");
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
+              "+proj=cart +ellps=WGS84 +step +proj=helmert +x=1 +y=2 +z=3 "
+              "+step +inv +proj=cart +ellps=WGS84 +step +proj=unitconvert "
+              "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 +step "
+              "+proj=pop +v_3");
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 3cb50352..d033508f 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -7659,15 +7659,14 @@ TEST(io, projparse_etmerc) {
 
     auto wkt1 = crs->exportToWKT(
         WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
-    EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\"") == std::string::npos)
-        << wkt1;
+    EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\"") == std::string::npos) << wkt1;
 }
 
 // ---------------------------------------------------------------------------
 
 TEST(io, projparse_tmerc_approx) {
-    auto obj =
-        PROJStringParser().createFromPROJString("+proj=tmerc +approx +type=crs");
+    auto obj = PROJStringParser().createFromPROJString(
+        "+proj=tmerc +approx +type=crs");
     auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
     ASSERT_TRUE(crs != nullptr);
     auto wkt2 = crs->exportToWKT(
@@ -8358,270 +8357,6 @@ TEST(io, projparse_ob_tran_longlat) {
 
 // ---------------------------------------------------------------------------
 
-TEST(io, projparse_helmert_translation) {
-    std::string projString("+proj=pipeline +step +proj=unitconvert +xy_in=deg "
-                           "+xy_out=rad +step +proj=cart +ellps=GRS80 +step "
-                           "+proj=helmert +x=1 +y=2 +z=3 +step +inv +proj=cart "
-                           "+ellps=GRS80 +step +proj=unitconvert +xy_in=rad "
-                           "+xy_out=deg");
-    auto obj = PROJStringParser().createFromPROJString(projString);
-    auto transf = nn_dynamic_pointer_cast<Transformation>(obj);
-    ASSERT_TRUE(transf != nullptr);
-    EXPECT_EQ(
-        transf->exportToPROJString(
-            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5)
-                .get()),
-        projString);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(io, projparse_helmert_translation_inv) {
-    std::string projString("+proj=pipeline +step +proj=unitconvert +xy_in=deg "
-                           "+xy_out=rad +step +proj=cart +ellps=GRS80 +step "
-                           "+inv +proj=helmert +x=1 +y=2 +z=3 +step +inv "
-                           "+proj=cart +ellps=GRS80 +step +proj=unitconvert "
-                           "+xy_in=rad +xy_out=deg");
-    auto obj = PROJStringParser().createFromPROJString(projString);
-    auto transf = nn_dynamic_pointer_cast<Transformation>(obj);
-    ASSERT_TRUE(transf != nullptr);
-    EXPECT_EQ(
-        transf->exportToPROJString(
-            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5)
-                .get()),
-        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-        "+proj=cart +ellps=GRS80 +step +proj=helmert +x=-1 +y=-2 +z=-3 +step "
-        "+inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad "
-        "+xy_out=deg");
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(io, projparse_helmert_position_vector) {
-    std::string projString(
-        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-        "+proj=cart +ellps=GRS80 +step +proj=helmert +x=1 +y=2 +z=3 +rx=4 "
-        "+ry=5 +rz=6 +s=7 +convention=position_vector +step +inv +proj=cart "
-        "+ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg");
-    auto obj = PROJStringParser().createFromPROJString(projString);
-    auto transf = nn_dynamic_pointer_cast<Transformation>(obj);
-    ASSERT_TRUE(transf != nullptr);
-    EXPECT_EQ(
-        transf->exportToPROJString(
-            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5)
-                .get()),
-        projString);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(io, projparse_helmert_position_vector_inv) {
-    std::string projString(
-        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-        "+proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=1 +y=2 +z=3 +rx=4 "
-        "+ry=5 +rz=6 +s=7 +convention=position_vector +step +inv +proj=cart "
-        "+ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg");
-    auto obj = PROJStringParser().createFromPROJString(projString);
-    auto transf = nn_dynamic_pointer_cast<CoordinateOperation>(obj);
-    ASSERT_TRUE(transf != nullptr);
-    EXPECT_EQ(
-        transf->exportToPROJString(
-            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5)
-                .get()),
-        projString);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(io, projparse_helmert_time_dependent_position_vector) {
-    std::string projString("+proj=pipeline +step +proj=unitconvert +xy_in=deg "
-                           "+xy_out=rad +step +proj=cart +ellps=GRS80 +step "
-                           "+proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 "
-                           "+s=7 +dx=0.1 +dy=0.2 +dz=0.3 +drx=0.4 +dry=0.5 "
-                           "+drz=0.6 +ds=0.7 +t_epoch=2018.5 "
-                           "+convention=position_vector +step +inv +proj=cart "
-                           "+ellps=GRS80 +step +proj=unitconvert +xy_in=rad "
-                           "+xy_out=deg");
-    auto obj = PROJStringParser().createFromPROJString(projString);
-    auto transf = nn_dynamic_pointer_cast<Transformation>(obj);
-    ASSERT_TRUE(transf != nullptr);
-    EXPECT_EQ(
-        transf->exportToPROJString(
-            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5)
-                .get()),
-        projString);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(io, projparse_helmert_coordinate_frame) {
-    std::string projString(
-        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-        "+proj=cart +ellps=GRS80 +step +proj=helmert +x=1 +y=2 +z=3 +rx=4 "
-        "+ry=5 +rz=6 +s=7 +convention=coordinate_frame +step +inv +proj=cart "
-        "+ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg");
-    auto obj = PROJStringParser().createFromPROJString(projString);
-    auto transf = nn_dynamic_pointer_cast<Transformation>(obj);
-    ASSERT_TRUE(transf != nullptr);
-    EXPECT_EQ(
-        transf->exportToPROJString(
-            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5)
-                .get()),
-        projString);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(io, projparse_helmert_time_dependent_coordinate_frame) {
-    std::string projString("+proj=pipeline +step +proj=unitconvert +xy_in=deg "
-                           "+xy_out=rad +step +proj=cart +ellps=GRS80 +step "
-                           "+proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 "
-                           "+s=7 +dx=0.1 +dy=0.2 +dz=0.3 +drx=0.4 +dry=0.5 "
-                           "+drz=0.6 +ds=0.7 +t_epoch=2018.5 "
-                           "+convention=coordinate_frame +step +inv +proj=cart "
-                           "+ellps=GRS80 +step +proj=unitconvert +xy_in=rad "
-                           "+xy_out=deg");
-    auto obj = PROJStringParser().createFromPROJString(projString);
-    auto transf = nn_dynamic_pointer_cast<Transformation>(obj);
-    ASSERT_TRUE(transf != nullptr);
-    EXPECT_EQ(
-        transf->exportToPROJString(
-            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5)
-                .get()),
-        projString);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(io, projparse_helmert_complex_pipeline) {
-    std::string projString(
-        "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-        "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-        "+ellps=WGS84 +step +proj=helmert +x=-1 +y=-2 +z=-3 +rx=-4 "
-        "+ry=-5 +rz=-6 +s=-7 +convention=position_vector +step +inv "
-        "+proj=cart +ellps=clrk80ign +step "
-        "+proj=longlat +ellps=clrk80ign +pm=paris +step "
-        "+proj=unitconvert +xy_in=rad +xy_out=grad +step +proj=axisswap "
-        "+order=2,1");
-    auto obj = PROJStringParser().createFromPROJString(projString);
-    auto transf = nn_dynamic_pointer_cast<Transformation>(obj);
-    ASSERT_TRUE(transf != nullptr);
-    EXPECT_EQ(
-        transf->exportToPROJString(
-            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5)
-                .get()),
-        projString);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(io, projparse_helmert_errors) {
-    // Missing convention
-    EXPECT_THROW(PROJStringParser().createFromPROJString("+proj=helmert +rx=4"),
-                 ParsingException);
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(io, projparse_molodensky) {
-    std::string projString(
-        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
-        "+proj=longlat +ellps=WGS84 +step +proj=molodensky +ellps=WGS84 "
-        "+dx=84.87 +dy=96.49 "
-        "+dz=116.95 +da=251 +df=1.41927e-05 +step +proj=longlat "
-        "+ellps=GRS80 +step +proj=unitconvert "
-        "+xy_in=rad +xy_out=deg");
-    auto obj = PROJStringParser().createFromPROJString(projString);
-    auto transf = nn_dynamic_pointer_cast<Transformation>(obj);
-    ASSERT_TRUE(transf != nullptr);
-
-    WKTFormatterNNPtr f(WKTFormatter::create());
-    f->simulCurNodeHasId();
-    f->setMultiLine(false);
-    transf->exportToWKT(f.get());
-    auto wkt = f->toString();
-    EXPECT_EQ(
-        wkt,
-        "COORDINATEOPERATION[\"unknown\",SOURCECRS[GEODCRS[\"unknown\",DATUM["
-        "\"Unknown based on WGS84 ellipsoid\",ELLIPSOID[\"WGS "
-        "84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]]],PRIMEM["
-        "\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],CS["
-        "ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\","
-        "0.0174532925199433]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT["
-        "\"degree\",0.0174532925199433]]]],TARGETCRS[GEODCRS[\"unknown\",DATUM["
-        "\"Unknown based on GRS80 ellipsoid\",ELLIPSOID[\"GRS "
-        "1980\",6378137,298.257222101,LENGTHUNIT[\"metre\",1]]],PRIMEM["
-        "\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],CS["
-        "ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\","
-        "0.0174532925199433]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT["
-        "\"degree\",0.0174532925199433]]]],METHOD[\"Molodensky\",ID[\"EPSG\","
-        "9604]],PARAMETER[\"X-axis "
-        "translation\",84.87,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8605]],"
-        "PARAMETER[\"Y-axis "
-        "translation\",96.49,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8606]],"
-        "PARAMETER[\"Z-axis "
-        "translation\",116.95,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8607]],"
-        "PARAMETER[\"Semi-major axis length "
-        "difference\",251,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8654]],PARAMETER["
-        "\"Flattening difference\",1.41927E-05,ID[\"EPSG\",8655]]]");
-
-    EXPECT_EQ(
-        transf->exportToPROJString(
-            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5)
-                .get()),
-        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-        "+proj=molodensky +ellps=WGS84 +dx=84.87 +dy=96.49 +dz=116.95 +da=251 "
-        "+df=1.41927e-05 +step +proj=unitconvert +xy_in=rad +xy_out=deg");
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(io, projparse_molodensky_inv) {
-    std::string projString(
-        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
-        "+proj=longlat +ellps=WGS84 +step +inv +proj=molodensky "
-        "+ellps=WGS84 +dx=84.87 +dy=96.49 "
-        "+dz=116.95 +da=251 +df=1.41927e-05 +step +proj=longlat "
-        "+ellps=GRS80 +step +proj=unitconvert "
-        "+xy_in=rad +xy_out=deg");
-    auto obj = PROJStringParser().createFromPROJString(projString);
-    auto transf = nn_dynamic_pointer_cast<CoordinateOperation>(obj);
-    ASSERT_TRUE(transf != nullptr);
-    EXPECT_EQ(
-        transf->exportToPROJString(
-            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5)
-                .get()),
-        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-        "+proj=molodensky +ellps=GRS80 +dx=-84.87 +dy=-96.49 +dz=-116.95 "
-        "+da=-251 +df=-1.41927e-05 +step +proj=unitconvert +xy_in=rad "
-        "+xy_out=deg");
-}
-
-// ---------------------------------------------------------------------------
-
-TEST(io, projparse_molodensky_abridged) {
-    std::string projString(
-        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
-        "+proj=longlat +ellps=WGS84 +step +proj=molodensky +ellps=WGS84 "
-        "+dx=84.87 +dy=96.49 "
-        "+dz=116.95 +da=251 +df=1.41927e-05 +abridged +step +proj=longlat "
-        "+ellps=GRS80 +step +proj=unitconvert "
-        "+xy_in=rad +xy_out=deg");
-    auto obj = PROJStringParser().createFromPROJString(projString);
-    auto transf = nn_dynamic_pointer_cast<Transformation>(obj);
-    ASSERT_TRUE(transf != nullptr);
-    EXPECT_EQ(
-        transf->exportToPROJString(
-            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_5)
-                .get()),
-        "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
-        "+proj=molodensky +ellps=WGS84 +dx=84.87 +dy=96.49 +dz=116.95 +da=251 "
-        "+df=1.41927e-05 +abridged +step +proj=unitconvert +xy_in=rad "
-        "+xy_out=deg");
-}
-
-// ---------------------------------------------------------------------------
-
 TEST(io, projparse_longlat_title) {
     std::string projString("+title=Ile d'Amsterdam 1963 +proj=longlat "
                            "+towgs84=109.7530,-528.1330,-362.2440 "
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 586226e2..e855a818 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -561,11 +561,12 @@ TEST(operation, transformation_createGeocentricTranslations) {
     EXPECT_EQ(inv_transf_as_transf->getTOWGS84Parameters(), expected_inv);
 
     EXPECT_EQ(transf->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-              "+ellps=GRS80 +step +proj=helmert +x=1 +y=2 +z=3 +step +inv "
-              "+proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
-              "+xy_out=deg +step +proj=axisswap +order=2,1");
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
+              "+proj=cart +ellps=GRS80 +step +proj=helmert +x=1 +y=2 +z=3 "
+              "+step +inv +proj=cart +ellps=WGS84 +step +proj=unitconvert "
+              "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 +step "
+              "+proj=pop +v_3");
 }
 
 // ---------------------------------------------------------------------------
@@ -668,13 +669,13 @@ TEST(operation, transformation_createPositionVector) {
     EXPECT_EQ(transf->getTOWGS84Parameters(), expected);
 
     EXPECT_EQ(transf->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-              "+ellps=GRS80 +step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 "
-              "+rz=6 +s=7 +convention=position_vector +step +inv +proj=cart "
-              "+ellps=WGS84 +step "
-              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
-              "+order=2,1");
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
+              "+proj=cart +ellps=GRS80 +step +proj=helmert +x=1 +y=2 +z=3 "
+              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
+              "+proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
+              "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop "
+              "+v_3");
 
     auto inv_transf = transf->inverse();
     ASSERT_EQ(inv_transf->coordinateOperationAccuracies().size(), 1);
@@ -691,12 +692,12 @@ TEST(operation, transformation_createPositionVector) {
 #else
     EXPECT_EQ(
         inv_transf->exportToPROJString(PROJStringFormatter::create().get()),
-        "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-        "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-        "+ellps=WGS84 +step +inv +proj=helmert +x=1 +y=2 +z=3 +rx=4 "
-        "+ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
-        "+proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad "
-        "+xy_out=deg +step +proj=axisswap +order=2,1");
+        "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap +order=2,1 "
+        "+step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
+        "+ellps=WGS84 +step +inv +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 "
+        "+rz=6 +s=7 +convention=position_vector +step +inv +proj=cart "
+        "+ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
+        "+proj=axisswap +order=2,1 +step +proj=pop +v_3");
 
     // In WKT, use approximate formula
     auto wkt = inv_transf->exportToWKT(WKTFormatter::create().get());
@@ -737,13 +738,13 @@ TEST(operation, transformation_createCoordinateFrameRotation) {
     EXPECT_EQ(params, expected);
 
     EXPECT_EQ(transf->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-              "+ellps=GRS80 +step +proj=helmert +x=1 +y=2 +z=3 +rx=-4 +ry=-5 "
-              "+rz=-6 +s=7 +convention=coordinate_frame +step +inv +proj=cart "
-              "+ellps=WGS84 +step "
-              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
-              "+order=2,1");
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
+              "+proj=cart +ellps=GRS80 +step +proj=helmert +x=1 +y=2 +z=3 "
+              "+rx=-4 +ry=-5 +rz=-6 +s=7 +convention=coordinate_frame +step "
+              "+inv +proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
+              "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop "
+              "+v_3");
 
     auto inv_transf = transf->inverse();
     ASSERT_EQ(inv_transf->coordinateOperationAccuracies().size(), 0);
@@ -760,12 +761,12 @@ TEST(operation, transformation_createCoordinateFrameRotation) {
 #else
     EXPECT_EQ(
         inv_transf->exportToPROJString(PROJStringFormatter::create().get()),
-        "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-        "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-        "+ellps=WGS84 +step +inv +proj=helmert +x=1 +y=2 +z=3 +rx=-4 "
-        "+ry=-5 +rz=-6 +s=7 +convention=coordinate_frame +step +inv "
-        "+proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad "
-        "+xy_out=deg +step +proj=axisswap +order=2,1");
+        "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap +order=2,1 "
+        "+step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
+        "+ellps=WGS84 +step +inv +proj=helmert +x=1 +y=2 +z=3 +rx=-4 +ry=-5 "
+        "+rz=-6 +s=7 +convention=coordinate_frame +step +inv +proj=cart "
+        "+ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
+        "+proj=axisswap +order=2,1 +step +proj=pop +v_3");
 
     // In WKT, use approximate formula
     auto wkt = inv_transf->exportToWKT(WKTFormatter::create().get());
@@ -4197,13 +4198,13 @@ TEST(operation, geogCRS_to_geogCRS_context_default) {
 
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
-            "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-            "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-            "+ellps=krass +step +proj=helmert +x=2.3287 +y=-147.0425 "
-            "+z=-92.0802 +rx=0.3092483 +ry=-0.32482185 +rz=-0.49729934 "
-            "+s=5.68906266 +convention=coordinate_frame +step +inv "
-            "+proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad "
-            "+xy_out=deg +step +proj=axisswap +order=2,1");
+            "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+            "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
+            "+proj=cart +ellps=krass +step +proj=helmert +x=2.3287 "
+            "+y=-147.0425 +z=-92.0802 +rx=0.3092483 +ry=-0.32482185 "
+            "+rz=-0.49729934 +s=5.68906266 +convention=coordinate_frame +step "
+            "+inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad "
+            "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop +v_3");
     }
 
     // Reverse case
@@ -4218,13 +4219,13 @@ TEST(operation, geogCRS_to_geogCRS_context_default) {
 
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
-            "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-            "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-            "+ellps=GRS80 +step +inv +proj=helmert +x=2.3287 "
+            "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+            "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
+            "+proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=2.3287 "
             "+y=-147.0425 +z=-92.0802 +rx=0.3092483 +ry=-0.32482185 "
-            "+rz=-0.49729934 +s=5.68906266 +convention=coordinate_frame "
-            "+step +inv +proj=cart +ellps=krass +step +proj=unitconvert "
-            "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1");
+            "+rz=-0.49729934 +s=5.68906266 +convention=coordinate_frame +step "
+            "+inv +proj=cart +ellps=krass +step +proj=unitconvert +xy_in=rad "
+            "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop +v_3");
     }
 }
 
@@ -4358,11 +4359,12 @@ TEST(operation, geogCRS_to_geogCRS_context_inverse_needed) {
         ASSERT_EQ(list.size(), 3);
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
-            "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-            "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-            "+ellps=clrk80ign +step +proj=helmert +x=-168 +y=-60 +z=320 "
-            "+step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert "
-            "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1");
+            "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+            "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
+            "+proj=cart +ellps=clrk80ign +step +proj=helmert +x=-168 +y=-60 "
+            "+z=320 +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert "
+            "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 +step "
+            "+proj=pop +v_3");
         EXPECT_EQ(list[1]->exportToPROJString(
                       PROJStringFormatter::create(
                           PROJStringFormatter::Convention::PROJ_5,
@@ -5299,15 +5301,15 @@ TEST(operation, boundCRS_of_geogCRS_to_geogCRS) {
     auto op = CoordinateOperationFactory::create()->createOperation(
         boundCRS, GeographicCRS::EPSG_4326);
     ASSERT_TRUE(op != nullptr);
-    EXPECT_EQ(
-        op->exportToPROJString(PROJStringFormatter::create().get()),
-        "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-        "+proj=unitconvert +xy_in=grad +xy_out=rad +step +inv +proj=longlat "
-        "+ellps=clrk80ign +pm=paris +step +proj=cart +ellps=clrk80ign "
-        "+step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
-        "+convention=position_vector +step "
-        "+inv +proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
-        "+xy_out=deg +step +proj=axisswap +order=2,1");
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=grad +xy_out=rad "
+              "+step +inv +proj=longlat +ellps=clrk80ign +pm=paris +step "
+              "+proj=cart +ellps=clrk80ign +step +proj=helmert +x=1 +y=2 +z=3 "
+              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
+              "+proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
+              "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop "
+              "+v_3");
 }
 
 // ---------------------------------------------------------------------------
@@ -5321,12 +5323,13 @@ TEST(operation, boundCRS_of_geogCRS_to_geogCRS_with_area) {
         boundCRS, authFactory->createCoordinateReferenceSystem("4326"));
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-              "+ellps=clrk66 +step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 "
-              "+rz=6 +s=7 +convention=position_vector +step +inv +proj=cart "
-              "+ellps=WGS84 +step +proj=unitconvert +xy_in=rad +xy_out=deg "
-              "+step +proj=axisswap +order=2,1");
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
+              "+proj=cart +ellps=clrk66 +step +proj=helmert +x=1 +y=2 +z=3 "
+              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
+              "+proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
+              "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop "
+              "+v_3");
 }
 
 // ---------------------------------------------------------------------------
@@ -5362,11 +5365,11 @@ TEST(operation, createOperation_boundCRS_identified_by_datum) {
         NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dest));
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad "
-              "+step +proj=cart +ellps=WGS84 "
-              "+step +proj=helmert +x=263 +y=-6 +z=-431 "
-              "+step +inv +proj=cart +ellps=clrk80ign "
-              "+step +proj=utm +zone=32 +ellps=clrk80ign");
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=unitconvert "
+              "+xy_in=deg +xy_out=rad +step +proj=cart +ellps=WGS84 +step "
+              "+proj=helmert +x=263 +y=-6 +z=-431 +step +inv +proj=cart "
+              "+ellps=clrk80ign +step +proj=pop +v_3 +step +proj=utm +zone=32 "
+              "+ellps=clrk80ign");
 
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), std::string());
@@ -5436,12 +5439,13 @@ TEST(operation, boundCRS_of_projCRS_to_geogCRS) {
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +inv +proj=utm +zone=31 +ellps=clrk80ign "
-              "+pm=paris "
-              "+step +proj=cart +ellps=clrk80ign +step +proj=helmert +x=1 +y=2 "
-              "+z=3 +rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step "
-              "+inv +proj=cart +ellps=WGS84 "
+              "+pm=paris +step +proj=longlat +ellps=clrk80ign +pm=paris +step "
+              "+proj=push +v_3 +step +inv +proj=longlat +ellps=clrk80ign "
+              "+pm=paris +step +proj=cart +ellps=clrk80ign +step +proj=helmert "
+              "+x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
+              "+convention=position_vector +step +inv +proj=cart +ellps=WGS84 "
               "+step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
-              "+proj=axisswap +order=2,1");
+              "+proj=axisswap +order=2,1 +step +proj=pop +v_3");
 }
 
 // ---------------------------------------------------------------------------
@@ -5456,15 +5460,14 @@ TEST(operation, boundCRS_of_geogCRS_to_projCRS) {
     auto op =
         CoordinateOperationFactory::create()->createOperation(boundCRS, utm31);
     ASSERT_TRUE(op != nullptr);
-    EXPECT_EQ(
-        op->exportToPROJString(PROJStringFormatter::create().get()),
-        "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-        "+proj=unitconvert +xy_in=grad +xy_out=rad +step +inv +proj=longlat "
-        "+ellps=clrk80ign +pm=paris +step +proj=cart +ellps=clrk80ign "
-        "+step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
-        "+convention=position_vector +step "
-        "+inv +proj=cart +ellps=WGS84 +step +proj=utm +zone=31 "
-        "+ellps=WGS84");
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=grad +xy_out=rad "
+              "+step +inv +proj=longlat +ellps=clrk80ign +pm=paris +step "
+              "+proj=cart +ellps=clrk80ign +step +proj=helmert +x=1 +y=2 +z=3 "
+              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
+              "+proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=utm "
+              "+zone=31 +ellps=WGS84");
 }
 
 // ---------------------------------------------------------------------------
@@ -5476,13 +5479,14 @@ TEST(operation, geogCRS_to_boundCRS_of_geogCRS) {
         GeographicCRS::EPSG_4326, boundCRS);
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-              "+ellps=WGS84 +step +inv +proj=helmert +x=1 +y=2 +z=3 +rx=4 "
-              "+ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
+              "+proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=1 +y=2 +z=3 "
+              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
               "+proj=cart +ellps=clrk80ign +step +proj=longlat "
               "+ellps=clrk80ign +pm=paris +step +proj=unitconvert +xy_in=rad "
-              "+xy_out=grad +step +proj=axisswap +order=2,1");
+              "+xy_out=grad +step +proj=axisswap +order=2,1 +step +proj=pop "
+              "+v_3");
 }
 
 // ---------------------------------------------------------------------------
@@ -5505,12 +5509,14 @@ TEST(operation, boundCRS_to_boundCRS) {
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +inv +proj=utm +zone=31 +ellps=clrk80ign "
-              "+pm=paris "
-              "+step +proj=cart +ellps=clrk80ign +step +proj=helmert +x=1 +y=2 "
-              "+z=3 +rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step "
-              "+inv +proj=helmert +x=8 +y=9 +z=10 +rx=11 +ry=12 +rz=13 +s=14 "
-              "+convention=position_vector +step +inv +proj=cart +ellps=GRS80 "
-              "+step +proj=utm +zone=32 +ellps=GRS80");
+              "+pm=paris +step +proj=longlat +ellps=clrk80ign +pm=paris +step "
+              "+proj=push +v_3 +step +inv +proj=longlat +ellps=clrk80ign "
+              "+pm=paris +step +proj=cart +ellps=clrk80ign +step +proj=helmert "
+              "+x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
+              "+convention=position_vector +step +inv +proj=helmert +x=8 +y=9 "
+              "+z=10 +rx=11 +ry=12 +rz=13 +s=14 +convention=position_vector "
+              "+step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step "
+              "+proj=utm +zone=32 +ellps=GRS80");
 }
 
 // ---------------------------------------------------------------------------
@@ -5523,13 +5529,13 @@ TEST(operation, boundCRS_to_boundCRS_noop_for_TOWGS84) {
     auto op = CoordinateOperationFactory::create()->createOperation(boundCRS1,
                                                                     boundCRS2);
     ASSERT_TRUE(op != nullptr);
-    EXPECT_EQ(
-        op->exportToPROJString(PROJStringFormatter::create().get()),
-        "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-        "+proj=unitconvert +xy_in=grad +xy_out=rad +step +inv +proj=longlat "
-        "+ellps=clrk80ign +pm=paris +step +proj=cart +ellps=clrk80ign "
-        "+step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert "
-        "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1");
+    EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=grad +xy_out=rad "
+              "+step +inv +proj=longlat +ellps=clrk80ign +pm=paris +step "
+              "+proj=cart +ellps=clrk80ign +step +inv +proj=cart +ellps=GRS80 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
+              "+proj=axisswap +order=2,1 +step +proj=pop +v_3");
 }
 
 // ---------------------------------------------------------------------------
@@ -5927,12 +5933,13 @@ TEST(operation, compoundCRS_with_boundGeogCRS_to_geogCRS) {
         compound, GeographicCRS::EPSG_4979);
     ASSERT_TRUE(op != nullptr);
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart "
-              "+ellps=WGS84 +step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 "
-              "+rz=6 +s=7 +convention=position_vector +step +inv +proj=cart "
-              "+ellps=WGS84 +step +proj=unitconvert +xy_in=rad +xy_out=deg "
-              "+step +proj=axisswap +order=2,1");
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step "
+              "+proj=cart +ellps=WGS84 +step +proj=helmert +x=1 +y=2 +z=3 "
+              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
+              "+proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad "
+              "+xy_out=deg +step +proj=axisswap +order=2,1 +step +proj=pop "
+              "+v_3");
 }
 
 // ---------------------------------------------------------------------------
@@ -5950,15 +5957,15 @@ TEST(operation, compoundCRS_with_boundGeogCRS_and_boundVerticalCRS_to_geogCRS) {
     // Not completely sure the order of horizontal and vertical operations
     // makes sense
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
-              "+proj=unitconvert +xy_in=grad +xy_out=rad +step +inv "
-              "+proj=longlat +ellps=clrk80ign +pm=paris +step +proj=cart "
-              "+ellps=clrk80ign +step +proj=helmert +x=1 +y=2 +z=3 +rx=4 +ry=5 "
-              "+rz=6 +s=7 +convention=position_vector +step +inv +proj=cart "
-              "+ellps=WGS84 +step +proj=vgridshift +grids=egm08_25.gtx "
-              "+multiplier=1 +step "
-              "+proj=unitconvert +xy_in=rad +xy_out=deg +step "
-              "+proj=axisswap +order=2,1");
+              "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
+              "+order=2,1 +step +proj=unitconvert +xy_in=grad +xy_out=rad "
+              "+step +inv +proj=longlat +ellps=clrk80ign +pm=paris +step "
+              "+proj=cart +ellps=clrk80ign +step +proj=helmert +x=1 +y=2 +z=3 "
+              "+rx=4 +ry=5 +rz=6 +s=7 +convention=position_vector +step +inv "
+              "+proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step "
+              "+proj=vgridshift +grids=egm08_25.gtx +multiplier=1 +step "
+              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
+              "+order=2,1");
 
     auto grids = op->gridsNeeded(DatabaseContext::create());
     EXPECT_EQ(grids.size(), 1);
@@ -5989,12 +5996,14 @@ TEST(operation, compoundCRS_with_boundProjCRS_and_boundVerticalCRS_to_geogCRS) {
     // makes sense
     EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +inv +proj=utm +zone=31 +ellps=clrk80ign "
+              "+pm=paris +step +proj=longlat +ellps=clrk80ign +pm=paris +step "
+              "+proj=push +v_3 +step +inv +proj=longlat +ellps=clrk80ign "
               "+pm=paris +step +proj=cart +ellps=clrk80ign +step +proj=helmert "
               "+x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 "
               "+convention=position_vector +step +inv +proj=cart +ellps=WGS84 "
-              "+step +proj=vgridshift +grids=egm08_25.gtx +multiplier=1 +step "
-              "+proj=unitconvert +xy_in=rad +xy_out=deg +step "
-              "+proj=axisswap +order=2,1");
+              "+step +proj=pop +v_3 +step +proj=vgridshift +grids=egm08_25.gtx "
+              "+multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
 
     auto opInverse = CoordinateOperationFactory::create()->createOperation(
         GeographicCRS::EPSG_4979, compound);
@@ -6062,21 +6071,23 @@ TEST(operation, compoundCRS_to_compoundCRS_with_vertical_transform) {
     {
         auto formatter = PROJStringFormatter::create();
         formatter->setUseApproxTMerc(true);
-        EXPECT_EQ(op->exportToPROJString(formatter.get()),
-                  "+proj=pipeline +step +inv +proj=tmerc +approx +lat_0=1 +lon_0=2 "
-                  "+k=3 +x_0=4 +y_0=5 +ellps=WGS84 +step "
-                  "+proj=vgridshift +grids=bla.gtx +multiplier=0.001 +step "
-                  "+proj=utm +approx +zone=32 "
-                  "+ellps=WGS84");
+        EXPECT_EQ(
+            op->exportToPROJString(formatter.get()),
+            "+proj=pipeline +step +inv +proj=tmerc +approx +lat_0=1 +lon_0=2 "
+            "+k=3 +x_0=4 +y_0=5 +ellps=WGS84 +step "
+            "+proj=vgridshift +grids=bla.gtx +multiplier=0.001 +step "
+            "+proj=utm +approx +zone=32 "
+            "+ellps=WGS84");
     }
     {
         auto formatter = PROJStringFormatter::create();
         formatter->setUseApproxTMerc(true);
-        EXPECT_EQ(op->inverse()->exportToPROJString(formatter.get()),
-                  "+proj=pipeline +step +inv +proj=utm +approx +zone=32 +ellps=WGS84 "
-                  "+step +inv +proj=vgridshift +grids=bla.gtx "
-                  "+multiplier=0.001 +step +proj=tmerc +approx +lat_0=1 +lon_0=2 "
-                  "+k=3 +x_0=4 +y_0=5 +ellps=WGS84");
+        EXPECT_EQ(
+            op->inverse()->exportToPROJString(formatter.get()),
+            "+proj=pipeline +step +inv +proj=utm +approx +zone=32 +ellps=WGS84 "
+            "+step +inv +proj=vgridshift +grids=bla.gtx "
+            "+multiplier=0.001 +step +proj=tmerc +approx +lat_0=1 +lon_0=2 "
+            "+k=3 +x_0=4 +y_0=5 +ellps=WGS84");
     }
 
     auto opInverse = CoordinateOperationFactory::create()->createOperation(
@@ -6198,10 +6209,11 @@ TEST(operation, IGNF_LAMB1_TO_EPSG_4326) {
     EXPECT_EQ(list[1]->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline +step +inv +proj=lcc +lat_1=49.5 +lat_0=49.5 "
               "+lon_0=0 +k_0=0.99987734 +x_0=600000 +y_0=200000 "
-              "+ellps=clrk80ign +pm=paris +step +proj=cart +ellps=clrk80ign "
-              "+step +proj=helmert +x=-168 +y=-60 +z=320 +step +inv +proj=cart "
-              "+ellps=WGS84 +step +proj=unitconvert +xy_in=rad +xy_out=deg "
-              "+step +proj=axisswap +order=2,1");
+              "+ellps=clrk80ign +pm=paris +step +proj=push +v_3 +step "
+              "+proj=cart +ellps=clrk80ign +step +proj=helmert +x=-168 +y=-60 "
+              "+z=320 +step +inv +proj=cart +ellps=WGS84 +step "
+              "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
+              "+order=2,1 +step +proj=pop +v_3");
 
     auto list2 = CoordinateOperationFactory::create()->createOperations(
         AuthorityFactory::create(DatabaseContext::create(), "EPSG")