Merge pull request #1733 from rouault/remove_preferred_hub

createOperations(): remove the concept of geodetic_datum_preferred_hub

6 files changed, 71 insertions, 215 deletions

diff --git a/data/sql/customizations.sql b/data/sql/customizations.sql
index 6301bda9..8003c921 100644
--- a/data/sql/customizations.sql
+++ b/data/sql/customizations.sql
@@ -66,15 +66,6 @@ INSERT INTO "axis" VALUES('PROJ','1','Easting','E','east','PROJ','ENh',1,'EPSG',
 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);
diff --git a/data/sql/proj_db_table_defs.sql b/data/sql/proj_db_table_defs.sql
index a97e75a8..ee91aa3d 100644
--- a/data/sql/proj_db_table_defs.sql
+++ b/data/sql/proj_db_table_defs.sql
@@ -126,19 +126,6 @@ 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 17f0fea5..52ffebad 100644
--- a/include/proj/io.hpp
+++ b/include/proj/io.hpp
@@ -1064,6 +1064,8 @@ class PROJ_GCC_DLL AuthorityFactory {
         const std::vector<std::pair<std::string, std::string>>
             &intermediateCRSAuthCodes,
         ObjectType allowedIntermediateObjectType = ObjectType::CRS,
+        const std::vector<std::string> &allowedAuthorities =
+            std::vector<std::string>(),
         const metadata::ExtentPtr &intersectingExtent1 = nullptr,
         const metadata::ExtentPtr &intersectingExtent2 = nullptr) const;
 
@@ -1114,10 +1116,6 @@ class PROJ_GCC_DLL AuthorityFactory {
     createCoordinateReferenceSystem(const std::string &code,
                                     bool allowCompound) const;
 
-    PROJ_INTERNAL std::list<datum::GeodeticReferenceFrameNNPtr>
-    getPreferredHubGeodeticReferenceFrames(
-        const std::string &geodeticReferenceFrameCode) const;
-
     PROJ_INTERNAL std::vector<operation::CoordinateOperationNNPtr>
     getTransformationsForGeoid(const std::string &geoidName,
                                bool usePROJAlternativeGridNames) const;
diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index 67d4837f..4a1f6b90 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/scripts/reference_exported_symbols.txt
@@ -300,7 +300,7 @@ osgeo::proj::io::AuthorityFactory::createEllipsoid(std::string const&) const
 osgeo::proj::io::AuthorityFactory::createExtent(std::string const&) const
 osgeo::proj::io::AuthorityFactory::createFromCoordinateReferenceSystemCodes(std::string const&, std::string const&) const
 osgeo::proj::io::AuthorityFactory::createFromCoordinateReferenceSystemCodes(std::string const&, std::string const&, std::string const&, std::string const&, bool, bool, bool, bool, bool, std::shared_ptr<osgeo::proj::metadata::Extent> const&, std::shared_ptr<osgeo::proj::metadata::Extent> const&) const
-osgeo::proj::io::AuthorityFactory::createFromCRSCodesWithIntermediates(std::string const&, std::string const&, std::string const&, std::string const&, bool, bool, bool, std::vector<std::pair<std::string, std::string>, std::allocator<std::pair<std::string, std::string> > > const&, osgeo::proj::io::AuthorityFactory::ObjectType, std::shared_ptr<osgeo::proj::metadata::Extent> const&, std::shared_ptr<osgeo::proj::metadata::Extent> const&) const
+osgeo::proj::io::AuthorityFactory::createFromCRSCodesWithIntermediates(std::string const&, std::string const&, std::string const&, std::string const&, bool, bool, bool, std::vector<std::pair<std::string, std::string>, std::allocator<std::pair<std::string, std::string> > > const&, osgeo::proj::io::AuthorityFactory::ObjectType, std::vector<std::string, std::allocator<std::string> > const&, std::shared_ptr<osgeo::proj::metadata::Extent> const&, std::shared_ptr<osgeo::proj::metadata::Extent> const&) const
 osgeo::proj::io::AuthorityFactory::createGeodeticCRS(std::string const&) const
 osgeo::proj::io::AuthorityFactory::createGeodeticDatum(std::string const&) const
 osgeo::proj::io::AuthorityFactory::createGeographicCRS(std::string const&) const
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 11c10e74..937ba34f 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -10428,7 +10428,6 @@ struct CoordinateOperationFactory::Private {
         const metadata::ExtentPtr &extent2;
         const CoordinateOperationContextNNPtr &context;
         bool inCreateOperationsWithDatumPivotAntiRecursion = false;
-        bool inCreateOperationsThroughPreferredHub = false;
         bool inCreateOperationsGeogToVertWithAlternativeGeog = false;
         bool inCreateOperationsGeogToVertWithIntermediateVert = false;
         bool skipHorizontalTransformation = false;
@@ -10579,12 +10578,6 @@ 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);
@@ -11587,46 +11580,51 @@ CoordinateOperationFactory::Private::findsOpsInRegistryWithIntermediate(
 
             const auto authorities(getCandidateAuthorities(
                 authFactory, srcAuthName, targetAuthName));
-            for (const auto &authority : authorities) {
-                const auto tmpAuthFactory = io::AuthorityFactory::create(
-                    authFactory->databaseContext(),
-                    authority == "any" ? std::string() : authority);
+            assert(!authorities.empty());
 
-                io::AuthorityFactory::ObjectType intermediateObjectType =
-                    io::AuthorityFactory::ObjectType::CRS;
-
-                // If doing GeogCRS --> GeogCRS, only use GeogCRS as
-                // intermediate CRS
-                // Avoid weird behaviour when doing NAD83 -> NAD83(2011)
-                // that would go through NAVD88 otherwise.
-                if (context.context->getIntermediateCRS().empty() &&
-                    dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get()) &&
-                    dynamic_cast<const crs::GeographicCRS *>(targetCRS.get())) {
-                    intermediateObjectType =
-                        io::AuthorityFactory::ObjectType::GEOGRAPHIC_CRS;
-                }
-                auto res = tmpAuthFactory->createFromCRSCodesWithIntermediates(
-                    srcAuthName, srcCode, targetAuthName, targetCode,
-                    context.context->getUsePROJAlternativeGridNames(),
-                    context.context->getGridAvailabilityUse() ==
-                        CoordinateOperationContext::GridAvailabilityUse::
-                            DISCARD_OPERATION_IF_MISSING_GRID,
-                    context.context->getDiscardSuperseded(),
-                    context.context->getIntermediateCRS(),
-                    intermediateObjectType, context.extent1, context.extent2);
-                if (!res.empty()) {
+            const auto tmpAuthFactory = io::AuthorityFactory::create(
+                authFactory->databaseContext(),
+                (authFactory->getAuthority() == "any" || authorities.size() > 1)
+                    ? std::string()
+                    : authorities.front());
+
+            io::AuthorityFactory::ObjectType intermediateObjectType =
+                io::AuthorityFactory::ObjectType::CRS;
+
+            // If doing GeogCRS --> GeogCRS, only use GeogCRS as
+            // intermediate CRS
+            // Avoid weird behaviour when doing NAD83 -> NAD83(2011)
+            // that would go through NAVD88 otherwise.
+            if (context.context->getIntermediateCRS().empty() &&
+                dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get()) &&
+                dynamic_cast<const crs::GeographicCRS *>(targetCRS.get())) {
+                intermediateObjectType =
+                    io::AuthorityFactory::ObjectType::GEOGRAPHIC_CRS;
+            }
+            auto res = tmpAuthFactory->createFromCRSCodesWithIntermediates(
+                srcAuthName, srcCode, targetAuthName, targetCode,
+                context.context->getUsePROJAlternativeGridNames(),
+                context.context->getGridAvailabilityUse() ==
+                    CoordinateOperationContext::GridAvailabilityUse::
+                        DISCARD_OPERATION_IF_MISSING_GRID,
+                context.context->getDiscardSuperseded(),
+                context.context->getIntermediateCRS(), intermediateObjectType,
+                authFactory->getAuthority() != "any" && authorities.size() > 1
+                    ? authorities
+                    : std::vector<std::string>(),
+                context.extent1, context.extent2);
+            if (!res.empty()) {
 
-                    auto resFiltered =
-                        FilterResults(res, context.context, context.extent1,
-                                      context.extent2, false)
-                            .getRes();
+                auto resFiltered =
+                    FilterResults(res, context.context, context.extent1,
+                                  context.extent2, false)
+                        .getRes();
 #ifdef TRACE_CREATE_OPERATIONS
-                    logTrace("filtering reduced from " +
-                             toString(static_cast<int>(res.size())) + " to " +
-                             toString(static_cast<int>(resFiltered.size())));
+                logTrace("filtering reduced from " +
+                         toString(static_cast<int>(res.size())) + " to " +
+                         toString(static_cast<int>(resFiltered.size())));
 #endif
-                    return resFiltered;
-                }
+                return resFiltered;
             }
         }
     }
@@ -12446,118 +12444,6 @@ 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 TRACE_CREATE_OPERATIONS
-    ENTER_BLOCK("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, geodDst, 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());
-        }
-    }
-
-    for (const auto &intermCRS : candidatesIntermCRS) {
-#ifdef TRACE_CREATE_OPERATIONS
-        ENTER_BLOCK("try " + objectAsStr(sourceCRS.get()) + "->" +
-                    objectAsStr(intermCRS.get()) + "->" +
-                    objectAsStr(targetCRS.get()) + ")");
-#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) {
@@ -12890,16 +12776,6 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
         res.insert(res.end(), resWithIntermediate.begin(),
                    resWithIntermediate.end());
         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();
-        }
     }
 
     if (doFilterAndCheckPerfectOp) {
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index d9917996..309e778a 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -3775,6 +3775,11 @@ static bool useIrrelevantPivot(const operation::CoordinateOperationNNPtr &op,
  * @param allowedIntermediateObjectType Restrict the type of the intermediate
  * object considered.
  * Only ObjectType::CRS and ObjectType::GEOGRAPHIC_CRS supported currently
+ * @param allowedAuthorities One or several authority name allowed for the two
+ * coordinate operations that are going to be searched. When this vector is
+ * no empty, it overrides the authority of this object. This is useful for
+ * example when the coordinate operations to chain belong to two different
+ * allowed authorities.
  * @param intersectingExtent1 Optional extent that the resulting operations
  * must intersect.
  * @param intersectingExtent2 Optional extent that the resulting operations
@@ -3793,6 +3798,7 @@ AuthorityFactory::createFromCRSCodesWithIntermediates(
     const std::vector<std::pair<std::string, std::string>>
         &intermediateCRSAuthCodes,
     ObjectType allowedIntermediateObjectType,
+    const std::vector<std::string> &allowedAuthorities,
     const metadata::ExtentPtr &intersectingExtent1,
     const metadata::ExtentPtr &intersectingExtent2) const {
 
@@ -3887,6 +3893,27 @@ AuthorityFactory::createFromCRSCodesWithIntermediates(
         "AND v1.deprecated = 0 AND v2.deprecated = 0 "
         "AND intersects_bbox(south_lat1, west_lon1, north_lat1, east_lon1, "
         "south_lat2, west_lon2, north_lat2, east_lon2) = 1 ");
+    if (!allowedAuthorities.empty()) {
+        additionalWhere += "AND v1.auth_name IN (";
+        for (size_t i = 0; i < allowedAuthorities.size(); i++) {
+            if (i > 0)
+                additionalWhere += ',';
+            additionalWhere += '?';
+        }
+        additionalWhere += ") AND v2.auth_name IN (";
+        for (size_t i = 0; i < allowedAuthorities.size(); i++) {
+            if (i > 0)
+                additionalWhere += ',';
+            additionalWhere += '?';
+        }
+        additionalWhere += ')';
+        for (const auto &allowedAuthority : allowedAuthorities) {
+            params.emplace_back(allowedAuthority);
+        }
+        for (const auto &allowedAuthority : allowedAuthorities) {
+            params.emplace_back(allowedAuthority);
+        }
+    }
     if (d->hasAuthorityRestriction()) {
         additionalWhere += "AND v1.auth_name = ? AND v2.auth_name = ? ";
         params.emplace_back(d->authority());
@@ -5357,29 +5384,6 @@ 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
 std::vector<operation::CoordinateOperationNNPtr>
 AuthorityFactory::getTransformationsForGeoid(
     const std::string &geoidName, bool usePROJAlternativeGridNames) const {