Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 1 commit into from
Nov 29, 2020
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions include/proj/io.hpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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