-
Notifications
You must be signed in to change notification settings - Fork 802
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Temporal gridshifting #1015
Temporal gridshifting #1015
Conversation
src/PJ_vgridshift.c
Outdated
|
||
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; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The risk is probably low. A TODO in the code + issue to remember should be good enough
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Created #1016
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You can probably significantly reduce (but not eliminate) the risk by saying:
struct tm cucumber;
...
cucumber = *localtime(...);
Essentially doing a copy of the the localtime
internal struct tm
, rather than storing a pointer to it.
src/PJ_hgridshift.c
Outdated
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) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
@@ -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), | |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
hgridshift.rst should also be similarly updated
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This review consists of just one comment about code flow and readability, but it applies 4 places :-)
src/PJ_hgridshift.c
Outdated
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; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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;
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Good suggestion. Updated the code.
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; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
cf. comment in PJ_hgridshift.c regarding inverted conditionals
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; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
cf. comment in PJ_hgridshift.c regarding inverted conditionals
9322ec2
to
1e03829
Compare
@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. |
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.