Generalize generalize_proj_crs_create_bound_vertical_crs_to_WGS84()

In recent commits, we added a generalize_proj_crs_create_bound_vertical_crs_to_WGS84()
function, but there are situations where more accurate results can be obtained, if
instead of specifying WGS84 as the hub CRS, the user can specify the exact hub CRS.
For example the GEOID2018 grid is against NAD83(2011).
So replace this function with proj_crs_create_bound_vertical_crs()

4 files changed, 111 insertions, 22 deletions

diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index ceb1ebdf..73a36a70 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/scripts/reference_exported_symbols.txt
@@ -904,7 +904,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_bound_vertical_crs
 proj_crs_create_projected_3D_crs_from_2D
 proj_crs_demote_to_2D
 proj_crs_get_coordinate_system
diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 4fedbe05..f5f7ba55 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -1810,7 +1810,8 @@ 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.
+/** \brief Returns a BoundCRS, with a transformation to a hub geographic 3D crs
+ * (use EPSG:4979 for WGS84 for example), using a grid.
  *
  * The returned object must be unreferenced with proj_destroy() after
  * use.
@@ -1818,33 +1819,43 @@ PJ *proj_crs_create_bound_crs_to_WGS84(PJ_CONTEXT *ctx, const PJ *crs,
  *
  * @param ctx PROJ context, or NULL for default context
  * @param vert_crs Object of type VerticalCRS (must not be NULL)
+ * @param hub_geographic_3D_crs Object of type Geographic 3D CRS (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) {
+PJ *proj_crs_create_bound_vertical_crs(PJ_CONTEXT *ctx, const PJ *vert_crs,
+                                       const PJ *hub_geographic_3D_crs,
+                                       const char *grid_name) {
     SANITIZE_CTX(ctx);
     assert(vert_crs);
+    assert(hub_geographic_3D_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");
+        proj_log_error(ctx, __FUNCTION__, "vert_crs is not a VerticalCRS");
+        return nullptr;
+    }
+    auto hub_crs =
+        std::dynamic_pointer_cast<CRS>(hub_geographic_3D_crs->iso_obj);
+    if (!hub_crs) {
+        proj_log_error(ctx, __FUNCTION__, "hub_geographic_3D_crs is not a CRS");
         return nullptr;
     }
     try {
         auto nnCRS = NN_NO_CHECK(l_crs);
+        auto nnHubCRS = NN_NO_CHECK(hub_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));
+                                  "unknown to " + hub_crs->nameStr() +
+                                      " ellipsoidal height"),
+                nnCRS, nnHubCRS, nullptr, std::string(grid_name),
+                std::vector<PositionalAccuracyNNPtr>());
+        return pj_obj_create(ctx,
+                             BoundCRS::create(nnCRS, nnHubCRS, transformation));
     } catch (const std::exception &e) {
         proj_log_error(ctx, __FUNCTION__, e.what());
         if (ctx->cpp_context) {
diff --git a/src/proj_experimental.h b/src/proj_experimental.h
index 63b858d5..c6d5bc45 100644
--- a/src/proj_experimental.h
+++ b/src/proj_experimental.h
@@ -317,9 +317,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);
+PJ PROJ_DLL *proj_crs_create_bound_vertical_crs(PJ_CONTEXT *ctx,
+                                                const PJ* vert_crs,
+                                                const PJ* hub_geographic_3D_crs,
+                                                const char* grid_name);
 
 /* BEGIN: Generated by scripts/create_c_api_projections.py*/
 PJ PROJ_DLL *proj_create_conversion_utm(
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index f8b447ca..3b6a02f2 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -4019,15 +4019,23 @@ TEST_F(CApi, proj_crs_create_projected_3D_crs_from_2D) {
 
 // ---------------------------------------------------------------------------
 
-TEST_F(CApi, proj_crs_create_bound_vertical_crs_to_WGS84) {
+TEST_F(CApi, proj_crs_create_bound_vertical_crs) {
 
     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");
+    auto crs4979 = proj_create_from_wkt(
+        m_ctxt,
+        GeographicCRS::EPSG_4979->exportToWKT(WKTFormatter::create().get())
+            .c_str(),
+        nullptr, nullptr, nullptr);
+    ObjectKeeper keeper_crs4979(crs4979);
+    ASSERT_NE(crs4979, nullptr);
+
+    auto bound_crs = proj_crs_create_bound_vertical_crs(m_ctxt, vert_crs,
+                                                        crs4979, "foo.gtx");
     ObjectKeeper keeper_bound_crs(bound_crs);
     ASSERT_NE(bound_crs, nullptr);
 
@@ -4078,7 +4086,7 @@ TEST_F(CApi, proj_create_crs_to_crs_with_only_ballpark_transformations) {
 
 TEST_F(
     CApi,
-    proj_create_crs_to_crs_from_custom_compound_crs_with_NAD83_2011_and_geoidgrid_to_WGS84_G1762) {
+    proj_create_crs_to_crs_from_custom_compound_crs_with_NAD83_2011_and_geoidgrid_ref_against_WGS84_to_WGS84_G1762) {
 
     if (strcmp(proj_grid_info("egm96_15.gtx").format, "missing") == 0) {
         return; // use GTEST_SKIP() if we upgrade gtest
@@ -4094,10 +4102,15 @@ TEST_F(
                                               "DummyDatum", "metre", 1.0);
     ASSERT_NE(inDummyCrs, nullptr);
 
-    PJ *inCrsV = proj_crs_create_bound_vertical_crs_to_WGS84(m_ctxt, inDummyCrs,
-                                                             "egm96_15.gtx");
+    auto crs4979 = proj_create_from_database(m_ctxt, "EPSG", "4979",
+                                           PJ_CATEGORY_CRS, false, nullptr);
+    ASSERT_NE(crs4979, nullptr);
+
+    PJ *inCrsV = proj_crs_create_bound_vertical_crs(m_ctxt, inDummyCrs, crs4979,
+                                                    "egm96_15.gtx");
     ASSERT_NE(inCrsV, nullptr);
     proj_destroy(inDummyCrs);
+    proj_destroy(crs4979);
 
     PJ *inCompound =
         proj_create_compound_crs(m_ctxt, "Compound", inCrsH, inCrsV);
@@ -4115,7 +4128,7 @@ TEST_F(
     // as ballpark. That one used to be eliminated because by
     // proj_create_crs_to_crs() because there were non Ballpark transformations
     // available. This resulted thus in an error when transforming outside of
-    // those few subzones.s
+    // those few subzones.
     P = proj_create_crs_to_crs_from_pj(m_ctxt, inCompound, outCrs, nullptr,
                                        nullptr);
     ASSERT_NE(P, nullptr);
@@ -4136,4 +4149,68 @@ TEST_F(
     EXPECT_NEAR(outcoord.xyzt.z, 118.059, 1e-3);
 }
 
+// ---------------------------------------------------------------------------
+
+TEST_F(
+    CApi,
+    proj_create_crs_to_crs_from_custom_compound_crs_with_NAD83_2011_and_geoidgrid_ref_against_NAD83_2011_to_WGS84_G1762) {
+
+    if (strcmp(proj_grid_info("egm96_15.gtx").format, "missing") == 0) {
+        return; // use GTEST_SKIP() if we upgrade gtest
+    }
+
+    PJ *P;
+
+    // NAD83(2011) 2D
+    PJ *inCrsH = proj_create_from_database(m_ctxt, "EPSG", "6318",
+                                           PJ_CATEGORY_CRS, false, nullptr);
+    ASSERT_NE(inCrsH, nullptr);
+
+    PJ *inDummyCrs = proj_create_vertical_crs(m_ctxt, "VerticalDummyCrs",
+                                              "DummyDatum", "metre", 1.0);
+    ASSERT_NE(inDummyCrs, nullptr);
+
+    // NAD83(2011) 3D
+    PJ *inGeog3DCRS = proj_create_from_database(
+        m_ctxt, "EPSG", "6319", PJ_CATEGORY_CRS, false, nullptr);
+    ASSERT_NE(inCrsH, nullptr);
+
+    // Note: this is actually a bad example since we tell here that egm96_15.gtx
+    // is referenced against NAD83(2011)
+    PJ *inCrsV = proj_crs_create_bound_vertical_crs(
+        m_ctxt, inDummyCrs, inGeog3DCRS, "egm96_15.gtx");
+    ASSERT_NE(inCrsV, nullptr);
+    proj_destroy(inDummyCrs);
+    proj_destroy(inGeog3DCRS);
+
+    PJ *inCompound =
+        proj_create_compound_crs(m_ctxt, "Compound", inCrsH, inCrsV);
+    ASSERT_NE(inCompound, nullptr);
+    proj_destroy(inCrsH);
+    proj_destroy(inCrsV);
+
+    // WGS84 (G1762)
+    PJ *outCrs = proj_create(m_ctxt, "EPSG:7665");
+    ASSERT_NE(outCrs, nullptr);
+
+    P = proj_create_crs_to_crs_from_pj(m_ctxt, inCompound, outCrs, nullptr,
+                                       nullptr);
+    ASSERT_NE(P, nullptr);
+    proj_destroy(inCompound);
+    proj_destroy(outCrs);
+
+    PJ_COORD in_coord;
+    in_coord.xyzt.x = 35;
+    in_coord.xyzt.y = -118;
+    in_coord.xyzt.z = 0;
+    in_coord.xyzt.t = 2010;
+
+    PJ_COORD outcoord = proj_trans(P, PJ_FWD, in_coord);
+    proj_destroy(P);
+
+    EXPECT_NEAR(outcoord.xyzt.x, 35.000003665064803, 1e-9);
+    EXPECT_NEAR(outcoord.xyzt.y, -118.00001414221214, 1e-9);
+    EXPECT_NEAR(outcoord.xyzt.z, -32.5823, 1e-3);
+}
+
 } // namespace