Map EPSG:1046 operation method to PROJ +proj=vertoffset (fixes #3194)

scripts/build_db.py

@@ -669,7 +669,8 @@ def fill_other_transformation(proj_db_cursor):
     # 9660: Geographic3D offsets
     # 1068: Height Depth Reversal
     # 1069: Change of Vertical Unit
-    proj_db_cursor.execute("SELECT coord_op_code, coord_op_name, coord_op_method_code, coord_op_method_name, source_crs_code, target_crs_code, coord_op_accuracy, coord_tfm_version, epsg_coordoperation.deprecated, epsg_coordoperation.remarks FROM epsg.epsg_coordoperation LEFT JOIN epsg.epsg_coordoperationmethod USING (coord_op_method_code) WHERE coord_op_method_code IN (9601, 9616, 9618, 9619, 9624, 9660, 1068, 1069)")
+    # 1046: Vertical Offset and Slope
+    proj_db_cursor.execute("SELECT coord_op_code, coord_op_name, coord_op_method_code, coord_op_method_name, source_crs_code, target_crs_code, coord_op_accuracy, coord_tfm_version, epsg_coordoperation.deprecated, epsg_coordoperation.remarks FROM epsg.epsg_coordoperation LEFT JOIN epsg.epsg_coordoperationmethod USING (coord_op_method_code) WHERE coord_op_method_code IN (9601, 9616, 9618, 9619, 9624, 9660, 1068, 1069, 1046)")
     for (code, name, method_code, method_name, source_crs_code, target_crs_code, coord_op_accuracy, coord_tfm_version, deprecated, remarks) in proj_db_cursor.fetchall():
 
         # 1068 and 1069 are Height Depth Reversal and Change of Vertical Unit
@@ -700,11 +701,19 @@ def fill_other_transformation(proj_db_cursor):
         param_value = [None for i in range(max_n_params)]
         param_uom_auth_name = [None for i in range(max_n_params)]
         param_uom_code = [None for i in range(max_n_params)]
+        interpolation_crs_auth_name = None
+        interpolation_crs_code = None
 
         iterator = proj_db_cursor.execute("SELECT sort_order, cop.parameter_code, parameter_name, parameter_value, uom_code from epsg_coordoperationparam cop LEFT JOIN epsg_coordoperationparamvalue copv LEFT JOIN epsg_coordoperationparamusage copu ON cop.parameter_code = copv.parameter_code AND copu.parameter_code = copv.parameter_code WHERE copu.coord_op_method_code = copv.coord_op_method_code AND coord_op_code = ? AND copv.coord_op_method_code = ? ORDER BY sort_order", (code, method_code))
         for (order, parameter_code, parameter_name, parameter_value, uom_code) in iterator:
             assert order <= max_n_params
             assert order == expected_order
+            if method_code == 1046 and order == 6: # Vertical offset and slope
+                assert parameter_code == 1037 # EPSG code for Horizontal CRS
+                interpolation_crs_auth_name = EPSG_AUTHORITY
+                interpolation_crs_code = str(int(parameter_value)) # needed to avoid codes like XXXX.0
+                break
+
             param_auth_name[order - 1] = EPSG_AUTHORITY
             param_code[order - 1] = parameter_code
             param_name[order - 1] = parameter_name
@@ -713,6 +722,7 @@ def fill_other_transformation(proj_db_cursor):
             param_uom_code[order - 1] = uom_code
             expected_order += 1
 
+
         arg = (EPSG_AUTHORITY, code, name,
                remarks, 
                EPSG_AUTHORITY, method_code, method_name,
@@ -733,7 +743,7 @@ def fill_other_transformation(proj_db_cursor):
                param_uom_auth_name[5], param_uom_code[5], param_auth_name[6],
                param_code[6], param_name[6], param_value[6],
                param_uom_auth_name[6], param_uom_code[6],
-               None, None, # interpolation CRS
+               interpolation_crs_auth_name, interpolation_crs_code,
                coord_tfm_version,
                deprecated)
 

src/iso19111/operation/parammappings.cpp

@@ -966,6 +966,7 @@ const struct MethodNameCode methodNameCodes[] = {
     METHOD_NAME_CODE(GEOGRAPHIC2D_WITH_HEIGHT_OFFSETS),
     METHOD_NAME_CODE(GEOGRAPHIC3D_OFFSETS),
     METHOD_NAME_CODE(VERTICAL_OFFSET),
+    METHOD_NAME_CODE(VERTICAL_OFFSET_AND_SLOPE),
     METHOD_NAME_CODE(NTV2),
     METHOD_NAME_CODE(NTV1),
     METHOD_NAME_CODE(NADCON),
@@ -1050,6 +1051,8 @@ const struct ParamNameCode paramNameCodes[] = {
     PARAM_NAME_CODE(ORDINATE_2_EVAL_POINT),
     PARAM_NAME_CODE(ORDINATE_3_EVAL_POINT),
     PARAM_NAME_CODE(GEOCENTRIC_TRANSLATION_FILE),
+    PARAM_NAME_CODE(INCLINATION_IN_LATITUDE),
+    PARAM_NAME_CODE(INCLINATION_IN_LONGITUDE),
 };
 
 const ParamNameCode *getParamNameCodes(size_t &nElts) {
@@ -1256,6 +1259,21 @@ static const ParamMapping *const paramsGeographic3DOffsets[] = {
 static const ParamMapping *const paramsVerticalOffsets[] = {
     &paramVerticalOffset, nullptr};
 
+static const ParamMapping paramInclinationInLatitude = {
+    EPSG_NAME_PARAMETER_INCLINATION_IN_LATITUDE,
+    EPSG_CODE_PARAMETER_INCLINATION_IN_LATITUDE, nullptr,
+    common::UnitOfMeasure::Type::ANGULAR, nullptr};
+
+static const ParamMapping paramInclinationInLongitude = {
+    EPSG_NAME_PARAMETER_INCLINATION_IN_LONGITUDE,
+    EPSG_CODE_PARAMETER_INCLINATION_IN_LONGITUDE, nullptr,
+    common::UnitOfMeasure::Type::ANGULAR, nullptr};
+
+static const ParamMapping *const paramsVerticalOffsetAndSlope[] = {
+    &paramOrdinate1EvalPoint,     &paramOrdinate2EvalPoint,
+    &paramVerticalOffset,         &paramInclinationInLatitude,
+    &paramInclinationInLongitude, nullptr};
+
 static const ParamMapping paramLatitudeLongitudeDifferenceFile = {
     EPSG_NAME_PARAMETER_LATITUDE_LONGITUDE_DIFFERENCE_FILE,
     EPSG_CODE_PARAMETER_LATITUDE_LONGITUDE_DIFFERENCE_FILE, nullptr,
@@ -1450,6 +1468,10 @@ static const MethodMapping otherMethodMappings[] = {
     {EPSG_NAME_METHOD_VERTICAL_OFFSET, EPSG_CODE_METHOD_VERTICAL_OFFSET,
      nullptr, nullptr, nullptr, paramsVerticalOffsets},
 
+    {EPSG_NAME_METHOD_VERTICAL_OFFSET_AND_SLOPE,
+     EPSG_CODE_METHOD_VERTICAL_OFFSET_AND_SLOPE, nullptr, nullptr, nullptr,
+     paramsVerticalOffsetAndSlope},
+
     {EPSG_NAME_METHOD_NTV2, EPSG_CODE_METHOD_NTV2, nullptr, nullptr, nullptr,
      paramsNTV2},
 

src/iso19111/operation/transformation.cpp

@@ -2991,6 +2991,66 @@ void Transformation::_exportToPROJString(
         return;
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_VERTICAL_OFFSET_AND_SLOPE) {
+
+        const crs::CRS *srcCRS = sourceCRS().get();
+        const crs::CRS *tgtCRS = targetCRS().get();
+
+        const auto sourceCRSCompound =
+            dynamic_cast<const crs::CompoundCRS *>(srcCRS);
+        const auto targetCRSCompound =
+            dynamic_cast<const crs::CompoundCRS *>(tgtCRS);
+        if (sourceCRSCompound && targetCRSCompound &&
+            sourceCRSCompound->componentReferenceSystems()[0]->_isEquivalentTo(
+                targetCRSCompound->componentReferenceSystems()[0].get(),
+                util::IComparable::Criterion::EQUIVALENT)) {
+            srcCRS = sourceCRSCompound->componentReferenceSystems()[1].get();
+            tgtCRS = targetCRSCompound->componentReferenceSystems()[1].get();
+        }
+
+        auto sourceCRSVert = dynamic_cast<const crs::VerticalCRS *>(srcCRS);
+        if (!sourceCRSVert) {
+            throw io::FormattingException(
+                "Can apply Vertical offset and slope only to VerticalCRS");
+        }
+
+        auto targetCRSVert = dynamic_cast<const crs::VerticalCRS *>(tgtCRS);
+        if (!targetCRSVert) {
+            throw io::FormattingException(
+                "Can apply Vertical offset and slope only to VerticalCRS");
+        }
+
+        const auto latitudeEvaluationPoint =
+            parameterValueNumeric(EPSG_CODE_PARAMETER_ORDINATE_1_EVAL_POINT,
+                                  common::UnitOfMeasure::DEGREE);
+        const auto longitudeEvaluationPoint =
+            parameterValueNumeric(EPSG_CODE_PARAMETER_ORDINATE_2_EVAL_POINT,
+                                  common::UnitOfMeasure::DEGREE);
+        const auto offsetHeight =
+            parameterValueNumericAsSI(EPSG_CODE_PARAMETER_VERTICAL_OFFSET);
+        const auto inclinationLatitude =
+            parameterValueNumeric(EPSG_CODE_PARAMETER_INCLINATION_IN_LATITUDE,
+                                  common::UnitOfMeasure::ARC_SECOND);
+        const auto inclinationLongitude =
+            parameterValueNumeric(EPSG_CODE_PARAMETER_INCLINATION_IN_LONGITUDE,
+                                  common::UnitOfMeasure::ARC_SECOND);
+
+        formatter->startInversion();
+        sourceCRSVert->addLinearUnitConvert(formatter);
+        formatter->stopInversion();
+
+        formatter->addStep("vertoffset");
+        formatter->addParam("lat_0", latitudeEvaluationPoint);
+        formatter->addParam("lon_0", longitudeEvaluationPoint);
+        formatter->addParam("dh", offsetHeight);
+        formatter->addParam("slope_lat", inclinationLatitude);
+        formatter->addParam("slope_lon", inclinationLongitude);
+
+        targetCRSVert->addLinearUnitConvert(formatter);
+
+        return;
+    }
+
     // Substitute grid names with PROJ friendly names.
     if (formatter->databaseContext()) {
         auto alternate = substitutePROJAlternativeGridNames(

src/proj_constants.h

@@ -620,6 +620,9 @@
 #define EPSG_CODE_METHOD_VERTICAL_OFFSET 9616
 #define EPSG_NAME_METHOD_VERTICAL_OFFSET "Vertical Offset"
 
+#define EPSG_CODE_METHOD_VERTICAL_OFFSET_AND_SLOPE 1046
+#define EPSG_NAME_METHOD_VERTICAL_OFFSET_AND_SLOPE "Vertical Offset and Slope"
+
 #define EPSG_CODE_METHOD_GEOGRAPHIC2D_OFFSETS 9619
 #define EPSG_NAME_METHOD_GEOGRAPHIC2D_OFFSETS "Geographic2D offsets"
 
@@ -646,6 +649,12 @@
 #define EPSG_NAME_PARAMETER_GEOID_UNDULATION "Geoid undulation"
 #define EPSG_CODE_PARAMETER_GEOID_UNDULATION 8604
 
+#define EPSG_NAME_PARAMETER_INCLINATION_IN_LATITUDE "Inclination in latitude"
+#define EPSG_CODE_PARAMETER_INCLINATION_IN_LATITUDE 8730
+
+#define EPSG_NAME_PARAMETER_INCLINATION_IN_LONGITUDE "Inclination in longitude"
+#define EPSG_CODE_PARAMETER_INCLINATION_IN_LONGITUDE 8731
+
 /* ------------------------------------------------------------------------ */
 
 #define EPSG_CODE_METHOD_AFFINE_PARAMETRIC_TRANSFORMATION 9624

test/unit/test_operationfactory.cpp

@@ -4789,22 +4789,37 @@ TEST(
             authFactory->createCoordinateReferenceSystem("8360"),
             // ETRS89 + EVRF2007 height
             authFactory->createCoordinateReferenceSystem("7423"), ctxt);
-        ASSERT_GE(list.size(), 1U);
+        ASSERT_GE(list.size(), 2U);
+
+        // For Czechia
         EXPECT_EQ(
             list[0]->exportToPROJString(PROJStringFormatter::create().get()),
             "+proj=pipeline "
             "+step +proj=axisswap +order=2,1 "
             "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+            "+step +proj=vertoffset +lat_0=49.9166666666667 "
+            "+lon_0=15.25 +dh=0.13 +slope_lat=0.026 +slope_lon=0 "
+            "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+            "+step +proj=axisswap +order=2,1");
+        EXPECT_EQ(list[0]->nameStr(),
+                  "Baltic 1957 height to EVRF2007 height (1)");
+
+        // For Slovakia
+        EXPECT_EQ(
+            list[1]->exportToPROJString(PROJStringFormatter::create().get()),
+            "+proj=pipeline "
+            "+step +proj=axisswap +order=2,1 "
+            "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
             "+step +proj=vgridshift "
             "+grids=sk_gku_Slovakia_ETRS89h_to_Baltic1957.tif +multiplier=1 "
             "+step +inv +proj=vgridshift "
             "+grids=sk_gku_Slovakia_ETRS89h_to_EVRF2007.tif +multiplier=1 "
             "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
             "+step +proj=axisswap +order=2,1");
         EXPECT_EQ(
-            list[0]->nameStr(),
+            list[1]->nameStr(),
             "ETRS89 + Baltic 1957 height to ETRS89 + EVRF2007 height (1)");
-        EXPECT_EQ(list[0]->inverse()->nameStr(), "Inverse of 'ETRS89 + Baltic "
+        EXPECT_EQ(list[1]->inverse()->nameStr(), "Inverse of 'ETRS89 + Baltic "
                                                  "1957 height to ETRS89 + "
                                                  "EVRF2007 height (1)'");
     }