From 0b5c3622e37dcf697db45e8db05fd4ebf045674b Mon Sep 17 00:00:00 2001 From: Kristian Evers 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_hgridshift.c | 68 +++++++++++++++++++++++++++++++++++++++++++++++------ src/PJ_vgridshift.c | 63 +++++++++++++++++++++++++++++++++++++++++++++---- 2 files changed, 120 insertions(+), 11 deletions(-) (limited to 'src') diff --git a/src/PJ_hgridshift.c b/src/PJ_hgridshift.c index 3932ef51..2919a85c 100644 --- a/src/PJ_hgridshift.c +++ b/src/PJ_hgridshift.c @@ -1,12 +1,20 @@ #define PJ_LIB__ +#include +#include #include +#include #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? */ 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 #include +#include +#include #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