Temporal gridshifting #1015

kbevers · 2018-05-24T07:02:33Z

Extending the gridshift operations into the fourth dimension. The idea is to be able to use the grid shift operations as step functions when making transformations in areas subject to deformation after seismic events. This is useful when transforming coordinate from one epoch to another, e.g. in dynamic reference frames. This work is closely related to the proposed deformation model grid format in #1001.

I will expand the docs later to today, for now I am just checking that the code builds properly.

rouault · 2018-05-24T09:42:39Z

src/PJ_vgridshift.c


   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);
    }

+   if (pj_param(P->ctx, P->params, "tt_final").i) {
+        Q->t_final = pj_param (P->ctx, P->params, "dt_final").f;


we should likely have a pj_eval_param_time(const char*) function to do the below processing, since it is repeated above.
localtime() is known on POSIX to be thread-unsafe. It is better to use localtime_r() instead, when available

we should likely have a pj_eval_param_time(const char*) function to do the below processing, since it is repeated above.

Good point. I intend to add similar functionality to deformation and possibly helmert as well. I will take care of it then.

localtime() is known on POSIX to be thread-unsafe. It is better to use localtime_r() instead, when available

When can localtime_r be expected to be available?

When can localtime_r be expected to be available?

We'd need a configure / cmake check to test for its availability. Or alternatively "port" GDAL's https://github.com/OSGeo/gdal/blob/master/gdal/port/cpl_time.cpp#L64 which as indicated in the file header comes from public domain code.

I think it is a bit out of scope for this PR. How big a risk is it to include the code as it is now? localtime() is only called once during creation of the PJ object.

The risk is probably low. A TODO in the code + issue to remember should be good enough

Created #1016

You can probably significantly reduce (but not eliminate) the risk by saying:

struct tm cucumber; ... cucumber = *localtime(...);

Essentially doing a copy of the the localtimeinternal struct tm, rather than storing a pointer to it.

rouault · 2018-05-24T09:48:50Z

src/PJ_hgridshift.c

+    struct pj_opaque_hgridshift *Q = (struct pj_opaque_hgridshift *) P->opaque;
+    PJ_COORD point = obs;
+    if (Q->t_final != 0.0  && Q->t_epoch != 0.0) {
+        if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)


shouldn't that be be <= and >= ?

Or so if t_final and t_epoch are specified, we apply the grid shift only if the observation time is below the epoch, and the epoch below the final epoch. That doesn't immediately seem obvious to me. Should probably be explained in the .rst files

I have just updated the .rst files. Do they cover the topic better now?

From a physical point of view I would say that the case where obs.lpzt.t == Q->t_epoch is undefined, since the shift doesn't happen instantaneously but rather over a timespan of minutes to days. To me it makes the most sense to say that the offset is applied only if the observation epoch is before the event that caused the deformation. Should you happen to have acquired a coordinate during the deformation event it will be partly to fully shifted and we dont' really have the means to do anything about that.

A similar point can be made for the second part of the condition, but I think from a practical point of view it make sense to change it to Q->t_final >= Q->t_epoch.

Yes the docs make things really clear (sorry I started looking at the code before reading your PR general comment where you mentioned that would be later detailed :-) )
And actually no strong opinion if that should be strict or equal comparisons.

rouault · 2018-05-24T09:49:13Z

docs/source/operations/transformations/vgridshift.rst

@@ -11,9 +11,11 @@ Change Vertical datum change by grid shift
 +-----------------+--------------------------------------------------------------------+
 | **Domain**      | 3D                                                                 |
 +-----------------+--------------------------------------------------------------------+
-| **Input type**  | Geodetic coordinates (horizontal), meters (vertical).              |
+| **Input type**  | Geodetic coordinates (horizontal), meters (vertical),              |


hgridshift.rst should also be similarly updated

busstoptaktik

This review consists of just one comment about code flow and readability, but it applies 4 places :-)

busstoptaktik · 2018-05-24T11:15:32Z

src/PJ_hgridshift.c

-    return obs;
+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;


I think the code below (and the similar in the forward_4d case) will be more readable if the conditionals are inverted, so the if/if/else-branches can be decoupled:

/* 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;

Good suggestion. Updated the code.

busstoptaktik · 2018-05-24T11:17:03Z

src/PJ_vgridshift.c

-    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;


cf. comment in PJ_hgridshift.c regarding inverted conditionals

busstoptaktik · 2018-05-24T11:20:27Z

src/PJ_vgridshift.c

-    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;


cf. comment in PJ_hgridshift.c regarding inverted conditionals

kbevers · 2018-05-24T13:20:22Z

@rouault and @busstoptaktik, thanks for the swift and helpful reviews. I am merging this now so that it can be included in 5.1.0.

kbevers added 3 commits May 24, 2018 08:41

First take on temporal gridshift

9d7db9d

Merge branch 'master' into temporal-gridshift

29c2b56

fix grids

ca453c1

rouault reviewed May 24, 2018

View reviewed changes

Improve docs for [h|v]gridshift

3029de1

busstoptaktik reviewed May 24, 2018

View reviewed changes

kbevers mentioned this pull request May 24, 2018

Centralize time functionality #1016

Closed

3 tasks

More readable code courtesy of Thomas

1e03829

kbevers force-pushed the temporal-gridshift branch from 9322ec2 to 1e03829 Compare May 24, 2018 12:03

Add tests

d5f6453

busstoptaktik assigned kbevers May 24, 2018

kbevers merged commit 0b5c362 into OSGeo:master May 24, 2018

kbevers added this to the 5.1.0 milestone May 24, 2018

kbevers deleted the temporal-gridshift branch July 9, 2018 08:46

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Temporal gridshifting #1015

Temporal gridshifting #1015

kbevers commented May 24, 2018

rouault May 24, 2018 •

edited

kbevers May 24, 2018

rouault May 24, 2018

kbevers May 24, 2018

rouault May 24, 2018

kbevers May 24, 2018

busstoptaktik May 24, 2018

rouault May 24, 2018

kbevers May 24, 2018

rouault May 24, 2018

rouault May 24, 2018

busstoptaktik left a comment

busstoptaktik May 24, 2018

kbevers May 24, 2018

busstoptaktik May 24, 2018

busstoptaktik May 24, 2018

kbevers commented May 24, 2018

Temporal gridshifting #1015

Temporal gridshifting #1015

Conversation

kbevers commented May 24, 2018

rouault May 24, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

busstoptaktik left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kbevers commented May 24, 2018

rouault May 24, 2018 •

edited