findsOpsInRegistryWithIntermediate(): restrict to using known source/…

…target CRS that have the same originating authority

Fixes https://lists.osgeo.org/pipermail/proj/2024-September/011531.html
i.e.
```
projinfo -s EPSG:4269 -t EPSG:6318 --3d --spatial-test intersects
```

The fix consists in making sure that we recognize that the 3d-promoted
object EPSG:4269 (NAD83(86)) is still linked to EPSG, and thus discard
ESRI 3D objects in findsOpsInRegistryWithIntermediate()

Fixes a "regression" of OSGeo#4244
(one could argue which of the results is the best, given that NAD83(86)
as a 3D geographic CRS has no solid foundation)

include/proj/crs.hpp

@@ -165,6 +165,9 @@ class PROJ_GCC_DLL CRS : public common::ObjectUsage,
     getResolvedCRS(const CRSNNPtr &crs,
                    const io::AuthorityFactoryPtr &authFactory,
                    metadata::ExtentPtr &extentOut);
+
+    PROJ_INTERNAL std::string getOriginatingAuthName() const;
+
     //! @endcond
 
   protected:

src/iso19111/crs.cpp

@@ -92,6 +92,8 @@ NS_PROJ_START
 
 namespace crs {
 
+constexpr const char *PROMOTED_TO_3D_PRELUDE = "Promoted to 3D from ";
+
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
@@ -1257,6 +1259,40 @@ CRS::getNonDeprecated(const io::DatabaseContextNNPtr &dbContext) const {
 
 // ---------------------------------------------------------------------------
 
+//! @cond Doxygen_Suppress
+
+/** \brief Return the authority name to which this object is registered, or
+ * has an indirect provenance.
+ *
+ * Typically this method called on EPSG:4269 (NAD83) promoted to 3D will return
+ * "EPSG".
+ *
+ * Returns empty string if more than an authority or no originating authority is
+ * found.
+ */
+std::string CRS::getOriginatingAuthName() const {
+    const auto &ids = identifiers();
+    if (ids.size() == 1) {
+        return *(ids[0]->codeSpace());
+    }
+    if (ids.size() > 1) {
+        return std::string();
+    }
+    const auto &l_remarks = remarks();
+    if (starts_with(l_remarks, PROMOTED_TO_3D_PRELUDE)) {
+        const auto pos = l_remarks.find(':');
+        if (pos != std::string::npos) {
+            return l_remarks.substr(strlen(PROMOTED_TO_3D_PRELUDE),
+                                    pos - strlen(PROMOTED_TO_3D_PRELUDE));
+        }
+    }
+    return std::string();
+}
+
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
 /** \brief Return a variant of this CRS "promoted" to a 3D one, if not already
  * the case.
  *
@@ -1314,7 +1350,7 @@ CRSNNPtr CRS::promoteTo3D(const std::string &newName,
         const auto &l_identifiers = identifiers();
         const auto &l_remarks = remarks();
         if (l_identifiers.size() == 1) {
-            std::string remarks("Promoted to 3D from ");
+            std::string remarks(PROMOTED_TO_3D_PRELUDE);
             remarks += *(l_identifiers[0]->codeSpace());
             remarks += ':';
             remarks += l_identifiers[0]->code();

src/iso19111/operation/coordinateoperationfactory.cpp

@@ -2001,12 +2001,21 @@ CoordinateOperationFactory::Private::findsOpsInRegistryWithIntermediate(
         auto geodSrc = dynamic_cast<crs::GeodeticCRS *>(sourceCRS.get());
         auto geodDst = dynamic_cast<crs::GeodeticCRS *>(targetCRS.get());
         if (geodSrc) {
+            const std::string originatingAuthSrc =
+                geodSrc->getOriginatingAuthName();
             const auto dbContext = authFactory->databaseContext().as_nullable();
             const auto candidatesSrcGeod(findCandidateGeodCRSForDatum(
                 authFactory, geodSrc, geodSrc->datumNonNull(dbContext)));
             std::vector<CoordinateOperationNNPtr> res;
             for (const auto &candidateSrcGeod : candidatesSrcGeod) {
-                if (candidateSrcGeod->coordinateSystem()->axisList().size() ==
+                // Restrict to using only objects that have the same authority
+                // as geodSrc.
+                const auto &candidateIds = candidateSrcGeod->identifiers();
+                if ((originatingAuthSrc.empty() ||
+                     (candidateIds.size() == 1 &&
+                      *(candidateIds[0]->codeSpace()) == originatingAuthSrc) ||
+                     candidateIds.size() > 1) &&
+                    candidateSrcGeod->coordinateSystem()->axisList().size() ==
                         geodSrc->coordinateSystem()->axisList().size() &&
                     ((dynamic_cast<crs::GeographicCRS *>(sourceCRS.get()) !=
                       nullptr) ==
@@ -2022,6 +2031,7 @@ CoordinateOperationFactory::Private::findsOpsInRegistryWithIntermediate(
                             continue;
                         }
                     }
+
                     const auto opsWithIntermediate =
                         findsOpsInRegistryWithIntermediate(
                             candidateSrcGeod, targetCRS, context,

test/unit/test_operationfactory.cpp

@@ -1533,6 +1533,51 @@ TEST(operation, geogCRS_without_id_to_geogCRS_3D_context) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, geogCRS_promoted_to_3D_to_geogCRS_3D_context) {
+    auto dbContext = DatabaseContext::create();
+    auto authFactory = AuthorityFactory::create(dbContext, std::string());
+    auto authFactoryEPSG = AuthorityFactory::create(dbContext, "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    {
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            authFactoryEPSG->createCoordinateReferenceSystem("4269")
+                ->promoteTo3D(std::string(),
+                              dbContext), // NAD83 (86) promoted to 3D
+            authFactoryEPSG->createCoordinateReferenceSystem(
+                "6319"), // NAD83 (2011) 3D
+            ctxt);
+        ASSERT_GE(list.size(), 1U);
+        // Check we don't get an ESRI operation
+        EXPECT_EQ(list[0]->nameStr(), "NAD83 to NAD83(HARN) (47) + "
+                                      "NAD83(HARN) to NAD83(FBN) (1) + "
+                                      "NAD83(FBN) to NAD83(NSRS2007) (1) + "
+                                      "NAD83(NSRS2007) to NAD83(2011) (1)");
+    }
+    {
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            authFactoryEPSG->createCoordinateReferenceSystem(
+                "6319"), // NAD83 (2011) 3D
+            authFactoryEPSG->createCoordinateReferenceSystem("4269")
+                ->promoteTo3D(std::string(),
+                              dbContext), // NAD83 (86) promoted to 3D
+            ctxt);
+        ASSERT_GE(list.size(), 1U);
+        // Check we don't get an ESRI operation
+        EXPECT_EQ(list[0]->nameStr(),
+                  "Inverse of NAD83(NSRS2007) to NAD83(2011) (1) + "
+                  "Inverse of NAD83(FBN) to NAD83(NSRS2007) (1) + "
+                  "Inverse of NAD83(HARN) to NAD83(FBN) (1) + "
+                  "Inverse of NAD83 to NAD83(HARN) (47)");
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 static GeodeticCRSNNPtr createGeocentricDatumWGS84() {
     PropertyMap propertiesCRS;
     propertiesCRS.set(Identifier::CODESPACE_KEY, "EPSG")