Merge pull request #3119 from rouault/compound_to_2D_crs

Transformation: no longer do vertical trasnformation when doing compound CRS to 2D CRS / add --3d to cs2cs

11 files changed, 85 insertions, 33 deletions

diff --git a/docs/source/apps/cs2cs.rst b/docs/source/apps/cs2cs.rst
index 841ce092..8f35c751 100644
--- a/docs/source/apps/cs2cs.rst
+++ b/docs/source/apps/cs2cs.rst
@@ -13,7 +13,7 @@ Synopsis
 
     | **cs2cs** [**-eEfIlrstvwW** [args]]
     |           [[--area <name_or_code>] | [--bbox <west_long,south_lat,east_long,north_lat>]]
-    |           [--authority <name>] [--no-ballpark] [--accuracy <accuracy>]
+    |           [--authority <name>] [--no-ballpark] [--accuracy <accuracy>] [--3d]
     |           ([*+opt[=arg]* ...] [+to *+opt[=arg]* ...] | {source_crs} {target_crs})
     |           file ...
 
@@ -194,6 +194,17 @@ The following control parameters can appear in any order:
 
     This option is mutually exclusive with :option:`--bbox`.
 
+.. option:: --3d
+
+    .. versionadded:: 9.1
+
+    "Promote" 2D CRS(s) to their 3D version, where the vertical axis is the
+    ellipsoidal height in metres, using the ellipsoid of the base geodetic CRS.
+    Depending on PROJ versions and the exact nature of the CRS involved,
+    especially before PROJ 9.1, a mix of 2D and 3D CRS could lead to 2D or 3D
+    transformations. Starting with PROJ 9.1, both CRS need to be 3D for vertical
+    transformation to possibly happen.
+
 .. only:: man
 
     The *+opt* run-line arguments are associated with cartographic
diff --git a/docs/source/apps/projinfo.rst b/docs/source/apps/projinfo.rst
index 04ee7578..a19afa5b 100644
--- a/docs/source/apps/projinfo.rst
+++ b/docs/source/apps/projinfo.rst
@@ -300,12 +300,12 @@ The following control parameters can appear in any order:
 
     .. versionadded:: 6.3
 
-    "Promote" the CRS(s) to their 3D version. In the context of researching
-    available coordinate transformations, explicitly specifying this option is
-    not necessary, because when one of the source or target CRS has a vertical
-    component but not the other one, the one that has no vertical component is
-    automatically promoted to a 3D version, where its vertical axis is the
+    "Promote" 2D CRS(s) to their 3D version, where the vertical axis is the
     ellipsoidal height in metres, using the ellipsoid of the base geodetic CRS.
+    Depending on PROJ versions and the exact nature of the CRS involved,
+    especially before PROJ 9.1, a mix of 2D and 3D CRS could lead to 2D or 3D
+    transformations. Starting with PROJ 9.1, both CRS need to be 3D for vertical
+    transformation to possibly happen.
 
 .. option:: --output-id=AUTH:NAME
 
diff --git a/src/apps/cs2cs.cpp b/src/apps/cs2cs.cpp
index a8e152a7..e21ecd8b 100644
--- a/src/apps/cs2cs.cpp
+++ b/src/apps/cs2cs.cpp
@@ -78,7 +78,7 @@ static const char *oterr = "*\t*"; /* output line for unprojectable input */
 static const char *usage =
     "%s\nusage: %s [-dDeEfIlrstvwW [args]]\n"
     "              [[--area name_or_code] | [--bbox west_long,south_lat,east_long,north_lat]]\n"
-    "              [--authority {name}] [--accuracy {accuracy}] [--no-ballpark]\n"
+    "              [--authority {name}] [--accuracy {accuracy}] [--no-ballpark] [--3d]\n"
     "              [+opt[=arg] ...] [+to +opt[=arg] ...] [file ...]\n";
 
 static double (*informat)(const char *,
@@ -384,6 +384,7 @@ int main(int argc, char **argv) {
     const char* authority = nullptr;
     double accuracy = -1;
     bool allowBallpark = true;
+    bool promoteTo3D = false;
 
     /* process run line arguments */
     while (--argc > 0) { /* collect run line arguments */
@@ -449,6 +450,9 @@ int main(int argc, char **argv) {
         else if (strcmp(*argv, "--no-ballpark") == 0 ) {
             allowBallpark = false;
         }
+        else if (strcmp(*argv, "--3d") == 0 ) {
+            promoteTo3D = true;
+        }
         else if (**argv == '-') {
             for (arg = *argv;;) {
                 switch (*++arg) {
@@ -785,8 +789,7 @@ int main(int argc, char **argv) {
     src = proj_create(nullptr, pj_add_type_crs_if_needed(fromStr).c_str());
     dst = proj_create(nullptr, pj_add_type_crs_if_needed(toStr).c_str());
 
-    if( proj_get_type(src) == PJ_TYPE_COMPOUND_CRS ||
-        proj_get_type(dst) == PJ_TYPE_COMPOUND_CRS ) {
+    if( promoteTo3D ) {
         auto src3D = proj_crs_promote_to_3D(nullptr, nullptr, src);
         if( src3D ) {
             proj_destroy(src);
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index d4c6aec1..e09aee93 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -10525,12 +10525,18 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
                     }
                     auto geogCRS = dynamic_cast<GeographicCRS *>(crs);
                     if (geogCRS) {
+                        const auto &cs = geogCRS->coordinateSystem();
                         // Override with longitude latitude in degrees
                         return GeographicCRS::create(
                             properties, geogCRS->datum(),
                             geogCRS->datumEnsemble(),
-                            EllipsoidalCS::createLongitudeLatitude(
-                                UnitOfMeasure::DEGREE));
+                            cs->axisList().size() == 2
+                                ? EllipsoidalCS::createLongitudeLatitude(
+                                      UnitOfMeasure::DEGREE)
+                                : EllipsoidalCS::
+                                      createLongitudeLatitudeEllipsoidalHeight(
+                                          UnitOfMeasure::DEGREE,
+                                          cs->axisList()[2]->unit()));
                     }
                     auto projCRS = dynamic_cast<ProjectedCRS *>(crs);
                     if (projCRS) {
diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp
index 20042f22..a6a2c986 100644
--- a/src/iso19111/operation/coordinateoperationfactory.cpp
+++ b/src/iso19111/operation/coordinateoperationfactory.cpp
@@ -4926,6 +4926,19 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
             }
         }
 
+        // Only do a vertical transformation if the target CRS is 3D.
+        const auto dstSingle = dynamic_cast<crs::SingleCRS *>(targetCRS.get());
+        if (dstSingle &&
+            dstSingle->coordinateSystem()->axisList().size() == 2) {
+            auto tmp = createOperations(componentsSrc[0], targetCRS, context);
+            for (const auto &op : tmp) {
+                auto opClone = op->shallowClone();
+                setCRSs(opClone.get(), sourceCRS, targetCRS);
+                res.emplace_back(opClone);
+            }
+            return;
+        }
+
         std::vector<CoordinateOperationNNPtr> horizTransforms;
         auto srcGeogCRS = componentsSrc[0]->extractGeographicCRS();
         if (srcGeogCRS) {
diff --git a/test/cli/td_out.dist b/test/cli/td_out.dist
index 5c9b18d3..5957ebc3 100644
--- a/test/cli/td_out.dist
+++ b/test/cli/td_out.dist
@@ -35,5 +35,11 @@ NAD27 -> NAD83: 1st through ntv1 or ntv2, 2nd through conus
 55d00'00.000"N 111d00'00.000"W 0.0	55.0001	-111.0009 0.0000
 39d00'00.000"N 111d00'00.000"W 0.0	39.0000	-111.0007 0.0000
 ##############################################################
-WGS84 -> WGS84+EGM96
+WGS84 (2D) -> WGS84+EGM96
+2dE 49dN 0	2.00	49.00 0.00
+##############################################################
+WGS84 (3D) -> WGS84+EGM96
+2dE 49dN 0	2.00	49.00 -45.06
+##############################################################
+WGS84 (2D), promoted to 3D -> WGS84+EGM96
 2dE 49dN 0	2.00	49.00 -45.06
diff --git a/test/cli/testdatumfile b/test/cli/testdatumfile
index e974d34c..3496299f 100755
--- a/test/cli/testdatumfile
+++ b/test/cli/testdatumfile
@@ -124,12 +124,28 @@ EOF
 
 #
 echo "##############################################################" >> ${OUT}
-echo "WGS84 -> WGS84+EGM96" >> ${OUT}
+echo "WGS84 (2D) -> WGS84+EGM96" >> ${OUT}
 #
 $EXE +init=epsg:4326 +to +init=epsg:4326 +geoidgrids=egm96_15.gtx -E >>${OUT} <<EOF
 2dE 49dN 0
 EOF
 
+#
+echo "##############################################################" >> ${OUT}
+echo "WGS84 (3D) -> WGS84+EGM96" >> ${OUT}
+#
+$EXE +init=epsg:4979 +to +init=epsg:4326 +geoidgrids=egm96_15.gtx -E >>${OUT} <<EOF
+2dE 49dN 0
+EOF
+
+#
+echo "##############################################################" >> ${OUT}
+echo "WGS84 (2D), promoted to 3D -> WGS84+EGM96" >> ${OUT}
+#
+$EXE --3d +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/cli/testprojinfo b/test/cli/testprojinfo
index 2f8c1de5..016850a7 100755
--- a/test/cli/testprojinfo
+++ b/test/cli/testprojinfo
@@ -148,8 +148,8 @@ echo "Testing -s EPSG:4936 -t EPSG:4978 --spatial-test intersects --summary wher
 $EXE -s EPSG:4936 -t EPSG:4978 --spatial-test intersects --summary >>${OUT} 2>&1
 echo "" >>${OUT}
 
-echo "Testing -s "+proj=longlat +datum=WGS84 +geoidgrids=@foo.gtx +type=crs" -t EPSG:4326 -o PROJ -q" >> ${OUT}
-$EXE -s "+proj=longlat +datum=WGS84 +geoidgrids=@foo.gtx +type=crs" -t EPSG:4326 -o PROJ -q >>${OUT} 2>&1
+echo "Testing -s "+proj=longlat +datum=WGS84 +geoidgrids=@foo.gtx +type=crs" -t EPSG:4979 -o PROJ -q" >> ${OUT}
+$EXE -s "+proj=longlat +datum=WGS84 +geoidgrids=@foo.gtx +type=crs" -t EPSG:4979 -o PROJ -q >>${OUT} 2>&1
 echo "" >>${OUT}
 
 echo 'Testing -s "AGD66" -t "WGS 84 (G1762)" --spatial-test intersects --summary. Should include a transformation through GDA2020' >> ${OUT}
diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist
index cb80bdd0..952c6b38 100644
--- a/test/cli/testprojinfo_out.dist
+++ b/test/cli/testprojinfo_out.dist
@@ -1243,10 +1243,8 @@ PROJCRS["WGS 84 / UTM zone 31N",
     REMARK["Promoted to 3D from EPSG:32631"]]
 
 Testing -s EPSG:32631 -t EPSG:4326+3855 --summary
-Candidate operations found: 3
-unknown id, Inverse of UTM zone 31N + WGS 84 to EGM2008 height (1), 1 m, World.
-unknown id, Inverse of UTM zone 31N + WGS 84 to EGM2008 height (2), 0.5 m, World.
-unknown id, Inverse of UTM zone 31N + Inverse of Transformation from EGM2008 height to WGS 84 (ballpark vertical transformation, without ellipsoid height to vertical height correction), unknown accuracy, World, has ballpark transformation
+Candidate operations found: 1
+unknown id, Inverse of UTM zone 31N + Inverse of Null geographic offset from WGS 84 to WGS 84, 0 m, World.
 
 Testing -s EPSG:32631 -t EPSG:4326+3855 --3d --summary
 Candidate operations found: 3
@@ -1262,7 +1260,7 @@ Candidate operations found: 2
 unknown id, Ballpark geocentric translation from ETRS89 to WGS 84, unknown accuracy, World, has ballpark transformation
 INVERSE(EPSG):9225, Inverse of WGS 84 to ETRS89 (2), 0.1 m, Germany - offshore North Sea. Netherlands - offshore east of 5E.
 
-Testing -s +proj=longlat +datum=WGS84 +geoidgrids=@foo.gtx +type=crs -t EPSG:4326 -o PROJ -q
+Testing -s +proj=longlat +datum=WGS84 +geoidgrids=@foo.gtx +type=crs -t EPSG:4979 -o PROJ -q
 +proj=pipeline
   +step +proj=unitconvert +xy_in=deg +xy_out=rad
   +step +proj=vgridshift +grids=@foo.gtx +multiplier=1
diff --git a/test/unit/test_network.cpp b/test/unit/test_network.cpp
index f8b70ae8..73265008 100644
--- a/test/unit/test_network.cpp
+++ b/test/unit/test_network.cpp
@@ -1146,7 +1146,7 @@ TEST(networking, curl_vgridshift) {
 
     // WGS84 to EGM2008 height. Using egm08_25.tif
     auto P =
-        proj_create_crs_to_crs(ctx, "EPSG:4326", "EPSG:4326+3855", nullptr);
+        proj_create_crs_to_crs(ctx, "EPSG:4979", "EPSG:4326+3855", nullptr);
     ASSERT_NE(P, nullptr);
 
     PJ_COORD c;
diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
index 0c3ecae4..bb69e347 100644
--- a/test/unit/test_operationfactory.cpp
+++ b/test/unit/test_operationfactory.cpp
@@ -5441,7 +5441,9 @@ TEST(operation, compoundCRS_to_geogCRS_2D_promote_to_3D_context) {
     auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
     ASSERT_TRUE(src != nullptr);
     auto nnSrc = NN_NO_CHECK(src);
-    auto dst = authFactory->createCoordinateReferenceSystem("4269"); // NAD83
+    auto dst =
+        authFactory->createCoordinateReferenceSystem("4269")->promoteTo3D(
+            std::string(), authFactory->databaseContext()); // NAD83
 
     auto listCompoundToGeog2D =
         CoordinateOperationFactory::create()->createOperations(nnSrc, dst,
@@ -5454,18 +5456,11 @@ TEST(operation, compoundCRS_to_geogCRS_2D_promote_to_3D_context) {
         CoordinateOperationFactory::create()->createOperations(dst, nnSrc,
                                                                ctxt);
     EXPECT_EQ(listGeog2DToCompound.size(), listCompoundToGeog2D.size());
-
-    auto listCompoundToGeog3D =
-        CoordinateOperationFactory::create()->createOperations(
-            nnSrc,
-            dst->promoteTo3D(std::string(), authFactory->databaseContext()),
-            ctxt);
-    EXPECT_EQ(listCompoundToGeog3D.size(), listCompoundToGeog2D.size());
 }
 
 // ---------------------------------------------------------------------------
 
-TEST(operation, compoundCRS_of_projCRS_to_geogCRS_2D_context) {
+TEST(operation, compoundCRS_of_projCRS_to_geogCRS_3D_context) {
     auto authFactory =
         AuthorityFactory::create(DatabaseContext::create(), "EPSG");
     auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
@@ -5480,7 +5475,9 @@ TEST(operation, compoundCRS_of_projCRS_to_geogCRS_2D_context) {
     auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
     ASSERT_TRUE(src != nullptr);
     auto nnSrc = NN_NO_CHECK(src);
-    auto dst = authFactory->createCoordinateReferenceSystem("4269"); // NAD83
+    auto dst =
+        authFactory->createCoordinateReferenceSystem("4269")->promoteTo3D(
+            std::string(), authFactory->databaseContext()); // NAD83
 
     auto list = CoordinateOperationFactory::create()->createOperations(
         nnSrc, dst, ctxt);
@@ -5945,7 +5942,9 @@ TEST(operation, compoundCRS_of_vertCRS_with_geoid_model_to_geogCRS) {
         createFromUserInput(wkt, authFactory->databaseContext(), false);
     auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
     ASSERT_TRUE(src != nullptr);
-    auto dst = authFactory->createCoordinateReferenceSystem("4269"); // NAD83
+    auto dst =
+        authFactory->createCoordinateReferenceSystem("4269")->promoteTo3D(
+            std::string(), authFactory->databaseContext()); // NAD83
 
     auto list = CoordinateOperationFactory::create()->createOperations(
         NN_NO_CHECK(src), dst, ctxt);
@@ -6153,7 +6152,7 @@ TEST(operation, compoundCRS_with_non_meter_horiz_and_vertical_to_geog) {
         AuthorityFactory::create(DatabaseContext::create(), "EPSG");
     auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
     auto list = CoordinateOperationFactory::create()->createOperations(
-        NN_NO_CHECK(src), authFactory->createCoordinateReferenceSystem("4326"),
+        NN_NO_CHECK(src), authFactory->createCoordinateReferenceSystem("4979"),
         ctxt);
     ASSERT_EQ(list.size(), 1U);
     // Check that vertical unit conversion is done just once