improve identify() for projected and bound CRS

3 files changed, 204 insertions, 20 deletions

diff --git a/src/crs.cpp b/src/crs.cpp
index 1cc5d31a..21df43b4 100644
--- a/src/crs.cpp
+++ b/src/crs.cpp
@@ -45,6 +45,7 @@
 #include "proj/internal/internal.hpp"
 #include "proj/internal/io_internal.hpp"
 
+#include <algorithm>
 #include <cassert>
 #include <cstring>
 #include <iostream>
@@ -88,6 +89,20 @@ namespace crs {
 struct CRS::Private {
     BoundCRSPtr canonicalBoundCRS_{};
     std::string extensionProj4_{};
+    bool implicitCS_ = false;
+
+    void setImplicitCS(const util::PropertyMap &properties) {
+        auto oIter = properties.find("IMPLICIT_CS");
+        if (oIter != properties.end()) {
+            if (auto genVal = util::nn_dynamic_pointer_cast<util::BoxedValue>(
+                    oIter->second)) {
+                if (genVal->type() == util::BoxedValue::Type::BOOLEAN &&
+                    genVal->booleanValue()) {
+                    implicitCS_ = true;
+                }
+            }
+        }
+    }
 };
 //! @endcond
 
@@ -860,6 +875,7 @@ GeodeticCRS::create(const util::PropertyMap &properties,
     crs->setProperties(properties);
     properties.getStringValue("EXTENSION_PROJ4",
                               crs->CRS::getPrivate()->extensionProj4_);
+
     return crs;
 }
 
@@ -1183,6 +1199,12 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
     std::list<Pair> res;
     const auto &thisName(nameStr());
 
+    const bool l_implicitCS = CRS::getPrivate()->implicitCS_;
+    const auto crsCriterion =
+        l_implicitCS
+            ? util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS
+            : util::IComparable::Criterion::EQUIVALENT;
+
     const GeographicCRSNNPtr candidatesCRS[] = {GeographicCRS::EPSG_4326,
                                                 GeographicCRS::EPSG_4267,
                                                 GeographicCRS::EPSG_4269};
@@ -1190,8 +1212,7 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
         const bool nameEquivalent = metadata::Identifier::isEquivalentName(
             thisName.c_str(), crs->nameStr().c_str());
         const bool nameEqual = thisName == crs->nameStr();
-        const bool isEq = _isEquivalentTo(
-            crs.get(), util::IComparable::Criterion::EQUIVALENT);
+        const bool isEq = _isEquivalentTo(crs.get(), crsCriterion);
         if (nameEquivalent && isEq && (!authorityFactory || nameEqual)) {
             res.emplace_back(util::nn_static_pointer_cast<GeodeticCRS>(crs),
                              nameEqual ? 100 : 90);
@@ -1226,15 +1247,13 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
         const auto &thisDatum(datum());
 
         auto searchByDatum = [this, &authorityFactory, &res, &thisDatum,
-                              &geodetic_crs_type]() {
+                              &geodetic_crs_type, crsCriterion]() {
             for (const auto &id : thisDatum->identifiers()) {
                 try {
                     auto tempRes = authorityFactory->createGeodeticCRSFromDatum(
                         *id->codeSpace(), id->code(), geodetic_crs_type);
                     for (const auto &crs : tempRes) {
-                        if (_isEquivalentTo(
-                                crs.get(),
-                                util::IComparable::Criterion::EQUIVALENT)) {
+                        if (_isEquivalentTo(crs.get(), crsCriterion)) {
                             res.emplace_back(crs, 70);
                         }
                     }
@@ -1245,7 +1264,8 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
 
         const auto &thisEllipsoid(ellipsoid());
         auto searchByEllipsoid = [this, &authorityFactory, &res, &thisDatum,
-                                  &thisEllipsoid, &geodetic_crs_type]() {
+                                  &thisEllipsoid, &geodetic_crs_type,
+                                  l_implicitCS]() {
             const auto ellipsoids =
                 thisEllipsoid->identifiers().empty()
                     ? authorityFactory->createEllipsoidFromExisting(
@@ -1267,11 +1287,11 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                                 crsDatum->primeMeridian()->_isEquivalentTo(
                                     thisDatum->primeMeridian().get(),
                                     util::IComparable::Criterion::EQUIVALENT) &&
-                                coordinateSystem()->_isEquivalentTo(
-                                    crs->coordinateSystem().get(),
-                                    util::IComparable::Criterion::EQUIVALENT)
-
-                                    ) {
+                                (!l_implicitCS ||
+                                 coordinateSystem()->_isEquivalentTo(
+                                     crs->coordinateSystem().get(),
+                                     util::IComparable::Criterion::
+                                         EQUIVALENT))) {
                                 res.emplace_back(crs, 60);
                             }
                         }
@@ -1302,8 +1322,7 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                                    authorityFactory->databaseContext(),
                                    *id->codeSpace())
                                    ->createGeodeticCRS(id->code());
-                    bool match = _isEquivalentTo(
-                        crs.get(), util::IComparable::Criterion::EQUIVALENT);
+                    bool match = _isEquivalentTo(crs.get(), crsCriterion);
                     res.emplace_back(crs, match ? 100 : 25);
                     return res;
                 } catch (const std::exception &) {
@@ -1319,9 +1338,7 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                     auto crs = util::nn_dynamic_pointer_cast<GeodeticCRS>(obj);
                     assert(crs);
                     auto crsNN = NN_NO_CHECK(crs);
-                    if (_isEquivalentTo(
-                            crs.get(),
-                            util::IComparable::Criterion::EQUIVALENT)) {
+                    if (_isEquivalentTo(crs.get(), crsCriterion)) {
                         if (crs->nameStr() == thisName) {
                             res.clear();
                             res.emplace_back(crsNN, 100);
@@ -1567,6 +1584,7 @@ GeographicCRS::create(const util::PropertyMap &properties,
     crs->setProperties(properties);
     properties.getStringValue("EXTENSION_PROJ4",
                               crs->CRS::getPrivate()->extensionProj4_);
+    crs->CRS::getPrivate()->setImplicitCS(properties);
     return crs;
 }
 
@@ -2582,6 +2600,7 @@ ProjectedCRS::create(const util::PropertyMap &properties,
     crs->setDerivingConversionCRS();
     properties.getStringValue("EXTENSION_PROJ4",
                               crs->CRS::getPrivate()->extensionProj4_);
+    crs->CRS::getPrivate()->setImplicitCS(properties);
     return crs;
 }
 
@@ -2836,7 +2855,8 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                                    *id->codeSpace())
                                    ->createProjectedCRS(id->code());
                     bool match = _isEquivalentTo(
-                        crs.get(), util::IComparable::Criterion::EQUIVALENT);
+                        crs.get(), util::IComparable::Criterion::
+                                       EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS);
                     res.emplace_back(crs, match ? 100 : 25);
                     return res;
                 } catch (const std::exception &) {
@@ -2857,13 +2877,27 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                     foundEquivalentName |= eqName;
                     if (_isEquivalentTo(
                             crs.get(),
-                            util::IComparable::Criterion::EQUIVALENT)) {
+                            util::IComparable::Criterion::
+                                EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS)) {
                         if (crs->nameStr() == thisName) {
                             res.clear();
                             res.emplace_back(crsNN, 100);
                             return res;
                         }
                         res.emplace_back(crsNN, eqName ? 90 : 70);
+                    } else if (crs->nameStr() == thisName &&
+                               CRS::getPrivate()->implicitCS_ &&
+                               l_baseCRS->_isEquivalentTo(
+                                   crs->baseCRS().get(),
+                                   util::IComparable::Criterion::
+                                       EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS) &&
+                               derivingConversionRef()->_isEquivalentTo(
+                                   crs->derivingConversionRef().get(),
+                                   util::IComparable::Criterion::EQUIVALENT) &&
+                               objects.size() == 1) {
+                        res.clear();
+                        res.emplace_back(crsNN, 100);
+                        return res;
                     } else {
                         res.emplace_back(crsNN, 25);
                     }
@@ -2925,7 +2959,8 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                 }
 
                 if (_isEquivalentTo(crs.get(),
-                                    util::IComparable::Criterion::EQUIVALENT)) {
+                                    util::IComparable::Criterion::
+                                        EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS)) {
                     res.emplace_back(crs, unsignificantName ? 90 : 70);
                 } else if (ellipsoid->_isEquivalentTo(
                                crs->baseCRS()->ellipsoid().get(),
@@ -3663,9 +3698,26 @@ BoundCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                     }
                 } catch (const std::exception &) {
                 }
+                bool foundOp = false;
                 for (const auto &op : ops) {
                     std::string opTransfPROJString;
                     bool opTransfPROJStringValid = false;
+                    if (op->nameStr().find("Null geographic") == 0) {
+                        if (isTOWGS84Compatible()) {
+                            auto params =
+                                transformation()->getTOWGS84Parameters();
+                            if (params ==
+                                std::vector<double>{0, 0, 0, 0, 0, 0, 0}) {
+                                res.emplace_back(create(candidateBaseCRS,
+                                                        d->hubCRS_,
+                                                        transformation()),
+                                                 pair.second);
+                                foundOp = true;
+                                break;
+                            }
+                        }
+                        continue;
+                    }
                     try {
                         opTransfPROJString = op->exportToPROJString(
                             io::PROJStringFormatter::create().get());
@@ -3682,9 +3734,15 @@ BoundCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                                    NN_NO_CHECK(util::nn_dynamic_pointer_cast<
                                                operation::Transformation>(op))),
                             pair.second);
+                        foundOp = true;
                         break;
                     }
                 }
+                if (!foundOp) {
+                    res.emplace_back(
+                        create(candidateBaseCRS, d->hubCRS_, transformation()),
+                        std::min(70, pair.second));
+                }
             }
         }
     }
diff --git a/src/io.cpp b/src/io.cpp
index 691d96e3..f95a8079 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -2517,6 +2517,11 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) {
     auto props = buildProperties(node);
     addExtensionProj4ToProp(nodeP, props);
 
+    // No explicit AXIS node ? (WKT1)
+    if (isNull(nodeP->lookForChild(WKTConstants::AXIS))) {
+        props.set("IMPLICIT_CS", true);
+    }
+
     auto datum =
         !isNull(datumNode)
             ? buildGeodeticReferenceFrame(datumNode, primeMeridian, dynamicNode)
@@ -3361,6 +3366,11 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
     auto cs = buildCS(csNode, node, UnitOfMeasure::NONE);
     auto cartesianCS = nn_dynamic_pointer_cast<CartesianCS>(cs);
 
+    // No explicit AXIS node ? (WKT1)
+    if (isNull(nodeP->lookForChild(WKTConstants::AXIS))) {
+        props.set("IMPLICIT_CS", true);
+    }
+
     if (isNull(csNode) && node->countChildrenOfName(WKTConstants::AXIS) == 0) {
         const auto methodCode = conversion->method()->getEPSGCode();
         // Krovak south oriented ?
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index bbbf2154..91ed145c 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -1216,6 +1216,23 @@ TEST(crs, geodeticcrs_identify_no_db) {
         EXPECT_EQ(res.front().first.get(), GeographicCRS::EPSG_4326.get());
         EXPECT_EQ(res.front().second, 25);
     }
+
+    {
+        // WKT1 identification
+        auto obj =
+            WKTParser()
+                .attachDatabaseContext(DatabaseContext::create())
+                .createFromWKT(
+                    "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\","
+                    "6378137,298.257223563]],PRIMEM[\"Greenwich\",0],"
+                    "UNIT[\"Degree\",0.0174532925199433]]");
+        auto crs = nn_dynamic_pointer_cast<GeodeticCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        auto res = crs->identify(nullptr);
+        ASSERT_EQ(res.size(), 1);
+        EXPECT_EQ(res.front().first.get(), GeographicCRS::EPSG_4326.get());
+        EXPECT_EQ(res.front().second, 100);
+    }
 }
 
 // ---------------------------------------------------------------------------
@@ -1893,6 +1910,71 @@ TEST(crs, projectedCRS_identify_db) {
         EXPECT_EQ(res.front().second, 70);
     }
     {
+        // Identify from a WKT1 string wit explicit correct axis order
+        auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(
+            "PROJCS[\"ETRS89 / UTM zone 32N (N-E)\",GEOGCS[\"ETRS89\","
+            "DATUM[\"European_Terrestrial_Reference_System_1989\","
+            "SPHEROID[\"GRS 1980\",6378137,298.257222101]],"
+            "PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],"
+            "PROJECTION[\"Transverse_Mercator\"],"
+            "PARAMETER[\"latitude_of_origin\",0],"
+            "PARAMETER[\"central_meridian\",9],"
+            "PARAMETER[\"scale_factor\",0.9996],"
+            "PARAMETER[\"false_easting\",500000],"
+            "PARAMETER[\"false_northing\",0],"
+            "UNIT[\"metre\",1],"
+            "AXIS[\"Northing\",NORTH],AXIS[\"Easting\",EAST]]");
+        auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        auto res = crs->identify(factoryEPSG);
+        ASSERT_EQ(res.size(), 1);
+        EXPECT_EQ(res.front().first->getEPSGCode(), 3044);
+        EXPECT_EQ(res.front().second, 100);
+    }
+    {
+        // Identify from a WKT1 string wit wrong axis order
+        auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(
+            "PROJCS[\"ETRS89 / UTM zone 32N (N-E)\",GEOGCS[\"ETRS89\","
+            "DATUM[\"European_Terrestrial_Reference_System_1989\","
+            "SPHEROID[\"GRS 1980\",6378137,298.257222101]],"
+            "PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],"
+            "PROJECTION[\"Transverse_Mercator\"],"
+            "PARAMETER[\"latitude_of_origin\",0],"
+            "PARAMETER[\"central_meridian\",9],"
+            "PARAMETER[\"scale_factor\",0.9996],"
+            "PARAMETER[\"false_easting\",500000],"
+            "PARAMETER[\"false_northing\",0],"
+            "UNIT[\"metre\",1],"
+            "AXIS[\"Easting\",EAST], AXIS[\"Northing\",NORTH]]");
+        auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        auto res = crs->identify(factoryEPSG);
+        ASSERT_EQ(res.size(), 1);
+        EXPECT_EQ(res.front().first->getEPSGCode(), 3044);
+        EXPECT_EQ(res.front().second, 25);
+    }
+    {
+        // Identify from a WKT1 string, without explicit axis
+        auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(
+            "PROJCS[\"ETRS89 / UTM zone 32N (N-E)\",GEOGCS[\"ETRS89\","
+            "DATUM[\"European_Terrestrial_Reference_System_1989\","
+            "SPHEROID[\"GRS 1980\",6378137,298.257222101]],"
+            "PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],"
+            "PROJECTION[\"Transverse_Mercator\"],"
+            "PARAMETER[\"latitude_of_origin\",0],"
+            "PARAMETER[\"central_meridian\",9],"
+            "PARAMETER[\"scale_factor\",0.9996],"
+            "PARAMETER[\"false_easting\",500000],"
+            "PARAMETER[\"false_northing\",0],"
+            "UNIT[\"metre\",1]]");
+        auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        auto res = crs->identify(factoryEPSG);
+        ASSERT_EQ(res.size(), 1);
+        EXPECT_EQ(res.front().first->getEPSGCode(), 3044);
+        EXPECT_EQ(res.front().second, 100);
+    }
+    {
         // No equivalent CRS to input one in result set
         auto obj =
             PROJStringParser().createFromPROJString("+proj=tmerc +datum=WGS84");
@@ -3410,6 +3492,40 @@ TEST(crs, boundCRS_identify_db) {
                   "NZGD2000 to WGS 84 (1)");
         EXPECT_EQ(res.front().second, 50);
     }
+
+    {
+        // WKT has EPSG code but the definition doesn't match with the official
+        // one (namely linear units are different)
+        // https://github.com/OSGeo/gdal/issues/990
+        // Also test that we can handle the synthetic Null geographic offset
+        // between NAD83 and WGS84
+        auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(
+            "PROJCS[\"NAD83 / Ohio North\",GEOGCS[\"NAD83\","
+            "DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\","
+            "6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],"
+            "TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6269\"]],"
+            "PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],"
+            "UNIT[\"degree\",0.0174532925199433, AUTHORITY[\"EPSG\",\"9122\"]],"
+            "AUTHORITY[\"EPSG\",\"4269\"]],"
+            "PROJECTION[\"Lambert_Conformal_Conic_2SP\"],"
+            "PARAMETER[\"standard_parallel_1\",41.7],"
+            "PARAMETER[\"standard_parallel_2\",40.43333333333333],"
+            "PARAMETER[\"latitude_of_origin\",39.66666666666666],"
+            "PARAMETER[\"central_meridian\",-82.5],"
+            "PARAMETER[\"false_easting\",1968503.937007874],"
+            "PARAMETER[\"false_northing\",0],"
+            "UNIT[\"International Foot\",0.3048,AUTHORITY[\"EPSG\",\"9002\"]],"
+            "AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"32122\"]"
+            "]");
+        auto crs = nn_dynamic_pointer_cast<BoundCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        auto res = crs->identify(factoryEPSG);
+        ASSERT_EQ(res.size(), 1);
+        EXPECT_EQ(res.front().second, 25);
+        auto wkt = crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get());
+        EXPECT_TRUE(wkt.find("32122") != std::string::npos) << wkt;
+    }
 }
 
 // ---------------------------------------------------------------------------