From aaa583a01f0c48b7e39511f70fea456f0d0809b1 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 24 Nov 2020 19:16:12 +0100 Subject: [PATCH] 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 | 10 +++++ include/proj/io.hpp | 4 ++ src/apps/projinfo.cpp | 7 ++++ src/iso19111/c_api.cpp | 12 ++++++ src/iso19111/crs.cpp | 43 +++++++++++++++++++++ src/iso19111/io.cpp | 25 ++++++++++++ test/unit/test_c_api.cpp | 15 +++++++- test/unit/test_crs.cpp | 72 +++++++++++++++++++++++++++++++++++ 8 files changed, 186 insertions(+), 2 deletions(-) diff --git a/docs/source/apps/projinfo.rst b/docs/source/apps/projinfo.rst index a07cf70951..803c0a6582 100644 --- a/docs/source/apps/projinfo.rst +++ b/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 diff --git a/include/proj/io.hpp b/include/proj/io.hpp index de8a90fca7..bc11c2204e 100644 --- a/include/proj/io.hpp +++ b/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 : diff --git a/src/apps/projinfo.cpp b/src/apps/projinfo.cpp index 3d57dcd6f0..da885fbb28 100644 --- a/src/apps/projinfo.cpp +++ b/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")) { diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index cbbdcaa8d3..d3c01e8fdd 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. + *
  • 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. */ @@ -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 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 713a471de1..b6a09bb674 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,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) { diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index 9885d221ff..c417371dee 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 e470a6210f..9c81d0cb6a 100644 --- a/test/unit/test_crs.cpp +++ b/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")