diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp index a71bf610cc..b607ff9d8f 100644 --- a/include/proj/crs.hpp +++ b/include/proj/crs.hpp @@ -154,6 +154,7 @@ class PROJ_GCC_DLL CRS : public common::ObjectUsage, const cs::CoordinateSystemAxisNNPtr &verticalAxisIfNotAlreadyPresent) const; + PROJ_INTERNAL bool hasImplicitCS() const; //! @endcond protected: diff --git a/src/apps/projinfo.cpp b/src/apps/projinfo.cpp index e8c9757615..37b76346dd 100644 --- a/src/apps/projinfo.cpp +++ b/src/apps/projinfo.cpp @@ -1198,26 +1198,48 @@ int main(int argc, char **argv) { std::cout << *ids[0]->codeSpace() << ":" << ids[0]->code() << ": " << pair.second << " %" << std::endl; - } else { - auto boundCRS = dynamic_cast( - identifiedCRS.get()); - if (boundCRS && - !boundCRS->baseCRS() - ->identifiers() - .empty()) { - const auto &idsBase = - boundCRS->baseCRS()->identifiers(); - std::cout << "BoundCRS of " - << *idsBase[0]->codeSpace() << ":" - << idsBase[0]->code() << ": " - << pair.second << " %" - << std::endl; - } else { - std::cout - << "un-identifier CRS: " << pair.second - << " %" << std::endl; + continue; + } + + auto boundCRS = + dynamic_cast(identifiedCRS.get()); + if (boundCRS && + !boundCRS->baseCRS()->identifiers().empty()) { + const auto &idsBase = + boundCRS->baseCRS()->identifiers(); + std::cout << "BoundCRS of " + << *idsBase[0]->codeSpace() << ":" + << idsBase[0]->code() << ": " + << pair.second << " %" << std::endl; + continue; + } + + auto compoundCRS = dynamic_cast( + identifiedCRS.get()); + if (compoundCRS) { + const auto &components = + compoundCRS->componentReferenceSystems(); + if (components.size() == 2 && + !components[0]->identifiers().empty() && + !components[1]->identifiers().empty()) { + const auto &idH = + components[0]->identifiers().front(); + const auto &idV = + components[1]->identifiers().front(); + if (*idH->codeSpace() == + *idV->codeSpace()) { + std::cout << *idH->codeSpace() << ":" + << idH->code() << '+' + << idV->code() << ": " + << pair.second << " %" + << std::endl; + continue; + } } } + + std::cout << "un-identified CRS: " << pair.second + << " %" << std::endl; } } catch (const std::exception &e) { std::cerr << "Identification failed: " << e.what() diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 051f9463f1..52b101199a 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -132,6 +132,16 @@ CRS::~CRS() = default; // --------------------------------------------------------------------------- +//! @cond Doxygen_Suppress + +/** \brief Return whether the CRS has an implicit coordinate system + * (e.g from ESRI WKT) */ +bool CRS::hasImplicitCS() const { return d->implicitCS_; } + +//! @endcond + +// --------------------------------------------------------------------------- + /** \brief Return the BoundCRS potentially attached to this CRS. * * In the case this method is called on a object returned by @@ -1971,7 +1981,7 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { io::DatabaseContextPtr dbContext = authorityFactory ? authorityFactory->databaseContext().as_nullable() : nullptr; - const bool l_implicitCS = CRS::getPrivate()->implicitCS_; + const bool l_implicitCS = hasImplicitCS(); const auto crsCriterion = l_implicitCS ? util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS @@ -2841,7 +2851,25 @@ void VerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const { ? io::WKTConstants::VERTCS : io::WKTConstants::VERT_CS, !identifiers().empty()); - formatter->addQuotedString(nameStr()); + + auto l_name = nameStr(); + if (formatter->useESRIDialect()) { + bool aliasFound = false; + const auto &dbContext = formatter->databaseContext(); + if (dbContext) { + auto l_alias = dbContext->getAliasFromOfficialName( + l_name, "vertical_crs", "ESRI"); + if (!l_alias.empty()) { + l_name = l_alias; + aliasFound = true; + } + } + if (!aliasFound) { + l_name = io::WKTFormatter::morphNameToESRI(l_name); + } + } + + formatter->addQuotedString(l_name); exportDatumOrDatumEnsembleToWkt(formatter); const auto &cs = SingleCRS::getPrivate()->coordinateSystem; const auto &axisList = cs->axisList(); @@ -4020,7 +4048,7 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { } } - const bool l_implicitCS = CRS::getPrivate()->implicitCS_; + const bool l_implicitCS = hasImplicitCS(); const auto addCRS = [&](const ProjectedCRSNNPtr &crs, const bool eqName) { const auto &l_unit = cs->axisList()[0]->unit(); if (_isEquivalentTo(crs.get(), util::IComparable::Criterion:: @@ -4617,6 +4645,13 @@ CompoundCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { const auto &thisName(nameStr()); + const auto &components = componentReferenceSystems(); + const bool l_implicitCS = components[0]->hasImplicitCS(); + const auto crsCriterion = + l_implicitCS + ? util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS + : util::IComparable::Criterion::EQUIVALENT; + if (authorityFactory) { const io::DatabaseContextNNPtr &dbContext = authorityFactory->databaseContext(); @@ -4635,9 +4670,8 @@ CompoundCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { auto crs = io::AuthorityFactory::create( dbContext, *id->codeSpace()) ->createCompoundCRS(id->code()); - bool match = _isEquivalentTo( - crs.get(), util::IComparable::Criterion::EQUIVALENT, - dbContext); + bool match = + _isEquivalentTo(crs.get(), crsCriterion, dbContext); res.emplace_back(crs, match ? 100 : 25); return res; } catch (const std::exception &) { @@ -4657,9 +4691,7 @@ CompoundCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { const bool eqName = metadata::Identifier::isEquivalentName( thisName.c_str(), crs->nameStr().c_str()); foundEquivalentName |= eqName; - if (_isEquivalentTo( - crs.get(), util::IComparable::Criterion::EQUIVALENT, - dbContext)) { + if (_isEquivalentTo(crs.get(), crsCriterion, dbContext)) { if (crs->nameStr() == thisName) { res.clear(); res.emplace_back(crsNN, 100); @@ -4667,6 +4699,7 @@ CompoundCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { } res.emplace_back(crsNN, eqName ? 90 : 70); } else { + res.emplace_back(crsNN, 25); } } @@ -4725,9 +4758,7 @@ CompoundCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { continue; } - if (_isEquivalentTo(crs.get(), - util::IComparable::Criterion::EQUIVALENT, - dbContext)) { + if (_isEquivalentTo(crs.get(), crsCriterion, dbContext)) { res.emplace_back(crs, unsignificantName ? 90 : 70); } else { res.emplace_back(crs, 25); @@ -4737,6 +4768,33 @@ CompoundCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { res.sort(lambdaSort); } + // If we didn't find a match for the CompoundCRS, check if the + // horizontal and vertical parts are not themselves well known. + if (identifiers().empty() && res.empty() && components.size() == 2) { + auto candidatesHorizCRS = components[0]->identify(authorityFactory); + auto candidatesVertCRS = components[1]->identify(authorityFactory); + if (candidatesHorizCRS.size() == 1 && + candidatesVertCRS.size() == 1 && + candidatesHorizCRS.front().second >= 70 && + candidatesVertCRS.front().second >= 70) { + auto newCRS = CompoundCRS::create( + util::PropertyMap().set( + common::IdentifiedObject::NAME_KEY, + candidatesHorizCRS.front().first->nameStr() + " + " + + candidatesVertCRS.front().first->nameStr()), + {candidatesHorizCRS.front().first, + candidatesVertCRS.front().first}); + const bool eqName = metadata::Identifier::isEquivalentName( + thisName.c_str(), newCRS->nameStr().c_str()); + res.emplace_back( + newCRS, + std::min(thisName == newCRS->nameStr() ? 100 : eqName ? 90 + : 70, + std::min(candidatesHorizCRS.front().second, + candidatesVertCRS.front().second))); + } + } + // Keep only results of the highest confidence if (res.size() >= 2) { const auto highestConfidence = res.front().second; diff --git a/src/iso19111/datum.cpp b/src/iso19111/datum.cpp index 7d42fd1389..5bc8074cbd 100644 --- a/src/iso19111/datum.cpp +++ b/src/iso19111/datum.cpp @@ -1933,8 +1933,23 @@ void VerticalReferenceFrame::_exportToWKT( ? io::WKTConstants::VDATUM : io::WKTConstants::VERT_DATUM, !identifiers().empty()); - const auto &l_name = nameStr(); + auto l_name = nameStr(); if (!l_name.empty()) { + if (!isWKT2 && formatter->useESRIDialect()) { + bool aliasFound = false; + const auto &dbContext = formatter->databaseContext(); + if (dbContext) { + auto l_alias = dbContext->getAliasFromOfficialName( + l_name, "vertical_datum", "ESRI"); + if (!l_alias.empty()) { + l_name = l_alias; + aliasFound = true; + } + } + if (!aliasFound) { + l_name = io::WKTFormatter::morphNameToESRI(l_name); + } + } formatter->addQuotedString(l_name); } else { formatter->addQuotedString("unnamed"); diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 7a107f96c1..99fc6046de 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -4006,6 +4006,22 @@ VerticalReferenceFrameNNPtr WKTParser::Private::buildVerticalReferenceFrame( const auto *nodeP = node->GP(); const std::string &name(nodeP->value()); auto &props = buildProperties(node); + + if (esriStyle_ && dbContext_) { + std::string outTableName; + std::string authNameFromAlias; + std::string codeFromAlias; + auto authFactory = + AuthorityFactory::create(NN_NO_CHECK(dbContext_), std::string()); + const std::string datumName = stripQuotes(nodeP->children()[0]); + auto officialName = authFactory->getOfficialNameFromAlias( + datumName, "vertical_datum", "ESRI", false, outTableName, + authNameFromAlias, codeFromAlias); + if (!officialName.empty()) { + props.set(IdentifiedObject::NAME_KEY, officialName); + } + } + if (ci_equal(name, WKTConstants::VERT_DATUM)) { const auto &children = nodeP->children(); if (children.size() >= 2) { @@ -4137,6 +4153,16 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) { throw ParsingException("Missing VDATUM or ENSEMBLE node"); } + for (const auto &childNode : nodeP->children()) { + const auto &childNodeChildren = childNode->GP()->children(); + if (childNodeChildren.size() == 2 && + ci_equal(childNode->GP()->value(), WKTConstants::PARAMETER) && + childNodeChildren[0]->GP()->value() == "\"Vertical_Shift\"") { + esriStyle_ = true; + break; + } + } + auto &dynamicNode = nodeP->lookForChild(WKTConstants::DYNAMIC); auto datum = !isNull(datumNode) @@ -4162,6 +4188,21 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) { auto &props = buildProperties(node); + if (esriStyle_ && dbContext_) { + std::string outTableName; + std::string authNameFromAlias; + std::string codeFromAlias; + auto authFactory = + AuthorityFactory::create(NN_NO_CHECK(dbContext_), std::string()); + const std::string vertCRSName = stripQuotes(nodeP->children()[0]); + auto officialName = authFactory->getOfficialNameFromAlias( + vertCRSName, "vertical_crs", "ESRI", false, outTableName, + authNameFromAlias, codeFromAlias); + if (!officialName.empty()) { + props.set(IdentifiedObject::NAME_KEY, officialName); + } + } + // Deal with Lidar WKT1 VertCRS that embeds geoid model in CRS name, // following conventions from // https://pubs.usgs.gov/tm/11b4/pdf/tm11-B4.pdf diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp index ec30580a0c..a0fee905b5 100644 --- a/test/unit/test_crs.cpp +++ b/test/unit/test_crs.cpp @@ -3593,7 +3593,7 @@ TEST(crs, verticalCRS_as_WKT1_GDAL) { TEST(crs, verticalCRS_as_WKT1_ESRI) { auto crs = createVerticalCRS(); - auto expected = "VERTCS[\"ODN height\",VDATUM[\"Ordnance Datum Newlyn\"]," + auto expected = "VERTCS[\"ODN_height\",VDATUM[\"Ordnance_Datum_Newlyn\"]," "PARAMETER[\"Vertical_Shift\",0.0]," "PARAMETER[\"Direction\",1.0]," "UNIT[\"Meter\",1.0]]"; @@ -3603,6 +3603,23 @@ TEST(crs, verticalCRS_as_WKT1_ESRI) { WKTFormatter::create(WKTFormatter::Convention::WKT1_ESRI).get()), expected); } + +// --------------------------------------------------------------------------- + +TEST(crs, verticalCRS_as_WKT1_ESRI_context) { + auto crs = createVerticalCRS(); + auto expected = "VERTCS[\"Newlyn\",VDATUM[\"Ordnance_Datum_Newlyn\"]," + "PARAMETER[\"Vertical_Shift\",0.0]," + "PARAMETER[\"Direction\",1.0]," + "UNIT[\"Meter\",1.0]]"; + + EXPECT_EQ(crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_ESRI, + DatabaseContext::create()) + .get()), + expected); +} + // --------------------------------------------------------------------------- TEST(crs, verticalCRS_down_as_WKT1_ESRI) { @@ -4058,6 +4075,42 @@ TEST(crs, compoundCRS_identify_db) { // Just check we don't get an exception crs->identify(factory); } + // Identify from ESRI WKT + { + auto sourceCRS = factory->createCompoundCRS("7405"); + auto wkt = sourceCRS->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_ESRI, dbContext) + .get()); + auto obj = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + auto res = crs->identify(factory); + ASSERT_EQ(res.size(), 1U); + EXPECT_EQ(res.front().first->getEPSGCode(), 7405); + EXPECT_EQ(res.front().second, 100); + } + // Identify a CompoundCRS whose horizontal and vertical parts are known + // but not the composition. + { + auto obj = createFromUserInput("EPSG:4326+3855", dbContext); + auto sourceCRS = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(sourceCRS != nullptr); + auto wkt = sourceCRS->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_ESRI, dbContext) + .get()); + auto obj2 = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj2); + ASSERT_TRUE(crs != nullptr); + auto res = crs->identify(factory); + ASSERT_EQ(res.size(), 1U); + EXPECT_EQ(res.front().first->getEPSGCode(), 0); + EXPECT_EQ(res.front().second, 100); + const auto &components = res.front().first->componentReferenceSystems(); + EXPECT_EQ(components[0]->getEPSGCode(), 4326); + EXPECT_EQ(components[1]->getEPSGCode(), 3855); + } } // --------------------------------------------------------------------------- diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 1be0f2aad3..055b1e1d1f 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -2106,6 +2106,31 @@ TEST(wkt_parse, VERTCS_WKT1_ESRI) { // --------------------------------------------------------------------------- +TEST(wkt_parse, VERTCS_WKT1_ESRI_context) { + auto wkt = "VERTCS[\"EGM2008_Geoid\",VDATUM[\"EGM2008_Geoid\"]," + "PARAMETER[\"Vertical_Shift\",0.0]," + "PARAMETER[\"Direction\",1.0],UNIT[\"Meter\",1.0]]"; + + auto obj = WKTParser() + .attachDatabaseContext(DatabaseContext::create()) + .createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ(crs->nameStr(), "EGM2008 height"); + + auto datum = crs->datum(); + EXPECT_EQ(datum->nameStr(), "EGM2008 geoid"); + + auto cs = crs->coordinateSystem(); + ASSERT_EQ(cs->axisList().size(), 1U); + EXPECT_EQ(cs->axisList()[0]->direction(), AxisDirection::UP); + + EXPECT_EQ(WKTParser().guessDialect(wkt), + WKTParser::WKTGuessedDialect::WKT1_ESRI); +} + +// --------------------------------------------------------------------------- + TEST(wkt_parse, VERTCS_WKT1_ESRI_down) { auto wkt = "VERTCS[\"Caspian\",VDATUM[\"Caspian_Sea\"]," "PARAMETER[\"Vertical_Shift\",0.0],"