diff --git a/include/proj/io.hpp b/include/proj/io.hpp
index 20dcedfe06..2e796a9d83 100644
--- a/include/proj/io.hpp
+++ b/include/proj/io.hpp
@@ -342,6 +342,10 @@ class PROJ_GCC_DLL WKTFormatter {
PROJ_INTERNAL void ingestWKTNode(const WKTNodeNNPtr &node);
+ PROJ_INTERNAL WKTFormatter &
+ setAllowEllipsoidalHeightAsVerticalCRS(bool allow) noexcept;
+ PROJ_INTERNAL bool isAllowedEllipsoidalHeightAsVerticalCRS() const noexcept;
+
//! @endcond
protected:
diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 90a414c68e..577dae6300 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -1418,6 +1418,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.
+ * 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.
+ * 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).
*
* @return a string, or NULL in case of error.
*/
@@ -1468,6 +1475,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 5d2c6d7bb2..573dd6db0f 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 075729ba30..9876417ba5 100644
--- a/src/iso19111/io.cpp
+++ b/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,32 @@ bool WKTFormatter::isStrict() const noexcept { return d->params_.strict_; }
// ---------------------------------------------------------------------------
+//! @cond Doxygen_Suppress
+
+/** \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_;
+}
+
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
/** Returns the WKT string from the formatter. */
const std::string &WKTFormatter::toString() const {
if (d->indentLevel_ > 0 || d->level_ > 0) {
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 0b45cf8663..7f6e5fb4ea 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 af172b4c0a..7318d1d1c0 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")