importFromWKT: better deal with axis of the baseCRS of a projected CRS

4 files changed, 109 insertions, 16 deletions

diff --git a/include/proj/util.hpp b/include/proj/util.hpp
index c2f2b7fe..c40595f2 100644
--- a/include/proj/util.hpp
+++ b/include/proj/util.hpp
@@ -506,6 +506,7 @@ class PropertyMap {
     bool getStringValue(const std::string &key, std::string &outVal) const;
     bool getStringValue(const std::string &key,
                         optional<std::string> &outVal) const;
+    void unset(const std::string &key);
 
     static PropertyMap createAndSetName(const char *name);
     static PropertyMap createAndSetName(const std::string &name);
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 6175c415..effb3968 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -1149,7 +1149,7 @@ struct WKTParser::Private {
     Private(const Private &) = delete;
     Private &operator=(const Private &) = delete;
 
-    void emitRecoverableAssertion(const std::string &errorMsg);
+    void emitRecoverableWarning(const std::string &errorMsg);
 
     BaseObjectNNPtr build(const WKTNodeNNPtr &node);
 
@@ -1331,7 +1331,7 @@ std::list<std::string> WKTParser::warningList() const {
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
-void WKTParser::Private::emitRecoverableAssertion(const std::string &errorMsg) {
+void WKTParser::Private::emitRecoverableWarning(const std::string &errorMsg) {
     if (strict_) {
         throw ParsingException(errorMsg);
     } else {
@@ -1688,8 +1688,8 @@ UnitOfMeasure WKTParser::Private::buildUnit(const WKTNodeNNPtr &node,
         auto &idNode =
             nodeP->lookForChild(WKTConstants::ID, WKTConstants::AUTHORITY);
         if (!isNull(idNode) && idNode->GP()->childrenSize() < 2) {
-            emitRecoverableAssertion("not enough children in " +
-                                     idNode->GP()->value() + " node");
+            emitRecoverableWarning("not enough children in " +
+                                   idNode->GP()->value() + " node");
         }
         const bool hasValidIdNode =
             !isNull(idNode) && idNode->GP()->childrenSize() >= 2;
@@ -2553,7 +2553,7 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) {
         // PRIMEM is required in WKT1
         if (ci_equal(nodeName, WKTConstants::GEOGCS) ||
             ci_equal(nodeName, WKTConstants::GEOCCS)) {
-            emitRecoverableAssertion(nodeName + " should have a PRIMEM node");
+            emitRecoverableWarning(nodeName + " should have a PRIMEM node");
         }
     }
 
@@ -2596,8 +2596,43 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) {
     if (ellipsoidalCS) {
         assert(!ci_equal(nodeName, WKTConstants::GEOCCS));
         try {
-            return GeographicCRS::create(props, datum, datumEnsemble,
-                                         NN_NO_CHECK(ellipsoidalCS));
+            auto crs = GeographicCRS::create(props, datum, datumEnsemble,
+                                             NN_NO_CHECK(ellipsoidalCS));
+            // In case of missing CS node, or to check it, query the coordinate
+            // system from the DB if possible (typically for the baseCRS of a
+            // ProjectedCRS)
+            if (!crs->identifiers().empty() && dbContext_) {
+                GeographicCRSPtr dbCRS;
+                try {
+                    const auto &id = crs->identifiers()[0];
+                    auto authFactory = AuthorityFactory::create(
+                        NN_NO_CHECK(dbContext_), *id->codeSpace());
+                    dbCRS = authFactory->createGeographicCRS(id->code())
+                                .as_nullable();
+                } catch (const util::Exception &) {
+                }
+                if (dbCRS &&
+                    (!isNull(csNode) ||
+                     node->countChildrenOfName(WKTConstants::AXIS) != 0) &&
+                    !ellipsoidalCS->_isEquivalentTo(
+                        dbCRS->coordinateSystem().get(),
+                        util::IComparable::Criterion::EQUIVALENT)) {
+                    emitRecoverableWarning(
+                        "Coordinate system of GeographicCRS in the WKT "
+                        "definition is different from the one of the "
+                        "authority. Unsetting the identifier to avoid "
+                        "confusion");
+                    props.unset(Identifier::CODESPACE_KEY);
+                    props.unset(Identifier::AUTHORITY_KEY);
+                    props.unset(IdentifiedObject::IDENTIFIERS_KEY);
+                    crs = GeographicCRS::create(props, datum, datumEnsemble,
+                                                NN_NO_CHECK(ellipsoidalCS));
+                } else if (dbCRS) {
+                    crs = GeographicCRS::create(props, datum, datumEnsemble,
+                                                dbCRS->coordinateSystem());
+                }
+            }
+            return crs;
         } catch (const util::Exception &e) {
             throw ParsingException(std::string("buildGeodeticCRS: ") +
                                    e.what());
@@ -4376,13 +4411,13 @@ BaseObjectNNPtr WKTParser::createFromWKT(const std::string &wkt) {
         dialect == WKTGuessedDialect::WKT1_ESRI) {
         auto errorMsg = pj_wkt1_parse(wkt);
         if (!errorMsg.empty()) {
-            d->emitRecoverableAssertion(errorMsg);
+            d->emitRecoverableWarning(errorMsg);
         }
     } else if (dialect == WKTGuessedDialect::WKT2_2015 ||
                dialect == WKTGuessedDialect::WKT2_2018) {
         auto errorMsg = pj_wkt2_parse(wkt);
         if (!errorMsg.empty()) {
-            d->emitRecoverableAssertion(errorMsg);
+            d->emitRecoverableWarning(errorMsg);
         }
     }
 
diff --git a/src/iso19111/util.cpp b/src/iso19111/util.cpp
index ac6357a2..b8c6c439 100644
--- a/src/iso19111/util.cpp
+++ b/src/iso19111/util.cpp
@@ -282,6 +282,17 @@ const BaseObjectNNPtr *PropertyMap::get(const std::string &key) const {
     }
     return nullptr;
 }
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
+void PropertyMap::unset(const std::string &key) {
+    for (auto iter = d->list_.begin(); iter != d->list_.end(); ++iter) {
+        if (iter->first == key) {
+            d->list_.erase(iter);
+            return;
+        }
+    }
+}
 //! @endcond
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 755f8258..71604ff8 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -350,12 +350,16 @@ static void checkEPSG_4326(GeographicCRSPtr crs, bool latLong = true,
     auto cs = crs->coordinateSystem();
     ASSERT_EQ(cs->axisList().size(), 2);
     if (latLong) {
-        EXPECT_EQ(cs->axisList()[0]->nameStr(), "Latitude");
-        EXPECT_EQ(cs->axisList()[0]->abbreviation(), "lat");
+        EXPECT_TRUE(cs->axisList()[0]->nameStr() == "Latitude" ||
+                    cs->axisList()[0]->nameStr() == "Geodetic latitude")
+            << cs->axisList()[0]->nameStr();
+        EXPECT_EQ(tolower(cs->axisList()[0]->abbreviation()), "lat");
         EXPECT_EQ(cs->axisList()[0]->direction(), AxisDirection::NORTH);
 
-        EXPECT_EQ(cs->axisList()[1]->nameStr(), "Longitude");
-        EXPECT_EQ(cs->axisList()[1]->abbreviation(), "lon");
+        EXPECT_TRUE(cs->axisList()[1]->nameStr() == "Longitude" ||
+                    cs->axisList()[1]->nameStr() == "Geodetic longitude")
+            << cs->axisList()[1]->nameStr();
+        EXPECT_EQ(tolower(cs->axisList()[1]->abbreviation()), "lon");
         EXPECT_EQ(cs->axisList()[1]->direction(), AxisDirection::EAST);
     } else {
         EXPECT_EQ(cs->axisList()[0]->nameStr(), "Longitude");
@@ -1005,8 +1009,6 @@ TEST(wkt_parse, wkt1_projected) {
                "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
                "        UNIT[\"degree\",0.0174532925199433,\n"
                "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
-               "        AXIS[\"latitude\",NORTH],\n"
-               "        AXIS[\"longitude\",EAST],\n"
                "        AUTHORITY[\"EPSG\",\"4326\"]],\n"
                "    PROJECTION[\"Transverse_Mercator\"],\n"
                "    PARAMETER[\"latitude_of_origin\",0],\n"
@@ -1019,10 +1021,14 @@ TEST(wkt_parse, wkt1_projected) {
                "    AXIS[\"(E)\",EAST],\n"
                "    AXIS[\"(N)\",NORTH],\n"
                "    AUTHORITY[\"EPSG\",\"32631\"]]";
-    auto obj = WKTParser().createFromWKT(wkt);
+    auto obj = WKTParser()
+                   .attachDatabaseContext(DatabaseContext::create())
+                   .createFromWKT(wkt);
     auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
     ASSERT_TRUE(crs != nullptr);
     checkProjected(crs);
+
+    EXPECT_TRUE(!crs->baseCRS()->identifiers().empty());
 }
 
 // ---------------------------------------------------------------------------
@@ -1058,6 +1064,46 @@ TEST(wkt_parse, wkt1_projected_no_axis) {
 
 // ---------------------------------------------------------------------------
 
+TEST(wkt_parse, wkt1_projected_wrong_axis_geogcs) {
+    auto wkt = "PROJCS[\"WGS 84 / UTM zone 31N\",\n"
+               "    GEOGCS[\"WGS 84\",\n"
+               "        DATUM[\"WGS_1984\",\n"
+               "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+               "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+               "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+               "        PRIMEM[\"Greenwich\",0,\n"
+               "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+               "        UNIT[\"degree\",0.0174532925199433,\n"
+               "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
+               "        AXIS[\"longitude\",EAST],\n"
+               "        AXIS[\"latitude\",NORTH],\n"
+               "        AUTHORITY[\"EPSG\",\"4326\"]],\n"
+               "    PROJECTION[\"Transverse_Mercator\"],\n"
+               "    PARAMETER[\"latitude_of_origin\",0],\n"
+               "    PARAMETER[\"central_meridian\",3],\n"
+               "    PARAMETER[\"scale_factor\",0.9996],\n"
+               "    PARAMETER[\"false_easting\",500000],\n"
+               "    PARAMETER[\"false_northing\",0],\n"
+               "    UNIT[\"metre\",1,\n"
+               "        AUTHORITY[\"EPSG\",\"9001\"]],\n"
+               "    AUTHORITY[\"EPSG\",\"32631\"]]";
+    WKTParser parser;
+    parser.setStrict(false).attachDatabaseContext(DatabaseContext::create());
+    auto obj = parser.createFromWKT(wkt);
+    EXPECT_TRUE(!parser.warningList().empty());
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+
+    EXPECT_TRUE(crs->baseCRS()->identifiers().empty());
+
+    auto cs = crs->baseCRS()->coordinateSystem();
+    ASSERT_EQ(cs->axisList().size(), 2);
+    EXPECT_EQ(cs->axisList()[0]->direction(), AxisDirection::EAST);
+    EXPECT_EQ(cs->axisList()[1]->direction(), AxisDirection::NORTH);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(wkt_parse, wkt1_projected_with_PROJ4_extension) {
     auto wkt = "PROJCS[\"unnamed\",\n"
                "    GEOGCS[\"WGS 84\",\n"