Add +proj=set operation to set component(s) of a coordinate to a fixed value

Fixes #1846

7 files changed, 173 insertions, 0 deletions

diff --git a/docs/source/operations/conversions/index.rst b/docs/source/operations/conversions/index.rst
index 42ad079d..b0bd0d37 100644
--- a/docs/source/operations/conversions/index.rst
+++ b/docs/source/operations/conversions/index.rst
@@ -18,4 +18,5 @@ conversions.
    noop
    pop
    push
+   set
    unitconvert
diff --git a/docs/source/operations/conversions/set.rst b/docs/source/operations/conversions/set.rst
new file mode 100644
index 00000000..635d08d3
--- /dev/null
+++ b/docs/source/operations/conversions/set.rst
@@ -0,0 +1,76 @@
+.. _set:
+
+================================================================================
+Set coordinate value
+================================================================================
+
+.. versionadded:: 7.0.0
+
+Set component(s) of a coordinate to a fixed value.
+
++---------------------+--------------------------------------------------------+
+| **Alias**           | set                                                    |
++---------------------+--------------------------------------------------------+
+| **Domain**          | 4D                                                     |
++---------------------+--------------------------------------------------------+
+| **Input type**      | Any                                                    |
++---------------------+--------------------------------------------------------+
+| **Output type**     | Any                                                    |
++---------------------+--------------------------------------------------------+
+
+This operations allows for components of coordinates to be set to a fixed value.
+This may be useful in :ref:`pipeline<pipeline>` when a step requires some
+component, typically an elevation or a date, to be set to a fixed value.
+
+Example
+################################################################################
+
+In the ETRS89 to Dutch RD with NAP height transformation, the used ellipsoidal
+height for the Helmert transformation is not the NAP height, but the height is
+set to 0 m. This is an unconventional trick to get the same results as when the
+effect of the Helmert transformation is included in the horizontal NTv2 grid.
+For the forward transformation from ETRS89 to RD with NAP height, we need to set
+the ellipsoidal ETRS89 height for the Helmert transformation to the equivalent
+of 0 m NAP. This is 43 m for the centre of the Netherlands and this value can
+be used as an approximation elsewhere (the effect of this approximation is
+below 1 mm for the horizontal coordinates, in an area up to hundreds of km
+outside the Netherlands).
+
+The ``+proj=set +v_3=0`` close to the end of the pipeline is to make it usable in
+the reverse direction.
+
+::
+
+   $ cct -t 0 -d 4 +proj=pipeline \
+        +step +proj=unitconvert +xy_in=deg +xy_out=rad \
+        +step +proj=axisswap +order=2,1 \
+        +step +proj=vgridshift +grids=nlgeo2018.gtx \
+        +step +proj=push +v_3 \
+        +step +proj=set +v_3=43 \
+        +step +proj=cart +ellps=GRS80 \
+        +step +proj=helmert +x=-565.7346 +y=-50.4058 +z=-465.2895 +rx=-0.395023 +ry=0.330776 +rz=-1.876073 +s=-4.07242 +convention=coordinate_frame +exact \
+        +step +proj=cart +inv +ellps=bessel \
+        +step +proj=hgridshift +inv +grids=rdcorr2018.gsb,null \
+        +step +proj=sterea +lat_0=52.156160556 +lon_0=5.387638889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel \
+        +step +proj=set +v_3=0 \
+        +step +proj=pop +v_3
+
+Parameters
+################################################################################
+
+.. option:: +v_1=value
+
+   Set the first coordinate component to the specified value
+
+.. option:: +v_2=value
+
+   Set the second coordinate component to the specified value
+
+.. option:: +v_3=value
+
+   Set the third coordinate component to the specified value
+
+.. option:: +v_4=value
+
+   Set the fourth coordinate component to the specified value
+
diff --git a/src/Makefile.am b/src/Makefile.am
index 8ab30a33..de35a754 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -174,6 +174,7 @@ libproj_la_SOURCES = \
 	conversions/geoc.cpp \
 	conversions/geocent.cpp \
 	conversions/noop.cpp \
+	conversions/set.cpp \
 	conversions/unitconvert.cpp \
 	\
 	transformations/affine.cpp \
diff --git a/src/conversions/set.cpp b/src/conversions/set.cpp
new file mode 100644
index 00000000..7628bf4f
--- /dev/null
+++ b/src/conversions/set.cpp
@@ -0,0 +1,75 @@
+#define PJ_LIB__
+
+#include "proj_internal.h"
+#include <errno.h>
+
+PROJ_HEAD(set, "Set coordinate value");
+
+/* Projection specific elements for the PJ object */
+namespace { // anonymous namespace
+struct Set {
+    bool v1;
+    bool v2;
+    bool v3;
+    bool v4;
+    double v1_val;
+    double v2_val;
+    double v3_val;
+    double v4_val;
+};
+} // anonymous namespace
+
+static PJ_COORD set_fwd_inv(PJ_COORD point, PJ *P) {
+
+    struct Set *set = static_cast<struct Set*>(P->opaque);
+
+    if (set->v1)
+        point.v[0] = set->v1_val;
+    if (set->v2)
+        point.v[1] = set->v2_val;
+    if (set->v3)
+        point.v[2] = set->v3_val;
+    if (set->v4)
+        point.v[3] = set->v4_val;
+
+    return point;
+}
+
+PJ *OPERATION(set, 0) {
+    P->inv4d = set_fwd_inv;
+    P->fwd4d = set_fwd_inv;
+
+    auto set = static_cast<struct Set*>(pj_calloc (1, sizeof(struct Set)));
+    P->opaque = set;
+    if (nullptr==P->opaque)
+        return pj_default_destructor(P, ENOMEM);
+
+    if (pj_param_exists(P->params, "v_1"))
+    {
+        set->v1 = true;
+        set->v1_val = pj_param (P->ctx, P->params, "dv_1").f;
+    }
+
+    if (pj_param_exists(P->params, "v_2"))
+    {
+        set->v2 = true;
+        set->v2_val = pj_param (P->ctx, P->params, "dv_2").f;
+    }
+
+    if (pj_param_exists(P->params, "v_3"))
+    {
+        set->v3 = true;
+        set->v3_val = pj_param (P->ctx, P->params, "dv_3").f;
+    }
+
+    if (pj_param_exists(P->params, "v_4"))
+    {
+        set->v4 = true;
+        set->v4_val = pj_param (P->ctx, P->params, "dv_4").f;
+    }
+
+    P->left  = PJ_IO_UNITS_WHATEVER;
+    P->right = PJ_IO_UNITS_WHATEVER;
+
+    return P;
+}
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 704ece3d..d89cfade 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -187,6 +187,7 @@ set(SRC_LIBPROJ_CONVERSIONS
   conversions/geoc.cpp
   conversions/geocent.cpp
   conversions/noop.cpp
+  conversions/set.cpp
   conversions/unitconvert.cpp
 )
 
diff --git a/src/pj_list.h b/src/pj_list.h
index 9798a36b..2e1bf759 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -134,6 +134,7 @@ PROJ_HEAD(robin, "Robinson")
 PROJ_HEAD(rouss, "Roussilhe Stereographic")
 PROJ_HEAD(rpoly, "Rectangular Polyconic")
 PROJ_HEAD(sch, "Spherical Cross-track Height")
+PROJ_HEAD(set, "Set coordinate value")
 PROJ_HEAD(sinu, "Sinusoidal (Sanson-Flamsteed)")
 PROJ_HEAD(somerc, "Swiss. Obl. Mercator")
 PROJ_HEAD(stere, "Stereographic")
diff --git a/test/gie/4D-API_cs2cs-style.gie b/test/gie/4D-API_cs2cs-style.gie
index 3e4b9d2c..b8512162 100644
--- a/test/gie/4D-API_cs2cs-style.gie
+++ b/test/gie/4D-API_cs2cs-style.gie
@@ -536,4 +536,22 @@ in the creation of a PJ object.
 operation   this is a bogus CRS meant to trigger a generic error in proj_create()
 expect      failure   errno generic error
 
+-------------------------------------------------------------------------------
+Test proj=set
+-------------------------------------------------------------------------------
+
+operation   +proj=set
+accept      1 2 3 4
+expect      1 2 3 4
+roundtrip   1
+
+operation   +proj=set +v_1=10 +v_2=20 +v_3=30 +v_4=40
+accept      1 2 3 4
+expect      10 20 30 40
+
+operation   +proj=set +v_1=10 +v_2=20 +v_3=30 +v_4=40
+direction   inverse
+accept      1 2 3 4
+expect      10 20 30 40
+
 </gie>