Merge pull request #2450 from rouault/setAllowEllipsoidalHeightAsVerticalCRS

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

3 files changed, 80 insertions, 0 deletions

diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 8fdebcbd..6bc1f166 100644
--- a/src/iso19111/c_api.cpp
+++ b/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;
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 5d2c6d7b..573dd6db 100644
--- a/src/iso19111/crs.cpp
+++ b/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");
     }
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 125105ec..8a42e7ee 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -139,6 +139,7 @@ struct WKTFormatter::Private {
         bool primeMeridianInDegree_ = false;
         bool use2019Keywords_ = false;
         bool useESRIDialect_ = false;
+        bool allowEllipsoidalHeightAsVerticalCRS_ = false;
         OutputAxisRule outputAxis_ = WKTFormatter::OutputAxisRule::YES;
     };
     Params params_{};
@@ -250,6 +251,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;
@@ -263,6 +266,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) {