Merge pull request #3578 from rouault/similarity_transformation

Implement EPSG:9621 'Similarity transformation' and import related tranformations from EPSG

docs/source/operations/transformations/affine.rst

@@ -22,6 +22,32 @@ x,y,z coordinates, and translation and scaling on the temporal coordinate.
 By default, the parameters are set for an identity transforms. The transformation
 is reversible unless the determinant of the sji matrix is 0, or `tscale` is 0
 
+This can be used to implement:
+
+- EPSG "Affine parametric transformation" of code 9624 by setting `A0`, `A1`, `A2`,
+  `B0`, `B1`, `B2` EPSG parameters to respectively `xoff`, `s11`, `s12`, `yoff`,
+  `s21`, `s21` PROJ parameters.
+
+- EPSG "Similarity transformation" of code 9621 by setting:
+
+    * xoff to :math:`X_{T0}`
+    * yoff to :math:`Y_{T0}`
+    * s11 to :math:`M \cos \theta`
+    * s12 to :math:`M \sin \theta`
+    * s21 to :math:`-M \sin \theta`
+    * s22 to :math:`M \cos \theta`
+
+  where:
+
+    * :math:`X_{T0}` is the first ordinate of the origin point of the source
+      CRS expressed in the target CRS.
+    * :math:`Y_{T0}` is the second ordinate of the origin point of the source
+      CRS expressed in the target CRS.
+    * :math:`M` is the multiplication factor applied to coordinates in the
+      source CRS to obtain the correct scale of the target CRS.
+    * :math:`\theta` is the angle about which the axes of the source CRS need to
+      be rotated to coincide with the axes of the target CRS, counter-clockwise
+      being positive
 
 Parameters
 ################################################################################

scripts/build_db.py

@@ -820,7 +820,8 @@ def fill_other_transformation(proj_db_cursor):
     # 1068: Height Depth Reversal
     # 1069: Change of Vertical Unit
     # 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)")
+    # 9621 : Similarity transformation
+    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, 9621)")
     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
@@ -843,6 +844,11 @@ def fill_other_transformation(proj_db_cursor):
             source_crs_code = source_codes[0][0]
             target_crs_code = target_codes[0][0]
 
+        # Engineering CRS
+        if source_crs_code in (5800, 5817):
+            print("Skipping transformation %s as source CRS (%d) is not handled" % (name, source_crs_code))
+            continue
+
         expected_order = 1
         max_n_params = 7
         param_auth_name = [None for i in range(max_n_params)]

src/iso19111/operation/conversion.cpp

@@ -3580,6 +3580,8 @@ void Conversion::_exportToPROJString(
         methodEPSGCode == EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT_NO_CONV_FACTOR;
     const bool isAffineParametric =
         methodEPSGCode == EPSG_CODE_METHOD_AFFINE_PARAMETRIC_TRANSFORMATION;
+    const bool isSimilarity =
+        methodEPSGCode == EPSG_CODE_METHOD_SIMILARITY_TRANSFORMATION;
     const bool isGeographicGeocentric =
         methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC_GEOCENTRIC;
     const bool isGeographicOffsets =
@@ -3589,7 +3591,7 @@ void Conversion::_exportToPROJString(
     const bool isHeightDepthReversal =
         methodEPSGCode == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL;
     const bool applySourceCRSModifiers =
-        !isZUnitConversion && !isAffineParametric &&
+        !isZUnitConversion && !isAffineParametric && !isSimilarity &&
         !isAxisOrderReversal(methodEPSGCode) && !isGeographicGeocentric &&
         !isGeographicOffsets && !isHeightDepthReversal;
     bool applyTargetCRSModifiers = applySourceCRSModifiers;

src/iso19111/operation/parammappings.cpp

@@ -945,6 +945,7 @@ const struct MethodNameCode methodNameCodes[] = {
     // Transformations
     METHOD_NAME_CODE(LONGITUDE_ROTATION),
     METHOD_NAME_CODE(AFFINE_PARAMETRIC_TRANSFORMATION),
+    METHOD_NAME_CODE(SIMILARITY_TRANSFORMATION),
     METHOD_NAME_CODE(COORDINATE_FRAME_GEOCENTRIC),
     METHOD_NAME_CODE(COORDINATE_FRAME_GEOGRAPHIC_2D),
     METHOD_NAME_CODE(COORDINATE_FRAME_GEOGRAPHIC_3D),
@@ -1112,6 +1113,31 @@ static const ParamMapping paramB2 = {
 static const ParamMapping *const paramsAffineParametricTransformation[] = {
     &paramA0, &paramA1, &paramA2, &paramB0, &paramB1, &paramB2, nullptr};
 
+static const ParamMapping paramOrdinate1EvalPointTargetCRS = {
+    EPSG_NAME_PARAMETER_ORDINATE_1_EVAL_POINT_TARGET_CRS,
+    EPSG_CODE_PARAMETER_ORDINATE_1_EVAL_POINT_TARGET_CRS, nullptr,
+    common::UnitOfMeasure::Type::UNKNOWN, nullptr};
+
+static const ParamMapping paramOrdinate2EvalPointTargetCRS = {
+    EPSG_NAME_PARAMETER_ORDINATE_2_EVAL_POINT_TARGET_CRS,
+    EPSG_CODE_PARAMETER_ORDINATE_2_EVAL_POINT_TARGET_CRS, nullptr,
+    common::UnitOfMeasure::Type::UNKNOWN, nullptr};
+
+static const ParamMapping paramScaleFactorForSourceCRSAxes = {
+    EPSG_NAME_PARAMETER_SCALE_FACTOR_FOR_SOURCE_CRS_AXES,
+    EPSG_CODE_PARAMETER_SCALE_FACTOR_FOR_SOURCE_CRS_AXES, nullptr,
+    common::UnitOfMeasure::Type::SCALE, nullptr};
+
+static const ParamMapping paramRotationAngleOfSourceCRSAxes = {
+    EPSG_NAME_PARAMETER_ROTATION_ANGLE_OF_SOURCE_CRS_AXES,
+    EPSG_CODE_PARAMETER_ROTATION_ANGLE_OF_SOURCE_CRS_AXES, nullptr,
+    common::UnitOfMeasure::Type::ANGULAR, nullptr};
+
+static const ParamMapping *const paramsSimilarityTransformation[] = {
+    &paramOrdinate1EvalPointTargetCRS, &paramOrdinate2EvalPointTargetCRS,
+    &paramScaleFactorForSourceCRSAxes, &paramRotationAngleOfSourceCRSAxes,
+    nullptr};
+
 static const ParamMapping paramXTranslation = {
     EPSG_NAME_PARAMETER_X_AXIS_TRANSLATION,
     EPSG_CODE_PARAMETER_X_AXIS_TRANSLATION, nullptr,
@@ -1393,6 +1419,10 @@ static const MethodMapping otherMethodMappings[] = {
      EPSG_CODE_METHOD_AFFINE_PARAMETRIC_TRANSFORMATION, nullptr, nullptr,
      nullptr, paramsAffineParametricTransformation},
 
+    {EPSG_NAME_METHOD_SIMILARITY_TRANSFORMATION,
+     EPSG_CODE_METHOD_SIMILARITY_TRANSFORMATION, nullptr, nullptr, nullptr,
+     paramsSimilarityTransformation},
+
     {PROJ_WKT2_NAME_METHOD_POLE_ROTATION_GRIB_CONVENTION, 0, nullptr, nullptr,
      nullptr, paramsPoleRotationGRIBConvention},
 

src/iso19111/operation/singleoperation.cpp

@@ -2883,6 +2883,36 @@ bool SingleOperation::exportToPROJStringGeneric(
         return true;
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_SIMILARITY_TRANSFORMATION) {
+        const double XT0 =
+            parameterValueMeasure(
+                EPSG_CODE_PARAMETER_ORDINATE_1_EVAL_POINT_TARGET_CRS)
+                .value();
+        const double YT0 =
+            parameterValueMeasure(
+                EPSG_CODE_PARAMETER_ORDINATE_2_EVAL_POINT_TARGET_CRS)
+                .value();
+        const double M =
+            parameterValueMeasure(
+                EPSG_CODE_PARAMETER_SCALE_FACTOR_FOR_SOURCE_CRS_AXES)
+                .value();
+        const double q = parameterValueNumeric(
+            EPSG_CODE_PARAMETER_ROTATION_ANGLE_OF_SOURCE_CRS_AXES,
+            common::UnitOfMeasure::RADIAN);
+
+        // Do not mess with axis unit and order for that transformation
+
+        formatter->addStep("affine");
+        formatter->addParam("xoff", XT0);
+        formatter->addParam("s11", M * cos(q));
+        formatter->addParam("s12", M * sin(q));
+        formatter->addParam("yoff", YT0);
+        formatter->addParam("s21", -M * sin(q));
+        formatter->addParam("s22", M * cos(q));
+
+        return true;
+    }
+
     if (isAxisOrderReversal(methodEPSGCode)) {
         formatter->addStep("axisswap");
         formatter->addParam("order", "2,1");

src/proj_constants.h

@@ -705,6 +705,28 @@
 
 /* ------------------------------------------------------------------------ */
 
+#define EPSG_CODE_METHOD_SIMILARITY_TRANSFORMATION 9621
+#define EPSG_NAME_METHOD_SIMILARITY_TRANSFORMATION \
+    "Similarity transformation"
+
+#define EPSG_NAME_PARAMETER_ORDINATE_1_EVAL_POINT_TARGET_CRS \
+    "Ordinate 1 of evaluation point in target CRS"
+#define EPSG_CODE_PARAMETER_ORDINATE_1_EVAL_POINT_TARGET_CRS 8621
+
+#define EPSG_NAME_PARAMETER_ORDINATE_2_EVAL_POINT_TARGET_CRS \
+    "Ordinate 2 of evaluation point in target CRS"
+#define EPSG_CODE_PARAMETER_ORDINATE_2_EVAL_POINT_TARGET_CRS 8622
+
+#define EPSG_NAME_PARAMETER_SCALE_FACTOR_FOR_SOURCE_CRS_AXES \
+    "Scale factor for source CRS axes"
+#define EPSG_CODE_PARAMETER_SCALE_FACTOR_FOR_SOURCE_CRS_AXES 1061
+
+#define EPSG_NAME_PARAMETER_ROTATION_ANGLE_OF_SOURCE_CRS_AXES \
+    "Rotation angle of source CRS axes"
+#define EPSG_CODE_PARAMETER_ROTATION_ANGLE_OF_SOURCE_CRS_AXES 8614
+
+/* ------------------------------------------------------------------------ */
+
 #define EPSG_CODE_METHOD_AXIS_ORDER_REVERSAL_2D 9843
 #define EPSG_NAME_METHOD_AXIS_ORDER_REVERSAL_2D "Axis Order Reversal (2D)"
 

test/cli/testvarious

@@ -1156,6 +1156,20 @@ $EXE --only-best --no-ballpark EPSG:4326+3855 EPSG:4979 -E >>${OUT} 2>&1 <<EOF
 49 2 0
 EOF
 
+echo "##############################################################" >> ${OUT}
+echo "Test Similarity Transformation (example from EPSG Guidance Note 7.2)" >> ${OUT}
+#
+$EXE -d 3 EPSG:23031 EPSG:25831 -E >>${OUT} 2>&1 <<EOF
+300000 4500000
+EOF
+
+echo "##############################################################" >> ${OUT}
+echo "Test inverse of Similarity Transformation (example from EPSG Guidance Note 7.2)" >> ${OUT}
+#
+$EXE -d 3 EPSG:25831 EPSG:23031 -E >>${OUT} 2>&1 <<EOF
+299905.060 4499796.515
+EOF
+
 
 # Done!
 # do 'diff' with distribution results

test/cli/tv_out.dist

@@ -568,3 +568,9 @@ Rel. 9.2.0, March 1st, 2023
 cannot initialize transformation
 cause: File not found or invalid
 program abnormally terminated
+##############################################################
+Test Similarity Transformation (example from EPSG Guidance Note 7.2)
+300000 4500000	299905.060	4499796.515 0.000
+##############################################################
+Test inverse of Similarity Transformation (example from EPSG Guidance Note 7.2)
+299905.060 4499796.515	300000.000	4500000.000 0.000