Yet another MEGA (Make ETRS89 Great Again) fix!

I've been made aware that ETRS89 to MGI was semi broken since the EPSG
update that added a ETRS89-AUT [2002] datum and migrated the 'MGI to
ETRS89 (8)' transformation using the last released grid AT_GIS_GRID_2021_09_28
to be 'MGI to ETRS89-AUT [2002]' and thus leaving user with the older
transformation 'MGI to ETRS89 (5)'.
This is just one manifestation of a problem that can be seen with ETRS89
in the Netherlands, Romania, etc.
Until EPSG perhaps figures out a plan to restore cleanly backward
compatibility (they have been made aware of those issues), I've devised
the following s/hack/strategy that seems to work (TM):

When looking for ETRS89 (EPSG:4258 the datum ensemble) to XXXX
transformations, also query the ``alias_name`` table for old names
of transformations that start with 'ETRS89 to XXXX' or
'XXXX to ETRS89', and are aliases to 'ETRS89-YYY to XXXX[something]'
or 'XXXX to ETRS89-YYY[something]' newly named transformation. And then
consider those 'migrated' transformations to be also used for plain
ETRS89, because they used to be used for it!

No new test cases, but actually reverting to old reference values before
commit 047e5db (update to EPSG 12.033)
that changed them.

Author: rouault

include/proj/io.hpp

@@ -1328,6 +1328,14 @@ class PROJ_GCC_DLL AuthorityFactory {
         const metadata::ExtentPtr &intersectingExtent1,
         const metadata::ExtentPtr &intersectingExtent2) const;
 
+    PROJ_INTERNAL std::vector<operation::CoordinateOperationNNPtr>
+    getOperationsFromAlias(const std::string &crs1Name,
+                           const std::string &crs2Name,
+                           bool usePROJAlternativeGridNames,
+                           bool discardIfMissingGrid,
+                           bool considerKnownGridsAsAvailable,
+                           bool discardSuperseded) const;
+
     typedef std::pair<common::IdentifiedObjectNNPtr, std::string>
         PairObjectName;
     PROJ_INTERNAL std::list<PairObjectName>

src/iso19111/factory.cpp

@@ -10431,6 +10431,62 @@ AuthorityFactory::getPointMotionOperationsFor(
     return res;
 }
 
+// ---------------------------------------------------------------------------
+
+std::vector<operation::CoordinateOperationNNPtr>
+AuthorityFactory::getOperationsFromAlias(const std::string &crs1Name,
+                                         const std::string &crs2Name,
+                                         bool usePROJAlternativeGridNames,
+                                         bool discardIfMissingGrid,
+                                         bool considerKnownGridsAsAvailable,
+                                         bool discardSuperseded) const {
+    std::string sql("SELECT alias.auth_name, alias.code FROM alias_name alias "
+                    "JOIN coordinate_operation_view cov "
+                    "ON alias.table_name = cov.table_name "
+                    "AND alias.auth_name = cov.auth_name "
+                    "AND alias.code = cov.code "
+                    "WHERE alias.table_name IN ('grid_transformation', "
+                    "'helmert_transformation', 'other_transformation', "
+                    "'concatenated_operation') "
+                    "AND (alt_name LIKE ? OR alt_name LIKE ?) "
+                    "AND NOT (alt_name LIKE ? OR alt_name LIKE ? OR alt_name "
+                    "LIKE ? OR alt_name LIKE ?) "
+                    "AND cov.deprecated = 0");
+    if (discardSuperseded) {
+        sql += " AND NOT EXISTS (SELECT 1 FROM supersession ss WHERE "
+               "ss.superseded_table_name = cov.table_name AND "
+               "ss.superseded_auth_name = cov.auth_name AND "
+               "ss.superseded_code = cov.code AND "
+               "ss.superseded_table_name = ss.replacement_table_name AND "
+               "ss.same_source_target_crs = 1)";
+    }
+    ListOfParams params{
+        std::string(crs1Name).append(" to ").append(crs2Name).append("%"),
+        std::string(crs2Name).append(" to ").append(crs1Name).append("%"),
+        //  If looking for "MGI to ETRS89", don't match "ETRS89 to ETRS89-XXX"
+        std::string(crs1Name).append(" to ").append(crs1Name).append("-%"),
+        std::string(crs2Name).append(" to ").append(crs2Name).append("-%"),
+        // If looking for "MGI to ETRS89", don't match "MGI to ETRS89-XXX" (not
+        // an actual example, but just in case)
+        std::string(crs1Name).append(" to ").append(crs2Name).append("-%"),
+        std::string(crs2Name).append(" to ").append(crs1Name).append("-%"),
+    };
+
+    std::vector<operation::CoordinateOperationNNPtr> res;
+    auto sqlRes = d->run(sql, params);
+    for (const auto &row : sqlRes) {
+        const auto &auth_name = row[0];
+        const auto &code = row[1];
+        auto op = d->createFactory(auth_name)->createCoordinateOperation(
+            code, usePROJAlternativeGridNames);
+        if (!discardIfMissingGrid ||
+            !d->rejectOpDueToMissingGrid(op, considerKnownGridsAsAvailable)) {
+            res.emplace_back(op);
+        }
+    }
+    return res;
+}
+
 //! @endcond
 
 // ---------------------------------------------------------------------------

src/iso19111/operation/coordinateoperationfactory.cpp

@@ -4270,6 +4270,80 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
         }
     }
 
+    // ETRS89 specific hack to deal with the fact that a lot of ETRS89-XXX
+    // datum/CRS have been created, with older transformations going from/to
+    // generic ETRS89 being re-attributed to that specific ETRS89-XXX datum/CRS
+    // which breaks backwards compatibility
+    // e.g we want ETRS89 to MGI to be able to use the previously named
+    // "MGI to ETRS89 (8)" operation, that is now "MGI to ETRS89-AUT [2002] (8)"
+    const bool srcIsETRS89 = sourceCRS->nameStr() == "ETRS89";
+    const bool dstIsETRS89 = targetCRS->nameStr() == "ETRS89";
+    if (geodSrc && geodDst && (srcIsETRS89 || dstIsETRS89) &&
+        !(srcIsETRS89 && dstIsETRS89)) {
+        const auto &authFactory = context.context->getAuthorityFactory();
+        const auto gridAvailabilityUse =
+            context.context->getGridAvailabilityUse();
+        const auto opFromAlias = authFactory->getOperationsFromAlias(
+            sourceCRS->nameStr(), targetCRS->nameStr(),
+            context.context->getUsePROJAlternativeGridNames(),
+            /* discardIfMissingGrid = */
+            gridAvailabilityUse ==
+                    CoordinateOperationContext::GridAvailabilityUse::
+                        DISCARD_OPERATION_IF_MISSING_GRID ||
+                gridAvailabilityUse == CoordinateOperationContext::
+                                           GridAvailabilityUse::KNOWN_AVAILABLE,
+            /* considerKnownGridsAsAvailable = */
+            gridAvailabilityUse == CoordinateOperationContext::
+                                       GridAvailabilityUse::KNOWN_AVAILABLE,
+            context.context->getDiscardSuperseded());
+        for (const auto &aliasOp : opFromAlias) {
+            const auto aliasOpSourceCRS = aliasOp->sourceCRS();
+            const auto aliasOpTargetCRS = aliasOp->targetCRS();
+            if (aliasOpSourceCRS && aliasOpTargetCRS &&
+                (aliasOpSourceCRS->nameStr().find("ETRS89-") !=
+                     std::string::npos ||
+                 aliasOpTargetCRS->nameStr().find("ETRS89-") !=
+                     std::string::npos)) {
+                const bool bSourceAliasIsEquivalentToSourceCRS = starts_with(
+                    aliasOpSourceCRS->nameStr(), sourceCRS->nameStr());
+                auto newOp = aliasOp->shallowClone();
+                // Patch the transformation to modify its source/target CRS
+                // to the expected one
+                if (bSourceAliasIsEquivalentToSourceCRS) {
+                    setCRSs(newOp.get(), sourceCRS, targetCRS);
+                } else {
+                    setCRSs(newOp.get(), targetCRS, sourceCRS);
+                    newOp = newOp->inverse();
+                }
+
+                const auto getNatureOf =
+                    [](const crs::GeodeticCRS &crs) -> std::string {
+                    if (dynamic_cast<const crs::GeographicCRS *>(&crs) !=
+                        nullptr) {
+                        if (crs.coordinateSystem()->axisList().size() == 2)
+                            return "geographic2D";
+                        else
+                            return "geographic3D";
+                    } else {
+                        return "geodetic";
+                    }
+                };
+                // TODO? Ideally we'd create 2D<-->3D<-->geocentric conversions
+                // when needed
+                const auto newOpSrc = dynamic_cast<const crs::GeodeticCRS *>(
+                    newOp->sourceCRS().get());
+                const auto newOpDst = dynamic_cast<const crs::GeodeticCRS *>(
+                    newOp->targetCRS().get());
+                if (newOpSrc && newOpDst &&
+                    getNatureOf(*newOpSrc) == getNatureOf(*geodSrc) &&
+                    getNatureOf(*newOpDst) == getNatureOf(*geodDst)) {
+                    res.emplace_back(newOp);
+                    doFilterAndCheckPerfectOp = true;
+                }
+            }
+        }
+    }
+
     if (doFilterAndCheckPerfectOp) {
         // If we get at least a result with perfect accuracy, do not bother
         // generating synthetic transforms.

test/unit/gie_self_tests.cpp

@@ -1297,15 +1297,15 @@ TEST(gie, proj_create_crs_to_crs_PULKOVO42_ETRS89) {
     c.xyzt.z = 0;
     c.xyzt.t = HUGE_VAL;
     c = proj_trans(P, PJ_FWD, c);
-    EXPECT_NEAR(c.xy.x, 45.000346961, 1e-9);
-    EXPECT_NEAR(c.xy.y, 25.001524075, 1e-9);
+    EXPECT_NEAR(c.xy.x, 44.999701238, 1e-9);
+    EXPECT_NEAR(c.xy.y, 24.998474948, 1e-9);
     EXPECT_EQ(std::string(proj_pj_info(P).definition),
               "proj=pipeline step proj=axisswap order=2,1 "
               "step proj=unitconvert xy_in=deg xy_out=rad "
               "step proj=push v_3 "
               "step proj=cart "
-              "ellps=krass step proj=helmert x=68.1564 y=32.7756 z=80.2249 "
-              "rx=2.20333014 ry=2.19256447 rz=-2.54166911 s=-0.14155333 "
+              "ellps=krass step proj=helmert x=2.3287 y=-147.0425 z=-92.0802 "
+              "rx=0.3092483 ry=-0.32482185 rz=-0.49729934 s=5.68906266 "
               "convention=coordinate_frame step inv proj=cart ellps=GRS80 "
               "step proj=pop v_3 "
               "step proj=unitconvert xy_in=rad xy_out=deg step proj=axisswap "
@@ -1322,8 +1322,8 @@ TEST(gie, proj_create_crs_to_crs_PULKOVO42_ETRS89) {
     proj_trans_generic(P, PJ_FWD, &(c.xyz.x), sizeof(double), 1, &(c.xyz.y),
                        sizeof(double), 1, &(c.xyz.z), sizeof(double), 1,
                        nullptr, 0, 0);
-    EXPECT_NEAR(c.xy.x, 45.000346961, 1e-9);
-    EXPECT_NEAR(c.xy.y, 25.001524075, 1e-9);
+    EXPECT_NEAR(c.xy.x, 44.999701238, 1e-9);
+    EXPECT_NEAR(c.xy.y, 24.998474948, 1e-9);
 
     // Poland
     c.xyz.x = 52; // Lat

test/unit/test_network.cpp

@@ -1386,7 +1386,8 @@ TEST(networking, network_endpoint_api_and_not_reachable_hgridshift) {
         EXPECT_EQ(c2.xyz.x, HUGE_VAL);
         EXPECT_EQ(proj_errno(P), PROJ_ERR_OTHER_NETWORK_ERROR);
         PJ *last_op = proj_trans_get_last_used_operation(P);
-        EXPECT_STREQ(proj_pj_info(last_op).description, "MGI to ETRS89 (5)");
+        EXPECT_STREQ(proj_pj_info(last_op).description,
+                     "MGI to ETRS89-AUT [2002] (8)");
         proj_destroy(last_op);
     }
 
@@ -1398,7 +1399,8 @@ TEST(networking, network_endpoint_api_and_not_reachable_hgridshift) {
         EXPECT_EQ(c2.xyz.x, HUGE_VAL);
         EXPECT_EQ(proj_errno(P), PROJ_ERR_OTHER_NETWORK_ERROR);
         PJ *last_op = proj_trans_get_last_used_operation(P);
-        EXPECT_STREQ(proj_pj_info(last_op).description, "MGI to ETRS89 (5)");
+        EXPECT_STREQ(proj_pj_info(last_op).description,
+                     "MGI to ETRS89-AUT [2002] (8)");
         proj_destroy(last_op);
     }
 
@@ -1410,7 +1412,8 @@ TEST(networking, network_endpoint_api_and_not_reachable_hgridshift) {
         EXPECT_EQ(c2.xyz.x, HUGE_VAL);
         EXPECT_EQ(proj_errno(P), PROJ_ERR_OTHER_NETWORK_ERROR);
         PJ *last_op = proj_trans_get_last_used_operation(P);
-        EXPECT_STREQ(proj_pj_info(last_op).description, "MGI to ETRS89 (5)");
+        EXPECT_STREQ(proj_pj_info(last_op).description,
+                     "MGI to ETRS89-AUT [2002] (8)");
         proj_destroy(last_op);
     }
 

test/unit/test_operationfactory.cpp

@@ -95,20 +95,24 @@ TEST(operation, geogCRS_to_geogCRS_context_default) {
             authFactory->createCoordinateReferenceSystem("4179"), // Pulkovo 42
             authFactory->createCoordinateReferenceSystem("4258"), // ETRS89
             ctxt);
-        ASSERT_EQ(list.size(), 3U);
+        ASSERT_EQ(list.size(), 4U);
         // Romania has a larger area than Poland (given our approx formula)
-        EXPECT_EQ(list[0]->getEPSGCode(), 15993); // Romania - 10m
-        EXPECT_EQ(list[1]->getEPSGCode(), 1644);  // Poland - 1m
-        EXPECT_EQ(list[2]->nameStr(),
+        EXPECT_EQ(
+            list[0]->nameStr(),
+            "Pulkovo 1942(58) to ETRS89-ROU [ETRF2000] (4)"); // Romania - 10m
+        EXPECT_EQ(list[0]->getEPSGCode(), 15994);             // Romania - 3m
+        EXPECT_EQ(list[1]->getEPSGCode(), 15993);             // Romania - 10m
+        EXPECT_EQ(list[2]->getEPSGCode(), 1644);              // Poland - 1m
+        EXPECT_EQ(list[3]->nameStr(),
                   "Ballpark geographic offset from Pulkovo 1942(58) to ETRS89");
 
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
             "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
             "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
-            "+step +proj=cart +ellps=krass +step +proj=helmert +x=68.1564 "
-            "+y=32.7756 +z=80.2249 +rx=2.20333014 +ry=2.19256447 "
-            "+rz=-2.54166911 +s=-0.14155333 +convention=coordinate_frame +step "
+            "+step +proj=cart +ellps=krass +step +proj=helmert +x=2.3287 "
+            "+y=-147.0425 +z=-92.0802 +rx=0.3092483 +ry=-0.32482185 "
+            "+rz=-0.49729934 +s=5.68906266 +convention=coordinate_frame +step "
             "+inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step "
             "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
             "+order=2,1");
@@ -119,18 +123,21 @@ TEST(operation, geogCRS_to_geogCRS_context_default) {
         auto list = CoordinateOperationFactory::create()->createOperations(
             authFactory->createCoordinateReferenceSystem("4258"),
             authFactory->createCoordinateReferenceSystem("4179"), ctxt);
-        ASSERT_EQ(list.size(), 3U);
+        ASSERT_EQ(list.size(), 4U);
         // Romania has a larger area than Poland (given our approx formula)
-        EXPECT_EQ(list[0]->nameStr(),
-                  "Inverse of Pulkovo 1942(58) to ETRS89 (3)"); // Romania - 10m
+        EXPECT_EQ(
+            list[0]->nameStr(),
+            "Inverse of Pulkovo 1942(58) to ETRS89-ROU [ETRF2000] (4)"); // Romania
+                                                                         // -
+                                                                         // 10m
 
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
             "+proj=pipeline +step +proj=axisswap +order=2,1 +step "
             "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 "
-            "+step +proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=68.1564 "
-            "+y=32.7756 +z=80.2249 +rx=2.20333014 +ry=2.19256447 "
-            "+rz=-2.54166911 +s=-0.14155333 +convention=coordinate_frame +step "
+            "+step +proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=2.3287 "
+            "+y=-147.0425 +z=-92.0802 +rx=0.3092483 +ry=-0.32482185 "
+            "+rz=-0.49729934 +s=5.68906266 +convention=coordinate_frame +step "
             "+inv +proj=cart +ellps=krass +step +proj=pop +v_3 +step "
             "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap "
             "+order=2,1");
@@ -208,8 +215,9 @@ TEST(operation, geogCRS_to_geogCRS_context_filter_bbox) {
         auto list = CoordinateOperationFactory::create()->createOperations(
             authFactory->createCoordinateReferenceSystem("4179"),
             authFactory->createCoordinateReferenceSystem("4258"), ctxt);
-        ASSERT_EQ(list.size(), 1U);
-        EXPECT_EQ(list[0]->getEPSGCode(), 15993); // Romania - 10m
+        ASSERT_EQ(list.size(), 2U);
+        EXPECT_EQ(list[0]->getEPSGCode(), 15994); // Romania - 3m
+        EXPECT_EQ(list[1]->getEPSGCode(), 15993); // Romania - 10m
     }
     {
         auto ctxt = CoordinateOperationContext::create(
@@ -220,8 +228,9 @@ TEST(operation, geogCRS_to_geogCRS_context_filter_bbox) {
         auto list = CoordinateOperationFactory::create()->createOperations(
             authFactory->createCoordinateReferenceSystem("4179"),
             authFactory->createCoordinateReferenceSystem("4258"), ctxt);
-        ASSERT_EQ(list.size(), 1U);
-        EXPECT_EQ(list[0]->getEPSGCode(), 15993); // Romania - 10m
+        ASSERT_EQ(list.size(), 2U);
+        EXPECT_EQ(list[0]->getEPSGCode(), 15994); // Romania - 3m
+        EXPECT_EQ(list[1]->getEPSGCode(), 15993); // Romania - 10m
     }
     {
         auto ctxt = CoordinateOperationContext::create(
@@ -1510,7 +1519,7 @@ TEST(operation, geogCRS_without_id_to_geogCRS_3D_context) {
     auto list =
         CoordinateOperationFactory::create()->createOperations(src, dst, ctxt);
     ASSERT_GE(list.size(), 1U);
-    auto wkt2 = "GEOGCRS[\"unnamed\",\n"
+    auto wkt2 = "GEOGCRS[\"Amersfoort\",\n"
                 "    DATUM[\"Amersfoort\",\n"
                 "        ELLIPSOID[\"Bessel 1841\",6377397.155,299.1528128,\n"
                 "            LENGTHUNIT[\"metre\",1]]],\n"
@@ -12294,10 +12303,9 @@ TEST(operation, createOperation_ETRS89_to_Amersfoort) {
             // Amersfoort
             authFactoryEPSG->createCoordinateReferenceSystem("4289"), ctxt);
         ASSERT_GE(list.size(), 1U);
-        // We check that we go through ETRS89-NLD [AGRS2010] to use the most
+        // We check that we use the most
         // precise "Amersfoort to ETRS89-NLD [AGRS2010] (9)" operation.
         EXPECT_EQ(list[0]->nameStr(),
-                  "ETRS89 to ETRS89-NLD [AGRS2010] + "
                   "Inverse of Amersfoort to ETRS89-NLD [AGRS2010] (9)");
     }
     {
@@ -12307,10 +12315,9 @@ TEST(operation, createOperation_ETRS89_to_Amersfoort) {
             // ETRS89
             authFactoryEPSG->createCoordinateReferenceSystem("4258"), ctxt);
         ASSERT_GE(list.size(), 1U);
-        // We check that we go through ETRS89-NLD [AGRS2010] to use the most
+        // We check that we use the most
         // precise "Amersfoort to ETRS89-NLD [AGRS2010] (9)" operation.
         EXPECT_EQ(list[0]->nameStr(),
-                  "Amersfoort to ETRS89-NLD [AGRS2010] (9) + "
-                  "Inverse of ETRS89 to ETRS89-NLD [AGRS2010]");
+                  "Amersfoort to ETRS89-NLD [AGRS2010] (9)");
     }
 }