vgridshift: add a +operation=add or +operation=substract parameter

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

docs/source/operations/transformations/vgridshift.rst

@@ -113,3 +113,11 @@ Optional
 .. include:: ../options/t_final.rst
 .. versionadded:: 5.1.0
 
+.. option:: +operation=add or +operation=substract
+
+    Specify the operation to apply to the grid value in the forward transformation
+    direction. If `add`, the grid value will be added to the input z coordinate.
+    If `substract`, the grid value will be substracted to the input z coordinate.
+    Note that the default is `substract` for historical reasons.
+
+.. versionadded:: 5.2.0

src/PJ_vgridshift.c

@@ -13,30 +13,33 @@ PROJ_HEAD(vgridshift, "Vertical grid shift");
 struct pj_opaque_vgridshift {
     double t_final;
     double t_epoch;
+    int 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;
 }
 
 
 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,19 @@ 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 foward direction substracts the grid offset. */
+   Q->forward_multiplier = -1.0;
+   if (pj_param(P->ctx, P->params, "toperation").i) {
+        const char* operation = pj_param(P->ctx, P->params, "soperation").s;
+        if (!strcmp("add", operation)) {
+            Q->forward_multiplier = 1.0;
+        } else if (!strcmp("substract", operation)) {
+            Q->forward_multiplier = -1.0;
+        } else {
+            proj_log_error(P, "vgridshift: +operation parameter supported values are 'add or 'substract'.");
+            return pj_default_destructor(P, PJD_ERR_INVALID_ARG);
+        }
+    }
 
     /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */
     proj_vgrid_init(P, "grids");

test/gie/more_builtins.gie

@@ -205,7 +205,24 @@ operation proj=vgridshift grids=nonexistinggrid.gtx
 expect    failure errno failed_to_load_grid
 -------------------------------------------------------------------------------
 
+-------------------------------------------------------------------------------
+operation  proj=vgridshift  grids=egm96_15.gtx  ellps=GRS80     operation=substract
+tolerance 15 cm
+ignore    pjd_err_failed_to_load_grid
+accept    12.5 55.5   0                  0
+expect    12.5 55.5 -36.021305084228516  0
+
+-------------------------------------------------------------------------------
+operation  proj=vgridshift  grids=egm96_15.gtx  ellps=GRS80     operation=add
+tolerance 15 cm
+ignore    pjd_err_failed_to_load_grid
+accept    12.5 55.5   0                  0
+expect    12.5 55.5  36.021305084228516  0
 
+-------------------------------------------------------------------------------
+operation  proj=vgridshift  grids=egm96_15.gtx  ellps=GRS80     operation=invalid
+ignore    pjd_err_failed_to_load_grid
+expect    failure pjd_err_invalid_arg
 
 
 -------------------------------------------------------------------------------