Merge pull request #2386 from rouault/improve_compoundcrs_identification

Improve CompoundCRS identification and name morphing in VerticalCRS with ESRI WKT1 

7 files changed, 247 insertions, 32 deletions

diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp
index a71bf610..b607ff9d 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 e8c97576..37b76346 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<BoundCRS *>(
-                                    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<BoundCRS *>(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<CompoundCRS *>(
+                                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 051f9463..52b10119 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 7d42fd13..5bc8074c 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 7a107f96..99fc6046 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 ec30580a..a0fee905 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<CompoundCRS>(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<CompoundCRS>(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<CompoundCRS>(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 1be0f2aa..055b1e1d 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<VerticalCRS>(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],"