Implement EPSG:9656 'Cartesian Grid Offsets' operation method, and im…

…port related records

include/proj/coordinateoperation.hpp

@@ -1683,6 +1683,12 @@ class PROJ_GCC_DLL Transformation : public SingleOperation {
         const common::Angle &offsetLong, const common::Length &offsetHeight,
         const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies);
 
+    PROJ_DLL static TransformationNNPtr createCartesianGridOffsets(
+        const util::PropertyMap &properties, const crs::CRSNNPtr &sourceCRSIn,
+        const crs::CRSNNPtr &targetCRSIn, const common::Length &eastingOffset,
+        const common::Length &northingOffset,
+        const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies);
+
     PROJ_DLL static TransformationNNPtr createVerticalOffset(
         const util::PropertyMap &properties, const crs::CRSNNPtr &sourceCRSIn,
         const crs::CRSNNPtr &targetCRSIn, const common::Length &offsetHeight,

scripts/build_db.py

@@ -829,8 +829,9 @@ def fill_other_transformation(proj_db_cursor):
     # 1068: Height Depth Reversal
     # 1069: Change of Vertical Unit
     # 1046: Vertical Offset and Slope
-    # 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)")
+    # 9621: Similarity transformation
+    # 9656: Cartesian Grid Offsets
+    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, 9656)")
     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
@@ -854,7 +855,7 @@ def fill_other_transformation(proj_db_cursor):
             target_crs_code = target_codes[0][0]
 
         # Engineering CRS
-        if source_crs_code in (5800, 5817):
+        if source_crs_code in (5800, 5817, 6715):
             print("Skipping transformation %s as source CRS (%d) is not handled" % (name, source_crs_code))
             continue
 

src/iso19111/operation/parammappings.cpp

@@ -1011,6 +1011,7 @@ const struct MethodNameCode methodNameCodesList[] = {
     METHOD_NAME_CODE(GEOGRAPHIC2D_OFFSETS),
     METHOD_NAME_CODE(GEOGRAPHIC2D_WITH_HEIGHT_OFFSETS),
     METHOD_NAME_CODE(GEOGRAPHIC3D_OFFSETS),
+    METHOD_NAME_CODE(CARTESIAN_GRID_OFFSETS),
     METHOD_NAME_CODE(VERTICAL_OFFSET),
     METHOD_NAME_CODE(VERTICAL_OFFSET_AND_SLOPE),
     METHOD_NAME_CODE(NTV2),
@@ -1336,6 +1337,17 @@ static const ParamMapping paramVerticalOffset = {
 static const ParamMapping *const paramsGeographic3DOffsets[] = {
     &paramLatitudeOffset, &paramLongitudeOffset, &paramVerticalOffset, nullptr};
 
+static const ParamMapping paramEastingOffset = {
+    EPSG_NAME_PARAMETER_EASTING_OFFSET, EPSG_CODE_PARAMETER_EASTING_OFFSET,
+    nullptr, common::UnitOfMeasure::Type::LINEAR, nullptr};
+
+static const ParamMapping paramNorthingOffset = {
+    EPSG_NAME_PARAMETER_NORTHING_OFFSET, EPSG_CODE_PARAMETER_NORTHING_OFFSET,
+    nullptr, common::UnitOfMeasure::Type::LINEAR, nullptr};
+
+static const ParamMapping *const paramsCartesianGridOffsets[] = {
+    &paramEastingOffset, &paramNorthingOffset, nullptr};
+
 static const ParamMapping *const paramsVerticalOffsets[] = {
     &paramVerticalOffset, nullptr};
 
@@ -1569,6 +1581,10 @@ static const MethodMapping otherMethodMappings[] = {
      EPSG_CODE_METHOD_GEOGRAPHIC3D_OFFSETS, nullptr, nullptr, nullptr,
      paramsGeographic3DOffsets},
 
+    {EPSG_NAME_METHOD_CARTESIAN_GRID_OFFSETS,
+     EPSG_CODE_METHOD_CARTESIAN_GRID_OFFSETS, nullptr, nullptr, nullptr,
+     paramsCartesianGridOffsets},
+
     {EPSG_NAME_METHOD_VERTICAL_OFFSET, EPSG_CODE_METHOD_VERTICAL_OFFSET,
      nullptr, nullptr, nullptr, paramsVerticalOffsets},
 

src/iso19111/operation/singleoperation.cpp

@@ -3829,6 +3829,43 @@ bool SingleOperation::exportToPROJStringGeneric(
         return true;
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_CARTESIAN_GRID_OFFSETS) {
+        double eastingOffset = parameterValueNumeric(
+            EPSG_CODE_PARAMETER_EASTING_OFFSET, common::UnitOfMeasure::METRE);
+        double northingOffset = parameterValueNumeric(
+            EPSG_CODE_PARAMETER_NORTHING_OFFSET, common::UnitOfMeasure::METRE);
+
+        auto l_sourceCRS = sourceCRS();
+        auto sourceCRSProj =
+            dynamic_cast<const crs::ProjectedCRS *>(l_sourceCRS.get());
+        if (!sourceCRSProj) {
+            throw io::FormattingException(
+                "Can apply Cartesian grid offsets only to ProjectedCRS");
+        }
+
+        auto l_targetCRS = targetCRS();
+        auto targetCRSProj =
+            dynamic_cast<const crs::ProjectedCRS *>(l_targetCRS.get());
+        if (!targetCRSProj) {
+            throw io::FormattingException(
+                "Can apply Cartesian grid offsets only to ProjectedCRS");
+        }
+
+        formatter->startInversion();
+        sourceCRSProj->addUnitConvertAndAxisSwap(formatter, false);
+        formatter->stopInversion();
+
+        if (eastingOffset != 0.0 || northingOffset != 0.0) {
+            formatter->addStep("affine");
+            formatter->addParam("xoff", eastingOffset);
+            formatter->addParam("yoff", northingOffset);
+        }
+
+        targetCRSProj->addUnitConvertAndAxisSwap(formatter, false);
+
+        return true;
+    }
+
     if (methodEPSGCode == EPSG_CODE_METHOD_VERTICAL_OFFSET) {
 
         const crs::CRS *srcCRS = sourceCRS().get();

src/iso19111/operation/transformation.cpp

@@ -1273,6 +1273,38 @@ TransformationNNPtr Transformation::createGeographic2DWithHeightOffsets(
 
 // ---------------------------------------------------------------------------
 
+/** \brief Instantiate a transformation with method Cartesian grid offsets
+ *
+ * This method is defined as
+ * <a href="https://epsg.org/coord-operation-method_9656/index.html">
+ * EPSG:9656</a>.
+ *
+ * @param properties See \ref general_properties of the Transformation.
+ * At minimum the name should be defined.
+ * @param sourceCRSIn Source CRS.
+ * @param targetCRSIn Target CRS.
+ * @param eastingOffset Easting offset to add.
+ * @param northingOffset Northing offset to add.
+ * @param accuracies Vector of positional accuracy (might be empty).
+ * @return new Transformation.
+ * @since PROJ 9.5.0
+ */
+TransformationNNPtr Transformation::createCartesianGridOffsets(
+    const util::PropertyMap &properties, const crs::CRSNNPtr &sourceCRSIn,
+    const crs::CRSNNPtr &targetCRSIn, const common::Length &eastingOffset,
+    const common::Length &northingOffset,
+    const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies) {
+    return create(
+        properties, sourceCRSIn, targetCRSIn, nullptr,
+        createMethodMapNameEPSGCode(EPSG_CODE_METHOD_CARTESIAN_GRID_OFFSETS),
+        VectorOfParameters{
+            createOpParamNameEPSGCode(EPSG_CODE_PARAMETER_EASTING_OFFSET),
+            createOpParamNameEPSGCode(EPSG_CODE_PARAMETER_NORTHING_OFFSET)},
+        VectorOfValues{eastingOffset, northingOffset}, accuracies);
+}
+
+// ---------------------------------------------------------------------------
+
 /** \brief Instantiate a transformation with method Vertical Offset.
  *
  * This method is defined as
@@ -1650,6 +1682,23 @@ TransformationNNPtr Transformation::inverseAsTransformation() const {
                       newOffsetHeight, coordinateOperationAccuracies()));
     }
 
+    if (methodEPSGCode == EPSG_CODE_METHOD_CARTESIAN_GRID_OFFSETS) {
+        const auto &eastingOffset =
+            parameterValueMeasure(EPSG_CODE_PARAMETER_EASTING_OFFSET);
+        const common::Length newEastingOffset(negate(eastingOffset.value()),
+                                              eastingOffset.unit());
+
+        const auto &northingOffset =
+            parameterValueMeasure(EPSG_CODE_PARAMETER_NORTHING_OFFSET);
+        const common::Length newNorthingOffset(negate(northingOffset.value()),
+                                               northingOffset.unit());
+        return Private::registerInv(
+            this, createCartesianGridOffsets(
+                      createPropertiesForInverse(this, false, false),
+                      l_targetCRS, l_sourceCRS, newEastingOffset,
+                      newNorthingOffset, coordinateOperationAccuracies()));
+    }
+
     if (methodEPSGCode == EPSG_CODE_METHOD_VERTICAL_OFFSET) {
 
         const auto &offsetHeight =

src/proj_constants.h

@@ -819,6 +819,17 @@
 
 /* ------------------------------------------------------------------------ */
 
+#define EPSG_CODE_METHOD_CARTESIAN_GRID_OFFSETS 9656
+#define EPSG_NAME_METHOD_CARTESIAN_GRID_OFFSETS "Cartesian Grid Offsets"
+
+#define EPSG_CODE_PARAMETER_EASTING_OFFSET 8728
+#define EPSG_NAME_PARAMETER_EASTING_OFFSET "Easting offset"
+
+#define EPSG_CODE_PARAMETER_NORTHING_OFFSET 8729
+#define EPSG_NAME_PARAMETER_NORTHING_OFFSET "Northing offset"
+
+/* ------------------------------------------------------------------------ */
+
 #define EPSG_NAME_METHOD_GEOCENTRIC_TOPOCENTRIC                                \
     "Geocentric/topocentric conversions"
 #define EPSG_CODE_METHOD_GEOCENTRIC_TOPOCENTRIC 9836

test/unit/test_factory.cpp

@@ -182,8 +182,6 @@ TEST(factory, AuthorityFactory_identifyBodyFromSemiMajorAxis) {
 TEST(factory, AuthorityFactory_createEllipsoid) {
     auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
     EXPECT_THROW(factory->createEllipsoid("-1"), NoSuchAuthorityCodeException);
-    EXPECT_TRUE(nn_dynamic_pointer_cast<Ellipsoid>(
-                    factory->createObject("7030")) != nullptr);
     auto ellipsoid = factory->createEllipsoid("7030");
     ASSERT_EQ(ellipsoid->identifiers().size(), 1U);
     EXPECT_EQ(ellipsoid->identifiers()[0]->code(), "7030");

test/unit/test_operationfactory.cpp

@@ -2578,6 +2578,66 @@ TEST(operation, projCRS_to_projCRS_context_incompatible_areas_ballpark) {
 
 // ---------------------------------------------------------------------------
 
+TEST(operation, projCRS_to_projCRS_context_grid_offsets) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    {
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            authFactory->createCoordinateReferenceSystem(
+                "3392"), // Karbala 1979 / UTM zone 38N
+            authFactory->createCoordinateReferenceSystem(
+                "3891"), // IGRS / UTM zone 38N
+            ctxt);
+        ASSERT_EQ(list.size(), 2U);
+        EXPECT_EQ(list[1]->nameStr(),
+                  "Karbala 1979 / UTM zone 38N to IGRS / UTM zone 38N (1)");
+        EXPECT_EQ(
+            list[1]->exportToPROJString(PROJStringFormatter::create().get()),
+            "+proj=affine +xoff=-287.54 +yoff=278.25");
+    }
+    {
+        auto list = CoordinateOperationFactory::create()->createOperations(
+            authFactory->createCoordinateReferenceSystem(
+                "3891"), // IGRS / UTM zone 38N
+            authFactory->createCoordinateReferenceSystem(
+                "3392"), // Karbala 1979 / UTM zone 38N
+            ctxt);
+        ASSERT_EQ(list.size(), 2U);
+        EXPECT_EQ(list[1]->nameStr(), "Inverse of Karbala 1979 / UTM zone 38N "
+                                      "to IGRS / UTM zone 38N (1)");
+        EXPECT_EQ(
+            list[1]->exportToPROJString(PROJStringFormatter::create().get()),
+            "+proj=affine +xoff=287.54 +yoff=-278.25");
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, projCRS_to_projCRS_context_grid_offsets_non_metre_unit_noop) {
+    auto authFactory =
+        AuthorityFactory::create(DatabaseContext::create(), "EPSG");
+    auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        authFactory->createCoordinateReferenceSystem(
+            "10516"), // NAD83(2011) / Adjusted Jackson (ftUS)
+        authFactory->createCoordinateReferenceSystem(
+            "8162"), // NAD83(HARN) / WISCRS Jackson (ftUS)
+        ctxt);
+    ASSERT_EQ(list.size(), 2U);
+    EXPECT_EQ(list[0]->nameStr(), "NAD83(2011) / Adjusted Jackson (ftUS) to "
+                                  "NAD83(HARN) / WISCRS Jackson (ftUS) (1)");
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=noop");
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(
     operation,
     projCRS_to_projCRS_context_incompatible_areas_crs_extent_use_intersection) {