Merge pull request #1394 from rouault/improve_createoperations_missing_ids

createOperations(): improve behaviour with input CRS from WKT that lacks intermediate IDs (fixes #1343)

2 files changed, 298 insertions, 20 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 24e5f8ab..15532a89 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -11311,16 +11311,37 @@ static bool hasIdentifiers(const CoordinateOperationNNPtr &op) {
 
 static std::vector<crs::CRSNNPtr>
 findCandidateGeodCRSForDatum(const io::AuthorityFactoryPtr &authFactory,
-                             const datum::GeodeticReferenceFramePtr &datum) {
+                             const datum::GeodeticReferenceFrame *datum) {
     std::vector<crs::CRSNNPtr> candidates;
-    for (const auto &id : datum->identifiers()) {
-        const auto &authName = *(id->codeSpace());
-        const auto &code = id->code();
-        if (!authName.empty()) {
-            auto l_candidates = authFactory->createGeodeticCRSFromDatum(
-                authName, code, std::string());
-            for (const auto &candidate : l_candidates) {
-                candidates.emplace_back(candidate);
+    assert(datum);
+    const auto &ids = datum->identifiers();
+    const auto &datumName = datum->nameStr();
+    if (!ids.empty()) {
+        for (const auto &id : ids) {
+            const auto &authName = *(id->codeSpace());
+            const auto &code = id->code();
+            if (!authName.empty()) {
+                auto l_candidates = authFactory->createGeodeticCRSFromDatum(
+                    authName, code, std::string());
+                for (const auto &candidate : l_candidates) {
+                    candidates.emplace_back(candidate);
+                }
+            }
+        }
+    } else if (datumName != "unknown" && datumName != "unnamed") {
+        auto matches = authFactory->createObjectsFromName(
+            datumName,
+            {io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME}, false,
+            2);
+        if (matches.size() == 1) {
+            const auto &match = matches.front();
+            if (datum->_isEquivalentTo(
+                    match.get(), util::IComparable::Criterion::EQUIVALENT) &&
+                !match->identifiers().empty()) {
+                return findCandidateGeodCRSForDatum(
+                    authFactory,
+                    dynamic_cast<const datum::GeodeticReferenceFrame *>(
+                        match.get()));
             }
         }
     }
@@ -11362,9 +11383,9 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
 
     const auto &authFactory = context.context->getAuthorityFactory();
     const auto candidatesSrcGeod(
-        findCandidateGeodCRSForDatum(authFactory, geodSrc->datum()));
+        findCandidateGeodCRSForDatum(authFactory, geodSrc->datum().get()));
     const auto candidatesDstGeod(
-        findCandidateGeodCRSForDatum(authFactory, geodDst->datum()));
+        findCandidateGeodCRSForDatum(authFactory, geodDst->datum().get()));
 
     auto createTransformations = [&](const crs::CRSNNPtr &candidateSrcGeod,
                                      const crs::CRSNNPtr &candidateDstGeod,
@@ -11605,12 +11626,8 @@ CoordinateOperationFactory::Private::createOperations(
             // but transformations are only available between their
             // corresponding geocentric CRS.
             const auto &srcDatum = geodSrc->datum();
-            const bool srcHasDatumWithId =
-                srcDatum && !srcDatum->identifiers().empty();
             const auto &dstDatum = geodDst->datum();
-            const bool dstHasDatumWithId =
-                dstDatum && !dstDatum->identifiers().empty();
-            if (srcHasDatumWithId && dstHasDatumWithId &&
+            if (srcDatum != nullptr && dstDatum != nullptr &&
                 !srcDatum->_isEquivalentTo(
                     dstDatum.get(), util::IComparable::Criterion::EQUIVALENT)) {
                 createOperationsWithDatumPivot(res, sourceCRS, targetCRS,
@@ -11955,6 +11972,31 @@ CoordinateOperationFactory::Private::createOperations(
     // A bit odd case as we are comparing apples to oranges, but in case
     // the vertical unit differ, do something useful.
     if (vertSrc && geogDst) {
+
+        if (vertSrc->identifiers().empty()) {
+            const auto &authFactory = context.context->getAuthorityFactory();
+            const auto &vertSrcName = vertSrc->nameStr();
+            if (authFactory != nullptr && vertSrcName != "unnamed" &&
+                vertSrcName != "unknown") {
+                auto matches = authFactory->createObjectsFromName(
+                    vertSrcName,
+                    {io::AuthorityFactory::ObjectType::VERTICAL_CRS}, false, 2);
+                if (matches.size() == 1) {
+                    const auto &match = matches.front();
+                    if (vertSrc->_isEquivalentTo(
+                            match.get(),
+                            util::IComparable::Criterion::EQUIVALENT) &&
+                        !match->identifiers().empty()) {
+                        return createOperations(
+                            NN_NO_CHECK(
+                                util::nn_dynamic_pointer_cast<crs::VerticalCRS>(
+                                    match)),
+                            targetCRS, context);
+                    }
+                }
+            }
+        }
+
         const double convSrc =
             vertSrc->coordinateSystem()->axisList()[0]->unit().conversionToSI();
         double convDst = 1.0;
@@ -12300,6 +12342,67 @@ CoordinateOperationFactory::Private::createOperations(
 
 // ---------------------------------------------------------------------------
 
+static crs::CRSNNPtr
+getResolvedCRS(const crs::CRSNNPtr &crs,
+               const CoordinateOperationContextNNPtr &context) {
+    const auto &authFactory = context->getAuthorityFactory();
+
+    auto projectedCrs = dynamic_cast<crs::ProjectedCRS *>(crs.get());
+    if (projectedCrs && authFactory) {
+        const auto &ids = projectedCrs->identifiers();
+        if (!ids.empty() && projectedCrs->baseCRS()->identifiers().empty()) {
+            const auto tmpAuthFactory = io::AuthorityFactory::create(
+                authFactory->databaseContext(), *ids.front()->codeSpace());
+            try {
+                crs::CRSNNPtr resolvedCrs(
+                    tmpAuthFactory->createProjectedCRS(ids.front()->code()));
+                if (resolvedCrs->_isEquivalentTo(
+                        crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
+                    return resolvedCrs;
+                }
+            } catch (const std::exception &) {
+            }
+        }
+    }
+
+    auto compoundCrs = dynamic_cast<crs::CompoundCRS *>(crs.get());
+    // If we get a CompoundCRS that has an EPSG code, but whose component CRS
+    // lack one, typically from WKT2, this might be an issue to get proper
+    // results in createOperations(), so import the CompoundCRS from the
+    // registry, and if equivalent to the original one, then use the version
+    // from the registry.
+    if (compoundCrs && authFactory) {
+        const auto &ids = compoundCrs->identifiers();
+        if (!ids.empty()) {
+            const auto &components = compoundCrs->componentReferenceSystems();
+            bool hasMissingId = false;
+            for (const auto &comp : components) {
+                if (comp->identifiers().empty()) {
+                    hasMissingId = true;
+                    break;
+                }
+            }
+            if (hasMissingId) {
+                const auto tmpAuthFactory = io::AuthorityFactory::create(
+                    authFactory->databaseContext(), *ids.front()->codeSpace());
+                try {
+                    crs::CRSNNPtr resolvedCrs(
+                        tmpAuthFactory->createCompoundCRS(ids.front()->code()));
+                    if (resolvedCrs->_isEquivalentTo(
+                            crs.get(),
+                            util::IComparable::Criterion::EQUIVALENT)) {
+                        return resolvedCrs;
+                    }
+                } catch (const std::exception &) {
+                }
+            }
+        }
+    }
+    return crs;
+}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Find a list of CoordinateOperation from sourceCRS to targetCRS.
  *
  * The operations are sorted with the most relevant ones first: by
@@ -12328,10 +12431,14 @@ CoordinateOperationFactory::createOperations(
     auto l_sourceCRS = srcBoundCRS ? NN_NO_CHECK(srcBoundCRS) : sourceCRS;
     auto l_targetCRS = targetBoundCRS ? NN_NO_CHECK(targetBoundCRS) : targetCRS;
 
-    Private::Context contextPrivate(sourceCRS, targetCRS, context);
-    return filterAndSort(
-        Private::createOperations(l_sourceCRS, l_targetCRS, contextPrivate),
-        context, l_sourceCRS, l_targetCRS);
+    auto l_resolvedSourceCRS = getResolvedCRS(l_sourceCRS, context);
+    auto l_resolvedTargetCRS = getResolvedCRS(l_targetCRS, context);
+    Private::Context contextPrivate(l_resolvedSourceCRS, l_resolvedTargetCRS,
+                                    context);
+    return filterAndSort(Private::createOperations(l_resolvedSourceCRS,
+                                                   l_resolvedTargetCRS,
+                                                   contextPrivate),
+                         context, l_resolvedSourceCRS, l_resolvedTargetCRS);
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 4e3ceb56..32acbf22 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -4746,6 +4746,52 @@ TEST(operation, geogCRS_to_geogCRS_3D) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, geogCRS_without_id_to_geogCRS_3D_context) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    auto src =
+        authFactory->createCoordinateReferenceSystem("4289"); // Amersfoort
+    auto dst =
+        authFactory->createCoordinateReferenceSystem("4937"); // ETRS89 3D
+    auto list =
+        CoordinateOperationFactory::create()->createOperations(src, dst, ctxt);
+    ASSERT_GE(list.size(), 1U);
+    auto wkt2 = "GEOGCRS[\"unnamed\",\n"
+                "    DATUM[\"Amersfoort\",\n"
+                "        ELLIPSOID[\"Bessel 1841\",6377397.155,299.1528128,\n"
+                "            LENGTHUNIT[\"metre\",1]]],\n"
+                "    PRIMEM[\"Greenwich\",0,\n"
+                "        ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+                "    CS[ellipsoidal,2],\n"
+                "        AXIS[\"geodetic latitude (Lat)\",north,\n"
+                "            ORDER[1],\n"
+                "            ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
+                "        AXIS[\"geodetic longitude (Lon)\",east,\n"
+                "            ORDER[2],\n"
+                "            ANGLEUNIT[\"degree\",0.0174532925199433]],"
+                "    USAGE[\n"
+                "        SCOPE[\"unknown\"],\n"
+                "        AREA[\"Netherlands - onshore\"],\n"
+                "        BBOX[50.75,3.2,53.7,7.22]]]\n";
+
+    auto obj = WKTParser().createFromWKT(wkt2);
+    auto src_from_wkt2 = nn_dynamic_pointer_cast<CRS>(obj);
+    ASSERT_TRUE(src_from_wkt2 != nullptr);
+    auto list2 = CoordinateOperationFactory::create()->createOperations(
+        NN_NO_CHECK(src_from_wkt2), dst, ctxt);
+    ASSERT_GE(list.size(), list2.size());
+    for (size_t i = 0; i < list.size(); i++) {
+        const auto &op = list[i];
+        const auto &op2 = list2[i];
+        EXPECT_TRUE(
+            op->isEquivalentTo(op2.get(), IComparable::Criterion::EQUIVALENT))
+            << op->nameStr() << " " << op2->nameStr();
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, geocentricCRS_to_geogCRS_same_datum) {
 
     auto op = CoordinateOperationFactory::create()->createOperation(
@@ -5172,6 +5218,50 @@ TEST(operation, projCRS_to_geogCRS) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, projCRS_no_id_to_geogCRS_context) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    auto src = authFactory->createCoordinateReferenceSystem(
+        "28992"); // Amersfoort / RD New
+    auto dst =
+        authFactory->createCoordinateReferenceSystem("4258"); // ETRS89 2D
+    auto list =
+        CoordinateOperationFactory::create()->createOperations(src, dst, ctxt);
+    ASSERT_GE(list.size(), 1U);
+    auto wkt2 =
+        "PROJCRS[\"unknown\",\n"
+        "    BASEGEOGCRS[\"Amersfoort\",\n"
+        "       DATUM[\"Amersfoort\",\n"
+        "            ELLIPSOID[\"Bessel 1841\",6377397.155,299.1528128]]],\n"
+        "    CONVERSION[\"unknown\",\n"
+        "        METHOD[\"Oblique Stereographic\"],\n"
+        "        PARAMETER[\"Latitude of natural origin\",52.1561605555556],\n"
+        "        PARAMETER[\"Longitude of natural origin\",5.38763888888889],\n"
+        "        PARAMETER[\"Scale factor at natural origin\",0.9999079],\n"
+        "        PARAMETER[\"False easting\",155000],\n"
+        "        PARAMETER[\"False northing\",463000]],\n"
+        "    CS[Cartesian,2],\n"
+        "        AXIS[\"(E)\",east],\n"
+        "        AXIS[\"(N)\",north],\n"
+        "    LENGTHUNIT[\"metre\",1],\n"
+        "    ID[\"EPSG\",28992]]";
+    auto obj = WKTParser().createFromWKT(wkt2);
+    auto src_from_wkt2 = nn_dynamic_pointer_cast<CRS>(obj);
+    ASSERT_TRUE(src_from_wkt2 != nullptr);
+    auto list2 = CoordinateOperationFactory::create()->createOperations(
+        NN_NO_CHECK(src_from_wkt2), dst, ctxt);
+    ASSERT_GE(list.size(), list2.size() - 1);
+    for (size_t i = 0; i < list.size(); i++) {
+        const auto &op = list[i];
+        const auto &op2 = list2[i];
+        EXPECT_TRUE(
+            op->isEquivalentTo(op2.get(), IComparable::Criterion::EQUIVALENT));
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, projCRS_to_projCRS) {
 
     auto op = CoordinateOperationFactory::create()->createOperation(
@@ -6487,6 +6577,87 @@ TEST(operation, compoundCRS_to_geogCRS_3D_context) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, compoundCRS_from_WKT2_to_geogCRS_3D_context) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    auto src = authFactory->createCoordinateReferenceSystem(
+        "7415"); // Amersfoort / RD New + NAP height
+    auto dst =
+        authFactory->createCoordinateReferenceSystem("4937"); // ETRS89 3D
+    auto list =
+        CoordinateOperationFactory::create()->createOperations(src, dst, ctxt);
+    ASSERT_GE(list.size(), 1U);
+    auto wkt2 = src->exportToWKT(
+        WKTFormatter::create(WKTFormatter::Convention::WKT2_2018).get());
+    auto obj = WKTParser().createFromWKT(wkt2);
+    auto src_from_wkt2 = nn_dynamic_pointer_cast<CRS>(obj);
+    ASSERT_TRUE(src_from_wkt2 != nullptr);
+    auto list2 = CoordinateOperationFactory::create()->createOperations(
+        NN_NO_CHECK(src_from_wkt2), dst, ctxt);
+    ASSERT_GE(list.size(), list2.size());
+    for (size_t i = 0; i < list.size(); i++) {
+        const auto &op = list[i];
+        const auto &op2 = list2[i];
+        EXPECT_TRUE(
+            op->isEquivalentTo(op2.get(), IComparable::Criterion::EQUIVALENT));
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, compoundCRS_from_WKT2_no_id_to_geogCRS_3D_context) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    auto src = authFactory->createCoordinateReferenceSystem(
+        "7415"); // Amersfoort / RD New + NAP height
+    auto dst =
+        authFactory->createCoordinateReferenceSystem("4937"); // ETRS89 3D
+    auto list =
+        CoordinateOperationFactory::create()->createOperations(src, dst, ctxt);
+    ASSERT_GE(list.size(), 1U);
+    auto wkt2 =
+        "COMPOUNDCRS[\"unknown\",\n"
+        "  PROJCRS[\"unknown\",\n"
+        "    BASEGEOGCRS[\"Amersfoort\",\n"
+        "            DATUM[\"Amersfoort\",\n"
+        "                ELLIPSOID[\"Bessel "
+        "1841\",6377397.155,299.1528128]]],\n"
+        "    CONVERSION[\"unknown\",\n"
+        "        METHOD[\"Oblique Stereographic\"],\n"
+        "        PARAMETER[\"Latitude of natural origin\",52.1561605555556],\n"
+        "        PARAMETER[\"Longitude of natural origin\",5.38763888888889],\n"
+        "        PARAMETER[\"Scale factor at natural origin\",0.9999079],\n"
+        "        PARAMETER[\"False easting\",155000],\n"
+        "        PARAMETER[\"False northing\",463000]],\n"
+        "    CS[Cartesian,2],\n"
+        "        AXIS[\"(E)\",east],\n"
+        "        AXIS[\"(N)\",north],\n"
+        "    LENGTHUNIT[\"metre\",1]],\n"
+        "  VERTCRS[\"NAP height\",\n"
+        "    VDATUM[\"Normaal Amsterdams Peil\"],\n"
+        "    CS[vertical,1],\n"
+        "        AXIS[\"gravity-related height (H)\",up,\n"
+        "            LENGTHUNIT[\"metre\",1]]]]";
+    auto obj = WKTParser().createFromWKT(wkt2);
+    auto src_from_wkt2 = nn_dynamic_pointer_cast<CRS>(obj);
+    ASSERT_TRUE(src_from_wkt2 != nullptr);
+    auto list2 = CoordinateOperationFactory::create()->createOperations(
+        NN_NO_CHECK(src_from_wkt2), dst, ctxt);
+    ASSERT_GE(list.size(), list2.size() - 1);
+    for (size_t i = 0; i < list.size(); i++) {
+        const auto &op = list[i];
+        const auto &op2 = list2[i];
+        EXPECT_TRUE(
+            op->isEquivalentTo(op2.get(), IComparable::Criterion::EQUIVALENT));
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, boundCRS_to_compoundCRS) {
     auto objSrc = PROJStringParser().createFromPROJString(
         "+proj=longlat +ellps=GRS67 +nadgrids=@foo.gsb +type=crs");