Add option to allow export of Geographic/Projected 3D CRS in WKT1_GDAL

as CompoundCRS with a VerticalCRS being an ellipsoidal height, which is
not conformant. But needed for LAS 1.4 that only supports WKT1

This is a partial backport of https://github.com/OSGeo/PROJ/pull/2450,
with only the new ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS=YES option
of proj_as_wkt()

2 files changed, 89 insertions, 2 deletions

diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 0b45cf86..7f6e5fb4 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -433,16 +433,27 @@ TEST_F(CApi, proj_as_wkt) {
     ObjectKeeper keeper_crs4979(crs4979);
     ASSERT_NE(crs4979, nullptr);
 
+    EXPECT_EQ(proj_as_wkt(m_ctxt, crs4979, PJ_WKT1_GDAL, nullptr), nullptr);
+
     // STRICT=NO
     {
-        EXPECT_EQ(proj_as_wkt(m_ctxt, crs4979, PJ_WKT1_GDAL, nullptr), nullptr);
-
         const char *const options[] = {"STRICT=NO", nullptr};
         auto wkt = proj_as_wkt(m_ctxt, crs4979, PJ_WKT1_GDAL, options);
         ASSERT_NE(wkt, nullptr);
         EXPECT_TRUE(std::string(wkt).find("GEOGCS[\"WGS 84\"") == 0) << wkt;
     }
 
+    // ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS=YES
+    {
+        const char *const options[] = {
+            "ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS=YES", nullptr};
+        auto wkt = proj_as_wkt(m_ctxt, crs4979, PJ_WKT1_GDAL, options);
+        ASSERT_NE(wkt, nullptr);
+        EXPECT_TRUE(std::string(wkt).find(
+                        "COMPD_CS[\"WGS 84 + Ellipsoid (metre)\"") == 0)
+            << wkt;
+    }
+
     // unsupported option
     {
         const char *const options[] = {"unsupported=yes", nullptr};
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index af172b4c..7318d1d1 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -527,6 +527,36 @@ TEST(crs, EPSG_4979_as_WKT1_GDAL) {
 
 // ---------------------------------------------------------------------------
 
+#ifdef notavailable_since_setAllowEllipsoidalHeightAsVerticalCRS_is_internal
+TEST(crs, EPSG_4979_as_WKT1_GDAL_with_ellipsoidal_height_as_vertical_crs) {
+    auto crs = GeographicCRS::EPSG_4979;
+    auto wkt = crs->exportToWKT(
+        &(WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL,
+                               DatabaseContext::create())
+              ->setAllowEllipsoidalHeightAsVerticalCRS(true)));
+
+    // For LAS 1.4 WKT1...
+    EXPECT_EQ(wkt, "COMPD_CS[\"WGS 84 + Ellipsoid (metre)\",\n"
+                   "    GEOGCS[\"WGS 84\",\n"
+                   "        DATUM[\"WGS_1984\",\n"
+                   "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+                   "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+                   "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+                   "        PRIMEM[\"Greenwich\",0,\n"
+                   "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+                   "        UNIT[\"degree\",0.0174532925199433,\n"
+                   "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
+                   "        AUTHORITY[\"EPSG\",\"4326\"]],\n"
+                   "    VERT_CS[\"Ellipsoid (metre)\",\n"
+                   "        VERT_DATUM[\"Ellipsoid\",2002],\n"
+                   "        UNIT[\"metre\",1,\n"
+                   "            AUTHORITY[\"EPSG\",\"9001\"]],\n"
+                   "        AXIS[\"Ellipsoidal height\",UP]]]");
+}
+#endif
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, EPSG_4979_as_WKT1_ESRI) {
     auto crs = GeographicCRS::EPSG_4979;
     WKTFormatterNNPtr f(
@@ -2051,6 +2081,52 @@ TEST(crs, projectedCRS_as_WKT1_ESRI) {
 
 // ---------------------------------------------------------------------------
 
+#ifdef notavailable_since_setAllowEllipsoidalHeightAsVerticalCRS_is_internal
+TEST(crs,
+     projectedCRS_3D_as_WKT1_GDAL_with_ellipsoidal_height_as_vertical_crs) {
+    auto dbContext = DatabaseContext::create();
+    auto crs = AuthorityFactory::create(dbContext, "EPSG")
+                   ->createProjectedCRS("32631")
+                   ->promoteTo3D(std::string(), dbContext);
+    auto wkt = crs->exportToWKT(
+        &(WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL, dbContext)
+              ->setAllowEllipsoidalHeightAsVerticalCRS(true)));
+
+    // For LAS 1.4 WKT1...
+    EXPECT_EQ(wkt,
+              "COMPD_CS[\"WGS 84 / UTM zone 31N + Ellipsoid (metre)\",\n"
+              "    PROJCS[\"WGS 84 / UTM zone 31N\",\n"
+              "        GEOGCS[\"WGS 84\",\n"
+              "            DATUM[\"WGS_1984\",\n"
+              "                SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+              "                    AUTHORITY[\"EPSG\",\"7030\"]],\n"
+              "                AUTHORITY[\"EPSG\",\"6326\"]],\n"
+              "            PRIMEM[\"Greenwich\",0,\n"
+              "                AUTHORITY[\"EPSG\",\"8901\"]],\n"
+              "            UNIT[\"degree\",0.0174532925199433,\n"
+              "                AUTHORITY[\"EPSG\",\"9122\"]],\n"
+              "            AUTHORITY[\"EPSG\",\"4326\"]],\n"
+              "        PROJECTION[\"Transverse_Mercator\"],\n"
+              "        PARAMETER[\"latitude_of_origin\",0],\n"
+              "        PARAMETER[\"central_meridian\",3],\n"
+              "        PARAMETER[\"scale_factor\",0.9996],\n"
+              "        PARAMETER[\"false_easting\",500000],\n"
+              "        PARAMETER[\"false_northing\",0],\n"
+              "        UNIT[\"metre\",1,\n"
+              "            AUTHORITY[\"EPSG\",\"9001\"]],\n"
+              "        AXIS[\"Easting\",EAST],\n"
+              "        AXIS[\"Northing\",NORTH],\n"
+              "        AUTHORITY[\"EPSG\",\"32631\"]],\n"
+              "    VERT_CS[\"Ellipsoid (metre)\",\n"
+              "        VERT_DATUM[\"Ellipsoid\",2002],\n"
+              "        UNIT[\"metre\",1,\n"
+              "            AUTHORITY[\"EPSG\",\"9001\"]],\n"
+              "        AXIS[\"Ellipsoidal height\",UP]]]");
+}
+#endif
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, projectedCRS_with_ESRI_code_as_WKT1_ESRI) {
     auto dbContext = DatabaseContext::create();
     auto crs = AuthorityFactory::create(dbContext, "ESRI")