Add C API to work with datum ensemble

Add:
- proj_crs_get_datum_ensemble()
- proj_crs_get_datum_forced()
- proj_datum_ensemble_get_member_count()
- proj_datum_ensemble_get_accuracy()
- proj_datum_ensemble_get_member()

Make proj_create_geographic_crs_from_datum() and
proj_create_geocentric_crs_from_datum() accept a datum ensemble.

5 files changed, 311 insertions, 20 deletions

diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index f37dde16..2ce13865 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/scripts/reference_exported_symbols.txt
@@ -971,6 +971,8 @@ proj_crs_demote_to_2D
 proj_crs_get_coordinate_system
 proj_crs_get_coordoperation
 proj_crs_get_datum
+proj_crs_get_datum_ensemble
+proj_crs_get_datum_forced
 proj_crs_get_geodetic_crs
 proj_crs_get_horizontal_datum
 proj_crs_get_sub_crs
@@ -979,6 +981,9 @@ proj_crs_promote_to_3D
 proj_cs_get_axis_count
 proj_cs_get_axis_info
 proj_cs_get_type
+proj_datum_ensemble_get_accuracy
+proj_datum_ensemble_get_member
+proj_datum_ensemble_get_member_count
 proj_degree_input
 proj_degree_output
 proj_destroy
diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 21a07729..fb0a3af7 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -3044,31 +3044,35 @@ PJ *proj_create_geographic_crs(PJ_CONTEXT *ctx, const char *crs_name,
  *
  * @param ctx PROJ context, or NULL for default context
  * @param crs_name Name of the GeographicCRS. Or NULL
- * @param datum Datum. Must not be NULL.
+ * @param datum_or_datum_ensemble Datum or DatumEnsemble (DatumEnsemble possible
+ * since 7.2). Must not be NULL.
  * @param ellipsoidal_cs Coordinate system. Must not be NULL.
  *
  * @return Object of type GeographicCRS that must be unreferenced with
  * proj_destroy(), or NULL in case of error.
  */
 PJ *proj_create_geographic_crs_from_datum(PJ_CONTEXT *ctx, const char *crs_name,
-                                          PJ *datum, PJ *ellipsoidal_cs) {
+                                          PJ *datum_or_datum_ensemble,
+                                          PJ *ellipsoidal_cs) {
 
     SANITIZE_CTX(ctx);
-    auto l_datum =
-        std::dynamic_pointer_cast<GeodeticReferenceFrame>(datum->iso_obj);
-    if (!l_datum) {
+    if (datum_or_datum_ensemble == nullptr) {
         proj_log_error(ctx, __FUNCTION__,
-                       "datum is not a GeodeticReferenceFrame");
+                       "Missing input datum_or_datum_ensemble");
         return nullptr;
     }
+    auto l_datum = std::dynamic_pointer_cast<GeodeticReferenceFrame>(
+        datum_or_datum_ensemble->iso_obj);
+    auto l_datum_ensemble = std::dynamic_pointer_cast<DatumEnsemble>(
+        datum_or_datum_ensemble->iso_obj);
     auto cs = std::dynamic_pointer_cast<EllipsoidalCS>(ellipsoidal_cs->iso_obj);
     if (!cs) {
         return nullptr;
     }
     try {
         auto geogCRS =
-            GeographicCRS::create(createPropertyMapName(crs_name),
-                                  NN_NO_CHECK(l_datum), NN_NO_CHECK(cs));
+            GeographicCRS::create(createPropertyMapName(crs_name), l_datum,
+                                  l_datum_ensemble, NN_NO_CHECK(cs));
         return pj_obj_create(ctx, geogCRS);
     } catch (const std::exception &e) {
         proj_log_error(ctx, __FUNCTION__, e.what());
@@ -3141,7 +3145,8 @@ PJ *proj_create_geocentric_crs(
  *
  * @param ctx PROJ context, or NULL for default context
  * @param crs_name Name of the GeographicCRS. Or NULL
- * @param datum Datum. Must not be NULL.
+ * @param datum_or_datum_ensemble Datum or DatumEnsemble (DatumEnsemble possible
+ * since 7.2). Must not be NULL.
  * @param linear_units Name of the linear units. Or NULL for Metre
  * @param linear_units_conv Conversion factor from the linear unit to metre. Or
  * 0 for Metre if linear_units == NULL. Otherwise should be not NULL
@@ -3150,22 +3155,24 @@ PJ *proj_create_geocentric_crs(
  * proj_destroy(), or NULL in case of error.
  */
 PJ *proj_create_geocentric_crs_from_datum(PJ_CONTEXT *ctx, const char *crs_name,
-                                          const PJ *datum,
+                                          const PJ *datum_or_datum_ensemble,
                                           const char *linear_units,
                                           double linear_units_conv) {
     SANITIZE_CTX(ctx);
+    if (datum_or_datum_ensemble == nullptr) {
+        proj_log_error(ctx, __FUNCTION__,
+                       "Missing input datum_or_datum_ensemble");
+        return nullptr;
+    }
+    auto l_datum = std::dynamic_pointer_cast<GeodeticReferenceFrame>(
+        datum_or_datum_ensemble->iso_obj);
+    auto l_datum_ensemble = std::dynamic_pointer_cast<DatumEnsemble>(
+        datum_or_datum_ensemble->iso_obj);
     try {
         const UnitOfMeasure linearUnit(
             createLinearUnit(linear_units, linear_units_conv));
-        auto l_datum =
-            std::dynamic_pointer_cast<GeodeticReferenceFrame>(datum->iso_obj);
-        if (!l_datum) {
-            proj_log_error(ctx, __FUNCTION__,
-                           "datum is not a GeodeticReferenceFrame");
-            return nullptr;
-        }
         auto geodCRS = GeodeticCRS::create(
-            createPropertyMapName(crs_name), NN_NO_CHECK(l_datum),
+            createPropertyMapName(crs_name), l_datum, l_datum_ensemble,
             cs::CartesianCS::createGeocentric(linearUnit));
         return pj_obj_create(ctx, geodCRS);
     } catch (const std::exception &e) {
@@ -7981,6 +7988,173 @@ PJ *proj_crs_get_datum(PJ_CONTEXT *ctx, const PJ *crs) {
 
 // ---------------------------------------------------------------------------
 
+/** \brief Returns the datum ensemble of a SingleCRS.
+ *
+ * The returned object must be unreferenced with proj_destroy() after
+ * use.
+ * It should be used by at most one thread at a time.
+ *
+ * @param ctx PROJ context, or NULL for default context
+ * @param crs Object of type SingleCRS (must not be NULL)
+ * @return Object that must be unreferenced with proj_destroy(), or NULL
+ * in case of error (or if there is no datum ensemble)
+ *
+ * @since 7.2
+ */
+PJ *proj_crs_get_datum_ensemble(PJ_CONTEXT *ctx, const PJ *crs) {
+    SANITIZE_CTX(ctx);
+    if (!crs) {
+        proj_log_error(ctx, __FUNCTION__, "missing required input");
+        return nullptr;
+    }
+    auto l_crs = dynamic_cast<const SingleCRS *>(crs->iso_obj.get());
+    if (!l_crs) {
+        proj_log_error(ctx, __FUNCTION__, "Object is not a SingleCRS");
+        return nullptr;
+    }
+    const auto &datumEnsemble = l_crs->datumEnsemble();
+    if (!datumEnsemble) {
+        return nullptr;
+    }
+    return pj_obj_create(ctx, NN_NO_CHECK(datumEnsemble));
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Returns the number of members of a datum ensemble.
+ *
+ * @param ctx PROJ context, or NULL for default context
+ * @param datum_ensemble Object of type DatumEnsemble (must not be NULL)
+ *
+ * @since 7.2
+ */
+int proj_datum_ensemble_get_member_count(PJ_CONTEXT *ctx,
+                                         const PJ *datum_ensemble) {
+    SANITIZE_CTX(ctx);
+    if (!datum_ensemble) {
+        proj_log_error(ctx, __FUNCTION__, "missing required input");
+        return 0;
+    }
+    auto l_datum_ensemble =
+        dynamic_cast<const DatumEnsemble *>(datum_ensemble->iso_obj.get());
+    if (!l_datum_ensemble) {
+        proj_log_error(ctx, __FUNCTION__, "Object is not a DatumEnsemble");
+        return 0;
+    }
+    return static_cast<int>(l_datum_ensemble->datums().size());
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Returns the positional accuracy of the datum ensemble.
+ *
+ * @param ctx PROJ context, or NULL for default context
+ * @param datum_ensemble Object of type DatumEnsemble (must not be NULL)
+ * @return the accuracy, or -1 in case of error.
+ *
+ * @since 7.2
+ */
+double proj_datum_ensemble_get_accuracy(PJ_CONTEXT *ctx,
+                                        const PJ *datum_ensemble) {
+    SANITIZE_CTX(ctx);
+    if (!datum_ensemble) {
+        proj_log_error(ctx, __FUNCTION__, "missing required input");
+        return -1;
+    }
+    auto l_datum_ensemble =
+        dynamic_cast<const DatumEnsemble *>(datum_ensemble->iso_obj.get());
+    if (!l_datum_ensemble) {
+        proj_log_error(ctx, __FUNCTION__, "Object is not a DatumEnsemble");
+        return -1;
+    }
+    const auto &accuracy = l_datum_ensemble->positionalAccuracy();
+    try {
+        return c_locale_stod(accuracy->value());
+    } catch (const std::exception &) {
+    }
+    return -1;
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Returns a member from a datum ensemble.
+ *
+ * The returned object must be unreferenced with proj_destroy() after
+ * use.
+ * It should be used by at most one thread at a time.
+ *
+ * @param ctx PROJ context, or NULL for default context
+ * @param datum_ensemble Object of type DatumEnsemble (must not be NULL)
+ * @param member_index Index of the datum member to extract (between 0 and
+ * proj_datum_ensemble_get_member_count()-1)
+ * @return Object that must be unreferenced with proj_destroy(), or NULL
+ * in case of error (or if there is no datum ensemble)
+ *
+ * @since 7.2
+ */
+PJ *proj_datum_ensemble_get_member(PJ_CONTEXT *ctx, const PJ *datum_ensemble,
+                                   int member_index) {
+    SANITIZE_CTX(ctx);
+    if (!datum_ensemble) {
+        proj_log_error(ctx, __FUNCTION__, "missing required input");
+        return nullptr;
+    }
+    auto l_datum_ensemble =
+        dynamic_cast<const DatumEnsemble *>(datum_ensemble->iso_obj.get());
+    if (!l_datum_ensemble) {
+        proj_log_error(ctx, __FUNCTION__, "Object is not a DatumEnsemble");
+        return nullptr;
+    }
+    if (member_index < 0 ||
+        member_index >= static_cast<int>(l_datum_ensemble->datums().size())) {
+        proj_log_error(ctx, __FUNCTION__, "Invalid member_index");
+        return nullptr;
+    }
+    return pj_obj_create(ctx, l_datum_ensemble->datums()[member_index]);
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Returns a datum for a SingleCRS.
+ *
+ * If the SingleCRS has a datum, then this datum is returned.
+ * Otherwise, the SingleCRS has a datum ensemble, and this datum ensemble is
+ * returned as a regular datum instead of a datum ensemble.
+ *
+ * The returned object must be unreferenced with proj_destroy() after
+ * use.
+ * It should be used by at most one thread at a time.
+ *
+ * @param ctx PROJ context, or NULL for default context
+ * @param crs Object of type SingleCRS (must not be NULL)
+ * @return Object that must be unreferenced with proj_destroy(), or NULL
+ * in case of error (or if there is no datum)
+ *
+ * @since 7.2
+ */
+PJ *proj_crs_get_datum_forced(PJ_CONTEXT *ctx, const PJ *crs) {
+    SANITIZE_CTX(ctx);
+    if (!crs) {
+        proj_log_error(ctx, __FUNCTION__, "missing required input");
+        return nullptr;
+    }
+    auto l_crs = dynamic_cast<const SingleCRS *>(crs->iso_obj.get());
+    if (!l_crs) {
+        proj_log_error(ctx, __FUNCTION__, "Object is not a SingleCRS");
+        return nullptr;
+    }
+    const auto &datum = l_crs->datum();
+    if (datum) {
+        return pj_obj_create(ctx, NN_NO_CHECK(datum));
+    }
+    const auto &datumEnsemble = l_crs->datumEnsemble();
+    assert(datumEnsemble);
+    auto dbContext = getDBcontextNoException(ctx, __FUNCTION__);
+    return pj_obj_create(ctx, datumEnsemble->asDatum(dbContext));
+}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Returns the coordinate system of a SingleCRS.
  *
  * The returned object must be unreferenced with proj_destroy() after
diff --git a/src/proj.h b/src/proj.h
index 512dcd84..92966526 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -1249,6 +1249,20 @@ PJ PROJ_DLL *proj_crs_get_sub_crs(PJ_CONTEXT *ctx, const PJ *crs, int index);
 
 PJ PROJ_DLL *proj_crs_get_datum(PJ_CONTEXT *ctx, const PJ *crs);
 
+PJ PROJ_DLL *proj_crs_get_datum_ensemble(PJ_CONTEXT *ctx, const PJ *crs);
+
+PJ PROJ_DLL *proj_crs_get_datum_forced(PJ_CONTEXT *ctx, const PJ *crs);
+
+int PROJ_DLL proj_datum_ensemble_get_member_count(PJ_CONTEXT *ctx,
+                                                  const PJ *datum_ensemble);
+
+double PROJ_DLL proj_datum_ensemble_get_accuracy(PJ_CONTEXT *ctx,
+                                                 const PJ *datum_ensemble);
+
+PJ PROJ_DLL *proj_datum_ensemble_get_member(PJ_CONTEXT *ctx,
+                                            const PJ *datum_ensemble,
+                                            int member_index);
+
 PJ PROJ_DLL *proj_crs_get_coordinate_system(PJ_CONTEXT *ctx, const PJ *crs);
 
 PJ_COORDINATE_SYSTEM_TYPE PROJ_DLL proj_cs_get_type(PJ_CONTEXT *ctx,
diff --git a/src/proj_experimental.h b/src/proj_experimental.h
index 50f1fb9f..05e49451 100644
--- a/src/proj_experimental.h
+++ b/src/proj_experimental.h
@@ -168,7 +168,7 @@ PJ PROJ_DLL *proj_create_geographic_crs(
 PJ PROJ_DLL *proj_create_geographic_crs_from_datum(
                             PJ_CONTEXT *ctx,
                             const char *crs_name,
-                            PJ* datum,
+                            PJ* datum_or_datum_ensemble,
                             PJ* ellipsoidal_cs);
 
 PJ PROJ_DLL *proj_create_geocentric_crs(
@@ -187,7 +187,7 @@ PJ PROJ_DLL *proj_create_geocentric_crs(
 PJ PROJ_DLL *proj_create_geocentric_crs_from_datum(
                             PJ_CONTEXT *ctx,
                             const char *crs_name,
-                            const PJ* datum,
+                            const PJ* datum_or_datum_ensemble,
                             const char *linear_units,
                             double linear_units_conv);
 
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 8e1dddfa..bd59e2bf 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -4917,4 +4917,102 @@ TEST_F(CApi, proj_is_equivalent_to_with_ctx) {
                                                PJ_COMP_EQUIVALENT));
 }
 
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, datum_ensemble) {
+    auto wkt =
+        "GEOGCRS[\"ETRS89\","
+        "    ENSEMBLE[\"European Terrestrial Reference System 1989 ensemble\","
+        "        MEMBER[\"European Terrestrial Reference Frame 1989\"],"
+        "        MEMBER[\"European Terrestrial Reference Frame 1990\"],"
+        "        MEMBER[\"European Terrestrial Reference Frame 1991\"],"
+        "        MEMBER[\"European Terrestrial Reference Frame 1992\"],"
+        "        MEMBER[\"European Terrestrial Reference Frame 1993\"],"
+        "        MEMBER[\"European Terrestrial Reference Frame 1994\"],"
+        "        MEMBER[\"European Terrestrial Reference Frame 1996\"],"
+        "        MEMBER[\"European Terrestrial Reference Frame 1997\"],"
+        "        MEMBER[\"European Terrestrial Reference Frame 2000\"],"
+        "        MEMBER[\"European Terrestrial Reference Frame 2005\"],"
+        "        MEMBER[\"European Terrestrial Reference Frame 2014\"],"
+        "        ELLIPSOID[\"GRS 1980\",6378137,298.257222101,"
+        "            LENGTHUNIT[\"metre\",1]],"
+        "        ENSEMBLEACCURACY[0.1]],"
+        "    PRIMEM[\"Greenwich\",0,"
+        "        ANGLEUNIT[\"degree\",0.0174532925199433]],"
+        "    CS[ellipsoidal,2],"
+        "        AXIS[\"geodetic latitude (Lat)\",north,"
+        "            ORDER[1],"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]],"
+        "        AXIS[\"geodetic longitude (Lon)\",east,"
+        "            ORDER[2],"
+        "            ANGLEUNIT[\"degree\",0.0174532925199433]]]";
+    auto from_wkt =
+        proj_create_from_wkt(m_ctxt, wkt, nullptr, nullptr, nullptr);
+    ObjectKeeper keeper_from_wkt(from_wkt);
+    EXPECT_NE(from_wkt, nullptr);
+
+    auto datum = proj_crs_get_datum(m_ctxt, from_wkt);
+    ObjectKeeper keeper_datum(datum);
+    ASSERT_EQ(datum, nullptr);
+
+    auto datum_ensemble = proj_crs_get_datum_ensemble(m_ctxt, from_wkt);
+    ObjectKeeper keeper_datum_ensemble(datum_ensemble);
+    ASSERT_NE(datum_ensemble, nullptr);
+
+    ASSERT_EQ(proj_datum_ensemble_get_member_count(m_ctxt, datum_ensemble), 11);
+    ASSERT_EQ(proj_datum_ensemble_get_member(m_ctxt, datum_ensemble, -1),
+              nullptr);
+    ASSERT_EQ(proj_datum_ensemble_get_member(m_ctxt, datum_ensemble, 11),
+              nullptr);
+
+    {
+        auto member = proj_datum_ensemble_get_member(m_ctxt, datum_ensemble, 0);
+        ObjectKeeper keeper_member(member);
+        ASSERT_NE(member, nullptr);
+
+        EXPECT_EQ(proj_get_name(member),
+                  std::string("European Terrestrial Reference Frame 1989"));
+    }
+
+    {
+        auto member =
+            proj_datum_ensemble_get_member(m_ctxt, datum_ensemble, 10);
+        ObjectKeeper keeper_member(member);
+        ASSERT_NE(member, nullptr);
+
+        EXPECT_EQ(proj_get_name(member),
+                  std::string("European Terrestrial Reference Frame 2014"));
+    }
+
+    ASSERT_EQ(proj_datum_ensemble_get_accuracy(m_ctxt, datum_ensemble), 0.1);
+
+    auto datum_forced = proj_crs_get_datum_forced(m_ctxt, from_wkt);
+    ObjectKeeper keeper_datum_forced(datum_forced);
+    ASSERT_NE(datum_forced, nullptr);
+
+    EXPECT_EQ(proj_get_name(datum_forced),
+              std::string("European Terrestrial Reference System 1989"));
+
+    auto cs = proj_crs_get_coordinate_system(m_ctxt, from_wkt);
+    ObjectKeeper keeper_cs(cs);
+    EXPECT_NE(cs, nullptr);
+
+    {
+        auto built_crs = proj_create_geographic_crs_from_datum(
+            m_ctxt, proj_get_name(from_wkt), datum_ensemble, cs);
+        ObjectKeeper keeper_built_crs(built_crs);
+        EXPECT_NE(built_crs, nullptr);
+
+        EXPECT_TRUE(proj_is_equivalent_to_with_ctx(m_ctxt, built_crs, from_wkt,
+                                                   PJ_COMP_EQUIVALENT));
+    }
+
+    {
+        auto built_crs = proj_create_geocentric_crs_from_datum(
+            m_ctxt, proj_get_name(from_wkt), datum_ensemble, "metre", 1.0);
+        ObjectKeeper keeper_built_crs(built_crs);
+        EXPECT_NE(built_crs, nullptr);
+    }
+}
+
 } // namespace