Skip to content
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.

Sign up for GitHub

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

Jump to bottom

Temporal gridshifting #1015

Merged
merged 6 commits into from May 24, 2018
Merged

Temporal gridshifting #1015

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
9d7db9d
First take on temporal gridshift
kbevers May 7, 2018
29c2b56
Merge branch 'master' into temporal-gridshift
kbevers May 24, 2018
ca453c1
fix grids
kbevers May 24, 2018
3029de1
Improve docs for [h|v]gridshift
kbevers May 24, 2018
1e03829
More readable code courtesy of Thomas
kbevers May 24, 2018
d5f6453
Add tests
kbevers May 24, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
10 changes: 8 additions & 2 deletions docs/source/operations/transformations/vgridshift.rst
View file
Open in desktop
Expand Up @@ -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), |
Copy link
Member

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

| | 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
Expand Down Expand Up @@ -50,3 +52,7 @@ Parameters
not available.

Grids are expected to be in GTX format.

.. option:: +t_final=<time>

.. option:: +t_epoch=<time>
59 changes: 50 additions & 9 deletions src/PJ_hgridshift.c
View file
Open in desktop
@@ -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;
Expand Down Expand Up @@ -34,21 +42,36 @@ static LPZ reverse_3d(XYZ xyz, PJ *P) {
return point.lpz;
}


static PJ_COORD forward_4d (PJ_COORD obs, PJ *P) {
obs.xyz = forward_3d (obs.lpz, P);
return obs;
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;
if (Q->t_final != 0.0 && Q->t_epoch != 0.0) {
if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
Copy link
Member

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

Copy link
Member Author

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.

Copy link
Member

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.

point.xyz = forward_3d (obs.lpz, P);
} else {
point.xyz = forward_3d (obs.lpz, P);
}
return point;
}


static PJ_COORD reverse_4d (PJ_COORD obs, PJ *P) {
obs.lpz = reverse_3d (obs.xyz, P);
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;
Copy link
Member

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;

Copy link
Member Author

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.

if (Q->t_final != 0.0 && Q->t_epoch != 0.0) {
if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
point.lpz = reverse_3d (obs.xyz, P);
} else {
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;
Expand All @@ -65,6 +88,24 @@ PJ *TRANSFORMATION(hgridshift,0) {
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;
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? */
Expand Down
50 changes: 46 additions & 4 deletions src/PJ_vgridshift.c
View file
Open in desktop
@@ -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}};
Expand Down Expand Up @@ -37,25 +44,60 @@ 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;
Copy link
Member

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

if (Q->t_final != 0.0 && Q->t_epoch != 0.0) {
if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
point.xyz = forward_3d (obs.lpz, P);
} else {
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;
Copy link
Member

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

if (Q->t_final != 0.0 && Q->t_epoch != 0.0) {
if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
point.lpz = reverse_3d (obs.xyz, P);
} else {
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);
}

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

@rouault rouault May 24, 2018 •
edited

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

Copy link
Member Author

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?

Copy link
Member

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.

Copy link
Member Author

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.

Copy link
Member

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

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Created #1016

Copy link
Member

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 localtimeinternal struct tm, rather than storing a pointer to it.

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

Expand Down
73 changes: 73 additions & 0 deletions test/gie/deformation.gie
View file
Open in desktop
Expand Up @@ -62,4 +62,77 @@ expect failure pjd_err_failed_to_load_grid

operation proj=deformation xy_grids=alaska z_grids=nonexisting ellps=GRS80
expect failure pjd_err_missing_args

-------------------------------------------------------------------------------
operation +proj=vgridshift +grids=egm96_15.gtx +t_epoch=2010.0 +t_final=2018.0
-------------------------------------------------------------------------------
tolerance 0.1 mm
ignore pjd_err_failed_to_load_grid

accept 12 56 0.0 2000.0
expect 12 56 -36.5966 2000.0

accept 12 56 0.0 2011.0
expect 12 56 0.0 2011.0

accept 12 56 0.0 2019.0
expect 12 56 0.0 2019.0

accept 12 56 0.0
expect 12 56 -36.5966


-------------------------------------------------------------------------------
operation +proj=vgridshift +grids=egm96_15.gtx +t_epoch=2010.0 +t_final=now
-------------------------------------------------------------------------------
tolerance 0.1 mm
ignore pjd_err_failed_to_load_grid

accept 12 56 0.0 2000.0
expect 12 56 -36.5966 2000.0

accept 12 56 0.0 2011.0
expect 12 56 0.0 2011.0

accept 12 56 0.0 3011.0
expect 12 56 0.0 3011.0


-------------------------------------------------------------------------------
operation +proj=hgridshift +grids=alaska +t_epoch=2010.0 +t_final=2018.0
-------------------------------------------------------------------------------
tolerance 0.1 mm
ignore pjd_err_failed_to_load_grid

accept -147.0 64.0 0.0 2000.0
expect -147.0023233121 63.9995792119 0.0 2000.0
roundtrip 100

accept -147.0 64.0 0.0 2011.0
expect -147.0 64.0 0.0 2011.0
roundtrip 100

accept -147.0 64.0 0.0 2011.0
expect -147.0 64.0 0.0 2020.0
roundtrip 100

-------------------------------------------------------------------------------
operation +proj=hgridshift +grids=alaska +t_epoch=2010.0 +t_final=now
-------------------------------------------------------------------------------
tolerance 0.1 mm
ignore pjd_err_failed_to_load_grid

accept -147.0 64.0 0.0 2000.0
expect -147.0023233121 63.9995792119 0.0 2000.0
roundtrip 100

accept -147.0 64.0 0.0 2011.0
expect -147.0 64.0 0.0 2011.0
roundtrip 100

accept -147.0 64.0 0.0 3011.0
expect -147.0 64.0 0.0 3011.0
roundtrip 100


</gie>