C API: add proj_crs_create_bound_vertical_crs_to_WGS84()

4 files changed, 82 insertions, 0 deletions

diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index 3206c3c8..b12907e5 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/scripts/reference_exported_symbols.txt
@@ -900,6 +900,7 @@ proj_crs_alter_geodetic_crs
 proj_crs_alter_parameters_linear_unit
 proj_crs_create_bound_crs
 proj_crs_create_bound_crs_to_WGS84
+proj_crs_create_bound_vertical_crs_to_WGS84
 proj_crs_create_projected_3D_crs_from_2D
 proj_crs_get_coordinate_system
 proj_crs_get_coordoperation
diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 3a1b66af..fd4090b2 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -1810,6 +1810,52 @@ PJ *proj_crs_create_bound_crs_to_WGS84(PJ_CONTEXT *ctx, const PJ *crs,
 
 // ---------------------------------------------------------------------------
 
+/** \brief Returns a BoundCRS, with a transformation to EPSG:4979 using a grid.
+ *
+ * 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 vert_crs Object of type VerticalCRS (must not be NULL)
+ * @param grid_name Grid name (typically a .gtx file)
+ * @return Object that must be unreferenced with proj_destroy(), or NULL
+ * in case of error.
+ * @since 7.0
+ */
+PJ *proj_crs_create_bound_vertical_crs_to_WGS84(PJ_CONTEXT *ctx,
+                                                const PJ *vert_crs,
+                                                const char *grid_name) {
+    SANITIZE_CTX(ctx);
+    assert(vert_crs);
+    assert(grid_name);
+    auto l_crs = std::dynamic_pointer_cast<VerticalCRS>(vert_crs->iso_obj);
+    if (!l_crs) {
+        proj_log_error(ctx, __FUNCTION__, "Object is not a VerticalCRS");
+        return nullptr;
+    }
+    try {
+        auto nnCRS = NN_NO_CHECK(l_crs);
+        auto transformation =
+            Transformation::createGravityRelatedHeightToGeographic3D(
+                PropertyMap().set(IdentifiedObject::NAME_KEY,
+                                  "unknown to WGS84 ellipsoidal height"),
+                nnCRS, GeographicCRS::EPSG_4979, nullptr,
+                std::string(grid_name), std::vector<PositionalAccuracyNNPtr>());
+        return pj_obj_create(
+            ctx,
+            BoundCRS::create(nnCRS, GeographicCRS::EPSG_4979, transformation));
+    } catch (const std::exception &e) {
+        proj_log_error(ctx, __FUNCTION__, e.what());
+        if (ctx->cpp_context) {
+            ctx->cpp_context->autoCloseDbIfNeeded();
+        }
+        return nullptr;
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Get the ellipsoid from a CRS or a GeodeticReferenceFrame.
  *
  * The returned object must be unreferenced with proj_destroy() after
diff --git a/src/proj_experimental.h b/src/proj_experimental.h
index 5e6bbb13..0452eb4b 100644
--- a/src/proj_experimental.h
+++ b/src/proj_experimental.h
@@ -313,6 +313,10 @@ PJ PROJ_DLL *proj_crs_create_bound_crs_to_WGS84(PJ_CONTEXT *ctx,
                                                         const PJ *crs,
                                                         const char *const *options);
 
+PJ PROJ_DLL *proj_crs_create_bound_vertical_crs_to_WGS84(PJ_CONTEXT *ctx,
+                                                         const PJ* vert_crs,
+                                                         const char* grid_name);
+
 /* BEGIN: Generated by scripts/create_c_api_projections.py*/
 PJ PROJ_DLL *proj_create_conversion_utm(
     PJ_CONTEXT *ctx,
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index b8dde430..3a1c5df5 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -3966,4 +3966,35 @@ TEST_F(CApi, proj_crs_create_projected_3D_crs_from_2D) {
     }
 }
 
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_crs_create_bound_vertical_crs_to_WGS84) {
+
+    auto vert_crs = proj_create_vertical_crs(m_ctxt, "myVertCRS", "myVertDatum",
+                                             nullptr, 0.0);
+    ObjectKeeper keeper_vert_crs(vert_crs);
+    ASSERT_NE(vert_crs, nullptr);
+
+    auto bound_crs = proj_crs_create_bound_vertical_crs_to_WGS84(
+        m_ctxt, vert_crs, "foo.gtx");
+    ObjectKeeper keeper_bound_crs(bound_crs);
+    ASSERT_NE(bound_crs, nullptr);
+
+    auto projCRS = proj_create_from_database(m_ctxt, "EPSG", "32631",
+                                             PJ_CATEGORY_CRS, false, nullptr);
+    ASSERT_NE(projCRS, nullptr);
+    ObjectKeeper keeper_projCRS(projCRS);
+
+    auto compound_crs =
+        proj_create_compound_crs(m_ctxt, "myCompoundCRS", projCRS, bound_crs);
+    ObjectKeeper keeper_compound_crss(compound_crs);
+    ASSERT_NE(compound_crs, nullptr);
+
+    auto proj_4 = proj_as_proj_string(m_ctxt, compound_crs, PJ_PROJ_4, nullptr);
+    ASSERT_NE(proj_4, nullptr);
+    EXPECT_EQ(std::string(proj_4),
+              "+proj=utm +zone=31 +datum=WGS84 +units=m +geoidgrids=foo.gtx "
+              "+vunits=m +no_defs +type=crs");
+}
+
 } // namespace