Make CRS identification work with CRS with DatumEnsemble

3 files changed, 345 insertions, 84 deletions

diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp
index bbdc9565..a028aceb 100644
--- a/include/proj/crs.hpp
+++ b/include/proj/crs.hpp
@@ -194,7 +194,10 @@ class PROJ_GCC_DLL SingleCRS : public CRS {
         PROJ_INTERNAL void
         exportDatumOrDatumEnsembleToWkt(io::WKTFormatter *formatter)
             const; // throw(io::FormattingException)
-                   //! @endcond
+
+    PROJ_INTERNAL const datum::DatumNNPtr
+    datumNonNull(const io::DatabaseContextPtr &dbContext) const;
+    //! @endcond
 
   protected:
     PROJ_INTERNAL SingleCRS(const datum::DatumPtr &datumIn,
@@ -289,6 +292,9 @@ class PROJ_GCC_DLL GeodeticCRS : virtual public SingleCRS,
         PROJ_INTERNAL void
         addDatumInfoToPROJString(io::PROJStringFormatter *formatter) const;
 
+    PROJ_INTERNAL const datum::GeodeticReferenceFrameNNPtr
+    datumNonNull(const io::DatabaseContextPtr &dbContext) const;
+
     PROJ_INTERNAL void addGeocentricUnitConversionIntoPROJString(
         io::PROJStringFormatter *formatter) const;
 
@@ -474,6 +480,9 @@ class PROJ_GCC_DLL VerticalCRS : virtual public SingleCRS,
         PROJ_INTERNAL void
         addLinearUnitConvert(io::PROJStringFormatter *formatter) const;
 
+    PROJ_INTERNAL const datum::VerticalReferenceFrameNNPtr
+    datumNonNull(const io::DatabaseContextPtr &dbContext) const;
+
     PROJ_INTERNAL void _exportToWKT(io::WKTFormatter *formatter)
         const override; // throw(io::FormattingException)
 
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index ecbd39e1..eb3918f7 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -1188,6 +1188,18 @@ const datum::DatumEnsemblePtr &SingleCRS::datumEnsemble() PROJ_PURE_DEFN {
 
 // ---------------------------------------------------------------------------
 
+//! @cond Doxygen_Suppress
+/** \brief Return the real datum or a synthetized one if a datumEnsemble.
+ */
+const datum::DatumNNPtr
+SingleCRS::datumNonNull(const io::DatabaseContextPtr &dbContext) const {
+    return d->datum ? NN_NO_CHECK(d->datum)
+                    : d->datumEnsemble->asDatum(dbContext);
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
 /** \brief Return the cs::CoordinateSystem associated with the CRS.
  *
  * @return a CoordinateSystem.
@@ -1207,20 +1219,41 @@ bool SingleCRS::baseIsEquivalentTo(
          !ObjectUsage::_isEquivalentTo(other, criterion, dbContext))) {
         return false;
     }
-    const auto &thisDatum = d->datum;
-    const auto &otherDatum = otherSingleCRS->d->datum;
-    if (thisDatum) {
-        if (!thisDatum->_isEquivalentTo(otherDatum.get(), criterion,
-                                        dbContext)) {
-            return false;
+
+    if (criterion == util::IComparable::Criterion::STRICT) {
+        const auto &thisDatum = d->datum;
+        const auto &otherDatum = otherSingleCRS->d->datum;
+        if (thisDatum) {
+            if (!thisDatum->_isEquivalentTo(otherDatum.get(), criterion,
+                                            dbContext)) {
+                return false;
+            }
+        } else {
+            if (otherDatum) {
+                return false;
+            }
+        }
+
+        const auto &thisDatumEnsemble = d->datumEnsemble;
+        const auto &otherDatumEnsemble = otherSingleCRS->d->datumEnsemble;
+        if (thisDatumEnsemble) {
+            if (!thisDatumEnsemble->_isEquivalentTo(otherDatumEnsemble.get(),
+                                                    criterion, dbContext)) {
+                return false;
+            }
+        } else {
+            if (otherDatumEnsemble) {
+                return false;
+            }
         }
     } else {
-        if (otherDatum) {
+        if (!datumNonNull(dbContext)->_isEquivalentTo(
+                otherSingleCRS->datumNonNull(dbContext).get(), criterion,
+                dbContext)) {
             return false;
         }
     }
 
-    // TODO test DatumEnsemble
     return d->coordinateSystem->_isEquivalentTo(
                otherSingleCRS->d->coordinateSystem.get(), criterion,
                dbContext) &&
@@ -1341,6 +1374,21 @@ const datum::GeodeticReferenceFramePtr &GeodeticCRS::datum() PROJ_PURE_DEFN {
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+/** \brief Return the real datum or a synthetized one if a datumEnsemble.
+ */
+const datum::GeodeticReferenceFrameNNPtr
+GeodeticCRS::datumNonNull(const io::DatabaseContextPtr &dbContext) const {
+    return NN_NO_CHECK(
+        d->datum_
+            ? d->datum_
+            : util::nn_dynamic_pointer_cast<datum::GeodeticReferenceFrame>(
+                  SingleCRS::getPrivate()->datumEnsemble->asDatum(dbContext)));
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 static datum::GeodeticReferenceFrame *oneDatum(const GeodeticCRS *crs) {
     const auto &l_datumEnsemble = crs->datumEnsemble();
     assert(l_datumEnsemble);
@@ -1704,8 +1752,8 @@ void GeodeticCRS::addDatumInfoToPROJString(
     const auto &TOWGS84Params = formatter->getTOWGS84Parameters();
     bool datumWritten = false;
     const auto &nadgrids = formatter->getHDatumExtension();
-    const auto &l_datum = datum();
-    if (formatter->getCRSExport() && l_datum && TOWGS84Params.empty() &&
+    const auto l_datum = datumNonNull(formatter->databaseContext());
+    if (formatter->getCRSExport() && TOWGS84Params.empty() &&
         nadgrids.empty()) {
         if (l_datum->_isEquivalentTo(
                 datum::GeodeticReferenceFrame::EPSG_6326.get(),
@@ -1953,11 +2001,12 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
 
     if (authorityFactory) {
 
-        const auto &thisDatum(datum());
+        const auto thisDatum(datumNonNull(dbContext));
 
-        auto searchByDatum = [this, &authorityFactory, &res, &thisDatum,
-                              &geodetic_crs_type, crsCriterion, &dbContext]() {
-            for (const auto &id : thisDatum->identifiers()) {
+        auto searchByDatumCode = [this, &authorityFactory, &res,
+                                  &geodetic_crs_type, crsCriterion, &dbContext](
+            const common::IdentifiedObjectNNPtr &l_datum) {
+            for (const auto &id : l_datum->identifiers()) {
                 try {
                     auto tempRes = authorityFactory->createGeodeticCRSFromDatum(
                         *id->codeSpace(), id->code(), geodetic_crs_type);
@@ -1972,10 +2021,10 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
             }
         };
 
-        const auto &thisEllipsoid(ellipsoid());
         auto searchByEllipsoid = [this, &authorityFactory, &res, &thisDatum,
-                                  &thisEllipsoid, &geodetic_crs_type,
-                                  l_implicitCS, &dbContext]() {
+                                  &geodetic_crs_type, l_implicitCS,
+                                  &dbContext]() {
+            const auto &thisEllipsoid = thisDatum->ellipsoid();
             const auto ellipsoids =
                 thisEllipsoid->identifiers().empty()
                     ? authorityFactory->createEllipsoidFromExisting(
@@ -1989,9 +2038,8 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                                 *id->codeSpace(), id->code(),
                                 geodetic_crs_type);
                         for (const auto &crs : tempRes) {
-                            const auto &crsDatum(crs->datum());
-                            if (crsDatum &&
-                                crsDatum->ellipsoid()->_isEquivalentTo(
+                            const auto crsDatum(crs->datumNonNull(dbContext));
+                            if (crsDatum->ellipsoid()->_isEquivalentTo(
                                     ellps.get(),
                                     util::IComparable::Criterion::EQUIVALENT,
                                     dbContext) &&
@@ -2013,18 +2061,32 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
             }
         };
 
+        const auto searchByDatumOrEllipsoid = [&authorityFactory, &res,
+                                               &thisDatum, searchByDatumCode,
+                                               searchByEllipsoid]() {
+            if (!thisDatum->identifiers().empty()) {
+                searchByDatumCode(thisDatum);
+            } else {
+                auto candidateDatums = authorityFactory->createObjectsFromName(
+                    thisDatum->nameStr(), {io::AuthorityFactory::ObjectType::
+                                               GEODETIC_REFERENCE_FRAME},
+                    false);
+                const size_t sizeBefore = res.size();
+                for (const auto &candidateDatum : candidateDatums) {
+                    searchByDatumCode(candidateDatum);
+                }
+                if (sizeBefore == res.size()) {
+                    searchByEllipsoid();
+                }
+            }
+        };
+
         const bool unsignificantName = thisName.empty() ||
                                        ci_equal(thisName, "unknown") ||
                                        ci_equal(thisName, "unnamed");
 
         if (unsignificantName) {
-            if (thisDatum) {
-                if (!thisDatum->identifiers().empty()) {
-                    searchByDatum();
-                } else {
-                    searchByEllipsoid();
-                }
-            }
+            searchByDatumOrEllipsoid();
         } else if (hasCodeCompatibleOfAuthorityFactory(this,
                                                        authorityFactory)) {
             // If the CRS has already an id, check in the database for the
@@ -2074,12 +2136,8 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                     break;
                 }
             }
-            if (!gotAbove25Pct && thisDatum) {
-                if (!thisDatum->identifiers().empty()) {
-                    searchByDatum();
-                } else {
-                    searchByEllipsoid();
-                }
+            if (!gotAbove25Pct) {
+                searchByDatumOrEllipsoid();
             }
         }
 
@@ -2106,22 +2164,20 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
             }
 
             // Then datum matching
-            const auto &aDatum(a.first->datum());
-            const auto &bDatum(b.first->datum());
-            if (thisDatum && aDatum && bDatum) {
-                const auto thisEquivADatum(thisDatum->_isEquivalentTo(
-                    aDatum.get(), util::IComparable::Criterion::EQUIVALENT,
-                    dbContext));
-                const auto thisEquivBDatum(thisDatum->_isEquivalentTo(
-                    bDatum.get(), util::IComparable::Criterion::EQUIVALENT,
-                    dbContext));
-
-                if (thisEquivADatum && !thisEquivBDatum) {
-                    return true;
-                }
-                if (!thisEquivADatum && thisEquivBDatum) {
-                    return false;
-                }
+            const auto aDatum(a.first->datumNonNull(dbContext));
+            const auto bDatum(b.first->datumNonNull(dbContext));
+            const auto thisEquivADatum(thisDatum->_isEquivalentTo(
+                aDatum.get(), util::IComparable::Criterion::EQUIVALENT,
+                dbContext));
+            const auto thisEquivBDatum(thisDatum->_isEquivalentTo(
+                bDatum.get(), util::IComparable::Criterion::EQUIVALENT,
+                dbContext));
+
+            if (thisEquivADatum && !thisEquivBDatum) {
+                return true;
+            }
+            if (!thisEquivADatum && thisEquivBDatum) {
+                return false;
             }
 
             // Then coordinate system matching
@@ -2153,23 +2209,21 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                 return false;
             }
 
-            if (aDatum && bDatum) {
-                // Favor the CRS whole ellipsoid names matches the ellipsoid
-                // name (WGS84...)
-                const bool aEllpsNameEqCRSName =
-                    metadata::Identifier::isEquivalentName(
-                        aDatum->ellipsoid()->nameStr().c_str(),
-                        a.first->nameStr().c_str());
-                const bool bEllpsNameEqCRSName =
-                    metadata::Identifier::isEquivalentName(
-                        bDatum->ellipsoid()->nameStr().c_str(),
-                        b.first->nameStr().c_str());
-                if (aEllpsNameEqCRSName && !bEllpsNameEqCRSName) {
-                    return true;
-                }
-                if (bEllpsNameEqCRSName && !aEllpsNameEqCRSName) {
-                    return false;
-                }
+            // Favor the CRS whole ellipsoid names matches the ellipsoid
+            // name (WGS84...)
+            const bool aEllpsNameEqCRSName =
+                metadata::Identifier::isEquivalentName(
+                    aDatum->ellipsoid()->nameStr().c_str(),
+                    a.first->nameStr().c_str());
+            const bool bEllpsNameEqCRSName =
+                metadata::Identifier::isEquivalentName(
+                    bDatum->ellipsoid()->nameStr().c_str(),
+                    b.first->nameStr().c_str());
+            if (aEllpsNameEqCRSName && !bEllpsNameEqCRSName) {
+                return true;
+            }
+            if (bEllpsNameEqCRSName && !aEllpsNameEqCRSName) {
+                return false;
             }
 
             // Arbitrary final sorting criterion
@@ -2341,6 +2395,12 @@ bool GeographicCRS::is2DPartOf3D(util::nn<const GeographicCRS *> other)
         return thisDatum->_isEquivalentTo(
             otherDatum.get(), util::IComparable::Criterion::EQUIVALENT);
     }
+    const auto &thisDatumEnsemble = datumEnsemble();
+    const auto &otherDatumEnsemble = other->datumEnsemble();
+    if (thisDatumEnsemble && otherDatumEnsemble) {
+        return thisDatumEnsemble->_isEquivalentTo(
+            otherDatumEnsemble.get(), util::IComparable::Criterion::EQUIVALENT);
+    }
     return false;
 }
 
@@ -2580,15 +2640,13 @@ void GeographicCRS::_exportToPROJString(
         if (formatter->getLegacyCRSToCRSContext() &&
             formatter->getHDatumExtension().empty() &&
             formatter->getTOWGS84Parameters().empty()) {
-            const auto &l_datum = datum();
-            if (l_datum &&
-                l_datum->_isEquivalentTo(
+            const auto l_datum = datumNonNull(formatter->databaseContext());
+            if (l_datum->_isEquivalentTo(
                     datum::GeodeticReferenceFrame::EPSG_6326.get(),
                     util::IComparable::Criterion::EQUIVALENT)) {
                 done = true;
                 formatter->addParam("ellps", "WGS84");
-            } else if (l_datum &&
-                       l_datum->_isEquivalentTo(
+            } else if (l_datum->_isEquivalentTo(
                            datum::GeodeticReferenceFrame::EPSG_6269.get(),
                            util::IComparable::Criterion::EQUIVALENT)) {
                 done = true;
@@ -2755,6 +2813,19 @@ const cs::VerticalCSNNPtr VerticalCRS::coordinateSystem() const {
 // ---------------------------------------------------------------------------
 
 //! @cond Doxygen_Suppress
+/** \brief Return the real datum or a synthetized one if a datumEnsemble.
+ */
+const datum::VerticalReferenceFrameNNPtr
+VerticalCRS::datumNonNull(const io::DatabaseContextPtr &dbContext) const {
+    return NN_NO_CHECK(
+        util::nn_dynamic_pointer_cast<datum::VerticalReferenceFrame>(
+            SingleCRS::datumNonNull(dbContext)));
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+//! @cond Doxygen_Suppress
 void VerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const {
     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
     formatter->startNode(isWKT2 ? io::WKTConstants::VERTCRS
@@ -3846,12 +3917,15 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
     std::list<Pair> res;
 
     const auto &thisName(nameStr());
+    io::DatabaseContextPtr dbContext =
+        authorityFactory ? authorityFactory->databaseContext().as_nullable()
+                         : nullptr;
 
     std::list<std::pair<GeodeticCRSNNPtr, int>> baseRes;
     const auto &l_baseCRS(baseCRS());
-    const auto l_datum = l_baseCRS->datum();
+    const auto l_datum = l_baseCRS->datumNonNull(dbContext);
     const bool significantNameForDatum =
-        l_datum != nullptr && !ci_starts_with(l_datum->nameStr(), "unknown") &&
+        !ci_starts_with(l_datum->nameStr(), "unknown") &&
         l_datum->nameStr() != "unnamed";
     const auto &ellipsoid = l_baseCRS->ellipsoid();
     auto geogCRS = dynamic_cast<const GeographicCRS *>(l_baseCRS.get());
@@ -3886,9 +3960,6 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
 
     const auto &conv = derivingConversionRef();
     const auto &cs = coordinateSystem();
-    io::DatabaseContextPtr dbContext =
-        authorityFactory ? authorityFactory->databaseContext().as_nullable()
-                         : nullptr;
 
     if (baseRes.size() == 1 && baseRes.front().second >= 70 &&
         conv->isUTM(zone, north) &&
@@ -3978,10 +4049,9 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
                                     util::IComparable::Criterion::EQUIVALENT,
                                     dbContext)) {
                 if (!significantNameForDatum ||
-                    (crs->baseCRS()->datum() &&
-                     l_datum->_isEquivalentTo(
-                         crs->baseCRS()->datum().get(),
-                         util::IComparable::Criterion::EQUIVALENT))) {
+                    l_datum->_isEquivalentTo(
+                        crs->baseCRS()->datumNonNull(dbContext).get(),
+                        util::IComparable::Criterion::EQUIVALENT)) {
                     res.emplace_back(crs, 70);
                 } else {
                     res.emplace_back(crs, 60);
@@ -4863,7 +4933,6 @@ BoundCRSNNPtr BoundCRS::createFromNadgrids(const CRSNNPtr &baseCRSIn,
             ? NN_NO_CHECK(std::static_pointer_cast<CRS>(sourceGeographicCRS))
             : baseCRSIn;
     if (sourceGeographicCRS != nullptr &&
-        sourceGeographicCRS->datum() != nullptr &&
         sourceGeographicCRS->primeMeridian()->longitude().getSIValue() != 0.0) {
         transformationSourceCRS = GeographicCRS::create(
             util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
@@ -4872,9 +4941,9 @@ BoundCRSNNPtr BoundCRS::createFromNadgrids(const CRSNNPtr &baseCRSIn,
             datum::GeodeticReferenceFrame::create(
                 util::PropertyMap().set(
                     common::IdentifiedObject::NAME_KEY,
-                    sourceGeographicCRS->datum()->nameStr() +
+                    sourceGeographicCRS->datumNonNull(nullptr)->nameStr() +
                         " (with Greenwich prime meridian)"),
-                sourceGeographicCRS->datum()->ellipsoid(),
+                sourceGeographicCRS->datumNonNull(nullptr)->ellipsoid(),
                 util::optional<std::string>(), datum::PrimeMeridian::GREENWICH),
             sourceGeographicCRS->coordinateSystem());
     }
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index e8758932..ec30580a 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -1680,6 +1680,90 @@ TEST(crs, geodeticcrs_identify_db) {
         EXPECT_EQ(res.front().first->identifiers()[0]->code(), "7844");
         EXPECT_EQ(res.front().second, 100);
     }
+    {
+        // Identify with DatumEnsemble
+        auto wkt =
+            "GEOGCRS[\"WGS 84\","
+            "  ENSEMBLE[\"World Geodetic System 1984 ensemble\","
+            "     MEMBER[\"World Geodetic System 1984 (Transit)\","
+            "       ID[\"EPSG\",1166]],"
+            "     MEMBER[\"World Geodetic System 1984 (G730)\","
+            "       ID[\"EPSG\",1152]],"
+            "     MEMBER[\"World Geodetic System 1984 (G873)\","
+            "       ID[\"EPSG\",1153]],"
+            "     MEMBER[\"World Geodetic System 1984 (G1150)\","
+            "       ID[\"EPSG\",1154]],"
+            "     MEMBER[\"World Geodetic System 1984 (G1674)\","
+            "       ID[\"EPSG\",1155]],"
+            "     MEMBER[\"World Geodetic System 1984 (G1762)\","
+            "       ID[\"EPSG\",1156]],"
+            "     ELLIPSOID[\"WGS 84\",6378137,298.257223563,"
+            "      LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]],"
+            "      ID[\"EPSG\",7030]],"
+            "     ENSEMBLEACCURACY[2]],"
+            "  PRIMEM[\"Greenwich\",0,"
+            "    ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9102]],"
+            "    ID[\"EPSG\",8901]],"
+            "  CS[ellipsoidal,2,"
+            "    ID[\"EPSG\",6422]],"
+            "  AXIS[\"Geodetic latitude (Lat)\",north,"
+            "    ORDER[1]],"
+            "  AXIS[\"Geodetic longitude (Lon)\",east,"
+            "    ORDER[2]],"
+            "  ANGLEUNIT[\"degree (supplier to define representation)\","
+            "0.0174532925199433,ID[\"EPSG\",9122]]]";
+        auto obj = WKTParser().createFromWKT(wkt);
+        auto crs = nn_dynamic_pointer_cast<GeodeticCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+
+        auto allFactory = AuthorityFactory::create(dbContext, std::string());
+        auto res = crs->identify(allFactory);
+        ASSERT_EQ(res.size(), 1U);
+        EXPECT_EQ(res.front().first->getEPSGCode(), 4326);
+        EXPECT_EQ(res.front().second, 100.0);
+    }
+    {
+        // Identify with DatumEnsemble and unknown CRS name
+        auto wkt =
+            "GEOGCRS[\"unknown\","
+            "  ENSEMBLE[\"World Geodetic System 1984 ensemble\","
+            "     MEMBER[\"World Geodetic System 1984 (Transit)\","
+            "       ID[\"EPSG\",1166]],"
+            "     MEMBER[\"World Geodetic System 1984 (G730)\","
+            "       ID[\"EPSG\",1152]],"
+            "     MEMBER[\"World Geodetic System 1984 (G873)\","
+            "       ID[\"EPSG\",1153]],"
+            "     MEMBER[\"World Geodetic System 1984 (G1150)\","
+            "       ID[\"EPSG\",1154]],"
+            "     MEMBER[\"World Geodetic System 1984 (G1674)\","
+            "       ID[\"EPSG\",1155]],"
+            "     MEMBER[\"World Geodetic System 1984 (G1762)\","
+            "       ID[\"EPSG\",1156]],"
+            "     ELLIPSOID[\"WGS 84\",6378137,298.257223563,"
+            "      LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]],"
+            "      ID[\"EPSG\",7030]],"
+            "     ENSEMBLEACCURACY[2]],"
+            "  PRIMEM[\"Greenwich\",0,"
+            "    ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9102]],"
+            "    ID[\"EPSG\",8901]],"
+            "  CS[ellipsoidal,2,"
+            "    ID[\"EPSG\",6422]],"
+            "  AXIS[\"Geodetic latitude (Lat)\",north,"
+            "    ORDER[1]],"
+            "  AXIS[\"Geodetic longitude (Lon)\",east,"
+            "    ORDER[2]],"
+            "  ANGLEUNIT[\"degree (supplier to define representation)\","
+            "0.0174532925199433,ID[\"EPSG\",9122]]]";
+        auto obj = WKTParser().createFromWKT(wkt);
+        auto crs = nn_dynamic_pointer_cast<GeodeticCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+
+        auto allFactory = AuthorityFactory::create(dbContext, std::string());
+        auto res = crs->identify(allFactory);
+        ASSERT_EQ(res.size(), 1U);
+        EXPECT_EQ(res.front().first->getEPSGCode(), 4326);
+        EXPECT_EQ(res.front().second, 70.0);
+    }
 }
 
 // ---------------------------------------------------------------------------
@@ -2694,6 +2778,105 @@ TEST(crs, projectedCRS_identify_db) {
         EXPECT_EQ(res.front().first->getEPSGCode(), 2154);
         EXPECT_EQ(res.front().second, 90);
     }
+    {
+        // Identify with DatumEnsemble
+        auto wkt =
+            "PROJCRS[\"WGS 84 / UTM zone 31N\","
+            "  BASEGEOGCRS[\"WGS 84\","
+            "    ENSEMBLE[\"World Geodetic System 1984 ensemble\","
+            "       MEMBER[\"World Geodetic System 1984 (Transit)\","
+            "         ID[\"EPSG\",1166]],"
+            "       MEMBER[\"World Geodetic System 1984 (G730)\","
+            "         ID[\"EPSG\",1152]],"
+            "       MEMBER[\"World Geodetic System 1984 (G873)\","
+            "         ID[\"EPSG\",1153]],"
+            "       MEMBER[\"World Geodetic System 1984 (G1150)\","
+            "         ID[\"EPSG\",1154]],"
+            "       MEMBER[\"World Geodetic System 1984 (G1674)\","
+            "         ID[\"EPSG\",1155]],"
+            "       MEMBER[\"World Geodetic System 1984 (G1762)\","
+            "         ID[\"EPSG\",1156]],"
+            "       ELLIPSOID[\"WGS 84\",6378137,298.257223563,"
+            "        LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],"
+            "       ENSEMBLEACCURACY[2]]],"
+            "  CONVERSION[\"UTM zone 31N\","
+            "    METHOD[\"Transverse Mercator\","
+            "      ID[\"EPSG\",9807]],"
+            "    PARAMETER[\"Latitude of natural origin\",0,"
+            "      ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9102]]],"
+            "    PARAMETER[\"Longitude of natural origin\",3,"
+            "      ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9102]]],"
+            "    PARAMETER[\"Scale factor at natural origin\",0.9996,"
+            "      SCALEUNIT[\"unity\",1,ID[\"EPSG\",9201]]],"
+            "    PARAMETER[\"False easting\",500000,"
+            "      LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],"
+            "    PARAMETER[\"False northing\",0,"
+            "      LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],"
+            "  CS[Cartesian,2],"
+            " AXIS[\"Easting (E)\",east,"
+            "   ORDER[1]],"
+            "  AXIS[\"Northing (N)\",north,"
+            "    ORDER[2]],"
+            "  LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]";
+        auto obj = WKTParser().createFromWKT(wkt);
+        auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+
+        auto allFactory = AuthorityFactory::create(dbContext, std::string());
+        auto res = crs->identify(allFactory);
+        ASSERT_EQ(res.size(), 1U);
+        EXPECT_EQ(res.front().first->getEPSGCode(), 32631);
+        EXPECT_EQ(res.front().second, 100.0);
+    }
+    {
+        // Identify with DatumEnsemble and unknown CRS name
+        auto wkt =
+            "PROJCRS[\"unknown\","
+            "  BASEGEOGCRS[\"unknown\","
+            "    ENSEMBLE[\"World Geodetic System 1984 ensemble\","
+            "       MEMBER[\"World Geodetic System 1984 (Transit)\","
+            "         ID[\"EPSG\",1166]],"
+            "       MEMBER[\"World Geodetic System 1984 (G730)\","
+            "         ID[\"EPSG\",1152]],"
+            "       MEMBER[\"World Geodetic System 1984 (G873)\","
+            "         ID[\"EPSG\",1153]],"
+            "       MEMBER[\"World Geodetic System 1984 (G1150)\","
+            "         ID[\"EPSG\",1154]],"
+            "       MEMBER[\"World Geodetic System 1984 (G1674)\","
+            "         ID[\"EPSG\",1155]],"
+            "       MEMBER[\"World Geodetic System 1984 (G1762)\","
+            "         ID[\"EPSG\",1156]],"
+            "       ELLIPSOID[\"WGS 84\",6378137,298.257223563,"
+            "        LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],"
+            "       ENSEMBLEACCURACY[2]]],"
+            "  CONVERSION[\"UTM zone 31N\","
+            "    METHOD[\"Transverse Mercator\","
+            "      ID[\"EPSG\",9807]],"
+            "    PARAMETER[\"Latitude of natural origin\",0,"
+            "      ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9102]]],"
+            "    PARAMETER[\"Longitude of natural origin\",3,"
+            "      ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9102]]],"
+            "    PARAMETER[\"Scale factor at natural origin\",0.9996,"
+            "      SCALEUNIT[\"unity\",1,ID[\"EPSG\",9201]]],"
+            "    PARAMETER[\"False easting\",500000,"
+            "      LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],"
+            "    PARAMETER[\"False northing\",0,"
+            "      LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],"
+            "  CS[Cartesian,2],"
+            " AXIS[\"Easting (E)\",east,"
+            "   ORDER[1]],"
+            "  AXIS[\"Northing (N)\",north,"
+            "    ORDER[2]],"
+            "  LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]";
+        auto obj = WKTParser().createFromWKT(wkt);
+        auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+        ASSERT_TRUE(crs != nullptr);
+
+        auto res = crs->identify(factoryEPSG);
+        ASSERT_EQ(res.size(), 1U);
+        EXPECT_EQ(res.front().first->getEPSGCode(), 32631);
+        EXPECT_GE(res.front().second, 70.0);
+    }
 }
 
 // ---------------------------------------------------------------------------