Merge pull request #1689 from rouault/generalize_proj_crs_create_bound_vertical_crs_to_WGS84

Generalize generalize_proj_crs_create_bound_vertical_crs_to_WGS84()

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