Vertical unit conversion: improve output for ft <--> us-ft

OSGeo · Nov 2, 2023 · f6765ed · f6765ed
1 parent d38713d
commit f6765ed
Show file tree

Hide file tree

Showing 4 changed files with 58 additions and 48 deletions.
diff --git a/include/proj/coordinateoperation.hpp b/include/proj/coordinateoperation.hpp
@@ -727,6 +727,10 @@ class PROJ_GCC_DLL SingleOperation : virtual public CoordinateOperation {
     createOperationParameterValueFromInterpolationCRS(int methodEPSGCode,
                                                       int crsEPSGCode);
 
+    PROJ_INTERNAL static void
+    exportToPROJStringChangeVerticalUnit(io::PROJStringFormatter *formatter,
+                                         double convFactor);
+
   private:
     PROJ_OPAQUE_PRIVATE_DATA
     SingleOperation &operator=(const SingleOperation &other) = delete;

diff --git a/src/iso19111/operation/conversion.cpp b/src/iso19111/operation/conversion.cpp
@@ -4144,29 +4144,7 @@ void Conversion::_exportToPROJString(
                     "requires an input and output vertical CRS");
             }
         }
-        auto uom = common::UnitOfMeasure(std::string(), convFactor,
-                                         common::UnitOfMeasure::Type::LINEAR)
-                       .exportToPROJString();
-        auto reverse_uom =
-            convFactor == 0.0
-                ? std::string()
-                : common::UnitOfMeasure(std::string(), 1.0 / convFactor,
-                                        common::UnitOfMeasure::Type::LINEAR)
-                      .exportToPROJString();
-        if (uom == "m") {
-            // do nothing
-        } else if (!uom.empty()) {
-            formatter->addStep("unitconvert");
-            formatter->addParam("z_in", uom);
-            formatter->addParam("z_out", "m");
-        } else if (!reverse_uom.empty()) {
-            formatter->addStep("unitconvert");
-            formatter->addParam("z_in", "m");
-            formatter->addParam("z_out", reverse_uom);
-        } else {
-            formatter->addStep("affine");
-            formatter->addParam("s33", convFactor);
-        }
+        exportToPROJStringChangeVerticalUnit(formatter, convFactor);
         bConversionDone = true;
         bEllipsoidParametersDone = true;
     } else if (methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC_TOPOCENTRIC) {

diff --git a/src/iso19111/operation/singleoperation.cpp b/src/iso19111/operation/singleoperation.cpp
@@ -3144,6 +3144,50 @@ static void setupPROJGeodeticTargetCRS(io::PROJStringFormatter *formatter,
 
 // ---------------------------------------------------------------------------
 
+/* static */
+void SingleOperation::exportToPROJStringChangeVerticalUnit(
+    io::PROJStringFormatter *formatter, double convFactor) {
+
+    const auto uom = common::UnitOfMeasure(std::string(), convFactor,
+                                           common::UnitOfMeasure::Type::LINEAR)
+                         .exportToPROJString();
+    const auto reverse_uom =
+        convFactor == 0.0
+            ? std::string()
+            : common::UnitOfMeasure(std::string(), 1.0 / convFactor,
+                                    common::UnitOfMeasure::Type::LINEAR)
+                  .exportToPROJString();
+    if (uom == "m") {
+        // do nothing
+    } else if (!uom.empty()) {
+        formatter->addStep("unitconvert");
+        formatter->addParam("z_in", uom);
+        formatter->addParam("z_out", "m");
+    } else if (!reverse_uom.empty()) {
+        formatter->addStep("unitconvert");
+        formatter->addParam("z_in", "m");
+        formatter->addParam("z_out", reverse_uom);
+    } else if (fabs(convFactor -
+                    common::UnitOfMeasure::FOOT.conversionToSI() /
+                        common::UnitOfMeasure::US_FOOT.conversionToSI()) <
+               1e-10) {
+        formatter->addStep("unitconvert");
+        formatter->addParam("z_in", "ft");
+        formatter->addParam("z_out", "us-ft");
+    } else if (fabs(convFactor -
+                    common::UnitOfMeasure::US_FOOT.conversionToSI() /
+                        common::UnitOfMeasure::FOOT.conversionToSI()) < 1e-10) {
+        formatter->addStep("unitconvert");
+        formatter->addParam("z_in", "us-ft");
+        formatter->addParam("z_out", "ft");
+    } else {
+        formatter->addStep("affine");
+        formatter->addParam("s33", convFactor);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
 bool SingleOperation::exportToPROJStringGeneric(
     io::PROJStringFormatter *formatter) const {
     const int methodEPSGCode = method()->getEPSGCode();
@@ -3267,30 +3311,7 @@ bool SingleOperation::exportToPROJStringGeneric(
     if (methodEPSGCode == EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT) {
         const double convFactor = parameterValueNumericAsSI(
             EPSG_CODE_PARAMETER_UNIT_CONVERSION_SCALAR);
-        const auto uom =
-            common::UnitOfMeasure(std::string(), convFactor,
-                                  common::UnitOfMeasure::Type::LINEAR)
-                .exportToPROJString();
-        const auto reverse_uom =
-            convFactor == 0.0
-                ? std::string()
-                : common::UnitOfMeasure(std::string(), 1.0 / convFactor,
-                                        common::UnitOfMeasure::Type::LINEAR)
-                      .exportToPROJString();
-        if (uom == "m") {
-            // do nothing
-        } else if (!uom.empty()) {
-            formatter->addStep("unitconvert");
-            formatter->addParam("z_in", uom);
-            formatter->addParam("z_out", "m");
-        } else if (!reverse_uom.empty()) {
-            formatter->addStep("unitconvert");
-            formatter->addParam("z_in", "m");
-            formatter->addParam("z_out", reverse_uom);
-        } else {
-            formatter->addStep("affine");
-            formatter->addParam("s33", convFactor);
-        }
+        exportToPROJStringChangeVerticalUnit(formatter, convFactor);
         return true;
     }
 

diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
@@ -6169,7 +6169,14 @@ TEST(operation, vertCRS_to_vertCRS) {
             NN_CHECK_ASSERT(vertcrs_ft), NN_CHECK_ASSERT(vertcrs_us_ft));
         ASSERT_TRUE(op != nullptr);
         EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
-                  "+proj=affine +s33=0.999998");
+                  "+proj=unitconvert +z_in=ft +z_out=us-ft");
+    }
+    {
+        auto op = CoordinateOperationFactory::create()->createOperation(
+            NN_CHECK_ASSERT(vertcrs_us_ft), NN_CHECK_ASSERT(vertcrs_ft));
+        ASSERT_TRUE(op != nullptr);
+        EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
+                  "+proj=unitconvert +z_in=us-ft +z_out=ft");
     }
 
     auto vertCRSMetreUp =