Merge pull request #3200 from rouault/vertoffset

Implement Vertical Offset and slope transformation method

docs/source/operations/transformations/index.rst

@@ -20,5 +20,6 @@ systems are based on different datums.
    molobadekas
    hgridshift
    tinshift
+   vertoffset
    vgridshift
    xyzgridshift

docs/source/operations/transformations/vertoffset.rst

@@ -0,0 +1,89 @@
+.. _vertoffset:
+
+================================================================================
+Vertical Offset And Slope
+================================================================================
+
+.. versionadded:: 9.1.0
+
+The Vertical Offset And Slope offsets transformation adds an offset to an
+orthometric height, taking account a constant offset and a inclinated plane,
+defined by its slope along latitude and longitude axis.
+
++---------------------+----------------------------------------------------------+
+| **Alias**           | vertoffset                                               |
++---------------------+----------------------------------------------------------+
+| **Domain**          | 3D                                                       |
++---------------------+----------------------------------------------------------+
+| **Input type**      | Geodetic coordinates (horizontal), meters (vertical)     |
++---------------------+----------------------------------------------------------+
+| **output type**     | Geodetic coordinates (horizontal), meters (vertical)     |
++---------------------+----------------------------------------------------------+
+
+It is documented as coordinate operation method code 1046 in the EPSG dataset (:cite:`IOGP2018`).
+It is typically used in Europe, to relate national vertical systems to
+pan-European vertical systems (EVRF200, EVRF2007).
+
+Examples
+###############################################################################
+
+Vertical offset from LN02 height to EVRF2000 height with horizontal coordinates in ETRS89::
+
+    +proj=vertoffset +lat_0=46.9166666666666666 +lon_0=8.183333333333334 \
+        +dh=-0.245 +slope_lat=-0.210 +slope_lon=-0.032 +ellps=GRS80
+
+Parameters
+################################################################################
+
+Required
+-------------------------------------------------------------------------------
+
+.. option:: +lat_0=<value>
+
+    Latitude of origin of the inclinated plane.
+
+.. option:: +lon_0=<value>
+
+    Longitude of origin of the inclinated plane
+
+.. option:: +dh=<value>
+
+    Offset in height, expressed in meter, to add.
+
+.. option:: +slope_lat=<value>
+
+    Slope parameter in the latitude domain, expressed in arc-seconds.
+
+.. option:: +slope_lon=<value>
+
+    Slope parameter in the longitude domain, expressed in arc-seconds.
+
+.. include:: ../options/ellps.rst
+
+Formula
+################################################################################
+
+The :math:`Z_{dest}` destination elevation is obtained from the
+:math:`Z_{src}` source elevation with:
+
+.. math::
+
+    \begin{align}
+        Z_{dest} = Z_{src} + \left( dh + slope_{lat} * {\rho}_0 * (\phi - {\phi}_0) + slope_{lon} * {\nu}_0 * (\lambda - {\lambda}_0) * cos(\phi) \right)
+    \end{align}
+
+where:
+
+* :math:`dh`, :math:`slope_{lat}` and :math:`slope_{lon}` are the above mentioned parameters
+* :math:`{\lambda}_0`, :math:`{\phi}_0` is the longitude, latitude of the point of origin of the inclinate plane (``+lon_0``, ``+lat_0``)
+* :math:`\lambda`, :math:`\phi` is the longitude, latitude of the point to evaluate
+* :math:`{\rho}_0` is the radius of curvature of the meridian at latitude :math:`{\phi}_0`
+* :math:`{\nu}_0` is the radius of curvature on the prime vertical (i.e. perpendicular to the meridian) at latitude :math:`{\phi}_0`
+
+The reverse formula is:
+
+.. math::
+
+    \begin{align}
+        Z_{src} = Z_{dest} - \left( dh + slope_{lat} * {\rho}_0 * (\phi - {\phi}_0) + slope_{lon} * {\nu}_0 * (\lambda - {\lambda}_0) * cos(\phi) \right)
+    \end{align}

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/lib_proj.cmake

@@ -175,6 +175,7 @@ set(SRC_LIBPROJ_TRANSFORMATIONS
   transformations/xyzgridshift.cpp
   transformations/defmodel.cpp
   transformations/tinshift.cpp
+  transformations/vertoffset.cpp
 )
 
 set(SRC_LIBPROJ_ISO19111

src/pj_list.h

@@ -170,6 +170,7 @@ PROJ_HEAD(vandg, "van der Grinten (I)")
 PROJ_HEAD(vandg2, "van der Grinten II")
 PROJ_HEAD(vandg3, "van der Grinten III")
 PROJ_HEAD(vandg4, "van der Grinten IV")
+PROJ_HEAD(vertoffset, "Vertical Offset and Slope")
 PROJ_HEAD(vitk1, "Vitkovsky I")
 PROJ_HEAD(vgridshift, "Vertical grid shift")
 PROJ_HEAD(wag1, "Wagner I (Kavrayskiy VI)")

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