[6.2 backport] cs2cs: autopromote CRS to 3D when there's a mix… (#1615)

[6.2 backport] cs2cs: autopromote CRS to 3D when there's a mix of 2D and 3D (fixes #1563)

5 files changed, 143 insertions, 1 deletions

diff --git a/test/cli/Makefile.am b/test/cli/Makefile.am
index 47cb1e7f..db47370c 100644
--- a/test/cli/Makefile.am
+++ b/test/cli/Makefile.am
@@ -39,7 +39,7 @@ testvarious-check:
 	PROJ_LIB=$(PROJ_LIB) $(TESTVARIOUS) $(CS2CSEXE)
 
 testdatumfile-check:
-	@if [ -f $(PROJ_LIB)/conus -a -f $(PROJ_LIB)/ntv1_can.dat -a -f $(PROJ_LIB)/MD -a -f $(PROJ_LIB)/ntf_r93.gsb ]; then \
+	@if [ -f $(PROJ_LIB)/conus -a -f $(PROJ_LIB)/ntv1_can.dat -a -f $(PROJ_LIB)/MD -a -f $(PROJ_LIB)/ntf_r93.gsb -a -f $(PROJ_LIB)/egm96_15.gtx ]; then \
 	  PROJ_LIB=$(PROJ_LIB) $(TESTDATUMFILE) $(CS2CSEXE) ; \
 	fi
 
diff --git a/test/cli/td_out.dist b/test/cli/td_out.dist
index cdbe0fa9..179ded11 100644
--- a/test/cli/td_out.dist
+++ b/test/cli/td_out.dist
@@ -27,3 +27,6 @@ edge or even a wee bit outside (#141).
 NAD27 -> NAD83: 1st through ntv1, 2nd through conus
 44d00'00.000"N 111d00'00.000"W 0.0	43d59'59.732"N	111d0'3.208"W 0.000
 39d00'00.000"N 111d00'00.000"W 0.0	38d59'59.912"N	111d0'2.604"W 0.000
+##############################################################
+WGS84 -> WGS84+EGM96
+2dE 49dN 0	2.00	49.00 -44.64
diff --git a/test/cli/testdatumfile b/test/cli/testdatumfile
index 5b56f077..45c29428 100755
--- a/test/cli/testdatumfile
+++ b/test/cli/testdatumfile
@@ -105,6 +105,14 @@ $EXE NAD27 NAD83 -E >>${OUT} <<EOF
 39d00'00.000"N 111d00'00.000"W 0.0
 EOF
 
+#
+echo "##############################################################" >> ${OUT}
+echo "WGS84 -> WGS84+EGM96" >> ${OUT}
+#
+$EXE +init=epsg:4326 +to +init=epsg:4326 +geoidgrids=egm96_15.gtx -E >>${OUT} <<EOF
+2dE 49dN 0
+EOF
+
 # Cleanup
 rm -rf "dir with \" space"
 
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 2a1c8577..db4b768b 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -3737,4 +3737,78 @@ TEST_F(CApi, proj_create_crs_to_crs_from_pj) {
               "+step +proj=utm +zone=31 +ellps=WGS84");
 }
 
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_crs_promote_to_3D) {
+
+    auto crs2D =
+        proj_create(m_ctxt, GeographicCRS::EPSG_4326
+                                ->exportToWKT(WKTFormatter::create().get())
+                                .c_str());
+    ObjectKeeper keeper_crs2D(crs2D);
+    EXPECT_NE(crs2D, nullptr);
+
+    auto crs3D = proj_crs_promote_to_3D(m_ctxt, nullptr, crs2D);
+    ObjectKeeper keeper_crs3D(crs3D);
+    EXPECT_NE(crs3D, nullptr);
+
+    auto cs = proj_crs_get_coordinate_system(m_ctxt, crs3D);
+    ASSERT_NE(cs, nullptr);
+    ObjectKeeper keeperCs(cs);
+    EXPECT_EQ(proj_cs_get_axis_count(m_ctxt, cs), 3);
+
+    auto code = proj_get_id_code(crs3D, 0);
+    ASSERT_TRUE(code != nullptr);
+    EXPECT_EQ(code, std::string("4979"));
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_crs_create_projected_3D_crs_from_2D) {
+
+    auto projCRS = proj_create_from_database(m_ctxt, "EPSG", "32631",
+                                             PJ_CATEGORY_CRS, false, nullptr);
+    ASSERT_NE(projCRS, nullptr);
+    ObjectKeeper keeper_projCRS(projCRS);
+
+    {
+        auto geog3DCRS = proj_create_from_database(
+            m_ctxt, "EPSG", "4979", PJ_CATEGORY_CRS, false, nullptr);
+        ASSERT_NE(geog3DCRS, nullptr);
+        ObjectKeeper keeper_geog3DCRS(geog3DCRS);
+
+        auto crs3D = proj_crs_create_projected_3D_crs_from_2D(
+            m_ctxt, nullptr, projCRS, geog3DCRS);
+        ObjectKeeper keeper_crs3D(crs3D);
+        EXPECT_NE(crs3D, nullptr);
+
+        EXPECT_EQ(proj_get_type(crs3D), PJ_TYPE_PROJECTED_CRS);
+
+        EXPECT_EQ(std::string(proj_get_name(crs3D)),
+                  std::string(proj_get_name(projCRS)));
+
+        auto cs = proj_crs_get_coordinate_system(m_ctxt, crs3D);
+        ASSERT_NE(cs, nullptr);
+        ObjectKeeper keeperCs(cs);
+        EXPECT_EQ(proj_cs_get_axis_count(m_ctxt, cs), 3);
+    }
+
+    {
+        auto crs3D = proj_crs_create_projected_3D_crs_from_2D(m_ctxt, nullptr,
+                                                              projCRS, nullptr);
+        ObjectKeeper keeper_crs3D(crs3D);
+        EXPECT_NE(crs3D, nullptr);
+
+        EXPECT_EQ(proj_get_type(crs3D), PJ_TYPE_PROJECTED_CRS);
+
+        EXPECT_EQ(std::string(proj_get_name(crs3D)),
+                  std::string(proj_get_name(projCRS)));
+
+        auto cs = proj_crs_get_coordinate_system(m_ctxt, crs3D);
+        ASSERT_NE(cs, nullptr);
+        ObjectKeeper keeperCs(cs);
+        EXPECT_EQ(proj_cs_get_axis_count(m_ctxt, cs), 3);
+    }
+}
+
 } // namespace
diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp
index 0f0304b1..0df128e6 100644
--- a/test/unit/test_crs.cpp
+++ b/test/unit/test_crs.cpp
@@ -5224,3 +5224,60 @@ TEST(crs, getNonDeprecated) {
         ASSERT_EQ(list.size(), 1U);
     }
 }
+
+// ---------------------------------------------------------------------------
+
+TEST(crs, promoteTo3D) {
+    auto dbContext = DatabaseContext::create();
+    {
+        auto crs = GeographicCRS::EPSG_4326;
+        auto crs3D = crs->promoteTo3D(std::string(), nullptr);
+        auto crs3DAsGeog = nn_dynamic_pointer_cast<GeographicCRS>(crs3D);
+        ASSERT_TRUE(crs3DAsGeog != nullptr);
+        EXPECT_EQ(crs3DAsGeog->coordinateSystem()->axisList().size(), 3U);
+        EXPECT_TRUE(crs3D->promoteTo3D(std::string(), nullptr)
+                        ->isEquivalentTo(crs3D.get()));
+    }
+    {
+        auto crs = GeographicCRS::EPSG_4326;
+        auto crs3D = crs->promoteTo3D(std::string(), dbContext);
+        auto crs3DAsGeog = nn_dynamic_pointer_cast<GeographicCRS>(crs3D);
+        ASSERT_TRUE(crs3DAsGeog != nullptr);
+        EXPECT_EQ(crs3DAsGeog->coordinateSystem()->axisList().size(), 3U);
+        EXPECT_TRUE(!crs3DAsGeog->identifiers().empty());
+    }
+    {
+        auto crs = createProjected();
+        auto crs3D = crs->promoteTo3D(std::string(), nullptr);
+        auto crs3DAsProjected = nn_dynamic_pointer_cast<ProjectedCRS>(crs3D);
+        ASSERT_TRUE(crs3DAsProjected != nullptr);
+        EXPECT_EQ(crs3DAsProjected->coordinateSystem()->axisList().size(), 3U);
+        EXPECT_EQ(
+            crs3DAsProjected->baseCRS()->coordinateSystem()->axisList().size(),
+            3U);
+        EXPECT_TRUE(crs3D->promoteTo3D(std::string(), nullptr)
+                        ->isEquivalentTo(crs3D.get()));
+    }
+    {
+        auto crs = createProjected();
+        auto crs3D = crs->promoteTo3D(std::string(), dbContext);
+        auto crs3DAsProjected = nn_dynamic_pointer_cast<ProjectedCRS>(crs3D);
+        ASSERT_TRUE(crs3DAsProjected != nullptr);
+        EXPECT_EQ(crs3DAsProjected->coordinateSystem()->axisList().size(), 3U);
+        EXPECT_EQ(
+            crs3DAsProjected->baseCRS()->coordinateSystem()->axisList().size(),
+            3U);
+        EXPECT_TRUE(!crs3DAsProjected->baseCRS()->identifiers().empty());
+    }
+    {
+        auto crs = BoundCRS::createFromTOWGS84(
+            createProjected(), std::vector<double>{1, 2, 3, 4, 5, 6, 7});
+        auto crs3D = crs->promoteTo3D(std::string(), dbContext);
+        auto crs3DAsBound = nn_dynamic_pointer_cast<BoundCRS>(crs3D);
+        ASSERT_TRUE(crs3DAsBound != nullptr);
+        auto baseCRS =
+            nn_dynamic_pointer_cast<ProjectedCRS>(crs3DAsBound->baseCRS());
+        ASSERT_TRUE(baseCRS != nullptr);
+        EXPECT_EQ(baseCRS->coordinateSystem()->axisList().size(), 3U);
+    }
+}