createOperations(): try to recover extent of CRS from the database when they are missing, especially for compound CRS. Helps having shorter/more relevant results

3 files changed, 211 insertions, 35 deletions

diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 6f99cd0c..9796221d 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -464,6 +464,31 @@ static const metadata::ExtentPtr &getExtent(const crs::CRSNNPtr &crs) {
     return nullExtent;
 }
 
+static const metadata::ExtentPtr
+getExtentPossiblySynthetized(const crs::CRSNNPtr &crs, bool &approxOut) {
+    const auto &rawExtent(getExtent(crs));
+    approxOut = false;
+    if (rawExtent)
+        return rawExtent;
+    const auto compoundCRS = dynamic_cast<const crs::CompoundCRS *>(crs.get());
+    if (compoundCRS) {
+        // For a compoundCRS, take the intersection of the extent of its
+        // components.
+        const auto &components = compoundCRS->componentReferenceSystems();
+        metadata::ExtentPtr extent;
+        approxOut = true;
+        for (const auto &component : components) {
+            const auto &componentExtent(getExtent(component));
+            if (extent && componentExtent)
+                extent = extent->intersection(NN_NO_CHECK(componentExtent));
+            else if (componentExtent)
+                extent = componentExtent;
+        }
+        return extent;
+    }
+    return rawExtent;
+}
+
 // ---------------------------------------------------------------------------
 
 static metadata::ExtentPtr
@@ -14033,57 +14058,129 @@ void CoordinateOperationFactory::Private::createOperationsBoundToCompound(
 
 static crs::CRSNNPtr
 getResolvedCRS(const crs::CRSNNPtr &crs,
-               const CoordinateOperationContextNNPtr &context) {
+               const CoordinateOperationContextNNPtr &context,
+               metadata::ExtentPtr &extentOut) {
     const auto &authFactory = context->getAuthorityFactory();
+    const auto &ids = crs->identifiers();
+    const auto &name = crs->nameStr();
+
+    bool approxExtent;
+    extentOut = getExtentPossiblySynthetized(crs, approxExtent);
+
+    // We try to "identify" the provided CRS with the ones of the database,
+    // but in a more restricted way that what identify() does.
+    // If we get a match from id in priority, and from name as a fallback, and
+    // that they are equivalent to the input CRS, then use the identified CRS.
+    // Even if they aren't equivalent, we update extentOut with the one of the
+    // identified CRS if our input one is absent/not reliable.
+
+    const auto tryToIdentifyByName = [&crs, &name, &authFactory, approxExtent,
+                                      &extentOut](
+        io::AuthorityFactory::ObjectType objectType) {
+        if (name != "unknown" && name != "unnamed") {
+            auto matches = authFactory->createObjectsFromName(
+                name, {objectType}, false, 2);
+            if (matches.size() == 1) {
+                const auto match =
+                    util::nn_static_pointer_cast<crs::CRS>(matches.front());
+                if (approxExtent || !extentOut) {
+                    extentOut = getExtent(match);
+                }
+                if (match->isEquivalentTo(
+                        crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
+                    return match;
+                }
+            }
+        }
+        return crs;
+    };
+
+    auto geogCRS = dynamic_cast<crs::GeographicCRS *>(crs.get());
+    if (geogCRS && authFactory) {
+        if (!ids.empty()) {
+            const auto tmpAuthFactory = io::AuthorityFactory::create(
+                authFactory->databaseContext(), *ids.front()->codeSpace());
+            try {
+                auto resolvedCrs(
+                    tmpAuthFactory->createGeographicCRS(ids.front()->code()));
+                if (approxExtent || !extentOut) {
+                    extentOut = getExtent(resolvedCrs);
+                }
+                if (resolvedCrs->isEquivalentTo(
+                        crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
+                    return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs);
+                }
+            } catch (const std::exception &) {
+            }
+        } else {
+            return tryToIdentifyByName(
+                geogCRS->coordinateSystem()->axisList().size() == 2
+                    ? io::AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS
+                    : io::AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS);
+        }
+    }
 
     auto projectedCrs = dynamic_cast<crs::ProjectedCRS *>(crs.get());
     if (projectedCrs && authFactory) {
-        const auto &ids = projectedCrs->identifiers();
-        if (!ids.empty() && projectedCrs->baseCRS()->identifiers().empty()) {
+        if (!ids.empty()) {
             const auto tmpAuthFactory = io::AuthorityFactory::create(
                 authFactory->databaseContext(), *ids.front()->codeSpace());
             try {
                 auto resolvedCrs(
                     tmpAuthFactory->createProjectedCRS(ids.front()->code()));
+                if (approxExtent || !extentOut) {
+                    extentOut = getExtent(resolvedCrs);
+                }
                 if (resolvedCrs->isEquivalentTo(
                         crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
                     return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs);
                 }
             } catch (const std::exception &) {
             }
+        } else {
+            return tryToIdentifyByName(
+                io::AuthorityFactory::ObjectType::PROJECTED_CRS);
         }
     }
 
     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;
+            const auto tmpAuthFactory = io::AuthorityFactory::create(
+                authFactory->databaseContext(), *ids.front()->codeSpace());
+            try {
+                auto resolvedCrs(
+                    tmpAuthFactory->createCompoundCRS(ids.front()->code()));
+                if (approxExtent || !extentOut) {
+                    extentOut = getExtent(resolvedCrs);
+                }
+                if (resolvedCrs->isEquivalentTo(
+                        crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
+                    return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs);
                 }
+            } catch (const std::exception &) {
             }
-            if (hasMissingId) {
-                const auto tmpAuthFactory = io::AuthorityFactory::create(
-                    authFactory->databaseContext(), *ids.front()->codeSpace());
-                try {
-                    auto resolvedCrs(
-                        tmpAuthFactory->createCompoundCRS(ids.front()->code()));
-                    if (resolvedCrs->isEquivalentTo(
-                            crs.get(),
-                            util::IComparable::Criterion::EQUIVALENT)) {
-                        return util::nn_static_pointer_cast<crs::CRS>(
-                            resolvedCrs);
-                    }
-                } catch (const std::exception &) {
+        } else {
+            auto outCrs = tryToIdentifyByName(
+                io::AuthorityFactory::ObjectType::COMPOUND_CRS);
+            if (outCrs.get() != crs.get()) {
+                return outCrs;
+            }
+            if (approxExtent || !extentOut) {
+                // If we still did not get a reliable extent, then try to
+                // resolve the components of the compoundCRS, and take the
+                // intersection of their extent.
+                extentOut = metadata::ExtentPtr();
+                const auto &components =
+                    compoundCrs->componentReferenceSystems();
+                for (const auto &component : components) {
+                    metadata::ExtentPtr componentExtent;
+                    getResolvedCRS(component, context, componentExtent);
+                    if (extentOut && componentExtent)
+                        extentOut = extentOut->intersection(
+                            NN_NO_CHECK(componentExtent));
+                    else if (componentExtent)
+                        extentOut = componentExtent;
                 }
             }
         }
@@ -14126,10 +14223,12 @@ CoordinateOperationFactory::createOperations(
     auto l_sourceCRS = srcBoundCRS ? NN_NO_CHECK(srcBoundCRS) : sourceCRS;
     auto l_targetCRS = targetBoundCRS ? NN_NO_CHECK(targetBoundCRS) : targetCRS;
 
-    auto l_resolvedSourceCRS = getResolvedCRS(l_sourceCRS, context);
-    auto l_resolvedTargetCRS = getResolvedCRS(l_targetCRS, context);
-    auto sourceCRSExtent(getExtent(l_resolvedSourceCRS));
-    auto targetCRSExtent(getExtent(l_resolvedTargetCRS));
+    metadata::ExtentPtr sourceCRSExtent;
+    auto l_resolvedSourceCRS =
+        getResolvedCRS(l_sourceCRS, context, sourceCRSExtent);
+    metadata::ExtentPtr targetCRSExtent;
+    auto l_resolvedTargetCRS =
+        getResolvedCRS(l_targetCRS, context, targetCRSExtent);
     Private::Context contextPrivate(sourceCRSExtent, targetCRSExtent, context);
 
     if (context->getSourceAndTargetCRSExtentUse() ==
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index f16ed7cf..26725e24 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -993,7 +993,8 @@ bool SingleCRS::baseIsEquivalentTo(
 
     // TODO test DatumEnsemble
     return d->coordinateSystem->_isEquivalentTo(
-        otherSingleCRS->d->coordinateSystem.get(), criterion);
+               otherSingleCRS->d->coordinateSystem.get(), criterion) &&
+           getExtensionProj4() == otherSingleCRS->getExtensionProj4();
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index d61c74c0..7d457dab 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -7196,8 +7196,7 @@ TEST(operation, compoundCRS_to_geogCRS_2D_promote_to_3D_context) {
             nnSrc,
             dst->promoteTo3D(std::string(), authFactory->databaseContext()),
             ctxt);
-    // It includes a trailing ballpart transformation
-    ASSERT_EQ(listCompoundToGeog3D.size(), listCompoundToGeog2D.size() + 1);
+    ASSERT_EQ(listCompoundToGeog3D.size(), listCompoundToGeog2D.size());
 }
 
 // ---------------------------------------------------------------------------
@@ -7255,6 +7254,83 @@ TEST(operation, compoundCRS_from_wkt_without_id_to_geogCRS) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation,
+     compoundCRS_of_projCRS_from_wkt_without_id_or_extent_to_geogCRS) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    auto wkt =
+        "COMPOUNDCRS[\"NAD83 / Pennsylvania South + NAVD88 height\",\n"
+        "    PROJCRS[\"NAD83 / Pennsylvania South\",\n"
+        "        BASEGEOGCRS[\"NAD83\",\n"
+        "            DATUM[\"North American Datum 1983\",\n"
+        "                ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n"
+        "                    LENGTHUNIT[\"metre\",1]]],\n"
+        "            PRIMEM[\"Greenwich\",0,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433]]],\n"
+        "        CONVERSION[\"SPCS83 Pennsylvania South zone (meters)\",\n"
+        "            METHOD[\"Lambert Conic Conformal (2SP)\",\n"
+        "                ID[\"EPSG\",9802]],\n"
+        "            PARAMETER[\"Latitude of false origin\",39.3333333333333,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8821]],\n"
+        "            PARAMETER[\"Longitude of false origin\",-77.75,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8822]],\n"
+        "            PARAMETER[\"Latitude of 1st standard "
+        "parallel\",40.9666666666667,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8823]],\n"
+        "            PARAMETER[\"Latitude of 2nd standard "
+        "parallel\",39.9333333333333,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8824]],\n"
+        "            PARAMETER[\"Easting at false origin\",600000,\n"
+        "                LENGTHUNIT[\"metre\",1],\n"
+        "                ID[\"EPSG\",8826]],\n"
+        "            PARAMETER[\"Northing at false origin\",0,\n"
+        "                LENGTHUNIT[\"metre\",1],\n"
+        "                ID[\"EPSG\",8827]]],\n"
+        "        CS[Cartesian,2],\n"
+        "            AXIS[\"easting (X)\",east,\n"
+        "                ORDER[1],\n"
+        "                LENGTHUNIT[\"metre\",1]],\n"
+        "            AXIS[\"northing (Y)\",north,\n"
+        "                ORDER[2],\n"
+        "                LENGTHUNIT[\"metre\",1]]],\n"
+        "    VERTCRS[\"NAVD88 height\",\n"
+        "        VDATUM[\"North American Vertical Datum 1988\"],\n"
+        "        CS[vertical,1],\n"
+        "            AXIS[\"gravity-related height (H)\",up,\n"
+        "                LENGTHUNIT[\"metre\",1]]]]";
+    auto srcObj =
+        createFromUserInput(wkt, authFactory->databaseContext(), false);
+    auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
+    ASSERT_TRUE(src != nullptr);
+    auto dst = authFactory->createCoordinateReferenceSystem("4269"); // NAD83
+
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        NN_NO_CHECK(src), dst, ctxt);
+    // NAD83 / Pennsylvania South + NAVD88 height
+    auto srcRefObj = createFromUserInput("EPSG:32129+5703",
+                                         authFactory->databaseContext(), false);
+    auto srcRef = nn_dynamic_pointer_cast<CRS>(srcRefObj);
+    ASSERT_TRUE(srcRef != nullptr);
+    ASSERT_TRUE(
+        src->isEquivalentTo(srcRef.get(), IComparable::Criterion::EQUIVALENT));
+    auto listRef = CoordinateOperationFactory::create()->createOperations(
+        NN_NO_CHECK(srcRef), dst, ctxt);
+
+    EXPECT_EQ(list.size(), listRef.size());
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(operation, compoundCRS_to_geogCRS_with_vertical_unit_change) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), "EPSG");