From 11012ebcc910c4d95d5acbff43f444715ae1ad32 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 14 May 2020 15:11:02 +0200 Subject: WKT/PROJJSON parsing: fix BoundCRS/PROJ4_GRIDS handling for vertical component when vertical CRS has non metre unit (fixes #2217) --- src/iso19111/io.cpp | 137 +++++++++++++++++++++++++++++++--------------------- 1 file changed, 83 insertions(+), 54 deletions(-) (limited to 'src/iso19111/io.cpp') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 97f741be..7fddd324 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -3987,6 +3987,79 @@ WKTParser::Private::buildParametricDatum(const WKTNodeNNPtr &node) { // --------------------------------------------------------------------------- +static CRSNNPtr +createBoundCRSSourceTransformationCRS(const crs::CRSPtr &sourceCRS, + const crs::CRSPtr &targetCRS) { + CRSPtr sourceTransformationCRS; + if (dynamic_cast(targetCRS.get())) { + GeographicCRSPtr sourceGeographicCRS = + sourceCRS->extractGeographicCRS(); + sourceTransformationCRS = sourceGeographicCRS; + if (sourceGeographicCRS) { + if (sourceGeographicCRS->datum() != nullptr && + sourceGeographicCRS->primeMeridian() + ->longitude() + .getSIValue() != 0.0) { + sourceTransformationCRS = + GeographicCRS::create( + util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + sourceGeographicCRS->nameStr() + + " (with Greenwich prime meridian)"), + datum::GeodeticReferenceFrame::create( + util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + sourceGeographicCRS->datum()->nameStr() + + " (with Greenwich prime meridian)"), + sourceGeographicCRS->datum()->ellipsoid(), + util::optional(), + datum::PrimeMeridian::GREENWICH), + sourceGeographicCRS->coordinateSystem()) + .as_nullable(); + } + } else { + auto vertSourceCRS = + std::dynamic_pointer_cast(sourceCRS); + if (!vertSourceCRS) { + throw ParsingException( + "Cannot find GeographicCRS or VerticalCRS in sourceCRS"); + } + const auto &axis = vertSourceCRS->coordinateSystem()->axisList()[0]; + if (axis->unit() == common::UnitOfMeasure::METRE && + &(axis->direction()) == &AxisDirection::UP) { + sourceTransformationCRS = sourceCRS; + } else { + std::string sourceTransformationCRSName( + vertSourceCRS->nameStr()); + if (ends_with(sourceTransformationCRSName, " (ftUS)")) { + sourceTransformationCRSName.resize( + sourceTransformationCRSName.size() - strlen(" (ftUS)")); + } + if (ends_with(sourceTransformationCRSName, " depth")) { + sourceTransformationCRSName.resize( + sourceTransformationCRSName.size() - strlen(" depth")); + } + if (!ends_with(sourceTransformationCRSName, " height")) { + sourceTransformationCRSName += " height"; + } + sourceTransformationCRS = + VerticalCRS::create( + PropertyMap().set(IdentifiedObject::NAME_KEY, + sourceTransformationCRSName), + vertSourceCRS->datum(), vertSourceCRS->datumEnsemble(), + VerticalCS::createGravityRelatedHeight( + common::UnitOfMeasure::METRE)) + .as_nullable(); + } + } + } else { + sourceTransformationCRS = sourceCRS; + } + return NN_NO_CHECK(sourceTransformationCRS); +} + +// --------------------------------------------------------------------------- + CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) { const auto *nodeP = node->GP(); auto &datumNode = @@ -4111,16 +4184,18 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) { gridName != "g2012a_conus.gtx,g2012a_alaska.gtx," "g2012a_guam.gtx,g2012a_hawaii.gtx," "g2012a_puertorico.gtx,g2012a_samoa.gtx") { - std::string transformationName(crs->nameStr()); - if (!ends_with(transformationName, " height")) { - transformationName += " height"; - } - transformationName += " to WGS84 ellipsoidal height"; + auto sourceTransformationCRS = + createBoundCRSSourceTransformationCRS( + crs.as_nullable(), + GeographicCRS::EPSG_4979.as_nullable()); auto transformation = Transformation:: createGravityRelatedHeightToGeographic3D( - PropertyMap().set(IdentifiedObject::NAME_KEY, - transformationName), - crs, GeographicCRS::EPSG_4979, nullptr, gridName, + PropertyMap().set( + IdentifiedObject::NAME_KEY, + sourceTransformationCRS->nameStr() + + " to WGS84 ellipsoidal height"), + sourceTransformationCRS, GeographicCRS::EPSG_4979, + nullptr, gridName, std::vector()); return nn_static_pointer_cast(BoundCRS::create( crs, GeographicCRS::EPSG_4979, transformation)); @@ -4190,52 +4265,6 @@ CRSNNPtr WKTParser::Private::buildCompoundCRS(const WKTNodeNNPtr &node) { // --------------------------------------------------------------------------- -static CRSNNPtr -createBoundCRSSourceTransformationCRS(const crs::CRSPtr &sourceCRS, - const crs::CRSPtr &targetCRS) { - CRSPtr sourceTransformationCRS; - if (dynamic_cast(targetCRS.get())) { - GeographicCRSPtr sourceGeographicCRS = - sourceCRS->extractGeographicCRS(); - sourceTransformationCRS = sourceGeographicCRS; - if (sourceGeographicCRS) { - if (sourceGeographicCRS->datum() != nullptr && - sourceGeographicCRS->primeMeridian() - ->longitude() - .getSIValue() != 0.0) { - sourceTransformationCRS = - GeographicCRS::create( - util::PropertyMap().set( - common::IdentifiedObject::NAME_KEY, - sourceGeographicCRS->nameStr() + - " (with Greenwich prime meridian)"), - datum::GeodeticReferenceFrame::create( - util::PropertyMap().set( - common::IdentifiedObject::NAME_KEY, - sourceGeographicCRS->datum()->nameStr() + - " (with Greenwich prime meridian)"), - sourceGeographicCRS->datum()->ellipsoid(), - util::optional(), - datum::PrimeMeridian::GREENWICH), - sourceGeographicCRS->coordinateSystem()) - .as_nullable(); - } - } else { - sourceTransformationCRS = - std::dynamic_pointer_cast(sourceCRS); - if (!sourceTransformationCRS) { - throw ParsingException( - "Cannot find GeographicCRS or VerticalCRS in sourceCRS"); - } - } - } else { - sourceTransformationCRS = sourceCRS; - } - return NN_NO_CHECK(sourceTransformationCRS); -} - -// --------------------------------------------------------------------------- - BoundCRSNNPtr WKTParser::Private::buildBoundCRS(const WKTNodeNNPtr &node) { const auto *nodeP = node->GP(); auto &abridgedNode = -- cgit v1.2.3 From 70424a9a17dbfae28dac287831a06f3aee8d2725 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 14 May 2020 16:08:06 +0200 Subject: WKT1 parsing: fix parsing of CompoundCRS that has a TOWGS84 for horizontal and PROJ4_GRIDS for vertical (refs #2217) --- src/iso19111/io.cpp | 35 +++++++++++++++++++++-------------- 1 file changed, 21 insertions(+), 14 deletions(-) (limited to 'src/iso19111/io.cpp') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 7fddd324..af957ada 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -4544,12 +4544,27 @@ CRSPtr WKTParser::Private::buildCRS(const WKTNodeNNPtr &node) { const auto *nodeP = node->GP(); const std::string &name(nodeP->value()); + const auto applyHorizontalBoundCRSParams = [&](const CRSNNPtr &crs) { + if (!toWGS84Parameters_.empty()) { + auto ret = BoundCRS::createFromTOWGS84(crs, toWGS84Parameters_); + toWGS84Parameters_.clear(); + return util::nn_static_pointer_cast(ret); + } else if (!datumPROJ4Grids_.empty()) { + auto ret = BoundCRS::createFromNadgrids(crs, datumPROJ4Grids_); + datumPROJ4Grids_.clear(); + return util::nn_static_pointer_cast(ret); + } + return crs; + }; + if (isGeodeticCRS(name)) { if (!isNull(nodeP->lookForChild(WKTConstants::BASEGEOGCRS, WKTConstants::BASEGEODCRS))) { - return buildDerivedGeodeticCRS(node); + return util::nn_static_pointer_cast( + applyHorizontalBoundCRSParams(buildDerivedGeodeticCRS(node))); } else { - return util::nn_static_pointer_cast(buildGeodeticCRS(node)); + return util::nn_static_pointer_cast( + applyHorizontalBoundCRSParams(buildGeodeticCRS(node))); } } @@ -4574,12 +4589,14 @@ CRSPtr WKTParser::Private::buildCRS(const WKTNodeNNPtr &node) { PROJStringParser().createFromPROJString(projString); auto crs = nn_dynamic_pointer_cast(projObj); if (crs) { - return crs; + return util::nn_static_pointer_cast( + applyHorizontalBoundCRSParams(NN_NO_CHECK(crs))); } } catch (const io::ParsingException &) { } } - return util::nn_static_pointer_cast(buildProjectedCRS(node)); + return util::nn_static_pointer_cast( + applyHorizontalBoundCRSParams(buildProjectedCRS(node))); } if (ci_equal(name, WKTConstants::VERT_CS) || @@ -4652,16 +4669,6 @@ BaseObjectNNPtr WKTParser::Private::build(const WKTNodeNNPtr &node) { auto crs = buildCRS(node); if (crs) { - if (!toWGS84Parameters_.empty()) { - return util::nn_static_pointer_cast( - BoundCRS::createFromTOWGS84(NN_NO_CHECK(crs), - toWGS84Parameters_)); - } - if (!datumPROJ4Grids_.empty()) { - return util::nn_static_pointer_cast( - BoundCRS::createFromNadgrids(NN_NO_CHECK(crs), - datumPROJ4Grids_)); - } return util::nn_static_pointer_cast(NN_NO_CHECK(crs)); } -- cgit v1.2.3 From 69ee50ca36316dd5429b212e7217def4e9ee1dcc Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 14 May 2020 17:04:26 +0200 Subject: Grid listing from a CoordinateOperation: take into account several grid names separated by comma (e.g. '+grids=foo1.gtx,foo2.gtx') --- src/iso19111/io.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) (limited to 'src/iso19111/io.cpp') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index af957ada..9d9fa4a3 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -7630,7 +7630,10 @@ std::set PROJStringFormatter::getUsedGridNames() const { for (const auto &step : d->steps_) { for (const auto ¶m : step.paramValues) { if (param.keyEquals("grids")) { - res.insert(param.value); + const auto gridNames = split(param.value, ","); + for (const auto &gridName : gridNames) { + res.insert(gridName); + } } } } -- cgit v1.2.3