Adjustments to deal with non-conformant WKT1 LAS COMPD_CS[] (refs #2171)

2 files changed, 60 insertions, 0 deletions

diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 5fbc4c4b..bad7deea 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -566,6 +566,12 @@ CRSNNPtr CRS::shallowClone() const { return _shallowClone(); }
 //! @cond Doxygen_Suppress
 
 CRSNNPtr CRS::allowNonConformantWKT1Export() const {
+    const auto boundCRS = dynamic_cast<const BoundCRS *>(this);
+    if (boundCRS) {
+        return BoundCRS::create(
+            boundCRS->baseCRS()->allowNonConformantWKT1Export(),
+            boundCRS->hubCRS(), boundCRS->transformation());
+    }
     auto crs(shallowClone());
     crs->d->allowNonConformantWKT1Export_ = true;
     return crs;
@@ -1422,7 +1428,10 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const {
             formatter->startNode(io::WKTConstants::COMPD_CS, false);
             formatter->addQuotedString(l_name + " + " + l_name);
             geogCRS2D->_exportToWKT(formatter);
+            const auto oldTOWGSParameters = formatter->getTOWGS84Parameters();
+            formatter->setTOWGS84Parameters({});
             geogCRS2D->_exportToWKT(formatter);
+            formatter->setTOWGS84Parameters(oldTOWGSParameters);
             formatter->endNode();
             return;
         }
@@ -4245,6 +4254,14 @@ CRSNNPtr CompoundCRS::createLax(const util::PropertyMap &properties,
         auto comp1 = components[1].get();
         auto comp0Geog = dynamic_cast<const GeographicCRS *>(comp0);
         auto comp0Proj = dynamic_cast<const ProjectedCRS *>(comp0);
+        if (comp0Geog == nullptr && comp0Proj == nullptr) {
+            auto comp0Bound = dynamic_cast<const BoundCRS *>(comp0);
+            if (comp0Bound) {
+                const auto *baseCRS = comp0Bound->baseCRS().get();
+                comp0Geog = dynamic_cast<const GeographicCRS *>(baseCRS);
+                comp0Proj = dynamic_cast<const ProjectedCRS *>(baseCRS);
+            }
+        }
         auto comp1Geog = dynamic_cast<const GeographicCRS *>(comp1);
         if ((comp0Geog != nullptr || comp0Proj != nullptr) &&
             comp1Geog != nullptr) {
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 50561953..d62e5c59 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -2485,6 +2485,49 @@ TEST(wkt_parse, COMPD_CS_non_conformant_horizontal_plus_horizontal_as_in_LAS) {
 // ---------------------------------------------------------------------------
 
 TEST(wkt_parse,
+     COMPD_CS_non_conformant_horizontal_TOWGS84_plus_horizontal_as_in_LAS) {
+
+    const auto wkt = "COMPD_CS[\"WGS 84 + WGS 84\",\n"
+                     "    GEOGCS[\"WGS 84\",\n"
+                     "        DATUM[\"WGS_1984\",\n"
+                     "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n"
+                     "                AUTHORITY[\"EPSG\",\"7030\"]],\n"
+                     "            TOWGS84[0,0,0,0,0,0,0],\n"
+                     "            AUTHORITY[\"EPSG\",\"6326\"]],\n"
+                     "        PRIMEM[\"Greenwich\",0,\n"
+                     "            AUTHORITY[\"EPSG\",\"8901\"]],\n"
+                     "        UNIT[\"degree\",0.0174532925199433,\n"
+                     "            AUTHORITY[\"EPSG\",\"9122\"]],\n"
+                     "        AUTHORITY[\"EPSG\",\"4326\"]],\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"
+                     "        AUTHORITY[\"EPSG\",\"4326\"]]]";
+    auto dbContext = DatabaseContext::create();
+    auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt);
+    auto crs = nn_dynamic_pointer_cast<BoundCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    auto baseCRS = nn_dynamic_pointer_cast<GeographicCRS>(crs->baseCRS());
+    ASSERT_TRUE(baseCRS != nullptr);
+    EXPECT_EQ(baseCRS->nameStr(), "WGS 84");
+    EXPECT_EQ(baseCRS->coordinateSystem()->axisList().size(), 3U);
+
+    EXPECT_EQ(
+        crs->exportToWKT(
+            WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL, dbContext)
+                .get()),
+        wkt);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(wkt_parse,
      COMPD_CS_horizontal_bound_geog_plus_vertical_ellipsoidal_height) {
     // See https://github.com/OSGeo/PROJ/issues/2228
     const char *wkt =