Parse '+proj=something_not_latlong +vunits=' without +geoidgrids as a Projected3D CRS and not a compound CRS with a unknown datum (fixes #2287) (#2288)

3 files changed, 180 insertions, 50 deletions

diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 5493684d..4fc9b70f 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -7877,8 +7877,7 @@ struct PROJStringParser::Private {
     GeodeticReferenceFrameNNPtr buildDatum(Step &step,
                                            const std::string &title);
     GeographicCRSNNPtr buildGeographicCRS(int iStep, int iUnitConvert,
-                                          int iAxisSwap, bool ignoreVUnits,
-                                          bool ignorePROJAxis);
+                                          int iAxisSwap, bool ignorePROJAxis);
     GeodeticCRSNNPtr buildGeocentricCRS(int iStep, int iUnitConvert);
     CRSNNPtr buildProjectedCRS(int iStep, GeographicCRSNNPtr geogCRS,
                                int iUnitConvert, int iAxisSwap);
@@ -7893,8 +7892,7 @@ struct PROJStringParser::Private {
                     AxisType axisType, bool ignorePROJAxis);
 
     EllipsoidalCSNNPtr buildEllipsoidalCS(int iStep, int iUnitConvert,
-                                          int iAxisSwap, bool ignoreVUnits,
-                                          bool ignorePROJAxis);
+                                          int iAxisSwap, bool ignorePROJAxis);
 };
 
 // ---------------------------------------------------------------------------
@@ -8623,10 +8621,8 @@ PROJStringParser::Private::processAxisSwap(Step &step,
 
 // ---------------------------------------------------------------------------
 
-EllipsoidalCSNNPtr
-PROJStringParser::Private::buildEllipsoidalCS(int iStep, int iUnitConvert,
-                                              int iAxisSwap, bool ignoreVUnits,
-                                              bool ignorePROJAxis) {
+EllipsoidalCSNNPtr PROJStringParser::Private::buildEllipsoidalCS(
+    int iStep, int iUnitConvert, int iAxisSwap, bool ignorePROJAxis) {
     auto &step = steps_[iStep];
     assert(iUnitConvert < 0 ||
            ci_equal(steps_[iUnitConvert].name, "unitconvert"));
@@ -8661,7 +8657,7 @@ PROJStringParser::Private::buildEllipsoidalCS(int iStep, int iUnitConvert,
         AxisAbbreviation::h, AxisDirection::UP,
         buildUnit(step, "vunits", "vto_meter"));
 
-    return (!ignoreVUnits && !hasParamValue(step, "geoidgrids") &&
+    return (!hasParamValue(step, "geoidgrids") &&
             (hasParamValue(step, "vunits") || hasParamValue(step, "vto_meter")))
                ? EllipsoidalCS::create(emptyPropertyMap, axis[0], axis[1], up)
                : EllipsoidalCS::create(emptyPropertyMap, axis[0], axis[1]);
@@ -8688,10 +8684,8 @@ namespace {
 template <class T> inline void ignoreRetVal(T) {}
 }
 
-GeographicCRSNNPtr
-PROJStringParser::Private::buildGeographicCRS(int iStep, int iUnitConvert,
-                                              int iAxisSwap, bool ignoreVUnits,
-                                              bool ignorePROJAxis) {
+GeographicCRSNNPtr PROJStringParser::Private::buildGeographicCRS(
+    int iStep, int iUnitConvert, int iAxisSwap, bool ignorePROJAxis) {
     auto &step = steps_[iStep];
 
     const bool l_isGeographicStep = isGeographicStep(step.name);
@@ -8705,8 +8699,8 @@ PROJStringParser::Private::buildGeographicCRS(int iStep, int iUnitConvert,
 
     auto props = PropertyMap().set(IdentifiedObject::NAME_KEY,
                                    title.empty() ? "unknown" : title);
-    auto cs = buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, ignoreVUnits,
-                                 ignorePROJAxis);
+    auto cs =
+        buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, ignorePROJAxis);
 
     if (l_isGeographicStep &&
         (hasUnusedParameters(step) ||
@@ -9262,15 +9256,18 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
             return DerivedGeographicCRS::create(
                 PropertyMap().set(IdentifiedObject::NAME_KEY, "unnamed"),
                 geogCRS, NN_NO_CHECK(conv),
-                buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, false,
-                                   false));
+                buildEllipsoidalCS(iStep, iUnitConvert, iAxisSwap, false));
         }
     }
 
     std::vector<CoordinateSystemAxisNNPtr> axis =
         processAxisSwap(step, unit, iAxisSwap, axisType, false);
 
-    auto cs = CartesianCS::create(emptyPropertyMap, axis[0], axis[1]);
+    auto csGeogCRS = geogCRS->coordinateSystem();
+    auto cs = csGeogCRS->axisList().size() == 2
+                  ? CartesianCS::create(emptyPropertyMap, axis[0], axis[1])
+                  : CartesianCS::create(emptyPropertyMap, axis[0], axis[1],
+                                        csGeogCRS->axisList()[2]);
 
     auto props = PropertyMap().set(IdentifiedObject::NAME_KEY,
                                    title.empty() ? "unknown" : title);
@@ -9286,19 +9283,6 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
                   props.set(IdentifiedObject::NAME_KEY, webMercatorName), cs)
             : ProjectedCRS::create(props, geogCRS, NN_NO_CHECK(conv), cs);
 
-    if (!hasParamValue(step, "geoidgrids") &&
-        (hasParamValue(step, "vunits") || hasParamValue(step, "vto_meter"))) {
-        auto vdatum = VerticalReferenceFrame::create(mapWithUnknownName);
-
-        const UnitOfMeasure vunit = buildUnit(step, "vunits", "vto_meter");
-
-        auto vcrs =
-            VerticalCRS::create(mapWithUnknownName, vdatum,
-                                VerticalCS::createGravityRelatedHeight(vunit));
-
-        crs = CompoundCRS::create(mapWithUnknownName,
-                                  std::vector<CRSNNPtr>{crs, vcrs});
-    }
     return crs;
 }
 
@@ -9766,13 +9750,11 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
             iSecondGeogStep < 0 && iProjStep < 0 &&
             (iFirstUnitConvert < 0 || iSecondUnitConvert < 0) &&
             (iFirstAxisSwap < 0 || iSecondAxisSwap < 0)) {
-            const bool ignoreVUnits = false;
             // First run is dry run to mark all recognized/unrecognized tokens
             for (int iter = 0; iter < 2; iter++) {
                 auto obj = d->buildBoundOrCompoundCRSIfNeeded(
-                    0,
-                    d->buildGeographicCRS(iFirstGeogStep, iFirstUnitConvert,
-                                          iFirstAxisSwap, ignoreVUnits, false));
+                    0, d->buildGeographicCRS(iFirstGeogStep, iFirstUnitConvert,
+                                             iFirstAxisSwap, false));
                 if (iter == 1) {
                     return nn_static_pointer_cast<BaseObject>(obj);
                 }
@@ -9783,7 +9765,6 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
             iSecondGeogStep < 0) {
             if (iFirstGeogStep < 0)
                 iFirstGeogStep = iProjStep;
-            const bool ignoreVUnits = true;
             // First run is dry run to mark all recognized/unrecognized tokens
             for (int iter = 0; iter < 2; iter++) {
                 auto obj = d->buildBoundOrCompoundCRSIfNeeded(
@@ -9796,7 +9777,7 @@ PROJStringParser::createFromPROJString(const std::string &projString) {
                                                 : -1,
                             iFirstAxisSwap < iFirstGeogStep ? iFirstAxisSwap
                                                             : -1,
-                            ignoreVUnits, true),
+                            true),
                         iFirstUnitConvert < iFirstGeogStep ? iSecondUnitConvert
                                                            : iFirstUnitConvert,
                         iFirstAxisSwap < iFirstGeogStep ? iSecondAxisSwap
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index 6c92b57d..e6e56d44 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -9296,18 +9296,11 @@ TEST(io, projparse_projected_to_meter_unknown) {
 TEST(io, projparse_projected_vunits) {
     auto obj = PROJStringParser().createFromPROJString(
         "+proj=tmerc +vunits=ft +type=crs");
-    auto crs = nn_dynamic_pointer_cast<CompoundCRS>(obj);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
     ASSERT_TRUE(crs != nullptr);
-    WKTFormatterNNPtr f(WKTFormatter::create());
-    f->simulCurNodeHasId();
-    f->setMultiLine(false);
-    crs->exportToWKT(f.get());
-    auto wkt = f->toString();
-    EXPECT_TRUE(wkt.find("CS[Cartesian,2]") != std::string::npos) << wkt;
-    EXPECT_TRUE(wkt.find("CS[vertical,1],AXIS[\"gravity-related height "
-                         "(H)\",up,LENGTHUNIT[\"foot\",0.3048]") !=
-                std::string::npos)
-        << wkt;
+    auto cs = crs->coordinateSystem();
+    ASSERT_EQ(cs->axisList().size(), 3U);
+    EXPECT_EQ(cs->axisList()[2]->unit().name(), "foot");
 }
 
 // ---------------------------------------------------------------------------
diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp
index 9ebb3794..6a39879a 100644
--- a/test/unit/test_operation.cpp
+++ b/test/unit/test_operation.cpp
@@ -8406,10 +8406,97 @@ TEST(operation, vertCRS_to_vertCRS_New_Zealand_context) {
 
 // ---------------------------------------------------------------------------
 
-TEST(operation, compoundCRS_to_geogCRS_3D) {
+TEST(operation, projCRS_3D_to_geogCRS_3D) {
 
     auto compoundcrs_ft_obj = PROJStringParser().createFromPROJString(
         "+proj=merc +vunits=ft +type=crs");
+    auto proj3DCRS_ft = nn_dynamic_pointer_cast<CRS>(compoundcrs_ft_obj);
+    ASSERT_TRUE(proj3DCRS_ft != nullptr);
+
+    auto geogcrs_m_obj = PROJStringParser().createFromPROJString(
+        "+proj=longlat +vunits=m +type=crs");
+    auto geogcrs_m = nn_dynamic_pointer_cast<CRS>(geogcrs_m_obj);
+    ASSERT_TRUE(geogcrs_m != nullptr);
+
+    {
+        auto op = CoordinateOperationFactory::create()->createOperation(
+            NN_CHECK_ASSERT(proj3DCRS_ft), NN_CHECK_ASSERT(geogcrs_m));
+        ASSERT_TRUE(op != nullptr);
+        EXPECT_FALSE(op->hasBallparkTransformation());
+        EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+                  "+proj=pipeline "
+                  "+step +proj=unitconvert +xy_in=m +z_in=ft "
+                  "+xy_out=m +z_out=m "
+                  "+step +inv +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 "
+                  "+ellps=WGS84 "
+                  "+step +proj=unitconvert +xy_in=rad +z_in=m "
+                  "+xy_out=deg +z_out=m");
+    }
+
+    {
+        auto op = CoordinateOperationFactory::create()->createOperation(
+            NN_CHECK_ASSERT(geogcrs_m), NN_CHECK_ASSERT(proj3DCRS_ft));
+        ASSERT_TRUE(op != nullptr);
+        EXPECT_FALSE(op->hasBallparkTransformation());
+        EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+                  "+proj=pipeline "
+                  "+step +proj=unitconvert +z_in=m +z_out=ft "
+                  "+step +proj=unitconvert +xy_in=deg +z_in=ft "
+                  "+xy_out=rad +z_out=m "
+                  "+step +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 "
+                  "+step +proj=unitconvert +xy_in=m +z_in=m "
+                  "+xy_out=m +z_out=ft");
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, compoundCRS_to_geogCRS_3D) {
+
+    auto compoundcrs_ft_obj = WKTParser().createFromWKT(
+        "COMPOUNDCRS[\"unknown\",\n"
+        "    PROJCRS[\"unknown\",\n"
+        "        BASEGEOGCRS[\"unknown\",\n"
+        "            DATUM[\"World Geodetic System 1984\",\n"
+        "                ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
+        "                    LENGTHUNIT[\"metre\",1]],\n"
+        "                ID[\"EPSG\",6326]],\n"
+        "            PRIMEM[\"Greenwich\",0,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8901]]],\n"
+        "        CONVERSION[\"unknown\",\n"
+        "            METHOD[\"Mercator (variant A)\",\n"
+        "                ID[\"EPSG\",9804]],\n"
+        "            PARAMETER[\"Latitude of natural origin\",0,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8801]],\n"
+        "            PARAMETER[\"Longitude of natural origin\",0,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8802]],\n"
+        "            PARAMETER[\"Scale factor at natural origin\",1,\n"
+        "                SCALEUNIT[\"unity\",1],\n"
+        "                ID[\"EPSG\",8805]],\n"
+        "            PARAMETER[\"False easting\",0,\n"
+        "                LENGTHUNIT[\"metre\",1],\n"
+        "                ID[\"EPSG\",8806]],\n"
+        "            PARAMETER[\"False northing\",0,\n"
+        "                LENGTHUNIT[\"metre\",1],\n"
+        "                ID[\"EPSG\",8807]]],\n"
+        "        CS[Cartesian,2],\n"
+        "            AXIS[\"(E)\",east,\n"
+        "                ORDER[1],\n"
+        "                LENGTHUNIT[\"metre\",1,\n"
+        "                    ID[\"EPSG\",9001]]],\n"
+        "            AXIS[\"(N)\",north,\n"
+        "                ORDER[2],\n"
+        "                LENGTHUNIT[\"metre\",1,\n"
+        "                    ID[\"EPSG\",9001]]]],\n"
+        "    VERTCRS[\"unknown\",\n"
+        "        VDATUM[\"unknown\"],\n"
+        "        CS[vertical,1],\n"
+        "            AXIS[\"gravity-related height (H)\",up,\n"
+        "                LENGTHUNIT[\"foot\",0.3048,\n"
+        "                    ID[\"EPSG\",9002]]]]]");
     auto compoundcrs_ft = nn_dynamic_pointer_cast<CRS>(compoundcrs_ft_obj);
     ASSERT_TRUE(compoundcrs_ft != nullptr);
 
@@ -9307,7 +9394,7 @@ TEST(operation, compoundCRS_from_WKT2_no_id_to_geogCRS_3D_context) {
 
 // ---------------------------------------------------------------------------
 
-TEST(operation, compoundCRS_with_non_meter_horiz_and_vertical_to_geog) {
+TEST(operation, proj3DCRS_with_non_meter_horiz_and_vertical_to_geog) {
     auto objSrc = PROJStringParser().createFromPROJString(
         "+proj=utm +zone=31 +datum=WGS84 +units=us-ft +vunits=us-ft +type=crs");
     auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
@@ -9322,6 +9409,75 @@ TEST(operation, compoundCRS_with_non_meter_horiz_and_vertical_to_geog) {
     // Check that vertical unit conversion is done just once
     EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=pipeline "
+              "+step +proj=unitconvert +xy_in=us-ft +z_in=us-ft "
+              "+xy_out=m +z_out=m "
+              "+step +inv +proj=utm +zone=31 +ellps=WGS84 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, compoundCRS_with_non_meter_horiz_and_vertical_to_geog) {
+    auto objSrc = WKTParser().createFromWKT(
+        "COMPOUNDCRS[\"unknown\",\n"
+        "    PROJCRS[\"unknown\",\n"
+        "        BASEGEOGCRS[\"unknown\",\n"
+        "            DATUM[\"World Geodetic System 1984\",\n"
+        "                ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
+        "                    LENGTHUNIT[\"metre\",1]],\n"
+        "                ID[\"EPSG\",6326]],\n"
+        "            PRIMEM[\"Greenwich\",0,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8901]]],\n"
+        "        CONVERSION[\"UTM zone 31N\",\n"
+        "            METHOD[\"Transverse Mercator\",\n"
+        "                ID[\"EPSG\",9807]],\n"
+        "            PARAMETER[\"Latitude of natural origin\",0,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8801]],\n"
+        "            PARAMETER[\"Longitude of natural origin\",3,\n"
+        "                ANGLEUNIT[\"degree\",0.0174532925199433],\n"
+        "                ID[\"EPSG\",8802]],\n"
+        "            PARAMETER[\"Scale factor at natural origin\",0.9996,\n"
+        "                SCALEUNIT[\"unity\",1],\n"
+        "                ID[\"EPSG\",8805]],\n"
+        "            PARAMETER[\"False easting\",500000,\n"
+        "                LENGTHUNIT[\"metre\",1],\n"
+        "                ID[\"EPSG\",8806]],\n"
+        "            PARAMETER[\"False northing\",0,\n"
+        "                LENGTHUNIT[\"metre\",1],\n"
+        "                ID[\"EPSG\",8807]],\n"
+        "            ID[\"EPSG\",16031]],\n"
+        "        CS[Cartesian,2],\n"
+        "            AXIS[\"(E)\",east,\n"
+        "                ORDER[1],\n"
+        "                LENGTHUNIT[\"US survey foot\",0.304800609601219,\n"
+        "                    ID[\"EPSG\",9003]]],\n"
+        "            AXIS[\"(N)\",north,\n"
+        "                ORDER[2],\n"
+        "                LENGTHUNIT[\"US survey foot\",0.304800609601219,\n"
+        "                    ID[\"EPSG\",9003]]]],\n"
+        "    VERTCRS[\"unknown\",\n"
+        "        VDATUM[\"unknown\"],\n"
+        "        CS[vertical,1],\n"
+        "            AXIS[\"gravity-related height (H)\",up,\n"
+        "                LENGTHUNIT[\"US survey foot\",0.304800609601219,\n"
+        "                    ID[\"EPSG\",9003]]]]]"
+
+        );
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+    auto authFactory =
+        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"),
+        ctxt);
+    ASSERT_EQ(list.size(), 1U);
+    // Check that vertical unit conversion is done just once
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
               "+step +proj=unitconvert +xy_in=us-ft +xy_out=m "
               "+step +inv +proj=utm +zone=31 +ellps=WGS84 "
               "+step +proj=unitconvert +xy_in=rad +z_in=us-ft "