Temporal gridshifting (#1015)

Temporal gridshifts allow [h|v]gridshift operations to be used as step functions
in a pipeline. This is useful in transformations dealing with deformations caused
by earthquakes.

See the included documentation for details.

docs/source/operations/options/t_epoch.rst

@@ -0,0 +1,3 @@
+.. option:: +t_epoch=<time>
+
+    Central epoch of the transformation.

docs/source/operations/options/t_final.rst

@@ -0,0 +1,7 @@
+.. option:: +t_final=<time>
+
+    Final epoch that the coordinate will be propagated to after transformation.
+    The special epoch *now* can be used instead of writing a specific period in
+    time. When *now* is used, it is replaced internally with the epoch of the
+    transformation. This means that the resulting coordinate will be slightly
+    different if carried out again at a later date.

docs/source/operations/pipeline.rst

@@ -83,6 +83,8 @@ does not.
 
 Will result in an error since :ref:`urm5` does not have an inverse operation defined.
 
+.. _global-pipeline-parameter:
+
 **4. Parameters added before the first `+step` are global and will be applied to all steps.**
 
 In the following the GRS80 ellipsoid will be applied to all steps.

docs/source/operations/transformations/hgridshift.rst

@@ -10,25 +10,27 @@ Change of horizontal datum by grid shift.
 
 
 +-----------------+--------------------------------------------------------------------+
-| **Domain**      | 2D                                                                 |
+| **Domain**      | 2D, 3D and 4D                                                      |
 +-----------------+--------------------------------------------------------------------+
-| **Input type**  | Geodetic coordinates.                                              |
+| **Input type**  | Geodetic coordinates (horizontal), meters (vertical),              |
+|                 | decimalyear (temporal)                                             |
 +-----------------+--------------------------------------------------------------------+
-| **Output type** | Geodetic coordinates.                                              |
+| **Output type** | Geodetic coordinates (horizontal), meters (vertical),              |
+|                 | decimalyear (temporal)                                             |
 +-----------------+--------------------------------------------------------------------+
 
 The horizontal grid shift is done by offsetting the planar input coordinates by
 a specific amount determined by the loaded grids. The simplest use case of the
 horizontal grid shift is applying a single grid::
 
-    +hgridshift +grids=nzgr2kgrid0005.gsb
+    +proj=hgridshift +grids=nzgr2kgrid0005.gsb
 
 
 More than one grid can be loaded at the same time, for instance in case the
 dataset needs to be transformed spans several countries. In this example grids
 of the continental US, Alaska and Canada is loaded at the same time::
 
-    +hgridshift +grids=@conus,@alaska,@ntv2_0.gsb,@ntv_can.dat
+    +proj=hgridshift +grids=@conus,@alaska,@ntv2_0.gsb,@ntv_can.dat
 
 The ``@`` in the above example states that the grid is optional, in case the grid
 is not found in the PROJ search path. The list of grids is prioritized so that
@@ -40,13 +42,75 @@ about all three formats can be found in the GDAL documentation. GDAL reads and
 writes all three formats. Using GDAL for construction of new grids is recommended.
 
 
+Temporal gridshifting
+################################################################################
+.. versionadded:: 5.1.0
+
+By initializing the horizontal gridshift operation with a central epoch, it can be
+used as a step function applying the grid offsets only if a coordinate is
+transformed from an epoch before grids central epoch to an epoch after. This is
+handy in transformations where it is necessary to handle deformations caused by
+seismic activity.
+
+The central epoch of the grid is controlled with :option:`+t_epoch` and the final
+epoch of the coordinate is set with :option:`+t_final`. The observation epoch of
+the coordinate is part of the coordinate tuple.
+
+Suppose we want to model the deformation of the 2008 earthquake in Iceland in
+a transformation of data from 2005 to 2009::
+
+    echo 63.992 -21.014 10.0 2005.0 | cct +proj=hgridshift +grids=iceland2008.gsb +t_epoch=2008.4071 +t_final=2009.0
+    63.9920021 -21.0140013 10.0 2005.0
+
+.. note::
+    The timestamp of the resulting coordinate is still 2005.0. The observation
+    time is always kept unchanged as it would otherwise be impossible to do the
+    inverse transformation.
+
+Temporal gridshifting is especially powerful in transformation pipelines where
+several gridshifts can be chained together, effectively acting as a series of
+step functions that can be applied to a coordinate that is propagated through
+time. In the following example we establish a pipeline that allows transformation
+of coordinates from any given epoch up until the current date, applying only
+those gridshifts that have central epochs between the observation epoch and
+the final epoch::
+
+    +proj=pipeline +t_final=now
+    +step +proj=hgridshift +grids=earthquake_1.gsb +t_epoch=2010.421
+    +step +proj=hgridshift +grids=earthquake_2.gsb +t_epoch=2013.853
+    +step +proj=hgridshift +grids=earthquake_3.gsb +t_epoch=2017.713
+
+.. note::
+    The special epoch *now* is used when specifying the final epoch with
+    :option:`+t_final`. This results in coordinates being transformed to the
+    current date. Additionally, :option:`+t_final` is used as a
+    :ref:`global pipeline parameter<global-pipeline-parameter>`, which means
+    that it is applied to all the steps in the pipeline.
+
+In the above transformation, a coordinate with observation epoch 2009.32 would
+be subject to all three gridshift steps in the pipeline. A coordinate with
+observation epoch 2014.12 would only by offset by the last step in the pipeline.
+
+
 Parameters
 ################################################################################
 
+Required
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
 .. option:: +grids=<list>
 
     Comma-separated list of grids to load. If a grid is prefixed by an `@` the
-    grid is consideres optional and PROJ will the not complain if the grid is
+    grid is considered optional and PROJ will the not complain if the grid is
     not available.
 
     Grids are expected to be in CTable2, NTv1 or NTv2 format.
+
+Optional
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+.. include:: ../options/t_epoch.rst
+.. versionadded:: 5.1.0
+.. include:: ../options/t_final.rst
+.. versionadded:: 5.1.0
+

docs/source/operations/transformations/vgridshift.rst

@@ -9,27 +9,29 @@ Vertical grid shift
 Change Vertical datum change by grid shift
 
 +-----------------+--------------------------------------------------------------------+
-| **Domain**      | 3D                                                                 |
+| **Domain**      | 3D and 4D                                                          |
 +-----------------+--------------------------------------------------------------------+
-| **Input type**  | Geodetic coordinates (horizontal), meters (vertical).              |
+| **Input type**  | Geodetic coordinates (horizontal), meters (vertical),              |
+|                 | decimalyear (temporal)                                             |
 +-----------------+--------------------------------------------------------------------+
-| **Output type** | Geodetic coordinates (horizontal), meters (vertical).              |
+| **Output type** | Geodetic coordinates (horizontal), meters (vertical),              |
+|                 | decimalyear (temporal)                                             |
 +-----------------+--------------------------------------------------------------------+
 
 The vertical grid shift is done by offsetting the vertical input coordinates by
 a specific amount determined by the loaded grids. The simplest use case of the
 horizontal grid shift is applying a single grid. Here we change the vertical
 reference from the ellipsoid to the global geoid model, EGM96::
 
-    +vgridshift +grids=egm96_16.gtx
+    +proj=vgridshift +grids=egm96_16.gtx
 
 
 More than one grid can be loaded at the same time, for instance in the case where
 a better geoid model than the global is available for a certain area. Here the
 gridshift is set up so that the local DVR90 geoid model takes precedence over
 the global model::
 
-    +vgridshift +grids=@dvr90.gtx,egm96_16.gtx
+    +proj=vgridshift +grids=@dvr90.gtx,egm96_16.gtx
 
 The ``@`` in the above example states that the grid is optional, in case the grid
 is not found in the PROJ search path. The list of grids is prioritized so that
@@ -40,13 +42,74 @@ PROJ supports the GTX file format for vertical grid corrections. Details
 about all the format can be found in the GDAL documentation. GDAL both reads and
 writes the format. Using GDAL for construction of new grids is recommended.
 
+Temporal gridshifting
+################################################################################
+.. versionadded:: 5.1.0
+
+By initializing the vertical gridshift operation with a central epoch, it can be
+used as a step function applying the grid offsets only if a coordinate is
+transformed from an epoch before grids central epoch to an epoch after. This is
+handy in transformations where it is necessary to handle deformations caused by
+seismic activity.
+
+The central epoch of the grid is controlled with :option:`+t_epoch` and the final
+epoch of the coordinate is set with :option:`+t_final`. The observation epoch of
+the coordinate is part of the coordinate tuple.
+
+Suppose we want to model the deformation of the 2008 earthquake in Iceland in
+a transformation of data from 2005 to 2009::
+
+    echo 63.992 -21.014 10.0 2005.0 | cct +proj=vgridshift +grids=iceland2008.gtx +t_epoch=2008.4071 +t_final=2009.0
+    63.992 -21.014 10.11 2005.0
+
+.. note::
+    The timestamp of the resulting coordinate is still 2005.0. The observation
+    time is always kept unchanged as it would otherwise be impossible to do the
+    inverse transformation.
+
+Temporal gridshifting is especially powerful in transformation pipelines where
+several gridshifts can be chained together, effectively acting as a series of
+step functions that can be applied to a coordinate that is propagated through
+time. In the following example we establish a pipeline that allows transformation
+of coordinates from any given epoch up until the current date, applying only
+those gridshifts that have central epochs between the observation epoch and
+the final epoch::
+
+    +proj=pipeline +t_final=now
+    +step +proj=vgridshift +grids=earthquake_1.gtx +t_epoch=2010.421
+    +step +proj=vgridshift +grids=earthquake_2.gtx +t_epoch=2013.853
+    +step +proj=vgridshift +grids=earthquake_3.gtx +t_epoch=2017.713
+
+.. note::
+    The special epoch *now* is used when specifying the final epoch with
+    :option:`+t_final`. This results in coordinates being transformed to the
+    current date. Additionally, :option:`+t_final` is used as a
+    :ref:`global pipeline parameter<global-pipeline-parameter>`, which means
+    that it is applied to all the steps in the pipeline.
+
+In the above transformation, a coordinate with observation epoch 2009.32 would
+be subject to all three gridshift steps in the pipeline. A coordinate with
+observation epoch 2014.12 would only by offset by the last step in the pipeline.
+
 Parameters
 ################################################################################
 
+Required
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
 .. option:: +grids=<list>
 
     Comma-separated list of grids to load. If a grid is prefixed by an `@` the
-    grid is consideres optional and PROJ will the not complain if the grid is
+    grid is considered optional and PROJ will the not complain if the grid is
     not available.
 
     Grids are expected to be in GTX format.
+
+Optional
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+.. include:: ../options/t_epoch.rst
+.. versionadded:: 5.1.0
+.. include:: ../options/t_final.rst
+.. versionadded:: 5.1.0
+

src/PJ_hgridshift.c

@@ -1,12 +1,20 @@
 #define PJ_LIB__
 
+#include <errno.h>
+#include <string.h>
 #include <stddef.h>
+#include <time.h>
 
 #include "proj_internal.h"
 #include "projects.h"
 
 PROJ_HEAD(hgridshift, "Horizontal grid shift");
 
+struct pj_opaque_hgridshift {
+    double t_final;
+    double t_epoch;
+};
+
 static XYZ forward_3d(LPZ lpz, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
     point.lpz = lpz;
@@ -34,21 +42,47 @@ static LPZ reverse_3d(XYZ xyz, PJ *P) {
     return point.lpz;
 }
 
+static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
+    struct pj_opaque_hgridshift *Q = (struct pj_opaque_hgridshift *) P->opaque;
+    PJ_COORD point = obs;
 
-static PJ_COORD forward_4d (PJ_COORD obs, PJ *P) {
-    obs.xyz = forward_3d (obs.lpz, P);
-    return obs;
-}
+    /* If transformation is not time restricted, we always call it */
+    if (Q->t_final==0 || Q->t_epoch==0) {
+        point.xyz = forward_3d (obs.lpz, P);
+        return point;
+    }
 
+    /* Time restricted - only apply transform if within time bracket */
+    if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
+        point.xyz = forward_3d (obs.lpz, P);
 
-static PJ_COORD reverse_4d (PJ_COORD obs, PJ *P) {
-    obs.lpz = reverse_3d (obs.xyz, P);
-    return obs;
+
+    return point;
 }
 
+static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
+    struct pj_opaque_hgridshift *Q = (struct pj_opaque_hgridshift *) P->opaque;
+    PJ_COORD point = obs;
+
+    /* If transformation is not time restricted, we always call it */
+    if (Q->t_final==0 || Q->t_epoch==0) {
+        point.lpz = reverse_3d (obs.xyz, P);
+        return point;
+    }
+
+    /* Time restricted - only apply transform if within time bracket */
+    if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
+        point.lpz = reverse_3d (obs.xyz, P);
+
+    return point;
+}
 
 
 PJ *TRANSFORMATION(hgridshift,0) {
+    struct pj_opaque_hgridshift *Q = pj_calloc (1, sizeof (struct pj_opaque_hgridshift));
+    if (0==Q)
+        return pj_default_destructor (P, ENOMEM);
+    P->opaque = (void *) Q;
 
     P->fwd4d  = forward_4d;
     P->inv4d  = reverse_4d;
@@ -65,6 +99,26 @@ PJ *TRANSFORMATION(hgridshift,0) {
         return pj_default_destructor (P, PJD_ERR_NO_ARGS);
     }
 
+   /* TODO: Refactor into shared function that can be used  */
+   /*       by both vgridshift and hgridshift               */
+   if (pj_param(P->ctx, P->params, "tt_final").i) {
+        Q->t_final = pj_param (P->ctx, P->params, "dt_final").f;
+        if (Q->t_final == 0) {
+            /* a number wasn't passed to +t_final, let's see if it was "now" */
+            /* and set the time accordingly.                                 */
+            if (!strcmp("now", pj_param(P->ctx, P->params, "st_final").s)) {
+                    time_t now;
+                    struct tm *date;
+                    time(&now);
+                    date = localtime(&now);
+                    Q->t_final = 1900.0 + date->tm_year + date->tm_yday/365.0;
+            }
+        }
+    }
+
+   if (pj_param(P->ctx, P->params, "tt_epoch").i)
+        Q->t_epoch = pj_param (P->ctx, P->params, "dt_epoch").f;
+
 
     proj_hgrid_init(P, "grids");
     /* Was gridlist compiled properly? */