Merge pull request #2259 from rouault/slovakian_geoids

Register Slovakian geoid transformations with needed code changes

8 files changed, 319 insertions, 28 deletions

diff --git a/data/sql/commit.sql b/data/sql/commit.sql
index d090ccb2..55df332f 100644
--- a/data/sql/commit.sql
+++ b/data/sql/commit.sql
@@ -32,9 +32,20 @@ FOR EACH ROW BEGIN
 
     -- test to check that our custom grid transformation overrides are really needed
     SELECT RAISE(ABORT, 'PROJ grid_transformation defined whereas EPSG has one')
-        WHERE EXISTS (SELECT 1 FROM grid_transformation g1, grid_transformation g2 WHERE
-            lower(g1.grid_name) = lower(g2.grid_name) AND
-            g1.auth_name = 'PROJ' AND g2.auth_name = 'EPSG');
+        WHERE EXISTS (SELECT 1 FROM grid_transformation g1
+                      JOIN grid_transformation g2
+                      ON g1.source_crs_auth_name = g2.source_crs_auth_name
+                      AND g1.source_crs_code = g2.source_crs_code
+                      AND g1.target_crs_auth_name = g2.target_crs_auth_name
+                      AND g1.target_crs_code = g2.target_crs_code
+                      WHERE g1.auth_name = 'PROJ' AND g2.auth_name = 'EPSG')
+        OR EXISTS (SELECT 1 FROM grid_transformation g1
+                      JOIN grid_transformation g2
+                      ON g1.source_crs_auth_name = g2.target_crs_auth_name
+                      AND g1.source_crs_code = g2.target_crs_code
+                      AND g1.target_crs_auth_name = g1.source_crs_auth_name
+                      AND g1.target_crs_code = g1.source_crs_code
+                      WHERE g1.auth_name = 'PROJ' AND g2.auth_name = 'EPSG');
 
     SELECT RAISE(ABORT, 'Arg! there is now a EPSG:102100 object. Hack in createFromUserInput() will no longer work')
         WHERE EXISTS(SELECT 1 FROM crs_view WHERE auth_name = 'EPSG' AND code = '102100');
diff --git a/data/sql/customizations.sql b/data/sql/customizations.sql
index efb0a4bc..815c371a 100644
--- a/data/sql/customizations.sql
+++ b/data/sql/customizations.sql
@@ -171,3 +171,42 @@ INSERT INTO "unit_of_measure" VALUES('PROJ','IN','inch','length',0.0254,'in',0);
 INSERT INTO "unit_of_measure" VALUES('PROJ','US_IN','US survey inch','length',0.025400050800101,'us-in',0);
 INSERT INTO "unit_of_measure" VALUES('PROJ','US_YD','US survey yard','length',0.914401828803658,'us-yd',0);
 INSERT INTO "unit_of_measure" VALUES('PROJ','IND_CH','Indian chain','length',20.11669506,'ind-ch',0);
+
+-- Deal with grid_transformation using EPSG:9635 'Geog3D to Geog2D+GravityRelatedHeight (US .gtx)' method
+-- We derive records using the more classic 'Geographic3D to GravityRelatedHeight' method
+-- We could probably do that at runtime too, but more simple and efficient to create records
+
+INSERT INTO "grid_transformation"
+SELECT
+    'PROJ' AS auth_name,
+    gt.auth_name || '_' || gt.code || '_RESTRICTED_TO_VERTCRS' AS code,
+    gcrs.name || ' to ' || vcrs.name || ' (from ' || gt.name || ')' AS name,
+    NULL AS description,
+    gt.scope,
+    'EPSG' AS method_auth_name,
+    '9665' AS method_code,
+    'Geographic3D to GravityRelatedHeight (gtx)' AS method_name,
+    gt.source_crs_auth_name,
+    gt.source_crs_code,
+    c.vertical_crs_auth_name AS target_crs_auth_name,
+    c.vertical_crs_code AS target_crs_code,
+    gt.area_of_use_auth_name,
+    gt.area_of_use_code,
+    gt.accuracy,
+    gt.grid_param_auth_name,
+    gt.grid_param_code,
+    gt.grid_param_name,
+    gt.grid_name,
+    gt.grid2_param_auth_name,
+    gt.grid2_param_code,
+    gt.grid2_param_name,
+    gt.grid2_name,
+    gt.interpolation_crs_auth_name,
+    gt.interpolation_crs_code,
+    gt.operation_version,
+    gt.deprecated
+FROM grid_transformation gt
+JOIN compound_crs c ON gt.target_crs_code = c.code AND gt.target_crs_auth_name = c.auth_name
+JOIN geodetic_crs gcrs ON gt.source_crs_auth_name = gcrs.auth_name AND gt.source_crs_code = gcrs.code
+JOIN vertical_crs vcrs on vcrs.auth_name = c.vertical_crs_auth_name AND vcrs.code = c.vertical_crs_code
+WHERE method_auth_name = 'EPSG' AND method_code = '9635' AND gt.deprecated = 0;
diff --git a/data/sql/grid_alternatives.sql b/data/sql/grid_alternatives.sql
index 3dd264ba..68b992c6 100644
--- a/data/sql/grid_alternatives.sql
+++ b/data/sql/grid_alternatives.sql
@@ -187,6 +187,8 @@ VALUES
 ('auckland-1946-to-nzvd2016-conversion.csv','nz_linz_auckht1946-nzvd2016.tif','auckht1946-nzvd2016.gtx','GTiff','vgridshift',0,NULL,'https://cdn.proj.org/nz_linz_auckht1946-nzvd2016.tif',1,1,NULL),
 -- Slovakia
 ('Slovakia_JTSK03_to_JTSK.LAS','sk_gku_JTSK03_to_JTSK.tif',NULL,'GTiff','hgridshift',0,NULL,'https://cdn.proj.org/sk_gku_JTSK03_to_JTSK.tif',1,1,NULL),
+('Slovakia_ETRS89h_to_Baltic1957.gtx','sk_gku_Slovakia_ETRS89h_to_Baltic1957.tif',NULL,'GTiff','geoid_like',0,NULL,'https://cdn.proj.org/sk_gku_Slovakia_ETRS89h_to_Baltic1957.tif',1,1,NULL),
+('Slovakia_ETRS89h_to_EVRF2007.gtx','sk_gku_Slovakia_ETRS89h_to_EVRF2007.tif',NULL,'GTiff','geoid_like',0,NULL,'https://cdn.proj.org/sk_gku_Slovakia_ETRS89h_to_EVRF2007.tif',1,1,NULL),
 -- Superseded entries
 ('bluff-1955-to-nzvd2016-conversion.csv','nz_linz_blufht1955-nzvd2016.tif','blufht1955-nzvd2016.gtx','GTiff','vgridshift',0,NULL,'https://cdn.proj.org/nz_linz_blufht1955-nzvd2016.tif',1,1,NULL),
 ('dunedin-1958-to-nzvd2016-conversion.csv','nz_linz_duneht1958-nzvd2016.tif','duneht1958-nzvd2016.gtx','GTiff','vgridshift',0,NULL,'https://cdn.proj.org/nz_linz_duneht1958-nzvd2016.tif',1,1,NULL),
diff --git a/include/proj/io.hpp b/include/proj/io.hpp
index d1f79aab..014d9987 100644
--- a/include/proj/io.hpp
+++ b/include/proj/io.hpp
@@ -439,6 +439,10 @@ class PROJ_GCC_DLL PROJStringFormatter {
     PROJ_INTERNAL void popOmitZUnitConversion();
     PROJ_INTERNAL bool omitZUnitConversion() const;
 
+    PROJ_INTERNAL void pushOmitHorizontalConversionInVertTransformation();
+    PROJ_INTERNAL void popOmitHorizontalConversionInVertTransformation();
+    PROJ_INTERNAL bool omitHorizontalConversionInVertTransformation() const;
+
     PROJ_INTERNAL void setLegacyCRSToCRSContext(bool legacyContext);
     PROJ_INTERNAL bool getLegacyCRSToCRSContext() const;
 
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index bdb2ad2e..ecc3df06 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -7822,9 +7822,34 @@ createPropertiesForInverse(const CoordinateOperation *op, bool derivedFrom,
     auto targetCRS = op->targetCRS();
     std::string name;
     if (!forwardName.empty()) {
-        if (starts_with(forwardName, INVERSE_OF) ||
-            forwardName.find(" + ") != std::string::npos) {
-            auto tokens = split(forwardName, " + ");
+        if (dynamic_cast<const Transformation *>(op) == nullptr &&
+            dynamic_cast<const ConcatenatedOperation *>(op) == nullptr &&
+            (starts_with(forwardName, INVERSE_OF) ||
+             forwardName.find(" + ") != std::string::npos)) {
+            std::vector<std::string> tokens;
+            std::string curToken;
+            bool inString = false;
+            for (size_t i = 0; i < forwardName.size(); ++i) {
+                if (inString) {
+                    curToken += forwardName[i];
+                    if (forwardName[i] == '\'') {
+                        inString = false;
+                    }
+                } else if (i + 3 < forwardName.size() &&
+                           memcmp(&forwardName[i], " + ", 3) == 0) {
+                    tokens.push_back(curToken);
+                    curToken.clear();
+                    i += 2;
+                } else if (forwardName[i] == '\'') {
+                    inString = true;
+                    curToken += forwardName[i];
+                } else {
+                    curToken += forwardName[i];
+                }
+            }
+            if (!curToken.empty()) {
+                tokens.push_back(curToken);
+            }
             for (size_t i = tokens.size(); i > 0;) {
                 i--;
                 if (!name.empty()) {
@@ -7841,7 +7866,11 @@ createPropertiesForInverse(const CoordinateOperation *op, bool derivedFrom,
             }
         } else if (!sourceCRS || !targetCRS ||
                    forwardName != buildOpName(opType, sourceCRS, targetCRS)) {
-            name = INVERSE_OF + forwardName;
+            if (forwardName.find(" + ") != std::string::npos) {
+                name = INVERSE_OF + '\'' + forwardName + '\'';
+            } else {
+                name = INVERSE_OF + forwardName;
+            }
         }
     }
     if (name.empty() && sourceCRS && targetCRS) {
@@ -8662,6 +8691,7 @@ isGeographic3DToGravityRelatedHeight(const OperationMethodNNPtr &method,
         "9663", // Geographic3D to GravityRelatedHeight (OSGM-GB)
         "9664", // Geographic3D to GravityRelatedHeight (IGN1997)
         "9665", // Geographic3D to GravityRelatedHeight (US .gtx)
+        "9635", // Geog3D to Geog2D+GravityRelatedHeight (US .gtx)
     };
 
     if (ci_find(methodName, "Geographic3D to GravityRelatedHeight") == 0) {
@@ -9697,6 +9727,21 @@ void Transformation::_exportToPROJString(
 
     const auto &heightFilename = _getHeightToGeographic3DFilename(this, true);
     if (!heightFilename.empty()) {
+        auto targetCRSGeog =
+            extractGeographicCRSIfGeographicCRSOrEquivalent(targetCRS());
+        if (!targetCRSGeog) {
+            throw io::FormattingException(
+                concat("Can apply ", methodName, " only to GeographicCRS"));
+        }
+
+        if (!formatter->omitHorizontalConversionInVertTransformation()) {
+            formatter->startInversion();
+            formatter->pushOmitZUnitConversion();
+            targetCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);
+            formatter->popOmitZUnitConversion();
+            formatter->stopInversion();
+        }
+
         if (isMethodInverseOf) {
             formatter->startInversion();
         }
@@ -9706,6 +9751,13 @@ void Transformation::_exportToPROJString(
         if (isMethodInverseOf) {
             formatter->stopInversion();
         }
+
+        if (!formatter->omitHorizontalConversionInVertTransformation()) {
+            formatter->pushOmitZUnitConversion();
+            targetCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);
+            formatter->popOmitZUnitConversion();
+        }
+
         return;
     }
 
@@ -9716,6 +9768,22 @@ void Transformation::_exportToPROJString(
         if (fileParameter &&
             fileParameter->type() == ParameterValue::Type::FILENAME) {
             auto filename = fileParameter->valueFile();
+
+            auto sourceCRSGeog =
+                extractGeographicCRSIfGeographicCRSOrEquivalent(sourceCRS());
+            if (!sourceCRSGeog) {
+                throw io::FormattingException(
+                    concat("Can apply ", methodName, " only to GeographicCRS"));
+            }
+
+            if (!formatter->omitHorizontalConversionInVertTransformation()) {
+                formatter->startInversion();
+                formatter->pushOmitZUnitConversion();
+                sourceCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);
+                formatter->popOmitZUnitConversion();
+                formatter->stopInversion();
+            }
+
             bool doInversion = isMethodInverseOf;
             // The EPSG Geog3DToHeight is the reverse convention of PROJ !
             doInversion = !doInversion;
@@ -9728,6 +9796,13 @@ void Transformation::_exportToPROJString(
             if (doInversion) {
                 formatter->stopInversion();
             }
+
+            if (!formatter->omitHorizontalConversionInVertTransformation()) {
+                formatter->pushOmitZUnitConversion();
+                sourceCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);
+                formatter->popOmitZUnitConversion();
+            }
+
             return;
         }
     }
@@ -12286,6 +12361,8 @@ createBallparkGeographicOffset(const crs::CRSNNPtr &sourceCRS,
 
 //! @cond Doxygen_Suppress
 
+// ---------------------------------------------------------------------------
+
 struct MyPROJStringExportableGeodToGeod final
     : public io::IPROJStringExportable {
     crs::GeodeticCRSPtr geodSrc{};
@@ -12332,14 +12409,18 @@ struct MyPROJStringExportableHorizVertical final
         _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();
 
         formatter->pushOmitZUnitConversion();
         geogDst->addAngularUnitConvertAndAxisSwap(formatter);
@@ -12385,7 +12466,9 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final
 
         formatter->popOmitZUnitConversion();
 
+        formatter->pushOmitHorizontalConversionInVertTransformation();
         verticalTransform->_exportToPROJString(formatter);
+        formatter->popOmitHorizontalConversionInVertTransformation();
 
         formatter->pushOmitZUnitConversion();
 
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index e7f8769c..38b407ee 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -6635,6 +6635,7 @@ struct PROJStringFormatter::Private {
     std::vector<InversionStackElt> inversionStack_{InversionStackElt()};
     bool omitProjLongLatIfPossible_ = false;
     std::vector<bool> omitZUnitConversion_{false};
+    std::vector<bool> omitHorizontalConversionInVertTransformation_{false};
     DatabaseContextPtr dbContext_{};
     bool useApproxTMerc_ = false;
     bool addNoDefs_ = true;
@@ -7739,6 +7740,25 @@ bool PROJStringFormatter::omitZUnitConversion() const {
 
 // ---------------------------------------------------------------------------
 
+void PROJStringFormatter::pushOmitHorizontalConversionInVertTransformation() {
+    d->omitHorizontalConversionInVertTransformation_.push_back(true);
+}
+
+// ---------------------------------------------------------------------------
+
+void PROJStringFormatter::popOmitHorizontalConversionInVertTransformation() {
+    assert(d->omitHorizontalConversionInVertTransformation_.size() > 1);
+    d->omitHorizontalConversionInVertTransformation_.pop_back();
+}
+
+// ---------------------------------------------------------------------------
+
+bool PROJStringFormatter::omitHorizontalConversionInVertTransformation() const {
+    return d->omitHorizontalConversionInVertTransformation_.back();
+}
+
+// ---------------------------------------------------------------------------
+
 void PROJStringFormatter::setLegacyCRSToCRSContext(bool legacyContext) {
     d->legacyCRSToCRSContext_ = legacyContext;
 }
diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist
index 9e03b003..f532a8ec 100644
--- a/test/cli/testprojinfo_out.dist
+++ b/test/cli/testprojinfo_out.dist
@@ -903,7 +903,7 @@ Operation No. 1:
 INVERSE(DERIVED_FROM(PROJ)):EPSG_4977_TO_EPSG_5613, Inverse of SWEREF99 to RH2000 height, unknown accuracy, Sweden - onshore
 
 PROJ string:
-+proj=vgridshift +grids=se_lantmateriet_SWEN17_RH2000.tif +multiplier=1
++proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=se_lantmateriet_SWEN17_RH2000.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 
 WKT2:2019 string:
 COORDINATEOPERATION["Inverse of SWEREF99 to RH2000 height",
@@ -967,7 +967,7 @@ Operation No. 1:
 INVERSE(DERIVED_FROM(EPSG)):8885, Inverse of RGF93 to NGF-IGN69 height (3), 0.01 m, France - mainland onshore
 
 PROJ string:
-+proj=vgridshift +grids=fr_ign_RAF18.tif +multiplier=1
++proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=fr_ign_RAF18.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 
 -------------------------------------
 Operation No. 2:
@@ -975,7 +975,7 @@ Operation No. 2:
 INVERSE(EPSG):10000, Inverse of RGF93 to NGF-IGN69 height (1), 0.5 m, France - mainland onshore
 
 PROJ string:
-+proj=vgridshift +grids=ggf97a.txt +multiplier=1
++proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=ggf97a.txt +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 
 Testing EPSG:32631 --3d
 PROJ.4 string:
@@ -1081,7 +1081,7 @@ Operation No. 1:
 DERIVED_FROM(EPSG):5656, GDA94 to AHD height (49), 0.03 m, Australia - mainland
 
 PROJ string:
-+proj=pipeline +step +inv +proj=vgridshift +grids=au_ga_AUSGeoid09_V1.01.tif +multiplier=1
++proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=au_ga_AUSGeoid09_V1.01.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 
 Testing -s "GDA2020" -t "AHD height" --grid-check none -o PROJ --spatial-test intersects
 Candidate operations found: 1
@@ -1091,7 +1091,7 @@ Operation No. 1:
 DERIVED_FROM(EPSG):8451, GDA2020 to AHD height (1), 0.03 m, Australia Christmas and Cocos - onshore
 
 PROJ string:
-+proj=pipeline +step +inv +proj=vgridshift +grids=au_ga_AUSGeoid2020_20180201.tif +multiplier=1
++proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=au_ga_AUSGeoid2020_20180201.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 
 Testing -k ellipsoid WGS84
 PROJ string:
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 51bbf50b..6ef6e686 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -4819,12 +4819,18 @@ TEST(operation, vertCRS_to_geogCRS_context) {
             authFactory->createCoordinateReferenceSystem("4979"), // WGS 84
             ctxt);
         ASSERT_EQ(list.size(), 2U);
-        EXPECT_EQ(list[1]->exportToPROJString(
-                      PROJStringFormatter::create(
-                          PROJStringFormatter::Convention::PROJ_5,
-                          authFactory->databaseContext())
-                          .get()),
-                  "+proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1");
+        EXPECT_EQ(
+            list[1]->exportToPROJString(
+                PROJStringFormatter::create(
+                    PROJStringFormatter::Convention::PROJ_5,
+                    authFactory->databaseContext())
+                    .get()),
+            "+proj=pipeline "
+            "+step +proj=axisswap +order=2,1 "
+            "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+            "+step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1 "
+            "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+            "+step +proj=axisswap +order=2,1");
     }
     {
         auto ctxt =
@@ -4837,7 +4843,12 @@ TEST(operation, vertCRS_to_geogCRS_context) {
         ASSERT_EQ(list.size(), 2U);
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
-            "+proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1");
+            "+proj=pipeline "
+            "+step +proj=axisswap +order=2,1 "
+            "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+            "+step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1 "
+            "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+            "+step +proj=axisswap +order=2,1");
     }
     {
         auto ctxt =
@@ -4850,9 +4861,13 @@ TEST(operation, vertCRS_to_geogCRS_context) {
         ASSERT_EQ(list.size(), 2U);
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
-            "+proj=pipeline +step +inv +proj=vgridshift "
-            "+grids=us_nga_egm08_25.tif "
-            "+multiplier=1");
+            "+proj=pipeline "
+            "+step +proj=axisswap +order=2,1 "
+            "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+            "+step +inv +proj=vgridshift +grids=us_nga_egm08_25.tif "
+            "+multiplier=1 "
+            "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+            "+step +proj=axisswap +order=2,1");
     }
     {
         auto ctxt =
@@ -4877,7 +4892,71 @@ TEST(operation, vertCRS_to_geogCRS_context) {
         ASSERT_EQ(list.size(), 1U);
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
-            "+proj=vgridshift +grids=nz_linz_nzgeoid2016.tif +multiplier=1");
+            "+proj=pipeline "
+            "+step +proj=axisswap +order=2,1 "
+            "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+            "+step +proj=vgridshift +grids=nz_linz_nzgeoid2016.tif "
+            "+multiplier=1 "
+            "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+            "+step +proj=axisswap +order=2,1");
+    }
+    {
+        // Test actually the database where we derive records using the more
+        // classic 'Geographic3D to GravityRelatedHeight' method from
+        // records using EPSG:9635
+        //'Geog3D to Geog2D+GravityRelatedHeight (US .gtx)' method
+        auto ctxt = CoordinateOperationContext::create(
+            AuthorityFactory::create(DatabaseContext::create(), std::string()),
+            nullptr, 0.0);
+        ctxt->setSpatialCriterion(
+            CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            // Baltic 1957 height
+            authFactory->createCoordinateReferenceSystem("8357"),
+            // ETRS89
+            authFactory->createCoordinateReferenceSystem("4937"), ctxt);
+        ASSERT_EQ(list.size(), 1U);
+        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=vgridshift "
+            "+grids=sk_gku_Slovakia_ETRS89h_to_Baltic1957.tif "
+            "+multiplier=1 "
+            "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+            "+step +proj=axisswap +order=2,1");
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, geog3DCRS_to_geog2DCRS_plus_vertCRS_context) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    {
+        auto ctxt =
+            CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+        ctxt->setSpatialCriterion(
+            CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            // ETRS89 (3D)
+            authFactory->createCoordinateReferenceSystem("4937"),
+            // ETRS89 + Baltic 1957 height
+            authFactory->createCoordinateReferenceSystem("8360"), ctxt);
+        ASSERT_GE(list.size(), 1U);
+        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 +inv +proj=vgridshift "
+            "+grids=sk_gku_Slovakia_ETRS89h_to_Baltic1957.tif +multiplier=1 "
+            "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+            "+step +proj=axisswap +order=2,1");
+
+        EXPECT_EQ(list[0]->inverse()->nameStr(),
+                  "Inverse of 'ETRS89 to ETRS89 + Baltic 1957 height (1)'");
     }
 }
 
@@ -6774,7 +6853,13 @@ static BoundCRSNNPtr createBoundVerticalCRS() {
 TEST(operation, transformation_height_to_PROJ_string) {
     auto transf = createBoundVerticalCRS()->transformation();
     EXPECT_EQ(transf->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1");
+              "+proj=pipeline "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +proj=vgridshift +grids=us_nga_egm08_25.tif "
+              "+multiplier=1 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
 
     auto grids = transf->gridsNeeded(DatabaseContext::create(), false);
     ASSERT_EQ(grids.size(), 1U);
@@ -6843,8 +6928,13 @@ TEST(operation, transformation_Geographic3D_to_GravityRelatedHeight_gtx) {
     // Check that we correctly inverse files in the case of
     // "Geographic3D to GravityRelatedHeight (US .gtx)"
     EXPECT_EQ(transf->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +inv +proj=vgridshift "
-              "+grids=naptrans2008.gtx +multiplier=1");
+              "+proj=pipeline "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +inv +proj=vgridshift "
+              "+grids=naptrans2008.gtx +multiplier=1 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
 }
 
 // ---------------------------------------------------------------------------
@@ -8111,6 +8201,44 @@ TEST(operation, compoundCRS_to_compoundCRS_WGS84_EGM96_to_WGS84_Belfast) {
 
 // ---------------------------------------------------------------------------
 
+TEST(
+    operation,
+    compoundCRS_to_compoundCRS_concatenated_operation_with_two_vert_transformation) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    {
+        auto ctxt =
+            CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+        ctxt->setSpatialCriterion(
+            CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            // ETRS89 + Baltic 1957 height
+            authFactory->createCoordinateReferenceSystem("8360"),
+            // ETRS89 + EVRF2007 height
+            authFactory->createCoordinateReferenceSystem("7423"), ctxt);
+        ASSERT_GE(list.size(), 1U);
+        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=vgridshift "
+            "+grids=sk_gku_Slovakia_ETRS89h_to_Baltic1957.tif +multiplier=1 "
+            "+step +inv +proj=vgridshift "
+            "+grids=sk_gku_Slovakia_ETRS89h_to_EVRF2007.tif +multiplier=1 "
+            "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+            "+step +proj=axisswap +order=2,1");
+        EXPECT_EQ(
+            list[0]->nameStr(),
+            "ETRS89 + Baltic 1957 height to ETRS89 + EVRF2007 height (1)");
+        EXPECT_EQ(list[0]->inverse()->nameStr(), "Inverse of 'ETRS89 + Baltic "
+                                                 "1957 height to ETRS89 + "
+                                                 "EVRF2007 height (1)'");
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, vertCRS_to_vertCRS) {
 
     auto vertcrs_m_obj = PROJStringParser().createFromPROJString("+vunits=m");
@@ -10710,6 +10838,10 @@ TEST(operation,
     // Test that even if the .gtx file is unknown, we export in the correct
     // direction
     EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()),
-              "+proj=pipeline +step +inv +proj=vgridshift +grids=foo.gtx "
-              "+multiplier=1");
+              "+proj=pipeline "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +inv +proj=vgridshift +grids=foo.gtx +multiplier=1 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
 }