Merge 8bc7c00 into 81bfa75

OSGeo · Jul 10, 2018 · ad0707e · ad0707e
2 parents 81bfa75 + 8bc7c00
commit ad0707e
Show file tree

Hide file tree

Showing 3 changed files with 32 additions and 4 deletions.
diff --git 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
@@ -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,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
@@ -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=1.0
+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
 
 -------------------------------------------------------------------------------
 Some tests from PJ_hgridshift.c