From 698e51f476772ebfdd8ba7b93c5a5beafcb90f64 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 22 Oct 2020 19:34:23 +0200 Subject: WKT parser: accept implicit compoundCRS from ESRI WKT, like "PROJCS[...],VERTCS[...]" --- src/iso19111/io.cpp | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) (limited to 'src/iso19111') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 0f4ffba0..503b6be5 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -6507,6 +6507,33 @@ BaseObjectNNPtr WKTParser::createFromWKT(const std::string &wkt) { } return d->buildGeodeticReferenceFrame(root, primeMeridian, null_node); + } else if (ci_equal(name, WKTConstants::GEOGCS) || + ci_equal(name, WKTConstants::PROJCS)) { + // Parse implicit compoundCRS from ESRI that is + // "PROJCS[...],VERTCS[...]" or "GEOGCS[...],VERTCS[...]" + if (indexEnd < wkt.size()) { + indexEnd = skipSpace(wkt, indexEnd); + if (indexEnd < wkt.size() && wkt[indexEnd] == ',') { + ++indexEnd; + indexEnd = skipSpace(wkt, indexEnd); + if (indexEnd < wkt.size() && + ci_starts_with(wkt.c_str() + indexEnd, + WKTConstants::VERTCS.c_str())) { + auto horizCRS = d->buildCRS(root); + if (horizCRS) { + auto vertCRS = + d->buildVerticalCRS(WKTNode::createFrom( + wkt, indexEnd, 0, indexEnd)); + return CompoundCRS::create( + util::PropertyMap().set( + IdentifiedObject::NAME_KEY, + horizCRS->nameStr() + " + " + + vertCRS->nameStr()), + {NN_NO_CHECK(horizCRS), vertCRS}); + } + } + } + } } return d->build(root); }; -- cgit v1.2.3 From f2cd2b3b99746dd9d08b6d388fa91055643f0747 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 22 Oct 2020 21:13:03 +0200 Subject: WKT parser: accept ESRI VERTCS[...,DATUM[...,SPHEROID[]] syntax to express ellipsoidal heights --- src/iso19111/io.cpp | 65 +++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 48 insertions(+), 17 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 503b6be5..b8e835d7 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -4145,11 +4145,19 @@ createBoundCRSSourceTransformationCRS(const crs::CRSPtr &sourceCRS, CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) { const auto *nodeP = node->GP(); - auto &datumNode = + const auto &nodeValue = nodeP->value(); + auto &vdatumNode = nodeP->lookForChild(WKTConstants::VDATUM, WKTConstants::VERT_DATUM, WKTConstants::VERTICALDATUM, WKTConstants::VRF); auto &ensembleNode = nodeP->lookForChild(WKTConstants::ENSEMBLE); - if (isNull(datumNode) && isNull(ensembleNode)) { + // like in ESRI VERTCS["WGS_1984",DATUM["D_WGS_1984", + // SPHEROID["WGS_1984",6378137.0,298.257223563]], + // PARAMETER["Vertical_Shift",0.0], + // PARAMETER["Direction",1.0],UNIT["Meter",1.0] + auto &geogDatumNode = ci_equal(nodeValue, WKTConstants::VERTCS) + ? nodeP->lookForChild(WKTConstants::DATUM) + : null_node; + if (isNull(vdatumNode) && isNull(geogDatumNode) && isNull(ensembleNode)) { throw ParsingException("Missing VDATUM or ENSEMBLE node"); } @@ -4164,28 +4172,50 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) { } auto &dynamicNode = nodeP->lookForChild(WKTConstants::DYNAMIC); - auto datum = - !isNull(datumNode) - ? buildVerticalReferenceFrame(datumNode, dynamicNode).as_nullable() - : nullptr; + auto vdatum = + !isNull(geogDatumNode) + ? VerticalReferenceFrame::create( + PropertyMap() + .set(IdentifiedObject::NAME_KEY, + buildGeodeticReferenceFrame(geogDatumNode, + PrimeMeridian::GREENWICH, + null_node) + ->nameStr()) + .set("VERT_DATUM_TYPE", "2002")) + .as_nullable() + : !isNull(vdatumNode) + ? buildVerticalReferenceFrame(vdatumNode, dynamicNode) + .as_nullable() + : nullptr; auto datumEnsemble = !isNull(ensembleNode) ? buildDatumEnsemble(ensembleNode, nullptr, false).as_nullable() : nullptr; auto &csNode = nodeP->lookForChild(WKTConstants::CS_); - const auto &nodeValue = nodeP->value(); if (isNull(csNode) && !ci_equal(nodeValue, WKTConstants::VERT_CS) && !ci_equal(nodeValue, WKTConstants::VERTCS) && !ci_equal(nodeValue, WKTConstants::BASEVERTCRS)) { ThrowMissing(WKTConstants::CS_); } - auto cs = buildCS(csNode, node, UnitOfMeasure::NONE); - auto verticalCS = nn_dynamic_pointer_cast(cs); + auto verticalCS = nn_dynamic_pointer_cast( + buildCS(csNode, node, UnitOfMeasure::NONE)); if (!verticalCS) { ThrowNotExpectedCSType("vertical"); } + if (vdatum && vdatum->getWKT1DatumType() == "2002" && + &(verticalCS->axisList()[0]->direction()) == &(AxisDirection::UP)) { + verticalCS = + VerticalCS::create( + util::PropertyMap(), + CoordinateSystemAxis::create( + util::PropertyMap().set(IdentifiedObject::NAME_KEY, + "ellipsoidal height"), + "h", AxisDirection::UP, verticalCS->axisList()[0]->unit())) + .as_nullable(); + } + auto &props = buildProperties(node); if (esriStyle_ && dbContext_) { @@ -4246,10 +4276,10 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) { "North American Vertical Datum 1988"); propsDatum.set(Identifier::CODE_KEY, 5103); propsDatum.set(Identifier::CODESPACE_KEY, Identifier::EPSG); - datum = + vdatum = VerticalReferenceFrame::create(propsDatum).as_nullable(); const auto dummyCRS = - VerticalCRS::create(PropertyMap(), datum, datumEnsemble, + VerticalCRS::create(PropertyMap(), vdatum, datumEnsemble, NN_NO_CHECK(verticalCS)); const auto model(Transformation::create( propsModel, dummyCRS, dummyCRS, nullptr, @@ -4265,7 +4295,7 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) { if (!isNull(geoidModelNode)) { auto &propsModel = buildProperties(geoidModelNode); const auto dummyCRS = VerticalCRS::create( - PropertyMap(), datum, datumEnsemble, NN_NO_CHECK(verticalCS)); + PropertyMap(), vdatum, datumEnsemble, NN_NO_CHECK(verticalCS)); const auto model(Transformation::create( propsModel, dummyCRS, dummyCRS, nullptr, OperationMethod::create(PropertyMap(), @@ -4275,10 +4305,10 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) { } auto crs = nn_static_pointer_cast(VerticalCRS::create( - props, datum, datumEnsemble, NN_NO_CHECK(verticalCS))); + props, vdatum, datumEnsemble, NN_NO_CHECK(verticalCS))); - if (!isNull(datumNode)) { - auto &extensionNode = datumNode->lookForChild(WKTConstants::EXTENSION); + if (!isNull(vdatumNode)) { + auto &extensionNode = vdatumNode->lookForChild(WKTConstants::EXTENSION); const auto &extensionChildren = extensionNode->GP()->children(); if (extensionChildren.size() == 2) { if (ci_equal(stripQuotes(extensionChildren[0]), "PROJ4_GRIDS")) { @@ -6524,12 +6554,13 @@ BaseObjectNNPtr WKTParser::createFromWKT(const std::string &wkt) { auto vertCRS = d->buildVerticalCRS(WKTNode::createFrom( wkt, indexEnd, 0, indexEnd)); - return CompoundCRS::create( + return CompoundCRS::createLax( util::PropertyMap().set( IdentifiedObject::NAME_KEY, horizCRS->nameStr() + " + " + vertCRS->nameStr()), - {NN_NO_CHECK(horizCRS), vertCRS}); + {NN_NO_CHECK(horizCRS), vertCRS}, + d->dbContext_); } } } -- cgit v1.2.3 From 82cd218a7ab44eafef92445f2caaec0fb6ec799a Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 22 Oct 2020 21:22:43 +0200 Subject: WKT1_ESRI export: export CompoundCRS as PROJCS[...],VERTCS[...] or GEOGCS[...],VERTCS[...] --- src/iso19111/crs.cpp | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 52b10119..1a17ac67 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -4528,15 +4528,21 @@ CRSNNPtr CompoundCRS::createLax(const util::PropertyMap &properties, //! @cond Doxygen_Suppress void CompoundCRS::_exportToWKT(io::WKTFormatter *formatter) const { const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2; - formatter->startNode(isWKT2 ? io::WKTConstants::COMPOUNDCRS - : io::WKTConstants::COMPD_CS, - !identifiers().empty()); - formatter->addQuotedString(nameStr()); - for (const auto &crs : componentReferenceSystems()) { - crs->_exportToWKT(formatter); + const auto &l_components = componentReferenceSystems(); + if (!isWKT2 && formatter->useESRIDialect() && l_components.size() == 2) { + l_components[0]->_exportToWKT(formatter); + l_components[1]->_exportToWKT(formatter); + } else { + formatter->startNode(isWKT2 ? io::WKTConstants::COMPOUNDCRS + : io::WKTConstants::COMPD_CS, + !identifiers().empty()); + formatter->addQuotedString(nameStr()); + for (const auto &crs : l_components) { + crs->_exportToWKT(formatter); + } + ObjectUsage::baseExportToWKT(formatter); + formatter->endNode(); } - ObjectUsage::baseExportToWKT(formatter); - formatter->endNode(); } //! @endcond -- cgit v1.2.3 From 9a112931555aa324ed9caa3c261a4fd11551d615 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 22 Oct 2020 23:16:58 +0200 Subject: WKT1_ESRI export: generate VERTCS[...,DATUM[...,SPHEROID[]] syntax when ellipsoidal height is found --- src/iso19111/crs.cpp | 31 +++++++++++++++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 1a17ac67..4f9098f2 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -2853,9 +2853,9 @@ void VerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const { !identifiers().empty()); auto l_name = nameStr(); + const auto &dbContext = formatter->databaseContext(); if (formatter->useESRIDialect()) { bool aliasFound = false; - const auto &dbContext = formatter->databaseContext(); if (dbContext) { auto l_alias = dbContext->getAliasFromOfficialName( l_name, "vertical_crs", "ESRI"); @@ -2870,7 +2870,34 @@ void VerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const { } formatter->addQuotedString(l_name); - exportDatumOrDatumEnsembleToWkt(formatter); + + const auto l_datum = datum(); + if (formatter->useESRIDialect() && l_datum && + l_datum->getWKT1DatumType() == "2002") { + bool foundMatch = false; + if (dbContext) { + auto authFactory = io::AuthorityFactory::create( + NN_NO_CHECK(dbContext), std::string()); + auto list = authFactory->createObjectsFromName( + l_datum->nameStr(), + {io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME}, + false /* approximate=false*/); + if (!list.empty()) { + auto gdatum = + util::nn_dynamic_pointer_cast(list.front()); + if (gdatum) { + gdatum->_exportToWKT(formatter); + foundMatch = true; + } + } + } + if (!foundMatch) { + // We should export a geodetic datum, but we cannot really do better + l_datum->_exportToWKT(formatter); + } + } else { + exportDatumOrDatumEnsembleToWkt(formatter); + } const auto &cs = SingleCRS::getPrivate()->coordinateSystem; const auto &axisList = cs->axisList(); -- cgit v1.2.3 From 4969076c15f73371401ee65f2e4617439239cd8b Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 22 Oct 2020 23:20:25 +0200 Subject: Database: import ESRI VERTCS that uses a (geodetic) datum to express ellipsoidal height --- src/iso19111/factory.cpp | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index ef5c6e02..211eb586 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -2100,6 +2100,9 @@ AuthorityFactory::createVerticalDatum(const std::string &code) const { if (!publication_date.empty()) { props.set("PUBLICATION_DATE", publication_date); } + if (d->authority() == "ESRI" && starts_with(code, "from_geogdatum_")) { + props.set("VERT_DATUM_TYPE", "2002"); + } auto anchor = util::optional(); if (frame_reference_epoch.empty()) { return datum::VerticalReferenceFrame::create(props, anchor); @@ -5719,9 +5722,12 @@ AuthorityFactory::createObjectsFromNameEx( "SELECT table_name, auth_name, code, name, deprecated, is_alias " "FROM ("); - const auto getTableAndTypeConstraints = [&allowedObjectTypes]() { + const auto getTableAndTypeConstraints = [&allowedObjectTypes, + &searchedName]() { typedef std::pair TableType; std::list res; + // Hide ESRI D_ vertical datums + const bool startsWithDUnderscore = starts_with(searchedName, "D_"); if (allowedObjectTypes.empty()) { for (const auto &tableName : {"prime_meridian", "ellipsoid", "geodetic_datum", @@ -5729,7 +5735,10 @@ AuthorityFactory::createObjectsFromNameEx( "vertical_crs", "compound_crs", "conversion", "helmert_transformation", "grid_transformation", "other_transformation", "concatenated_operation"}) { - res.emplace_back(TableType(tableName, std::string())); + if (!(startsWithDUnderscore && + strcmp(tableName, "vertical_datum") == 0)) { + res.emplace_back(TableType(tableName, std::string())); + } } } else { for (const auto type : allowedObjectTypes) { @@ -5744,8 +5753,10 @@ AuthorityFactory::createObjectsFromNameEx( case ObjectType::DATUM: res.emplace_back( TableType("geodetic_datum", std::string())); - res.emplace_back( - TableType("vertical_datum", std::string())); + if (!startsWithDUnderscore) { + res.emplace_back( + TableType("vertical_datum", std::string())); + } break; case ObjectType::GEODETIC_REFERENCE_FRAME: res.emplace_back( -- cgit v1.2.3 From 6351422cc2108072162f7b8cdff12916723ccc20 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 23 Oct 2020 00:03:04 +0200 Subject: WKT1_ESRI export: try to export Geographic3D and Projected3D CRS when we can find a corresponding ellipsoidal vertical datum --- src/iso19111/crs.cpp | 69 +++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 66 insertions(+), 3 deletions(-) (limited to 'src/iso19111') diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 4f9098f2..edc8a71f 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -1577,6 +1577,49 @@ GeodeticCRS::create(const util::PropertyMap &properties, // --------------------------------------------------------------------------- +//! @cond Doxygen_Suppress + +// Try to format a Geographic/ProjectedCRS 3D CRS as a +// GEOGCS[]/PROJCS[],VERTCS[...,DATUM[],...] if we find corresponding objects +static bool exportAsESRIWktCompoundCRSWithEllipsoidalHeight( + const CRS *self, const GeodeticCRS *geodCRS, io::WKTFormatter *formatter) { + const auto &dbContext = formatter->databaseContext(); + if (!dbContext) { + return false; + } + const auto l_datum = geodCRS->datumNonNull(formatter->databaseContext()); + auto l_alias = dbContext->getAliasFromOfficialName( + l_datum->nameStr(), "geodetic_datum", "ESRI"); + if (l_alias.empty()) { + return false; + } + auto authFactory = + io::AuthorityFactory::create(NN_NO_CHECK(dbContext), std::string()); + auto list = authFactory->createObjectsFromName( + l_alias, {io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME}, + false /* approximate=false*/); + if (list.empty()) { + return false; + } + auto gdatum = util::nn_dynamic_pointer_cast(list.front()); + if (gdatum == nullptr || gdatum->identifiers().empty()) { + return false; + } + const auto &gdatum_ids = gdatum->identifiers(); + auto vertCRSList = authFactory->createVerticalCRSFromDatum( + "ESRI", "from_geogdatum_" + *gdatum_ids[0]->codeSpace() + '_' + + gdatum_ids[0]->code()); + if (vertCRSList.size() != 1) { + return false; + } + self->demoteTo2D(std::string(), dbContext)->_exportToWKT(formatter); + vertCRSList.front()->_exportToWKT(formatter); + return true; +} +//! @endcond + +// --------------------------------------------------------------------------- + //! @cond Doxygen_Suppress void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const { const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2; @@ -1589,11 +1632,21 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const { auto l_name = nameStr(); const auto &dbContext = formatter->databaseContext(); - if (formatter->useESRIDialect()) { - if (axisList.size() != 2) { + if (!isWKT2 && formatter->useESRIDialect() && axisList.size() == 3) { + if (!isGeographic) { io::FormattingException::Throw( - "Only export of Geographic 2D CRS is supported in WKT1_ESRI"); + "Geocentric CRS not supported in WKT1_ESRI"); } + // Try to format the Geographic 3D CRS as a GEOGCS[],VERTCS[...,DATUM[]] + // if we find corresponding objects + if (dbContext) { + if (exportAsESRIWktCompoundCRSWithEllipsoidalHeight(this, this, + formatter)) { + return; + } + } + io::FormattingException::Throw( + "Cannot export this Geographic 3D CRS in WKT1_ESRI"); } if (!isWKT2 && formatter->isStrict() && isGeographic && @@ -3524,6 +3577,16 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const { } } + if (formatter->useESRIDialect() && dbContext) { + // Try to format the ProjecteD 3D CRS as a + // PROJCS[],VERTCS[...,DATUM[]] + // if we find corresponding objects + if (exportAsESRIWktCompoundCRSWithEllipsoidalHeight( + this, baseCRS().as_nullable().get(), formatter)) { + return; + } + } + if (!formatter->useESRIDialect() && CRS::getPrivate()->allowNonConformantWKT1Export_) { formatter->startNode(io::WKTConstants::COMPD_CS, false); -- cgit v1.2.3