vgridshift: add a +multiplier={value} parameter

As mentionned in #1071, it is often unclear how the offset of a vertical grid
is applied.

3 files changed, 32 insertions, 4 deletions

diff --git a/docs/source/operations/transformations/vgridshift.rst b/docs/source/operations/transformations/vgridshift.rst
index e671779b..484ad4bf 100644
--- a/docs/source/operations/transformations/vgridshift.rst
+++ b/docs/source/operations/transformations/vgridshift.rst
@@ -113,3 +113,19 @@ Optional
 .. include:: ../options/t_final.rst
 .. versionadded:: 5.1.0
 
+.. option:: +multiplier=<value>
+
+    Specify the multiplier to apply to the grid value in the forward transformation
+    direction, such that:
+
+    .. math::
+        :label: multiplier_formula
+
+        Z_{target} = Z_{source} + multiplier \times gridvalue
+
+    The multiplier can be used to control whether the gridvalue should be added
+    or sustracted, and if unit conversion must be done (the multiplied gridvalue
+    must be expressed in metre).
+
+    Note that the default is `-1.0` for historical reasons.
+.. versionadded:: 5.2.0
diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c
index 1cd0880e..f88527f7 100644
--- a/src/PJ_vgridshift.c
+++ b/src/PJ_vgridshift.c
@@ -13,16 +13,18 @@ PROJ_HEAD(vgridshift, "Vertical grid shift");
 struct pj_opaque_vgridshift {
     double t_final;
     double t_epoch;
+    double forward_multiplier;
 };
 
 static XYZ forward_3d(LPZ lpz, PJ *P) {
+    struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque;
     PJ_COORD point = {{0,0,0,0}};
     point.lpz = lpz;
 
     if (P->vgridlist_geoid != NULL) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.xyz.z -= proj_vgrid_value(P, point.lp);
+        point.xyz.z += Q->forward_multiplier * proj_vgrid_value(P, point.lp);
     }
 
     return point.xyz;
@@ -30,13 +32,14 @@ static XYZ forward_3d(LPZ lpz, PJ *P) {
 
 
 static LPZ reverse_3d(XYZ xyz, PJ *P) {
+    struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) P->opaque;
     PJ_COORD point = {{0,0,0,0}};
     point.xyz = xyz;
 
     if (P->vgridlist_geoid != NULL) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.xyz.z += proj_vgrid_value(P, point.lp);
+        point.xyz.z -= Q->forward_multiplier * proj_vgrid_value(P, point.lp);
     }
 
     return point.lpz;
@@ -110,6 +113,11 @@ PJ *TRANSFORMATION(vgridshift,0) {
    if (pj_param(P->ctx, P->params, "tt_epoch").i)
         Q->t_epoch = pj_param (P->ctx, P->params, "dt_epoch").f;
 
+    /* historical: the forward direction substracts the grid offset. */
+    Q->forward_multiplier = -1.0;
+    if (pj_param(P->ctx, P->params, "tmultiplier").i) {
+        Q->forward_multiplier = pj_param(P->ctx, P->params, "dmultiplier").f;
+    }
 
     /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */
     proj_vgrid_init(P, "grids");
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index 7e205e0a..3a6c1650 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -205,8 +205,12 @@ operation proj=vgridshift grids=nonexistinggrid.gtx
 expect    failure errno failed_to_load_grid
 -------------------------------------------------------------------------------
 
-
-
+-------------------------------------------------------------------------------
+operation  proj=vgridshift  grids=egm96_15.gtx  ellps=GRS80     multiplier=0.1
+tolerance 15 cm
+ignore    pjd_err_failed_to_load_grid
+accept    12.5 55.5   0                  0
+expect    12.5 55.5  3.6021305084228516  0
 
 -------------------------------------------------------------------------------
 Some tests from PJ_hgridshift.c