Coordinate transformation: improve transformations from/to WGS84 (Gxxxx)

Currently very few transformations from/to WGS84 (Gxxxx) are registered
in the EPSG database, and there isn't even transformations between WGS84 EPSG:4326
and those ones. Consequently transformations to those realizations often
ended up as no-operation, whereas going through WGS84 EPSG:4326 will bring
more meaningful results. So register those EPSG:4326<-->WGS 84 (Gxxx)
null transformations, and when having WGS 84 (Gxxx) as source/target,
consider EPSG:4326 as an intermediate.
This change has no effect on the existing direct transformations
from/to WGS 84 (Gxxx).

6 files changed, 259 insertions, 1 deletions

diff --git a/data/sql/customizations.sql b/data/sql/customizations.sql
index 5e038de4..98633045 100644
--- a/data/sql/customizations.sql
+++ b/data/sql/customizations.sql
@@ -65,3 +65,20 @@ INSERT INTO "coordinate_system" VALUES('PROJ','ENh','Cartesian',3);
 INSERT INTO "axis" VALUES('PROJ','1','Easting','E','east','PROJ','ENh',1,'EPSG','9001');
 INSERT INTO "axis" VALUES('PROJ','2','Northing','N','north','PROJ','ENh',2,'EPSG','9001');
 INSERT INTO "axis" VALUES('PROJ','3','Ellipsoidal height','h','up','PROJ','ENh',2,'EPSG','9001');
+
+---- Preferred hub for geodetic datum -----
+
+INSERT INTO "geodetic_datum_preferred_hub" VALUES('EPSG','1152','EPSG','6326'); -- WGS84 (G730) to WGS84
+INSERT INTO "geodetic_datum_preferred_hub" VALUES('EPSG','1153','EPSG','6326'); -- WGS84 (G873) to WGS84
+INSERT INTO "geodetic_datum_preferred_hub" VALUES('EPSG','1154','EPSG','6326'); -- WGS84 (G1150) to WGS84
+INSERT INTO "geodetic_datum_preferred_hub" VALUES('EPSG','1155','EPSG','6326'); -- WGS84 (G1674) to WGS84
+INSERT INTO "geodetic_datum_preferred_hub" VALUES('EPSG','1156','EPSG','6326'); -- WGS84 (G1762) to WGS84
+INSERT INTO "geodetic_datum_preferred_hub" VALUES('EPSG','1166','EPSG','6326'); -- WGS84 (Transit) to WGS84
+
+-- Consider all WGS84 related CRS are equivalent with an accuracy of 2m
+INSERT INTO "helmert_transformation" VALUES('PROJ','WGS84_TO_WGS84_G730','WGS 84 to WGS 84 (G730)','','Accuracy 2m','EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4326','EPSG','9053','EPSG','1262',2.0,0,0,0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0);
+INSERT INTO "helmert_transformation" VALUES('PROJ','WGS84_TO_WGS84_G873','WGS 84 to WGS 84 (G873)','','Accuracy 2m','EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4326','EPSG','9054','EPSG','1262',2.0,0,0,0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0);
+INSERT INTO "helmert_transformation" VALUES('PROJ','WGS84_TO_WGS84_G1150','WGS 84 to WGS 84 (G1150)','','Accuracy 2m','EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4326','EPSG','9055','EPSG','1262',2.0,0,0,0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0);
+INSERT INTO "helmert_transformation" VALUES('PROJ','WGS84_TO_WGS84_G1674','WGS 84 to WGS 84 (G1674)','','Accuracy 2m','EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4326','EPSG','9056','EPSG','1262',2.0,0,0,0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0);
+INSERT INTO "helmert_transformation" VALUES('PROJ','WGS84_TO_WGS84_G1762','WGS 84 to WGS 84 (G1762)','','Accuracy 2m','EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4326','EPSG','9057','EPSG','1262',2.0,0,0,0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0);
+INSERT INTO "helmert_transformation" VALUES('PROJ','WGS84_TO_WGS84_TRANSIT','WGS 84 to WGS 84 (Transit)','','Accuracy 2m','EPSG','9603','Geocentric translations (geog2D domain)','EPSG','4326','EPSG','8888','EPSG','1262',2.0,0,0,0,'EPSG','9001',NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'',0);
diff --git a/data/sql/proj_db_table_defs.sql b/data/sql/proj_db_table_defs.sql
index 255de827..6b38e736 100644
--- a/data/sql/proj_db_table_defs.sql
+++ b/data/sql/proj_db_table_defs.sql
@@ -126,6 +126,19 @@ FOR EACH ROW BEGIN
         WHERE EXISTS(SELECT 1 FROM area WHERE area.auth_name = NEW.area_of_use_auth_name AND area.code = NEW.area_of_use_code AND area.deprecated != 0) AND NEW.deprecated = 0;
 END;
 
+-- indicates that if there is no transformation from/into (src_auth_name, src_code),
+-- a research going through (hub_auth_name, hub_code) should be made
+CREATE TABLE geodetic_datum_preferred_hub(
+    src_auth_name TEXT NOT NULL CHECK (length(src_auth_name) >= 1),
+    src_code TEXT NOT NULL CHECK (length(src_code) >= 1),
+    hub_auth_name TEXT NOT NULL CHECK (length(hub_auth_name) >= 1),
+    hub_code TEXT NOT NULL CHECK (length(hub_code) >= 1),
+
+    CONSTRAINT unique_geodetic_datum_preferred_hub UNIQUE (src_auth_name, src_code, hub_auth_name, hub_code),
+    CONSTRAINT fk_geodetic_datum_preferred_hub_src FOREIGN KEY (src_auth_name, src_code) REFERENCES geodetic_datum(auth_name, code),
+    CONSTRAINT fk_geodetic_datum_preferred_hub_src FOREIGN KEY (hub_auth_name, hub_code) REFERENCES geodetic_datum(auth_name, code)
+);
+
 CREATE TABLE vertical_datum (
     auth_name TEXT NOT NULL CHECK (length(auth_name) >= 1),
     code TEXT NOT NULL CHECK (length(code) >= 1),
diff --git a/include/proj/io.hpp b/include/proj/io.hpp
index e439b9ef..12b3b111 100644
--- a/include/proj/io.hpp
+++ b/include/proj/io.hpp
@@ -1103,6 +1103,11 @@ class PROJ_GCC_DLL AuthorityFactory {
     PROJ_INTERNAL crs::CRSNNPtr
     createCoordinateReferenceSystem(const std::string &code,
                                     bool allowCompound) const;
+
+    PROJ_INTERNAL std::list<datum::GeodeticReferenceFrameNNPtr>
+    getPreferredHubGeodeticReferenceFrames(
+        const std::string &geodeticReferenceFrameCode) const;
+
     //! @endcond
 
   protected:
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index a1612339..f7ea385d 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -10171,6 +10171,7 @@ struct CoordinateOperationFactory::Private {
         const crs::CRSNNPtr &targetCRS;
         const CoordinateOperationContextNNPtr &context;
         bool inCreateOperationsWithDatumPivotAntiRecursion = false;
+        bool inCreateOperationsThroughPreferredHub = false;
 
         Context(const crs::CRSNNPtr &sourceCRSIn,
                 const crs::CRSNNPtr &targetCRSIn,
@@ -10195,6 +10196,12 @@ struct CoordinateOperationFactory::Private {
         const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst,
         Context &context);
 
+    static void createOperationsThroughPreferredHub(
+        std::vector<CoordinateOperationNNPtr> &res,
+        const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
+        const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst,
+        Context &context);
+
     static bool
     hasPerfectAccuracyResult(const std::vector<CoordinateOperationNNPtr> &res,
                              const Context &context);
@@ -12016,6 +12023,122 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
 
 // ---------------------------------------------------------------------------
 
+void CoordinateOperationFactory::Private::createOperationsThroughPreferredHub(
+    std::vector<CoordinateOperationNNPtr> &res, const crs::CRSNNPtr &sourceCRS,
+    const crs::CRSNNPtr &targetCRS, const crs::GeodeticCRS *geodSrc,
+    const crs::GeodeticCRS *geodDst, Private::Context &context) {
+
+    const auto &srcDatum = geodSrc->datum();
+    const auto &dstDatum = geodDst->datum();
+
+    if (!srcDatum || !dstDatum)
+        return;
+    const auto &srcDatumIds = srcDatum->identifiers();
+    const auto &dstDatumIds = dstDatum->identifiers();
+    if (srcDatumIds.empty() || dstDatumIds.empty())
+        return;
+
+    const auto &authFactory = context.context->getAuthorityFactory();
+
+    const auto srcAuthFactory = io::AuthorityFactory::create(
+        authFactory->databaseContext(), *(srcDatumIds.front()->codeSpace()));
+    const auto srcPreferredHubs =
+        srcAuthFactory->getPreferredHubGeodeticReferenceFrames(
+            srcDatumIds.front()->code());
+
+    const auto dstAuthFactory = io::AuthorityFactory::create(
+        authFactory->databaseContext(), *(dstDatumIds.front()->codeSpace()));
+    const auto dstPreferredHubs =
+        dstAuthFactory->getPreferredHubGeodeticReferenceFrames(
+            dstDatumIds.front()->code());
+    if (srcPreferredHubs.empty() && dstPreferredHubs.empty())
+        return;
+
+    // Currently if we have prefered hubs for both source and target, we
+    // will use only the one for target, arbitrarily... We could use boths
+    // but that would complicate things.
+    if (!srcPreferredHubs.empty() && dstPreferredHubs.empty()) {
+        std::vector<CoordinateOperationNNPtr> resTmp;
+        createOperationsThroughPreferredHub(resTmp, targetCRS, sourceCRS,
+                                            geodDst, geodSrc, context);
+        if (!resTmp.empty()) {
+            resTmp = applyInverse(resTmp);
+            res.insert(res.end(), resTmp.begin(), resTmp.end());
+        }
+        return;
+    }
+    assert(!dstPreferredHubs.empty());
+
+#ifdef DEBUG
+    EnterDebugLevel enterFunction;
+    debugTrace("createOperationsThroughPreferredHub(" +
+               objectAsStr(sourceCRS.get()) + "," +
+               objectAsStr(targetCRS.get()) + ")");
+#endif
+
+    struct AntiRecursionGuard {
+        Context &context;
+
+        explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) {
+            assert(!context.inCreateOperationsThroughPreferredHub);
+            context.inCreateOperationsThroughPreferredHub = true;
+        }
+
+        ~AntiRecursionGuard() {
+            context.inCreateOperationsThroughPreferredHub = false;
+        }
+    };
+    AntiRecursionGuard guard(context);
+
+    std::vector<crs::CRSNNPtr> candidatesIntermCRS;
+    for (const auto &datumHub : dstPreferredHubs) {
+        auto candidates =
+            findCandidateGeodCRSForDatum(authFactory, datumHub.get());
+        bool addedGeog2D = false;
+        for (const auto &intermCRS : candidates) {
+            auto geogCRS = dynamic_cast<crs::GeographicCRS *>(intermCRS.get());
+            if (geogCRS &&
+                geogCRS->coordinateSystem()->axisList().size() == 2) {
+                candidatesIntermCRS.emplace_back(intermCRS);
+                addedGeog2D = true;
+                break;
+            }
+        }
+        if (!addedGeog2D) {
+            candidatesIntermCRS.insert(candidatesIntermCRS.end(),
+                                       candidates.begin(), candidates.end());
+        }
+    }
+
+    const bool allowEmptyIntersection = true;
+    for (const auto &intermCRS : candidatesIntermCRS) {
+#ifdef DEBUG
+        EnterDebugLevel loopLevel;
+        debugTrace("try " + objectAsStr(sourceCRS.get()) + "->" +
+                   objectAsStr(intermCRS.get()) + "->" +
+                   objectAsStr(targetCRS.get()) + ")");
+        EnterDebugLevel loopLevel2;
+#endif
+        const auto opsFirst = createOperations(sourceCRS, intermCRS, context);
+        const auto opsLast = createOperations(intermCRS, targetCRS, context);
+        for (const auto &opFirst : opsFirst) {
+            for (const auto &opLast : opsLast) {
+                if (!opFirst->hasBallparkTransformation() ||
+                    !opLast->hasBallparkTransformation()) {
+                    try {
+                        res.emplace_back(
+                            ConcatenatedOperation::createComputeMetadata(
+                                {opFirst, opLast}, !allowEmptyIntersection));
+                    } catch (const InvalidOperationEmptyIntersection &) {
+                    }
+                }
+            }
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 static CoordinateOperationNNPtr
 createBallparkGeocentricTranslation(const crs::CRSNNPtr &sourceCRS,
                                     const crs::CRSNNPtr &targetCRS) {
@@ -12262,7 +12385,17 @@ CoordinateOperationFactory::Private::createOperations(
                     sourceCRS, targetCRS, context.context);
                 res.insert(res.end(), resWithIntermediate.begin(),
                            resWithIntermediate.end());
-                doFilterAndCheckPerfectOp = true;
+                doFilterAndCheckPerfectOp = !res.empty();
+
+                // If transforming from/to WGS84 (Gxxxx), try through 'neutral'
+                // WGS84
+                if (res.empty() && geodSrc && geodDst &&
+                    !context.inCreateOperationsThroughPreferredHub &&
+                    !context.inCreateOperationsWithDatumPivotAntiRecursion) {
+                    createOperationsThroughPreferredHub(
+                        res, sourceCRS, targetCRS, geodSrc, geodDst, context);
+                    doFilterAndCheckPerfectOp = !res.empty();
+                }
             }
         }
 
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index 1f40f1f0..4bf5427d 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -4982,6 +4982,29 @@ AuthorityFactory::createCompoundCRSFromExisting(
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+std::list<datum::GeodeticReferenceFrameNNPtr>
+AuthorityFactory::getPreferredHubGeodeticReferenceFrames(
+    const std::string &geodeticReferenceFrameCode) const {
+    std::list<datum::GeodeticReferenceFrameNNPtr> res;
+
+    const std::string sql("SELECT hub_auth_name, hub_code FROM "
+                          "geodetic_datum_preferred_hub WHERE "
+                          "src_auth_name = ? AND src_code = ?");
+    auto sqlRes = d->run(sql, {d->authority(), geodeticReferenceFrameCode});
+    for (const auto &row : sqlRes) {
+        const auto &auth_name = row[0];
+        const auto &code = row[1];
+        res.emplace_back(
+            d->createFactory(auth_name)->createGeodeticDatum(code));
+    }
+
+    return res;
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 FactoryException::FactoryException(const char *message) : Exception(message) {}
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index e9dec3f7..a3b49ba9 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -4495,6 +4495,73 @@ TEST(operation, geogCRS_to_geogCRS_context_NAD27_to_WGS84) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, geogCRS_to_geogCRS_context_NAD27_to_WGS84_G1762) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), std::string());
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+
+    auto authFactoryEPSG =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        // NAD27
+        authFactoryEPSG->createCoordinateReferenceSystem("4267"),
+        // WGS84 (G1762)
+        authFactoryEPSG->createCoordinateReferenceSystem("9057"), ctxt);
+    ASSERT_GE(list.size(), 78U);
+    EXPECT_EQ(list[0]->nameStr(),
+              "NAD27 to WGS 84 (33) + WGS 84 to WGS 84 (G1762)");
+    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=hgridshift +grids=ntv2_0.gsb "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
+    EXPECT_EQ(list[1]->nameStr(),
+              "NAD27 to WGS 84 (3) + WGS 84 to WGS 84 (G1762)");
+    EXPECT_EQ(list[2]->nameStr(),
+              "NAD27 to WGS 84 (79) + WGS 84 to WGS 84 (G1762)");
+    EXPECT_EQ(list[3]->nameStr(),
+              "NAD27 to WGS 84 (4) + WGS 84 to WGS 84 (G1762)");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, geogCRS_to_geogCRS_context_WGS84_G1674_to_WGS84_G1762) {
+    // Check that particular behaviour with WGS 84 (Gxxx) related to
+    // 'geodetic_datum_preferred_hub' table and custom no-op transformations
+    // between WGS 84 and WGS 84 (Gxxx) doesn't affect direct transformations
+    // to those realizations.
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), std::string());
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+
+    auto authFactoryEPSG =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        // WGS84 (G1674)
+        authFactoryEPSG->createCoordinateReferenceSystem("9056"),
+        // WGS84 (G1762)
+        authFactoryEPSG->createCoordinateReferenceSystem("9057"), ctxt);
+    ASSERT_EQ(list.size(), 1U);
+    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=cart +ellps=WGS84 "
+              "+step +proj=helmert +x=-0.004 +y=0.003 +z=0.004 +rx=0.00027 "
+              "+ry=-0.00027 +rz=0.00038 +s=-0.0069 "
+              "+convention=coordinate_frame "
+              "+step +inv +proj=cart +ellps=WGS84 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, vertCRS_to_geogCRS_context) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), "EPSG");