Skip to content

Commit

Permalink
Merge pull request #2470 from rouault/backport_7_2_pr2450
Browse files Browse the repository at this point in the history 
Add option to allow export of Geographic/Projected 3D CRS in WKT1_GDAL
  • Loading branch information
@rouault
rouault committed Nov 29, 2020
2 parents f7d17d4 + 9bb6e21 commit 0d3c815
Show file tree
Hide file tree
Showing 6 changed files with 177 additions and 2 deletions.
4 changes: 4 additions & 0 deletions include/proj/io.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
12 changes: 12 additions & 0 deletions src/iso19111/c_api.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.</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.
*/
Expand Down Expand Up @@ -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;
Expand Down
43 changes: 43 additions & 0 deletions src/iso19111/crs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

// ---------------------------------------------------------------------------
Expand Down Expand Up @@ -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.");
}
Expand Down Expand Up @@ -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");
}
Expand Down
29 changes: 29 additions & 0 deletions src/iso19111/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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_{};
Expand Down Expand Up @@ -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;
Expand All @@ -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) {
Expand Down
15 changes: 13 additions & 2 deletions test/unit/test_c_api.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down
76 changes: 76 additions & 0 deletions test/unit/test_crs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -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")
Expand Down

0 comments on commit 0d3c815

Please sign in to comment.