From db10e86c2c4c4606ebf90961483d630ea81f748b Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Thu, 10 Jun 2021 19:36:16 +0200
Subject: [PATCH] BoundCRS::identify(): avoid incompatible transformation for
 WKT1 / TOWGS84 export (fixes OSGeo/gdal#3958)

---
 src/iso19111/crs.cpp   | 10 +++++++++-
 test/unit/test_crs.cpp | 26 +++++++++++++++++++++++---
 2 files changed, 32 insertions(+), 4 deletions(-)

diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 51317c18d0..b0f8a30f79 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -46,6 +46,8 @@
 #include "proj/internal/internal.hpp"
 #include "proj/internal/io_internal.hpp"
 
+#include "operation/oputils.hpp"
+
 #include "proj_constants.h"
 #include "proj_json_streaming_writer.hpp"
 
@@ -5519,7 +5521,13 @@ BoundCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                     auto opNormalized = op->normalizeForVisualization();
                     std::string opTransfPROJString;
                     bool opTransfPROJStringValid = false;
-                    if (op->nameStr().find("Ballpark geographic") == 0) {
+                    const auto &opName = op->nameStr();
+                    if (starts_with(opName,
+                                    osgeo::proj::operation::
+                                        BALLPARK_GEOCENTRIC_TRANSLATION) ||
+                        starts_with(
+                            opName,
+                            osgeo::proj::operation::NULL_GEOGRAPHIC_OFFSET)) {
                         if (refIsNullTransform) {
                             res.emplace_back(create(candidateBaseCRS,
                                                     d->hubCRS_,
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index 614580f213..10b21ee580 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -4646,9 +4646,10 @@ TEST(crs, boundCRS_identify_db) {
         auto res = crs->identify(factoryEPSG);
         ASSERT_EQ(res.size(), 1U);
         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;
+        auto boundCRS = dynamic_cast<const BoundCRS *>(res.front().first.get());
+        ASSERT_TRUE(boundCRS != nullptr);
+        EXPECT_EQ(boundCRS->baseCRS()->getEPSGCode(), 32122);
+        EXPECT_EQ(boundCRS->transformation()->method()->getEPSGCode(), 9603);
     }
 
     {
@@ -4665,6 +4666,25 @@ TEST(crs, boundCRS_identify_db) {
         EXPECT_EQ(boundCRS->baseCRS()->getEPSGCode(), 3148);
         EXPECT_EQ(res.front().second, 70);
     }
+
+    {
+        // Identify a WKT with datum WGS84 and TOWGS84
+        auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(
+            "GEOGCS[\"WGS84 Coordinate System\",DATUM[\"WGS 1984\","
+            "SPHEROID[\"WGS 1984\",6378137,298.257223563],"
+            "TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6326\"]],"
+            "PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],"
+            "AUTHORITY[\"EPSG\",\"4326\"]]");
+        auto crs = nn_dynamic_pointer_cast<BoundCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        auto res = crs->identify(factoryEPSG);
+        ASSERT_EQ(res.size(), 1U);
+        EXPECT_EQ(res.front().second, 100);
+        auto boundCRS = dynamic_cast<const BoundCRS *>(res.front().first.get());
+        ASSERT_TRUE(boundCRS != nullptr);
+        EXPECT_EQ(boundCRS->baseCRS()->getEPSGCode(), 4326);
+        EXPECT_EQ(boundCRS->transformation()->method()->getEPSGCode(), 9606);
+    }
 }
 
 // ---------------------------------------------------------------------------