Reverse direction of deformation operations/transformations/deformation

Changed the direction of dt-calculation to follow the same convention as
helmert. Changed from dt = t_c - t_obs to dt = t_obs - t_c, which
effectively reverses the direction of the operation. Existing projstrings
using deformation can simply reverse the direction of the operation to
achieve the same results as before this commit.

4 files changed, 40 insertions, 16 deletions

diff --git a/docs/source/operations/transformations/deformation.rst b/docs/source/operations/transformations/deformation.rst
index 87b08453..0a772547 100644
--- a/docs/source/operations/transformations/deformation.rst
+++ b/docs/source/operations/transformations/deformation.rst
@@ -127,7 +127,7 @@ Mathematically speaking, application of a deformation model is simple. The defor
 represented as a grid of velocities in three dimensions. Coordinate corrections are
 applied in cartesian space. For a given coordinate, :math:`(X, Y, Z)`, velocities
 :math:`(V_X, V_Y, V_Z)` can be interpolated from the gridded model. The time span
-between :math:`t_c` and :math:`t_{obs}` determine the magnitude of the coordinate
+between :math:`t_{obs}` and :math:`t_c` determine the magnitude of the coordinate
 correcton as seen in eq. :eq:`apply_velocity` below.
 
 .. math::
@@ -144,7 +144,7 @@ correcton as seen in eq. :eq:`apply_velocity` below.
             Y \\
             Z \\
         \end{pmatrix}_A +
-        (t_c - t_{obs})
+        (t_{obs} - t_c)
         \begin{pmatrix}
             V_X \\
             V_Y \\
@@ -154,14 +154,14 @@ correcton as seen in eq. :eq:`apply_velocity` below.
 
 Corrections are done in cartesian space.
 
-Coordinates of the gridded model are in ENU (east, north, up) space because it would
-otherwise require an enormous 3 dimensional grid to handle the corrections in cartesian
-space. Keeping the correction in lat/long space reduces the complexity of the grid
-significantly. Consequently though, the input coordinates needs to be converted to
-lat/long space when searching for corrections in the grid. This is done with the *cart*
-operation. The converted grid corrections can then be applied to the input coordinates
-in cartesian space. The conversion from ENU space to cartesian space is done in the
-following way:
+Coordinates of the gridded model are in ENU (east, north, up) space because it
+would otherwise require an enormous 3 dimensional grid to handle the corrections
+in cartesian space. Keeping the correction in lat/long space reduces the
+complexity of the grid significantly. Consequently though, the input coordinates
+needs to be converted to lat/long space when searching for corrections in the
+grid. This is done with the :ref:`cart<cart>` operation. The converted grid
+corrections can then be applied to the input coordinates in cartesian space. The
+conversion from ENU space to cartesian space is done in the following way:
 
 .. math::
     :label: enu2xyz
diff --git a/docs/source/usage/differences.rst b/docs/source/usage/differences.rst
index 5694e4d0..caab865c 100644
--- a/docs/source/usage/differences.rst
+++ b/docs/source/usage/differences.rst
@@ -69,3 +69,25 @@ Before PROJ 6, the proj_def.dat was used to provide general and per-projection
 parameters, when +no_defs was not specified. It has now been removed. In case,
 no ellipsoid or datum specification is provided in the PROJ string, the
 default ellipsoid is GRS80 (was WGS84 in previous PROJ versions).
+
+Changes to :ref:`deformation<deformation>`
+------------------------------------------------------------------
+
+In the initial version of the of :ref:`deformation<deformation>` operation the time span between :math:`t_{obs}` and :math:`t_c` was determine by the expression
+
+.. math::
+
+    dt = t_c - t_{obs}
+
+With version 6.0.0 this has been reversed in order to behave similarly to
+the :ref:`Helmert operation<helmert>`, which determines time differences as
+
+.. math::
+
+    dt = t_{obs} - t_c
+
+Effectively this means that the direction of the operation has been reversed,
+so that what in PROJ 5 was a forward operation is now an inverse operation and
+vice versa.
+
+Pipelines written for PROJ 5 can be migrated to PROJ 6 by adding :option:`+inv` to forward steps involving the deformation operation. Similarly :option:`+inv` should be removed from inverse steps to be compatible with PROJ 6.
diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp
index 9221d39e..14f34e4e 100644
--- a/src/transformations/deformation.cpp
+++ b/src/transformations/deformation.cpp
@@ -172,7 +172,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
     out = in;
 
     if (Q->t_obs != HUGE_VAL) {
-        dt = Q->t_epoch - Q->t_obs;
+        dt = Q->t_obs - Q->t_epoch;
     } else {
         out = proj_coord_error(); /* in the 3D case +t_obs must be specified */
         proj_log_debug(P, "deformation: +t_obs must be specified");
@@ -196,9 +196,9 @@ static PJ_COORD forward_4d(PJ_COORD in, PJ *P) {
     PJ_COORD out = in;
 
     if (Q->t_obs != HUGE_VAL) {
-            dt = Q->t_epoch - Q->t_obs;
+            dt = Q->t_obs - Q->t_epoch;
         } else {
-            dt = Q->t_epoch - in.xyzt.t;
+            dt = in.xyzt.t - Q->t_epoch ;
     }
 
     shift = get_grid_shift(P, in.xyz);
@@ -219,7 +219,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ in, PJ *P) {
     out.xyz = in;
 
     if (Q->t_obs != HUGE_VAL) {
-        dt = Q->t_epoch - Q->t_obs;
+        dt = Q->t_obs - Q->t_epoch;
     } else {
         out = proj_coord_error(); /* in the 3D case +t_obs must be specified */
         proj_log_debug(P, "deformation: +t_obs must be specified");
@@ -238,9 +238,9 @@ static PJ_COORD reverse_4d(PJ_COORD in, PJ *P) {
 
 
     if (Q->t_obs != HUGE_VAL) {
-            dt = Q->t_epoch - Q->t_obs;
+            dt = Q->t_obs - Q->t_epoch;
         } else {
-            dt = Q->t_epoch - in.xyzt.t;
+            dt = in.xyzt.t - Q->t_epoch;
     }
 
     out.xyz = reverse_shift(P, in.xyz, dt);
diff --git a/test/gie/deformation.gie b/test/gie/deformation.gie
index 03071d1e..e4d87fe2 100644
--- a/test/gie/deformation.gie
+++ b/test/gie/deformation.gie
@@ -18,6 +18,7 @@ operation   +proj=deformation +xy_grids=alaska +z_grids=egm96_15.gtx
             +t_epoch=2016.0 +t_obs=2000.0 +ellps=GRS80
 -------------------------------------------------------------------------------
 tolerance   0.1 mm
+direction   inverse
 ignore      pjd_err_failed_to_load_grid
 accept      -3004295.5882503074 -1093474.1690603832 5500477.1338251457
 expect      -3004295.7025       -1093474.2106       5500477.3444
@@ -38,6 +39,7 @@ operation   +proj=deformation
             +xy_grids=alaska +z_grids=egm96_15.gtx +t_epoch=2016.0 +ellps=GRS80
 -------------------------------------------------------------------------------
 tolerance   0.1 mm
+direction   inverse
 ignore      pjd_err_failed_to_load_grid
 accept      -3004295.5882503074 -1093474.1690603832 5500477.1338251457 2000.0
 expect      -3004295.7025       -1093474.2106       5500477.3444       2000.0