Merge pull request #1710 from rouault/geoid_model

Add support for GEOIDMODEL

3 files changed, 317 insertions, 43 deletions

diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index b8310ce5..bdadc8b8 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -3602,7 +3602,7 @@ TEST_F(CApi, proj_as_projjson) {
         EXPECT_EQ(std::string(projjson),
                   "{\n"
                   "  \"$schema\": "
-                  "\"https://proj.org/schemas/v0.1/projjson.schema.json\",\n"
+                  "\"https://proj.org/schemas/v0.2/projjson.schema.json\",\n"
                   "  \"type\": \"Ellipsoid\",\n"
                   "  \"name\": \"WGS 84\",\n"
                   "  \"semi_major_axis\": 6378137,\n"
@@ -4215,4 +4215,111 @@ TEST_F(
     EXPECT_NEAR(outcoord.xyzt.z, -32.5823, 1e-3);
 }
 
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_create_vertical_crs_ex) {
+
+    // NAD83(2011) / UTM zone 11N
+    auto horiz_crs = proj_create_from_database(m_ctxt, "EPSG", "6340",
+                                               PJ_CATEGORY_CRS, false, nullptr);
+    ObjectKeeper keeper_horiz_crs(horiz_crs);
+    ASSERT_NE(horiz_crs, nullptr);
+
+    auto vert_crs = proj_create_vertical_crs_ex(
+        m_ctxt, "myVertCRS (ftUS)", "myVertDatum", nullptr, nullptr,
+        "US survey foot", 0.304800609601219, "PROJ @foo.gtx", nullptr, nullptr,
+        nullptr, nullptr);
+    ObjectKeeper keeper_vert_crs(vert_crs);
+    ASSERT_NE(vert_crs, nullptr);
+
+    auto compound =
+        proj_create_compound_crs(m_ctxt, "Compound", horiz_crs, vert_crs);
+    ObjectKeeper keeper_compound(compound);
+    ASSERT_NE(compound, nullptr);
+
+    // NAD83(2011) 3D
+    PJ *geog_crs = proj_create(m_ctxt, "EPSG:6319");
+    ObjectKeeper keeper_geog_crs(geog_crs);
+    ASSERT_NE(geog_crs, nullptr);
+
+    auto P = proj_create_crs_to_crs_from_pj(m_ctxt, compound, geog_crs, nullptr,
+                                            nullptr);
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+
+    auto name = proj_get_name(P);
+    ASSERT_TRUE(name != nullptr);
+    EXPECT_EQ(name,
+              std::string("Inverse of UTM zone 11N + "
+                          "Transformation from myVertCRS (ftUS) to myVertCRS + "
+                          "Transformation from myVertCRS to NAD83(2011)"));
+
+    auto proj_5 = proj_as_proj_string(m_ctxt, P, PJ_PROJ_5, nullptr);
+    ASSERT_NE(proj_5, nullptr);
+    EXPECT_EQ(std::string(proj_5),
+              "+proj=pipeline "
+              "+step +inv +proj=utm +zone=11 +ellps=GRS80 "
+              "+step +proj=unitconvert +z_in=us-ft +z_out=m "
+              "+step +proj=vgridshift +grids=@foo.gtx +multiplier=1 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_create_vertical_crs_ex_with_geog_crs) {
+
+    // NAD83(2011) / UTM zone 11N
+    auto horiz_crs = proj_create_from_database(m_ctxt, "EPSG", "6340",
+                                               PJ_CATEGORY_CRS, false, nullptr);
+    ObjectKeeper keeper_horiz_crs(horiz_crs);
+    ASSERT_NE(horiz_crs, nullptr);
+
+    // WGS84
+    PJ *wgs84 = proj_create(m_ctxt, "EPSG:4979");
+    ObjectKeeper keeper_wgs84(wgs84);
+    ASSERT_NE(wgs84, nullptr);
+
+    auto vert_crs = proj_create_vertical_crs_ex(
+        m_ctxt, "myVertCRS", "myVertDatum", nullptr, nullptr, "US survey foot",
+        0.304800609601219, "PROJ @foo.gtx", nullptr, nullptr, wgs84, nullptr);
+    ObjectKeeper keeper_vert_crs(vert_crs);
+    ASSERT_NE(vert_crs, nullptr);
+
+    auto compound =
+        proj_create_compound_crs(m_ctxt, "Compound", horiz_crs, vert_crs);
+    ObjectKeeper keeper_compound(compound);
+    ASSERT_NE(compound, nullptr);
+
+    // NAD83(2011) 3D
+    PJ *geog_crs = proj_create(m_ctxt, "EPSG:6319");
+    ObjectKeeper keeper_geog_crs(geog_crs);
+    ASSERT_NE(geog_crs, nullptr);
+
+    auto P = proj_create_crs_to_crs_from_pj(m_ctxt, compound, geog_crs, nullptr,
+                                            nullptr);
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+
+    auto name = proj_get_name(P);
+    ASSERT_TRUE(name != nullptr);
+    EXPECT_EQ(
+        name,
+        std::string("Inverse of UTM zone 11N + "
+                    "Ballpark geographic offset from NAD83(2011) to WGS 84 + "
+                    "Transformation from myVertCRS to myVertCRS (metre) + "
+                    "Transformation from myVertCRS (metre) to WGS 84 + "
+                    "Ballpark geographic offset from WGS 84 to NAD83(2011)"));
+
+    auto proj_5 = proj_as_proj_string(m_ctxt, P, PJ_PROJ_5, nullptr);
+    ASSERT_NE(proj_5, nullptr);
+    EXPECT_EQ(std::string(proj_5),
+              "+proj=pipeline "
+              "+step +inv +proj=utm +zone=11 +ellps=GRS80 "
+              "+step +proj=unitconvert +z_in=us-ft +z_out=m "
+              "+step +proj=vgridshift +grids=@foo.gtx +multiplier=1 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
+}
+
 } // namespace
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 8aff0908..a71b63bb 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -1994,16 +1994,21 @@ TEST(wkt_parse, vertcrs_VRF_WKT2) {
 // ---------------------------------------------------------------------------
 
 TEST(wkt_parse, vertcrs_with_GEOIDMODEL) {
-    auto wkt = "VERTCRS[\"CGVD2013\","
-               "    VRF[\"Canadian Geodetic Vertical Datum of 2013\"],"
-               "    CS[vertical,1],"
-               "        AXIS[\"gravity-related height (H)\",up],"
-               "        LENGTHUNIT[\"metre\",1.0],"
-               "    GEOIDMODEL[\"CGG2013\",ID[\"EPSG\",6648]]]";
+    auto wkt = "VERTCRS[\"CGVD2013\",\n"
+               "    VDATUM[\"Canadian Geodetic Vertical Datum of 2013\"],\n"
+               "    CS[vertical,1],\n"
+               "        AXIS[\"gravity-related height (H)\",up,\n"
+               "            LENGTHUNIT[\"metre\",1]],\n"
+               "    GEOIDMODEL[\"CGG2013\",\n"
+               "        ID[\"EPSG\",6648]]]";
 
     auto obj = WKTParser().createFromWKT(wkt);
     auto crs = nn_dynamic_pointer_cast<VerticalCRS>(obj);
     ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(
+        crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT2_2019).get()),
+        wkt);
 }
 
 // ---------------------------------------------------------------------------
@@ -2060,6 +2065,71 @@ TEST(wkt_parse, vertcrs_WKT1_GDAL_minimum) {
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, vertcrs_WKT1_LAS_ftUS) {
+    auto wkt = "VERT_CS[\"NAVD88 - Geoid03 (Feet)\","
+               "    VERT_DATUM[\"unknown\",2005],"
+               "    UNIT[\"US survey foot\",0.3048006096012192,"
+               "        AUTHORITY[\"EPSG\",\"9003\"]],"
+               "    AXIS[\"Up\",UP]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<VerticalCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->nameStr(), "NAVD88 height (ftUS)");
+    ASSERT_EQ(crs->identifiers().size(), 1U);
+    EXPECT_EQ(crs->identifiers()[0]->code(), "6360");
+    EXPECT_EQ(*(crs->identifiers()[0]->codeSpace()), "EPSG");
+
+    const auto &geoidModel = crs->geoidModel();
+    ASSERT_TRUE(!geoidModel.empty());
+    EXPECT_EQ(geoidModel[0]->nameStr(), "GEOID03");
+
+    auto datum = crs->datum();
+    EXPECT_EQ(datum->nameStr(), "North American Vertical Datum 1988");
+    ASSERT_EQ(datum->identifiers().size(), 1U);
+    EXPECT_EQ(datum->identifiers()[0]->code(), "5103");
+    EXPECT_EQ(*(datum->identifiers()[0]->codeSpace()), "EPSG");
+
+    const auto &axis = crs->coordinateSystem()->axisList()[0];
+    EXPECT_EQ(axis->direction(), AxisDirection::UP);
+    EXPECT_EQ(axis->unit().name(), "US survey foot");
+    EXPECT_NEAR(axis->unit().conversionToSI(), 0.3048006096012192, 1e-16);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse, vertcrs_WKT1_LAS_metre) {
+    auto wkt = "VERT_CS[\"NAVD88 via Geoid09\","
+               "    VERT_DATUM[\"unknown\",2005],"
+               "    UNIT[\"metre\",1.0,"
+               "        AUTHORITY[\"EPSG\",\"9001\"]],"
+               "    AXIS[\"Up\",UP]]";
+
+    auto obj = WKTParser().createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<VerticalCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->nameStr(), "NAVD88 height");
+    ASSERT_EQ(crs->identifiers().size(), 1U);
+    EXPECT_EQ(crs->identifiers()[0]->code(), "5703");
+    EXPECT_EQ(*(crs->identifiers()[0]->codeSpace()), "EPSG");
+
+    const auto &geoidModel = crs->geoidModel();
+    ASSERT_TRUE(!geoidModel.empty());
+    EXPECT_EQ(geoidModel[0]->nameStr(), "GEOID09");
+
+    auto datum = crs->datum();
+    EXPECT_EQ(datum->nameStr(), "North American Vertical Datum 1988");
+    ASSERT_EQ(datum->identifiers().size(), 1U);
+    EXPECT_EQ(datum->identifiers()[0]->code(), "5103");
+    EXPECT_EQ(*(datum->identifiers()[0]->codeSpace()), "EPSG");
+
+    const auto &axis = crs->coordinateSystem()->axisList()[0];
+    EXPECT_EQ(axis->direction(), AxisDirection::UP);
+    EXPECT_EQ(axis->unit(), UnitOfMeasure::METRE);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(wkt_parse, dynamic_vertical_reference_frame) {
     auto obj = WKTParser().createFromWKT(
         "VERTCRS[\"RH2000\","
@@ -10842,6 +10912,45 @@ TEST(json_import, vertical_crs_with_datum_ensemble) {
 
 // ---------------------------------------------------------------------------
 
+TEST(json_import, vertical_crs_with_geoid_model) {
+    auto json = "{\n"
+                "  \"$schema\": \"foo\",\n"
+                "  \"type\": \"VerticalCRS\",\n"
+                "  \"name\": \"CGVD2013\",\n"
+                "  \"datum\": {\n"
+                "    \"type\": \"VerticalReferenceFrame\",\n"
+                "    \"name\": \"Canadian Geodetic Vertical Datum of 2013\"\n"
+                "  },\n"
+                "  \"coordinate_system\": {\n"
+                "    \"subtype\": \"vertical\",\n"
+                "    \"axis\": [\n"
+                "      {\n"
+                "        \"name\": \"Gravity-related height\",\n"
+                "        \"abbreviation\": \"H\",\n"
+                "        \"direction\": \"up\",\n"
+                "        \"unit\": \"metre\"\n"
+                "      }\n"
+                "    ]\n"
+                "  },\n"
+                "  \"geoid_model\": {\n"
+                "    \"name\": \"CGG2013\",\n"
+                "    \"id\": {\n"
+                "      \"authority\": \"EPSG\",\n"
+                "      \"code\": 6648\n"
+                "    }\n"
+                "  }\n"
+                "}";
+
+    // No database
+    auto obj = createFromUserInput(json, nullptr);
+    auto vcrs = nn_dynamic_pointer_cast<VerticalCRS>(obj);
+    ASSERT_TRUE(vcrs != nullptr);
+    EXPECT_EQ(vcrs->exportToJSON(&(JSONFormatter::create()->setSchema("foo"))),
+              json);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(json_import, parametric_crs) {
     auto json = "{\n"
                 "  \"$schema\": \"foo\",\n"
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 8df785b1..07b0daea 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -7072,17 +7072,13 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
                   "+proj=pipeline "
                   "+step +proj=axisswap +order=2,1 "
                   "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
-                  // Inv here since the grid is not known
-                  "+step +inv +proj=vgridshift +grids=geoid09_conus.bin "
+                  "+step +proj=vgridshift +grids=geoid09_conus.gtx "
                   "+multiplier=1 "
                   "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
                   "+step +proj=axisswap +order=2,1");
     }
 
-    // CompoundCRS to Geog3DCRS, with same vertical unit, and with
-    // ellipsoid height <--> vertical height correction that requires a
-    // horizontal adjustment before and after (which is empty in practice here
-    // as NAD83 to NAD83(2011) is no-op)
+    // NAD83 + NAVD88 height --> WGS 84
     {
         auto ctxt =
             CoordinateOperationContext::create(authFactory, nullptr, 0.0);
@@ -7105,10 +7101,10 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
         ASSERT_GE(list.size(), 2U);
 
         EXPECT_EQ(list[0]->nameStr(),
-                  "NAD83 to NAD83(2011) (1) + "
-                  "Inverse of NAD83(2011) to NAVD88 height (1) + "
-                  "Inverse of NAD83 to NAD83(2011) (1) + "
-                  "NAD83 to WGS 84 (1)");
+                  "NAD83 to WGS 84 (1) + "
+                  "Inverse of NAD83(NSRS2007) to WGS 84 (1) + "
+                  "Inverse of NAD83(NSRS2007) to NAVD88 height (1) + "
+                  "NAD83(NSRS2007) to WGS 84 (1)");
         EXPECT_EQ(list[0]->exportToPROJString(
                       PROJStringFormatter::create(
                           PROJStringFormatter::Convention::PROJ_5,
@@ -7117,30 +7113,10 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
                   "+proj=pipeline "
                   "+step +proj=axisswap +order=2,1 "
                   "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
-                  "+step +proj=vgridshift +grids=g2012bu0.gtx "
+                  "+step +proj=vgridshift +grids=geoid09_conus.gtx "
                   "+multiplier=1 "
                   "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
                   "+step +proj=axisswap +order=2,1");
-
-        // Shows vertical step, and then horizontal step
-        EXPECT_EQ(list[1]->nameStr(),
-                  "NAD83 to NAD83(2011) (1) + "
-                  "Inverse of NAD83(2011) to NAVD88 height (1) + "
-                  "Inverse of NAD83 to NAD83(2011) (1) + "
-                  "NAD83 to WGS 84 (18)");
-        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=g2012bu0.gtx "
-                  "+multiplier=1 "
-                  "+step +proj=hgridshift +grids=FL "
-                  "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
-                  "+step +proj=axisswap +order=2,1");
     }
 
     // Another variation, but post horizontal adjustment is in two steps
@@ -7166,7 +7142,7 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
         ASSERT_GE(list.size(), 2U);
 
         EXPECT_EQ(list[0]->nameStr(),
-                  "Inverse of NAD83(2011) to NAVD88 height (1) + "
+                  "Inverse of NAD83(2011) to NAVD88 height (3) + "
                   "Inverse of NAD83 to NAD83(2011) (1) + "
                   "NAD83 to WGS 84 (1)");
         EXPECT_EQ(list[0]->exportToPROJString(
@@ -7177,14 +7153,14 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
                   "+proj=pipeline "
                   "+step +proj=axisswap +order=2,1 "
                   "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
-                  "+step +proj=vgridshift +grids=g2012bu0.gtx "
+                  "+step +proj=vgridshift +grids=g2018u0.gtx "
                   "+multiplier=1 "
                   "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
                   "+step +proj=axisswap +order=2,1");
 
         // Shows vertical step, and then horizontal step
         EXPECT_EQ(list[1]->nameStr(),
-                  "Inverse of NAD83(2011) to NAVD88 height (1) + "
+                  "Inverse of NAD83(2011) to NAVD88 height (3) + "
                   "Inverse of NAD83 to NAD83(2011) (1) + "
                   "NAD83 to WGS 84 (18)");
         EXPECT_EQ(list[1]->exportToPROJString(
@@ -7195,7 +7171,7 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
                   "+proj=pipeline "
                   "+step +proj=axisswap +order=2,1 "
                   "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
-                  "+step +proj=vgridshift +grids=g2012bu0.gtx "
+                  "+step +proj=vgridshift +grids=g2018u0.gtx "
                   "+multiplier=1 "
                   "+step +proj=hgridshift +grids=FL "
                   "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
@@ -7227,7 +7203,7 @@ TEST(operation, compoundCRS_to_geogCRS_2D_promote_to_3D_context) {
                                                                ctxt);
     // The checked value is not that important, but in case this changes,
     // likely due to a EPSG upgrade, worth checking
-    ASSERT_EQ(listCompoundToGeog2D.size(), 469U);
+    ASSERT_EQ(listCompoundToGeog2D.size(), 467U);
 
     auto listGeog2DToCompound =
         CoordinateOperationFactory::create()->createOperations(dst, nnSrc,
@@ -7623,6 +7599,88 @@ TEST(
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, compoundCRS_of_vertCRS_with_geoid_model_to_geogCRS) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    auto wkt =
+        "COMPOUNDCRS[\"NAD83 / Pennsylvania South + NAVD88 height\",\n"
+        "    PROJCRS[\"NAD83 / Pennsylvania South\",\n"
+        "        BASEGEOGCRS[\"NAD83\",\n"
+        "            DATUM[\"North American Datum 1983\",\n"
+        "                ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n"
+        "                    LENGTHUNIT[\"metre\",1]]],\n"
+        "            PRIMEM[\"Greenwich\",0,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433]]],\n"
+        "        CONVERSION[\"SPCS83 Pennsylvania South zone (meters)\",\n"
+        "            METHOD[\"Lambert Conic Conformal (2SP)\",\n"
+        "                ID[\"EPSG\",9802]],\n"
+        "            PARAMETER[\"Latitude of false origin\",39.3333333333333,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8821]],\n"
+        "            PARAMETER[\"Longitude of false origin\",-77.75,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8822]],\n"
+        "            PARAMETER[\"Latitude of 1st standard "
+        "parallel\",40.9666666666667,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8823]],\n"
+        "            PARAMETER[\"Latitude of 2nd standard "
+        "parallel\",39.9333333333333,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8824]],\n"
+        "            PARAMETER[\"Easting at false origin\",600000,\n"
+        "                LENGTHUNIT[\"metre\",1],\n"
+        "                ID[\"EPSG\",8826]],\n"
+        "            PARAMETER[\"Northing at false origin\",0,\n"
+        "                LENGTHUNIT[\"metre\",1],\n"
+        "                ID[\"EPSG\",8827]]],\n"
+        "        CS[Cartesian,2],\n"
+        "            AXIS[\"easting (X)\",east,\n"
+        "                ORDER[1],\n"
+        "                LENGTHUNIT[\"metre\",1]],\n"
+        "            AXIS[\"northing (Y)\",north,\n"
+        "                ORDER[2],\n"
+        "                LENGTHUNIT[\"metre\",1]]],\n"
+        "    VERTCRS[\"NAVD88 height\",\n"
+        "        VDATUM[\"North American Vertical Datum 1988\"],\n"
+        "        CS[vertical,1],\n"
+        "            AXIS[\"gravity-related height (H)\",up,\n"
+        "                LENGTHUNIT[\"metre\",1]],\n"
+        "        GEOIDMODEL[\"GEOID12B\"]]]";
+    auto srcObj =
+        createFromUserInput(wkt, authFactory->databaseContext(), false);
+    auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
+    ASSERT_TRUE(src != nullptr);
+    auto dst = authFactory->createCoordinateReferenceSystem("4269"); // NAD83
+
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        NN_NO_CHECK(src), dst, ctxt);
+    ASSERT_TRUE(!list.empty());
+    EXPECT_EQ(list[0]->nameStr(),
+              "Inverse of SPCS83 Pennsylvania South zone (meters) + "
+              "Ballpark geographic offset from NAD83 to NAD83(2011) + "
+              "Inverse of NAD83(2011) to NAVD88 height (1) + "
+              "Ballpark geographic offset from NAD83(2011) to NAD83");
+    auto op_proj =
+        list[0]->exportToPROJString(PROJStringFormatter::create().get());
+    EXPECT_EQ(op_proj,
+              "+proj=pipeline "
+              "+step +inv +proj=lcc +lat_0=39.3333333333333 +lon_0=-77.75 "
+              "+lat_1=40.9666666666667 +lat_2=39.9333333333333 +x_0=600000 "
+              "+y_0=0 +ellps=GRS80 "
+              "+step +proj=vgridshift +grids=g2012bu0.gtx +multiplier=1 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, compoundCRS_from_WKT2_to_geogCRS_3D_context) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), "EPSG");