ProjectedCRS identification: deal with switched 1st/2nd std parallels for LCC_2SP (#2150)

As switching the 2 standard parallels for Lambert Conformal Conic 2 standard
parallels result in a equivalently numerical projection definition, take
that into account when querying the database.

2 files changed, 87 insertions, 2 deletions

diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index c6c25f36..32ed6abb 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -6051,8 +6051,15 @@ AuthorityFactory::createProjectedCRSFromExisting(
         params.emplace_back(d->authority());
     }
 
-    int iParam = 1;
+    int iParam = 0;
+    bool hasLat1stStd = false;
+    double lat1stStd = 0;
+    int iParamLat1stStd = 0;
+    bool hasLat2ndStd = false;
+    double lat2ndStd = 0;
+    int iParamLat2ndStd = 0;
     for (const auto &genOpParamvalue : conv->parameterValues()) {
+        iParam++;
         auto opParamvalue =
             dynamic_cast<const operation::OperationParameterValue *>(
                 genOpParamvalue.get());
@@ -6070,6 +6077,23 @@ AuthorityFactory::createProjectedCRSFromExisting(
         const auto &unit = measure.unit();
         if (unit == common::UnitOfMeasure::DEGREE &&
             geogCRS->coordinateSystem()->axisList()[0]->unit() == unit) {
+            if (methodEPSGCode ==
+                EPSG_CODE_METHOD_LAMBERT_CONIC_CONFORMAL_2SP) {
+                // Special case for standard parallels of LCC_2SP. See below
+                if (paramEPSGCode ==
+                    EPSG_CODE_PARAMETER_LATITUDE_1ST_STD_PARALLEL) {
+                    hasLat1stStd = true;
+                    lat1stStd = measure.value();
+                    iParamLat1stStd = iParam;
+                    continue;
+                } else if (paramEPSGCode ==
+                           EPSG_CODE_PARAMETER_LATITUDE_2ND_STD_PARALLEL) {
+                    hasLat2ndStd = true;
+                    lat2ndStd = measure.value();
+                    iParamLat2ndStd = iParam;
+                    continue;
+                }
+            }
             const auto iParamAsStr(toString(iParam));
             sql += " AND conv.param";
             sql += iParamAsStr;
@@ -6084,7 +6108,44 @@ AuthorityFactory::createProjectedCRSFromExisting(
             params.emplace_back(measure.value() - 1);
             params.emplace_back(measure.value() + 1);
         }
-        iParam++;
+    }
+
+    // Special case for standard parallels of LCC_2SP: they can be switched
+    if (methodEPSGCode == EPSG_CODE_METHOD_LAMBERT_CONIC_CONFORMAL_2SP &&
+        hasLat1stStd && hasLat2ndStd) {
+        const auto iParam1AsStr(toString(iParamLat1stStd));
+        const auto iParam2AsStr(toString(iParamLat2ndStd));
+        sql += " AND conv.param";
+        sql += iParam1AsStr;
+        sql += "_code = ? AND conv.param";
+        sql += iParam1AsStr;
+        sql += "_auth_name = 'EPSG' AND conv.param";
+        sql += iParam2AsStr;
+        sql += "_code = ? AND conv.param";
+        sql += iParam2AsStr;
+        sql += "_auth_name = 'EPSG' AND ((";
+        params.emplace_back(
+            toString(EPSG_CODE_PARAMETER_LATITUDE_1ST_STD_PARALLEL));
+        params.emplace_back(
+            toString(EPSG_CODE_PARAMETER_LATITUDE_2ND_STD_PARALLEL));
+        double val1 = lat1stStd;
+        double val2 = lat2ndStd;
+        for (int i = 0; i < 2; i++) {
+            if (i == 1) {
+                sql += ") OR (";
+                std::swap(val1, val2);
+            }
+            sql += "conv.param";
+            sql += iParam1AsStr;
+            sql += "_value BETWEEN ? AND ? AND conv.param";
+            sql += iParam2AsStr;
+            sql += "_value BETWEEN ? AND ?";
+            params.emplace_back(val1 - 1);
+            params.emplace_back(val1 + 1);
+            params.emplace_back(val2 - 1);
+            params.emplace_back(val2 + 1);
+        }
+        sql += "))";
     }
     auto sqlRes = d->run(sql, params);
 
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index 2bde8ee3..ad32d98b 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -2550,6 +2550,30 @@ TEST(crs, projectedCRS_identify_db) {
         EXPECT_EQ(res.front().first->getEPSGCode(), 2954);
         EXPECT_EQ(res.front().second, 100);
     }
+    {
+        // Test identification of LCC_2SP with switched standard parallels.
+        auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(
+            "PROJCS[\"foo\",\n"
+            "    GEOGCS[\"RGF93\",\n"
+            "        DATUM[\"Reseau_Geodesique_Francais_1993\",\n"
+            "            SPHEROID[\"GRS 1980\",6378137,298.257222101]],\n"
+            "        PRIMEM[\"Greenwich\",0],\n"
+            "        UNIT[\"degree\",0.0174532925199433]],\n"
+            "    PROJECTION[\"Lambert_Conformal_Conic_2SP\"],\n"
+            "    PARAMETER[\"latitude_of_origin\",46.5],\n"
+            "    PARAMETER[\"central_meridian\",3],\n"
+            "    PARAMETER[\"standard_parallel_1\",44],\n"
+            "    PARAMETER[\"standard_parallel_2\",49],\n"
+            "    PARAMETER[\"false_easting\",700000],\n"
+            "    PARAMETER[\"false_northing\",6600000],\n"
+            "    UNIT[\"metre\",1]]");
+        auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+        auto res = crs->identify(factoryEPSG);
+        ASSERT_EQ(res.size(), 1U);
+        EXPECT_EQ(res.front().first->getEPSGCode(), 2154);
+        EXPECT_EQ(res.front().second, 70);
+    }
 }
 
 // ---------------------------------------------------------------------------