From 0b5c3622e37dcf697db45e8db05fd4ebf045674b Mon Sep 17 00:00:00 2001
From: Kristian Evers <kristianevers@gmail.com>
Date: Thu, 24 May 2018 15:22:43 +0200
Subject: 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.
---
 src/PJ_vgridshift.c | 63 +++++++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 59 insertions(+), 4 deletions(-)

(limited to 'src/PJ_vgridshift.c')

diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c
index 205806b2..1cd0880e 100644
--- a/src/PJ_vgridshift.c
+++ b/src/PJ_vgridshift.c
@@ -1,12 +1,19 @@
 #define PJ_LIB__
 
+#include <errno.h>
 #include <stddef.h>
+#include <string.h>
+#include <time.h>
 
 #include "proj_internal.h"
 #include "projects.h"
 
 PROJ_HEAD(vgridshift, "Vertical grid shift");
 
+struct pj_opaque_vgridshift {
+    double t_final;
+    double t_epoch;
+};
 
 static XYZ forward_3d(LPZ lpz, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
@@ -37,25 +44,73 @@ static LPZ reverse_3d(XYZ xyz, PJ *P) {
 
 
 static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
-    PJ_COORD point = {{0,0,0,0}};
-    point.xyz = forward_3d (obs.lpz, P);
+    struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) 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.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);
+
+
     return point;
 }
 
 static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
-    PJ_COORD point;
-    point.lpz = reverse_3d (obs.xyz, P);
+    struct pj_opaque_vgridshift *Q = (struct pj_opaque_vgridshift *) 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(vgridshift,0) {
+    struct pj_opaque_vgridshift *Q = pj_calloc (1, sizeof (struct pj_opaque_vgridshift));
+    if (0==Q)
+        return pj_default_destructor (P, ENOMEM);
+    P->opaque = (void *) Q;
 
    if (!pj_param(P->ctx, P->params, "tgrids").i) {
         proj_log_error(P, "vgridshift: +grids parameter missing.");
         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;
+
+
     /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */
     proj_vgrid_init(P, "grids");
 
-- 
cgit v1.2.3