Vertical CRS transformation: synthetize a vertical unit change transformation when needed, and also sort Null geographic offset transformation in last

5 files changed, 216 insertions, 78 deletions

diff --git a/include/proj/coordinateoperation.hpp b/include/proj/coordinateoperation.hpp
index 92b655f9..2f7c9bbf 100644
--- a/include/proj/coordinateoperation.hpp
+++ b/include/proj/coordinateoperation.hpp
@@ -1493,6 +1493,11 @@ class PROJ_GCC_DLL Transformation : public SingleOperation {
     PROJ_DLL TransformationNNPtr substitutePROJAlternativeGridNames(
         io::DatabaseContextNNPtr databaseContext) const;
 
+    PROJ_DLL static TransformationNNPtr createChangeVerticalUnit(
+        const util::PropertyMap &properties, const crs::CRSNNPtr &sourceCRSIn,
+        const crs::CRSNNPtr &targetCRSIn, const common::Scale &factor,
+        const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies);
+
     PROJ_PRIVATE :
         //! @cond Doxygen_Suppress
         PROJ_INTERNAL const std::string &
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index ed98832f..d3446460 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -106,6 +106,7 @@ constexpr double UTM_SOUTH_FALSE_NORTHING = 10000000.0;
 static const std::string INVERSE_OF = "Inverse of ";
 static const char *NULL_GEOCENTRIC_TRANSLATION = "Null geocentric translation";
 static const char *NULL_GEOGRAPHIC_OFFSET = "Null geographic offset";
+static const char *APPROXIMATE_TRANSFORMATION = " (approximate transformation)";
 //! @endcond
 
 //! @cond Doxygen_Suppress
@@ -7054,6 +7055,39 @@ TransformationNNPtr Transformation::createVerticalOffset(
 
 // ---------------------------------------------------------------------------
 
+/** \brief Instantiate a transformation based on the Change of Vertical Unit
+ * method.
+ *
+ * This method is defined as [EPSG:1069]
+ * (https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:method:EPSG::1069)
+ *
+ * @param properties See \ref general_properties of the conversion. If the name
+ * is not provided, it is automatically set.
+ * @param sourceCRSIn Source CRS.
+ * @param targetCRSIn Target CRS.
+ * @param factor Conversion factor
+ * @param accuracies Vector of positional accuracy (might be empty).
+ * @return a new Transformation.
+ */
+TransformationNNPtr Transformation::createChangeVerticalUnit(
+    const util::PropertyMap &properties, const crs::CRSNNPtr &sourceCRSIn,
+    const crs::CRSNNPtr &targetCRSIn, const common::Scale &factor,
+    const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies) {
+    return create(
+        properties, sourceCRSIn, targetCRSIn, nullptr,
+        createMethodMapNameEPSGCode(EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT),
+        VectorOfParameters{
+            createOpParamNameEPSGCode(
+                EPSG_CODE_PARAMETER_UNIT_CONVERSION_SCALAR),
+        },
+        VectorOfValues{
+            factor,
+        },
+        accuracies);
+}
+
+// ---------------------------------------------------------------------------
+
 static const char *getCRSQualifierStr(const crs::CRSPtr &crs) {
     auto geod = dynamic_cast<crs::GeodeticCRS *>(crs.get());
     if (geod) {
@@ -7481,6 +7515,17 @@ TransformationNNPtr Transformation::inverseAsTransformation() const {
                                  coordinateOperationAccuracies()));
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT) {
+        const double convFactor = parameterValueNumericAsSI(
+            EPSG_CODE_PARAMETER_UNIT_CONVERSION_SCALAR);
+        return d->registerInv(
+            shared_from_this(),
+            createChangeVerticalUnit(
+                createPropertiesForInverse(this, false, false), l_targetCRS,
+                l_sourceCRS, common::Scale(1.0 / convFactor),
+                coordinateOperationAccuracies()));
+    }
+
     return InverseTransformation::create(NN_NO_CHECK(
         util::nn_dynamic_pointer_cast<Transformation>(shared_from_this())));
 }
@@ -8835,6 +8880,31 @@ bool SingleOperation::exportToPROJStringGeneric(
                                       "conversion");
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT) {
+        double convFactor = parameterValueNumericAsSI(
+            EPSG_CODE_PARAMETER_UNIT_CONVERSION_SCALAR);
+        auto uom = common::UnitOfMeasure(std::string(), convFactor,
+                                         common::UnitOfMeasure::Type::LINEAR)
+                       .exportToPROJString();
+        auto reverse_uom =
+            common::UnitOfMeasure(std::string(), 1.0 / convFactor,
+                                  common::UnitOfMeasure::Type::LINEAR)
+                .exportToPROJString();
+        if (!uom.empty()) {
+            formatter->addStep("unitconvert");
+            formatter->addParam("z_in", uom);
+            formatter->addParam("z_out", "m");
+        } else if (!reverse_uom.empty()) {
+            formatter->addStep("unitconvert");
+            formatter->addParam("z_in", "m");
+            formatter->addParam("z_out", reverse_uom);
+        } else {
+            formatter->addStep("affine");
+            formatter->addParam("s33", convFactor);
+        }
+        return true;
+    }
+
     return false;
 }
 
@@ -9667,14 +9737,18 @@ struct PrecomputedOpCharacteristics {
     bool gridsAvailable_ = false;
     bool gridsKnown_ = false;
     size_t stepCount_ = 0;
+    bool isApprox_ = false;
+    bool isNullTransformation_ = false;
 
     PrecomputedOpCharacteristics() = default;
     PrecomputedOpCharacteristics(double area, double accuracy, bool hasGrids,
                                  bool gridsAvailable, bool gridsKnown,
-                                 size_t stepCount)
+                                 size_t stepCount, bool isApprox,
+                                 bool isNullTransformation)
         : area_(area), accuracy_(accuracy), hasGrids_(hasGrids),
           gridsAvailable_(gridsAvailable), gridsKnown_(gridsKnown),
-          stepCount_(stepCount) {}
+          stepCount_(stepCount), isApprox_(isApprox),
+          isNullTransformation_(isNullTransformation) {}
 };
 
 // ---------------------------------------------------------------------------
@@ -9701,6 +9775,22 @@ struct SortFunction {
         // CAUTION: the order of the comparisons is extremely important
         // to get the intended result.
 
+        if (!iterA->second.isApprox_ && iterB->second.isApprox_) {
+            return true;
+        }
+        if (iterA->second.isApprox_ && !iterB->second.isApprox_) {
+            return false;
+        }
+
+        if (!iterA->second.isNullTransformation_ &&
+            iterB->second.isNullTransformation_) {
+            return true;
+        }
+        if (iterA->second.isNullTransformation_ &&
+            !iterB->second.isNullTransformation_) {
+            return false;
+        }
+
         if (iterA->second.hasGrids_ && iterB->second.hasGrids_) {
             // Operations where grids are all available go before other
             if (iterA->second.gridsAvailable_ &&
@@ -9926,6 +10016,8 @@ struct FilterResults {
                     if (name.find(NULL_GEOGRAPHIC_OFFSET) ==
                             std::string::npos &&
                         name.find(NULL_GEOCENTRIC_TRANSLATION) ==
+                            std::string::npos &&
+                        name.find(APPROXIMATE_TRANSFORMATION) ==
                             std::string::npos) {
                         hasOpThatContainsAreaOfInterest = true;
                     }
@@ -9961,6 +10053,8 @@ struct FilterResults {
                     if (name.find(NULL_GEOGRAPHIC_OFFSET) ==
                             std::string::npos &&
                         name.find(NULL_GEOCENTRIC_TRANSLATION) ==
+                            std::string::npos &&
+                        name.find(APPROXIMATE_TRANSFORMATION) ==
                             std::string::npos) {
                         hasOpThatContainsAreaOfInterest = true;
                     }
@@ -10067,9 +10161,17 @@ struct FilterResults {
 
             const auto stepCount = getStepCount(op);
 
+            const bool isApprox =
+                op->nameStr().find(APPROXIMATE_TRANSFORMATION) !=
+                std::string::npos;
+            const bool isNullTransformation =
+                op->nameStr().find(NULL_GEOGRAPHIC_OFFSET) !=
+                    std::string::npos ||
+                op->nameStr().find(NULL_GEOCENTRIC_TRANSLATION) !=
+                    std::string::npos;
             map[op.get()] = PrecomputedOpCharacteristics(
                 area, getAccuracy(op), hasGrids, gridsAvailable, gridsKnown,
-                stepCount);
+                stepCount, isApprox, isNullTransformation);
         }
 
         // Sort !
@@ -10082,7 +10184,8 @@ struct FilterResults {
 
         // If we have more than one result, and than the last result is the
         // default "Null geographic offset" or "Null geocentric translation"
-        // operations we have synthetized, remove it as
+        // operations we have synthetized, and that at least one operation
+        // has the desired area of interest, remove it as
         // all previous results are necessarily better
         if (hasOpThatContainsAreaOfInterest && res.size() > 1) {
             const std::string &name = res.back()->nameStr();
@@ -11571,29 +11674,35 @@ CoordinateOperationFactory::Private::createOperations(
     if (vertSrc && vertDst) {
         const auto &srcDatum = vertSrc->datum();
         const auto &dstDatum = vertDst->datum();
-        if (srcDatum && dstDatum &&
-            srcDatum->_isEquivalentTo(
-                dstDatum.get(), util::IComparable::Criterion::EQUIVALENT)) {
-            const double convSrc = vertSrc->coordinateSystem()
-                                       ->axisList()[0]
-                                       ->unit()
-                                       .conversionToSI();
-            const double convDst = vertDst->coordinateSystem()
-                                       ->axisList()[0]
-                                       ->unit()
-                                       .conversionToSI();
-            if (convSrc != convDst) {
-                const double factor = convSrc / convDst;
+        const bool equivalentVDatum =
+            (srcDatum && dstDatum &&
+             srcDatum->_isEquivalentTo(
+                 dstDatum.get(), util::IComparable::Criterion::EQUIVALENT));
+
+        const double convSrc =
+            vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI();
+        const double convDst =
+            vertDst->coordinateSystem()->axisList()[0]->unit().conversionToSI();
+        if (convSrc != convDst) {
+            const double factor = convSrc / convDst;
+            auto name =
+                buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr());
+            if (!equivalentVDatum) {
+                name += APPROXIMATE_TRANSFORMATION;
+                auto conv = Transformation::createChangeVerticalUnit(
+                    util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                            name),
+                    sourceCRS, targetCRS, common::Scale(factor), {});
+                res.push_back(conv);
+            } else {
                 auto conv = Conversion::createChangeVerticalUnit(
-                    util::PropertyMap().set(
-                        common::IdentifiedObject::NAME_KEY,
-                        buildTransfName(sourceCRS->nameStr(),
-                                        targetCRS->nameStr())),
+                    util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                            name),
                     common::Scale(factor));
                 conv->setCRSs(sourceCRS, targetCRS, nullptr);
                 res.push_back(conv);
-                return res;
             }
+            return res;
         }
     }
 
diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist
index ebc59c40..15372803 100644
--- a/test/cli/testprojinfo_out.dist
+++ b/test/cli/testprojinfo_out.dist
@@ -168,9 +168,9 @@ DERIVED_FROM(EPSG):1312, NAD27 to NAD83 (3), 1.0 m, Canada
 DERIVED_FROM(EPSG):1313, NAD27 to NAD83 (4), 1.5 m, Canada - NAD27
 DERIVED_FROM(EPSG):1241, NAD27 to NAD83 (1), 0.15 m, USA - CONUS including EEZ
 DERIVED_FROM(EPSG):1243, NAD27 to NAD83 (2), 0.5 m, USA - Alaska including EEZ
-unknown id, Null geographic offset from NAD27 to NAD83, unknown accuracy, World
 EPSG:1462, NAD27 to NAD83 (5), 1.0 m, Canada - Quebec
 EPSG:1573, NAD27 to NAD83 (6), 1.5 m, Canada - Quebec
+unknown id, Null geographic offset from NAD27 to NAD83, unknown accuracy, World
 
 Testing projinfo -s NAD27 -t NAD83 --grid-check none --spatial-test intersects
 -------------------------------------
@@ -366,13 +366,13 @@ COORDINATEOPERATION["NAD27 to NAD83 (2)",
 -------------------------------------
 Operation n°5:
 
-unknown id, Null geographic offset from NAD27 to NAD83, unknown accuracy, World
+EPSG:1462, NAD27 to NAD83 (5), 1.0 m, Canada - Quebec
 
 PROJ string:
-
++proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=GS2783v1.QUE +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 
 WKT2_2018 string:
-COORDINATEOPERATION["Null geographic offset from NAD27 to NAD83",
+COORDINATEOPERATION["NAD27 to NAD83 (5)",
     SOURCECRS[
         GEOGCRS["NAD27",
             DATUM["North American Datum 1927",
@@ -386,12 +386,7 @@ COORDINATEOPERATION["Null geographic offset from NAD27 to NAD83",
                     ANGLEUNIT["degree",0.0174532925199433]],
                 AXIS["geodetic longitude (Lon)",east,
                     ORDER[2],
-                    ANGLEUNIT["degree",0.0174532925199433]],
-            USAGE[
-                SCOPE["unknown"],
-                AREA["North America - NAD27"],
-                BBOX[7.15,167.65,83.17,-47.74]],
-            ID["EPSG",4267]]],
+                    ANGLEUNIT["degree",0.0174532925199433]]]],
     TARGETCRS[
         GEOGCRS["NAD83",
             DATUM["North American Datum 1983",
@@ -405,35 +400,27 @@ COORDINATEOPERATION["Null geographic offset from NAD27 to NAD83",
                     ANGLEUNIT["degree",0.0174532925199433]],
                 AXIS["geodetic longitude (Lon)",east,
                     ORDER[2],
-                    ANGLEUNIT["degree",0.0174532925199433]],
-            USAGE[
-                SCOPE["unknown"],
-                AREA["North America - NAD83"],
-                BBOX[14.92,167.65,86.46,-47.74]],
-            ID["EPSG",4269]]],
-    METHOD["Geographic2D offsets",
-        ID["EPSG",9619]],
-    PARAMETER["Latitude offset",0,
-        ANGLEUNIT["degree",0.0174532925199433],
-        ID["EPSG",8601]],
-    PARAMETER["Longitude offset",0,
-        ANGLEUNIT["degree",0.0174532925199433],
-        ID["EPSG",8602]],
+                    ANGLEUNIT["degree",0.0174532925199433]]]],
+    METHOD["NTv1",
+        ID["EPSG",9614]],
+    PARAMETERFILE["Latitude and longitude difference file","GS2783v1.QUE"],
+    OPERATIONACCURACY[1.0],
     USAGE[
         SCOPE["unknown"],
-        AREA["World"],
-        BBOX[-90,-180,90,180]]]
+        AREA["Canada - Quebec"],
+        BBOX[44.99,-79.85,62.62,-57.1]],
+    ID["EPSG",1462]]
 
 -------------------------------------
 Operation n°6:
 
-EPSG:1462, NAD27 to NAD83 (5), 1.0 m, Canada - Quebec
+EPSG:1573, NAD27 to NAD83 (6), 1.5 m, Canada - Quebec
 
 PROJ string:
-+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=GS2783v1.QUE +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
++proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=QUE27-83.gsb +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
 
 WKT2_2018 string:
-COORDINATEOPERATION["NAD27 to NAD83 (5)",
+COORDINATEOPERATION["NAD27 to NAD83 (6)",
     SOURCECRS[
         GEOGCRS["NAD27",
             DATUM["North American Datum 1927",
@@ -462,26 +449,26 @@ COORDINATEOPERATION["NAD27 to NAD83 (5)",
                 AXIS["geodetic longitude (Lon)",east,
                     ORDER[2],
                     ANGLEUNIT["degree",0.0174532925199433]]]],
-    METHOD["NTv1",
-        ID["EPSG",9614]],
-    PARAMETERFILE["Latitude and longitude difference file","GS2783v1.QUE"],
-    OPERATIONACCURACY[1.0],
+    METHOD["NTv2",
+        ID["EPSG",9615]],
+    PARAMETERFILE["Latitude and longitude difference file","QUE27-83.gsb"],
+    OPERATIONACCURACY[1.5],
     USAGE[
         SCOPE["unknown"],
         AREA["Canada - Quebec"],
         BBOX[44.99,-79.85,62.62,-57.1]],
-    ID["EPSG",1462]]
+    ID["EPSG",1573]]
 
 -------------------------------------
 Operation n°7:
 
-EPSG:1573, NAD27 to NAD83 (6), 1.5 m, Canada - Quebec
+unknown id, Null geographic offset from NAD27 to NAD83, unknown accuracy, World
 
 PROJ string:
-+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=QUE27-83.gsb +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
+
 
 WKT2_2018 string:
-COORDINATEOPERATION["NAD27 to NAD83 (6)",
+COORDINATEOPERATION["Null geographic offset from NAD27 to NAD83",
     SOURCECRS[
         GEOGCRS["NAD27",
             DATUM["North American Datum 1927",
@@ -495,7 +482,12 @@ COORDINATEOPERATION["NAD27 to NAD83 (6)",
                     ANGLEUNIT["degree",0.0174532925199433]],
                 AXIS["geodetic longitude (Lon)",east,
                     ORDER[2],
-                    ANGLEUNIT["degree",0.0174532925199433]]]],
+                    ANGLEUNIT["degree",0.0174532925199433]],
+            USAGE[
+                SCOPE["unknown"],
+                AREA["North America - NAD27"],
+                BBOX[7.15,167.65,83.17,-47.74]],
+            ID["EPSG",4267]]],
     TARGETCRS[
         GEOGCRS["NAD83",
             DATUM["North American Datum 1983",
@@ -509,16 +501,24 @@ COORDINATEOPERATION["NAD27 to NAD83 (6)",
                     ANGLEUNIT["degree",0.0174532925199433]],
                 AXIS["geodetic longitude (Lon)",east,
                     ORDER[2],
-                    ANGLEUNIT["degree",0.0174532925199433]]]],
-    METHOD["NTv2",
-        ID["EPSG",9615]],
-    PARAMETERFILE["Latitude and longitude difference file","QUE27-83.gsb"],
-    OPERATIONACCURACY[1.5],
+                    ANGLEUNIT["degree",0.0174532925199433]],
+            USAGE[
+                SCOPE["unknown"],
+                AREA["North America - NAD83"],
+                BBOX[14.92,167.65,86.46,-47.74]],
+            ID["EPSG",4269]]],
+    METHOD["Geographic2D offsets",
+        ID["EPSG",9619]],
+    PARAMETER["Latitude offset",0,
+        ANGLEUNIT["degree",0.0174532925199433],
+        ID["EPSG",8601]],
+    PARAMETER["Longitude offset",0,
+        ANGLEUNIT["degree",0.0174532925199433],
+        ID["EPSG",8602]],
     USAGE[
         SCOPE["unknown"],
-        AREA["Canada - Quebec"],
-        BBOX[44.99,-79.85,62.62,-57.1]],
-    ID["EPSG",1573]]
+        AREA["World"],
+        BBOX[-90,-180,90,180]]]
 
 Testing projinfo -s EPSG:4230 -t EPSG:4258 --bbox 8,54.51,15.24,57.8 --summary
 Candidate operations found: 1
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 6205a9b8..8cfed9ad 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -1421,7 +1421,7 @@ TEST_F(CApi, proj_create_operations_with_pivot) {
         ASSERT_NE(res, nullptr);
         ObjListKeeper keeper_res(res);
         EXPECT_EQ(proj_list_get_count(res), 7);
-        auto op = proj_list_get(m_ctxt, res, 1);
+        auto op = proj_list_get(m_ctxt, res, 0);
         ASSERT_NE(op, nullptr);
         ObjectKeeper keeper_op(op);
 
@@ -1451,7 +1451,7 @@ TEST_F(CApi, proj_create_operations_with_pivot) {
         ASSERT_NE(res, nullptr);
         ObjListKeeper keeper_res(res);
         // includes results from ESRI
-        EXPECT_EQ(proj_list_get_count(res), 5);
+        EXPECT_EQ(proj_list_get_count(res), 4);
         auto op = proj_list_get(m_ctxt, res, 0);
         ASSERT_NE(op, nullptr);
         ObjectKeeper keeper_op(op);
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 9a968378..42f8fe76 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -4362,7 +4362,7 @@ TEST(operation, geogCRS_to_geogCRS_context_inverse_needed) {
             authFactory->createCoordinateReferenceSystem("4275"), // NTF
             authFactory->createCoordinateReferenceSystem("4258"), // ETRS89
             ctxt);
-        ASSERT_EQ(list.size(), 3U);
+        ASSERT_EQ(list.size(), 2U);
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
             "+proj=pipeline +step +proj=push +v_3 +step +proj=axisswap "
@@ -4376,12 +4376,6 @@ TEST(operation, geogCRS_to_geogCRS_context_inverse_needed) {
                           PROJStringFormatter::Convention::PROJ_5,
                           authFactory->databaseContext())
                           .get()),
-                  "");
-        EXPECT_EQ(list[2]->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=hgridshift +grids=ntf_r93.gsb +step +proj=unitconvert "
@@ -6136,7 +6130,8 @@ TEST(operation, compoundCRS_to_compoundCRS_context) {
         authFactory->createCoordinateReferenceSystem("7406"),
         // NAD83(NSRS2007) + NAVD88 height
         authFactory->createCoordinateReferenceSystem("5500"), ctxt);
-    ASSERT_EQ(list.size(), 88U);
+    // 152 or 155 depending if the VERTCON grids are there
+    ASSERT_GE(list.size(), 152U);
     EXPECT_EQ(list[0]->nameStr(), "NGVD29 height (ftUS) to NAVD88 height (3) + "
                                   "NAD27 to WGS 84 (79) + Inverse of "
                                   "NAD83(NSRS2007) to WGS 84 (1)");
@@ -6147,6 +6142,35 @@ TEST(operation, compoundCRS_to_compoundCRS_context) {
               "+step +proj=hgridshift +grids=conus +step +proj=push +v_3 +step "
               "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
               "+order=2,1 +step +proj=pop +v_3");
+
+    bool foundApprox = false;
+    for (size_t i = 0; i < list.size(); i++) {
+        auto projString =
+            list[i]->exportToPROJString(PROJStringFormatter::create().get());
+        EXPECT_TRUE(
+            projString.find("+proj=pipeline +step +proj=axisswap +order=2,1 "
+                            "+step +proj=unitconvert +xy_in=deg +z_in=us-ft "
+                            "+xy_out=rad +z_out=m") == 0)
+            << list[i]->nameStr();
+        if (list[i]->nameStr().find("Transformation from NGVD29 height (ftUS) "
+                                    "to NAVD88 height (approximate "
+                                    "transformation)") == 0) {
+            EXPECT_EQ(list[i]->nameStr(),
+                      "Transformation from NGVD29 height (ftUS) to NAVD88 "
+                      "height (approximate transformation) + NAD27 to WGS 84 "
+                      "(79) + Inverse of NAD83(NSRS2007) to WGS 84 (1)");
+            EXPECT_EQ(projString,
+                      "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
+                      "+proj=unitconvert +xy_in=deg +z_in=us-ft +xy_out=rad "
+                      "+z_out=m +step +proj=hgridshift +grids=conus +step "
+                      "+proj=push +v_3 +step +proj=unitconvert +xy_in=rad "
+                      "+xy_out=deg +step +proj=axisswap +order=2,1 +step "
+                      "+proj=pop +v_3");
+            foundApprox = true;
+            break;
+        }
+    }
+    EXPECT_TRUE(foundApprox);
 }
 
 // ---------------------------------------------------------------------------