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

Previously, when computing transformation between a compound CRS and a
geographic/projected 2D CRS, the behaviour was similar to implicitly
promoting the 2D CRS to 3D CRS in the pipeline computation logic, hence
a geoid model could be applied. But note that when doing a geographic 3D
to geographic/projected 2D CRS transformation, we *did* not do this implicit
promotion and if a Helmert transformation existed between the datums, it
was done only in 2D. So this is a bit inconsistent and triggered the
comment in https://github.com/OSGeo/PROJ/issues/2318#issuecomment-1068924513

With this commit, we no longer do any vertical transformation when doing
compound CRS to the 2D CRS, but just take the transformation of the
horizontal part of the compound CRS to the 2D CRS.
Said otherwise, NAD83+NAVD88 to NAD83 will no longer lead to the
application of the geoid model. Unless you explicitly ask for the
promotion NAD83 to 3D.

Also related, in https://github.com/OSGeo/PROJ/issues/1563 that went to 6.3.0,
I changed cs2cs to automatically promote to 3D the CRS as soon as one of
them was compound, for the sake of being consistent with the past
behaviour. But it then becomes difficult to predict PROJ behaviour
depending on which level of it you consider...
This commit undoes that and adds an explicit --3d switch to cs2cs, similarly to
projinfo, to ask for promotion.

Other bug fix found in the process, when using legacy syntax, +init=epsg:4979,
(WGS 84 3D), the resulting CRS was 2D and not 3D.

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