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

docs/source/apps/projinfo.rst

@@ -23,6 +23,7 @@ Synopsis
     |    [--grid-check none|discard_missing|sort|known_available]
     |    [--pivot-crs always|if_no_direct_transformation|never|{auth:code[,auth:code]*}]
     |    [--show-superseded] [--hide-ballpark]
+    |    [--allow-ellipsoidal-height-as-vertical-crs]
     |    [--boundcrs-to-wgs84]
     |    [--main-db-path path] [--aux-db-path path]*
     |    [--identify] [--3d]
@@ -212,6 +213,15 @@ The following control parameters can appear in any order:
 
     .. note:: only used for coordinate operation computation
 
+.. option:: --allow-ellipsoidal-height-as-vertical-crs
+
+    .. versionadded:: 8.0
+
+    Allow to export a geographic or projected 3D CRS as a compound CRS whose
+    vertical CRS represents the ellipsoidal height.
+
+    .. note:: only used for CRS, and with WKT1:GDAL output format
+
 .. option:: --boundcrs-to-wgs84
 
     When specified, this option researches a coordinate operation from the

include/proj/io.hpp

@@ -243,6 +243,10 @@ class PROJ_GCC_DLL WKTFormatter {
     PROJ_DLL WKTFormatter &setStrict(bool strict) noexcept;
     PROJ_DLL bool isStrict() const noexcept;
 
+    PROJ_DLL WKTFormatter &
+    setAllowEllipsoidalHeightAsVerticalCRS(bool allow) noexcept;
+    PROJ_DLL bool isAllowedEllipsoidalHeightAsVerticalCRS() const noexcept;
+
     PROJ_DLL const std::string &toString() const;
 
     PROJ_PRIVATE :

src/apps/projinfo.cpp

@@ -72,6 +72,7 @@ struct OutputOptions {
     bool singleLine = false;
     bool strict = true;
     bool ballparkAllowed = true;
+    bool allowEllipsoidalHeightAsVerticalCRS = false;
 };
 } // anonymous namespace
 
@@ -95,6 +96,8 @@ static void usage() {
         << "                [--pivot-crs always|if_no_direct_transformation|"
         << "never|{auth:code[,auth:code]*}]" << std::endl
         << "                [--show-superseded] [--hide-ballpark]" << std::endl
+        << "                [--allow-ellipsoidal-height-as-vertical-crs]"
+        << std::endl
         << "                [--boundcrs-to-wgs84]" << std::endl
         << "                [--main-db-path path] [--aux-db-path path]*"
         << std::endl
@@ -487,6 +490,8 @@ static void outputObject(
                     formatter->setMultiLine(false);
                 }
                 formatter->setStrict(outputOpt.strict);
+                formatter->setAllowEllipsoidalHeightAsVerticalCRS(
+                    outputOpt.allowEllipsoidalHeightAsVerticalCRS);
                 auto wkt = wktExportable->exportToWKT(formatter.get());
                 if (outputOpt.c_ify) {
                     wkt = c_ify_string(wkt);
@@ -1070,6 +1075,8 @@ int main(int argc, char **argv) {
             showSuperseded = true;
         } else if (arg == "--lax") {
             outputOpt.strict = false;
+        } else if (arg == "--allow-ellipsoidal-height-as-vertical-crs") {
+            outputOpt.allowEllipsoidalHeightAsVerticalCRS = true;
         } else if (arg == "--hide-ballpark") {
             outputOpt.ballparkAllowed = false;
         } else if (ci_equal(arg, "--3d")) {

src/iso19111/c_api.cpp

@@ -1421,6 +1421,13 @@ const char *proj_get_id_code(const PJ *obj, int index) {
  * variants, for WKT1_GDAL for ProjectedCRS with easting/northing ordering
  * (otherwise stripped), but not for WKT1_ESRI. Setting to YES will output
  * them unconditionally, and to NO will omit them unconditionally.</li>
+ * <li>STRICT=YES/NO. Default is YES. If NO, a Geographic 3D CRS can be for
+ * example exported as WKT1_GDAL with 3 axes, whereas this is normally not
+ * allowed.</li>
+ * <li>ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS=YES/NO. Default is NO. If set
+ * to YES and type == PJ_WKT1_GDAL, a Geographic 3D CRS or a Projected 3D CRS
+ * will be exported as a compound CRS whose vertical part represents an
+ * ellipsoidal height (for example for use with LAS 1.4 WKT1).</li>
  * </ul>
  * @return a string, or NULL in case of error.
  */
@@ -1471,6 +1478,11 @@ const char *proj_as_wkt(PJ_CONTEXT *ctx, const PJ *obj, PJ_WKT_TYPE type,
                 }
             } else if ((value = getOptionValue(*iter, "STRICT="))) {
                 formatter->setStrict(ci_equal(value, "YES"));
+            } else if ((value = getOptionValue(
+                            *iter,
+                            "ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS="))) {
+                formatter->setAllowEllipsoidalHeightAsVerticalCRS(
+                    ci_equal(value, "YES"));
             } else {
                 std::string msg("Unknown option :");
                 msg += *iter;

src/iso19111/crs.cpp

@@ -1620,6 +1620,34 @@ static bool exportAsESRIWktCompoundCRSWithEllipsoidalHeight(
     vertCRSList.front()->_exportToWKT(formatter);
     return true;
 }
+
+// ---------------------------------------------------------------------------
+
+// Try to format a Geographic/ProjectedCRS 3D CRS as a
+// GEOGCS[]/PROJCS[],VERTCS["Ellipsoid (metre)",DATUM["Ellipsoid",2002],...]
+static bool exportAsWKT1CompoundCRSWithEllipsoidalHeight(
+    const CRSNNPtr &base2DCRS,
+    const cs::CoordinateSystemAxisNNPtr &verticalAxis,
+    io::WKTFormatter *formatter) {
+    std::string verticalCRSName = "Ellipsoid (";
+    verticalCRSName += verticalAxis->unit().name();
+    verticalCRSName += ')';
+    auto vertDatum = datum::VerticalReferenceFrame::create(
+        util::PropertyMap()
+            .set(common::IdentifiedObject::NAME_KEY, "Ellipsoid")
+            .set("VERT_DATUM_TYPE", "2002"));
+    auto vertCRS = VerticalCRS::create(
+        util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
+                                verticalCRSName),
+        vertDatum.as_nullable(), nullptr,
+        cs::VerticalCS::create(util::PropertyMap(), verticalAxis));
+    formatter->startNode(io::WKTConstants::COMPD_CS, false);
+    formatter->addQuotedString(base2DCRS->nameStr() + " + " + verticalCRSName);
+    base2DCRS->_exportToWKT(formatter);
+    vertCRS->_exportToWKT(formatter);
+    formatter->endNode();
+    return true;
+}
 //! @endcond
 
 // ---------------------------------------------------------------------------
@@ -1687,6 +1715,13 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const {
             return;
         }
 
+        if (formatter->isAllowedEllipsoidalHeightAsVerticalCRS()) {
+            if (exportAsWKT1CompoundCRSWithEllipsoidalHeight(
+                    geogCRS2D, axisList[2], formatter)) {
+                return;
+            }
+        }
+
         io::FormattingException::Throw(
             "WKT1 does not support Geographic 3D CRS.");
     }
@@ -3640,6 +3675,14 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const {
             return;
         }
 
+        if (!formatter->useESRIDialect() &&
+            formatter->isAllowedEllipsoidalHeightAsVerticalCRS()) {
+            if (exportAsWKT1CompoundCRSWithEllipsoidalHeight(
+                    projCRS2D, axisList[2], formatter)) {
+                return;
+            }
+        }
+
         io::FormattingException::Throw(
             "Projected 3D CRS can only be exported since WKT2:2019");
     }

src/iso19111/io.cpp

@@ -140,6 +140,7 @@ struct WKTFormatter::Private {
         bool primeMeridianInDegree_ = false;
         bool use2019Keywords_ = false;
         bool useESRIDialect_ = false;
+        bool allowEllipsoidalHeightAsVerticalCRS_ = false;
         OutputAxisRule outputAxis_ = WKTFormatter::OutputAxisRule::YES;
     };
     Params params_{};
@@ -251,6 +252,8 @@ WKTFormatter::setOutputAxis(OutputAxisRule outputAxisIn) noexcept {
  *
  * The default is strict mode, in which case a FormattingException can be
  * thrown.
+ * In non-strict mode, a Geographic 3D CRS can be for example exported as
+ * WKT1_GDAL with 3 axes, whereas this is normally not allowed.
  */
 WKTFormatter &WKTFormatter::setStrict(bool strictIn) noexcept {
     d->params_.strict_ = strictIn;
@@ -264,6 +267,28 @@ bool WKTFormatter::isStrict() const noexcept { return d->params_.strict_; }
 
 // ---------------------------------------------------------------------------
 
+/** \brief Set whether the formatter should export, in WKT1, a Geographic or
+ * Projected 3D CRS as a compound CRS whose vertical part represents an
+ * ellipsoidal height.
+ */
+WKTFormatter &
+WKTFormatter::setAllowEllipsoidalHeightAsVerticalCRS(bool allow) noexcept {
+    d->params_.allowEllipsoidalHeightAsVerticalCRS_ = allow;
+    return *this;
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Return whether the formatter should export, in WKT1, a Geographic or
+ * Projected 3D CRS as a compound CRS whose vertical part represents an
+ * ellipsoidal height.
+ */
+bool WKTFormatter::isAllowedEllipsoidalHeightAsVerticalCRS() const noexcept {
+    return d->params_.allowEllipsoidalHeightAsVerticalCRS_;
+}
+
+// ---------------------------------------------------------------------------
+
 /** Returns the WKT string from the formatter. */
 const std::string &WKTFormatter::toString() const {
     if (d->indentLevel_ > 0 || d->level_ > 0) {

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};

test/unit/test_crs.cpp

@@ -513,6 +513,34 @@ TEST(crs, EPSG_4979_as_WKT1_GDAL) {
 
 // ---------------------------------------------------------------------------
 
+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]]]");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, EPSG_4979_as_WKT1_ESRI) {
     auto crs = GeographicCRS::EPSG_4979;
     WKTFormatterNNPtr f(
@@ -2037,6 +2065,50 @@ TEST(crs, projectedCRS_as_WKT1_ESRI) {
 
 // ---------------------------------------------------------------------------
 
+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]]]");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(crs, projectedCRS_with_ESRI_code_as_WKT1_ESRI) {
     auto dbContext = DatabaseContext::create();
     auto crs = AuthorityFactory::create(dbContext, "ESRI")