diff options
| -rw-r--r-- | docs/source/apps/gie.rst | 21 | ||||
| -rw-r--r-- | docs/source/index.rst | 4 | ||||
| -rw-r--r-- | docs/source/operations/options/t_epoch.rst | 3 | ||||
| -rw-r--r-- | docs/source/operations/options/t_final.rst | 7 | ||||
| -rw-r--r-- | docs/source/operations/pipeline.rst | 2 | ||||
| -rw-r--r-- | docs/source/operations/transformations/helmert.rst | 2 | ||||
| -rw-r--r-- | docs/source/operations/transformations/hgridshift.rst | 76 | ||||
| -rw-r--r-- | docs/source/operations/transformations/vgridshift.rst | 75 | ||||
| -rw-r--r-- | docs/source/resource_files.rst | 6 | ||||
| -rw-r--r-- | docs/source/usage/differences.rst | 47 | ||||
| -rw-r--r-- | docs/source/usage/index.rst | 2 | ||||
| -rw-r--r-- | docs/source/usage/transformation.rst | 2 | ||||
| -rw-r--r-- | src/Makefile.am | 4 | ||||
| -rw-r--r-- | src/PJ_hgridshift.c | 68 | ||||
| -rw-r--r-- | src/PJ_vgridshift.c | 63 | ||||
| -rw-r--r-- | src/bin_cct.cmake | 6 | ||||
| -rw-r--r-- | src/bin_gie.cmake | 6 | ||||
| -rw-r--r-- | src/cct.c | 4 | ||||
| -rw-r--r-- | src/gie.c | 31 | ||||
| -rw-r--r-- | src/proj_strtod.c | 6 | ||||
| -rw-r--r-- | src/proj_strtod.h | 4 | ||||
| -rw-r--r-- | test/gie/deformation.gie | 80 |
22 files changed, 466 insertions, 53 deletions
diff --git a/docs/source/apps/gie.rst b/docs/source/apps/gie.rst index b6c1c166..9d3fa0c2 100644 --- a/docs/source/apps/gie.rst +++ b/docs/source/apps/gie.rst @@ -207,17 +207,32 @@ gie command language to function. The accepted coordinate is passed to the operation first in it's forward mode, then the output from the forward operation is passed back to the inverse operation. This procedure is done ``n`` times. If the - resulting coordinate is within the set tolerance of the initial coordinate + resulting coordinate is within the set tolerance of the initial coordinate, the test is passed. - Example: + Example with the default 100 iterations and the default tolerance: .. code-block:: console operation proj=merc accept 12 55 - roundtrip 10000 5 mm + roundtrip + + Example with count and default tolerance: + + .. code-block:: console + + operation proj=merc + accept 12 55 + roundtrip 10000 + Example with count and tolerance: + + .. code-block:: console + + operation proj=merc + accept 12 55 + roundtrip 10000 5 mm .. option:: direction <direction> diff --git a/docs/source/index.rst b/docs/source/index.rst index 86393291..9e097f3b 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -50,6 +50,10 @@ PROJ removed from PROJ in version 7.0.0 scheduled for release February 1st 2020. + .. attention:: + + With the introduction of PROJ 5, behavioural changes has been made to + existing functionality. Consult :ref:`differences` for the details. diff --git a/docs/source/operations/options/t_epoch.rst b/docs/source/operations/options/t_epoch.rst new file mode 100644 index 00000000..76a4473a --- /dev/null +++ b/docs/source/operations/options/t_epoch.rst @@ -0,0 +1,3 @@ +.. option:: +t_epoch=<time> + + Central epoch of the transformation. diff --git a/docs/source/operations/options/t_final.rst b/docs/source/operations/options/t_final.rst new file mode 100644 index 00000000..37908209 --- /dev/null +++ b/docs/source/operations/options/t_final.rst @@ -0,0 +1,7 @@ +.. option:: +t_final=<time> + + Final epoch that the coordinate will be propagated to after transformation. + The special epoch *now* can be used instead of writing a specific period in + time. When *now* is used, it is replaced internally with the epoch of the + transformation. This means that the resulting coordinate will be slightly + different if carried out again at a later date. diff --git a/docs/source/operations/pipeline.rst b/docs/source/operations/pipeline.rst index 9223202d..33112536 100644 --- a/docs/source/operations/pipeline.rst +++ b/docs/source/operations/pipeline.rst @@ -83,6 +83,8 @@ does not. Will result in an error since :ref:`urm5` does not have an inverse operation defined. +.. _global-pipeline-parameter: + **4. Parameters added before the first `+step` are global and will be applied to all steps.** In the following the GRS80 ellipsoid will be applied to all steps. diff --git a/docs/source/operations/transformations/helmert.rst b/docs/source/operations/transformations/helmert.rst index 6a3751df..c2557321 100644 --- a/docs/source/operations/transformations/helmert.rst +++ b/docs/source/operations/transformations/helmert.rst @@ -72,7 +72,7 @@ Transformation from `ITRF2000@2017.0` to `ITRF93@2017.0` using 15 parameters: dx=-0.0029 dy=-0.0002 dz=-0.0006 ds=0.00001 rx=-0.00039 ry=0.00080 rz=-0.00114 drx=-0.00011 dry=-0.00019 drz=0.00007 - epoch=1988.0 tobs=2017.0 transpose + t_epoch=1988.0 t_obs=2017.0 transpose Parameters +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ diff --git a/docs/source/operations/transformations/hgridshift.rst b/docs/source/operations/transformations/hgridshift.rst index 35c84919..1d864a9e 100644 --- a/docs/source/operations/transformations/hgridshift.rst +++ b/docs/source/operations/transformations/hgridshift.rst @@ -10,25 +10,27 @@ Change of horizontal datum by grid shift. +-----------------+--------------------------------------------------------------------+ -| **Domain** | 2D | +| **Domain** | 2D, 3D and 4D | +-----------------+--------------------------------------------------------------------+ -| **Input type** | Geodetic coordinates. | +| **Input type** | Geodetic coordinates (horizontal), meters (vertical), | +| | decimalyear (temporal) | +-----------------+--------------------------------------------------------------------+ -| **Output type** | Geodetic coordinates. | +| **Output type** | Geodetic coordinates (horizontal), meters (vertical), | +| | decimalyear (temporal) | +-----------------+--------------------------------------------------------------------+ The horizontal grid shift is done by offsetting the planar input coordinates by a specific amount determined by the loaded grids. The simplest use case of the horizontal grid shift is applying a single grid:: - +hgridshift +grids=nzgr2kgrid0005.gsb + +proj=hgridshift +grids=nzgr2kgrid0005.gsb More than one grid can be loaded at the same time, for instance in case the dataset needs to be transformed spans several countries. In this example grids of the continental US, Alaska and Canada is loaded at the same time:: - +hgridshift +grids=@conus,@alaska,@ntv2_0.gsb,@ntv_can.dat + +proj=hgridshift +grids=@conus,@alaska,@ntv2_0.gsb,@ntv_can.dat The ``@`` in the above example states that the grid is optional, in case the grid is not found in the PROJ search path. The list of grids is prioritized so that @@ -40,13 +42,75 @@ about all three formats can be found in the GDAL documentation. GDAL reads and writes all three formats. Using GDAL for construction of new grids is recommended. +Temporal gridshifting +################################################################################ +.. versionadded:: 5.1.0 + +By initializing the horizontal gridshift operation with a central epoch, it can be +used as a step function applying the grid offsets only if a coordinate is +transformed from an epoch before grids central epoch to an epoch after. This is +handy in transformations where it is necessary to handle deformations caused by +seismic activity. + +The central epoch of the grid is controlled with :option:`+t_epoch` and the final +epoch of the coordinate is set with :option:`+t_final`. The observation epoch of +the coordinate is part of the coordinate tuple. + +Suppose we want to model the deformation of the 2008 earthquake in Iceland in +a transformation of data from 2005 to 2009:: + + echo 63.992 -21.014 10.0 2005.0 | cct +proj=hgridshift +grids=iceland2008.gsb +t_epoch=2008.4071 +t_final=2009.0 + 63.9920021 -21.0140013 10.0 2005.0 + +.. note:: + The timestamp of the resulting coordinate is still 2005.0. The observation + time is always kept unchanged as it would otherwise be impossible to do the + inverse transformation. + +Temporal gridshifting is especially powerful in transformation pipelines where +several gridshifts can be chained together, effectively acting as a series of +step functions that can be applied to a coordinate that is propagated through +time. In the following example we establish a pipeline that allows transformation +of coordinates from any given epoch up until the current date, applying only +those gridshifts that have central epochs between the observation epoch and +the final epoch:: + + +proj=pipeline +t_final=now + +step +proj=hgridshift +grids=earthquake_1.gsb +t_epoch=2010.421 + +step +proj=hgridshift +grids=earthquake_2.gsb +t_epoch=2013.853 + +step +proj=hgridshift +grids=earthquake_3.gsb +t_epoch=2017.713 + +.. note:: + The special epoch *now* is used when specifying the final epoch with + :option:`+t_final`. This results in coordinates being transformed to the + current date. Additionally, :option:`+t_final` is used as a + :ref:`global pipeline parameter<global-pipeline-parameter>`, which means + that it is applied to all the steps in the pipeline. + +In the above transformation, a coordinate with observation epoch 2009.32 would +be subject to all three gridshift steps in the pipeline. A coordinate with +observation epoch 2014.12 would only by offset by the last step in the pipeline. + + Parameters ################################################################################ +Required ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + .. option:: +grids=<list> Comma-separated list of grids to load. If a grid is prefixed by an `@` the - grid is consideres optional and PROJ will the not complain if the grid is + grid is considered optional and PROJ will the not complain if the grid is not available. Grids are expected to be in CTable2, NTv1 or NTv2 format. + +Optional ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. include:: ../options/t_epoch.rst +.. versionadded:: 5.1.0 +.. include:: ../options/t_final.rst +.. versionadded:: 5.1.0 + diff --git a/docs/source/operations/transformations/vgridshift.rst b/docs/source/operations/transformations/vgridshift.rst index ca2b421c..e671779b 100644 --- a/docs/source/operations/transformations/vgridshift.rst +++ b/docs/source/operations/transformations/vgridshift.rst @@ -9,11 +9,13 @@ Vertical grid shift Change Vertical datum change by grid shift +-----------------+--------------------------------------------------------------------+ -| **Domain** | 3D | +| **Domain** | 3D and 4D | +-----------------+--------------------------------------------------------------------+ -| **Input type** | Geodetic coordinates (horizontal), meters (vertical). | +| **Input type** | Geodetic coordinates (horizontal), meters (vertical), | +| | 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 @@ -21,7 +23,7 @@ a specific amount determined by the loaded grids. The simplest use case of the horizontal grid shift is applying a single grid. Here we change the vertical reference from the ellipsoid to the global geoid model, EGM96:: - +vgridshift +grids=egm96_16.gtx + +proj=vgridshift +grids=egm96_16.gtx More than one grid can be loaded at the same time, for instance in the case where @@ -29,7 +31,7 @@ a better geoid model than the global is available for a certain area. Here the gridshift is set up so that the local DVR90 geoid model takes precedence over the global model:: - +vgridshift +grids=@dvr90.gtx,egm96_16.gtx + +proj=vgridshift +grids=@dvr90.gtx,egm96_16.gtx The ``@`` in the above example states that the grid is optional, in case the grid is not found in the PROJ search path. The list of grids is prioritized so that @@ -40,13 +42,74 @@ PROJ supports the GTX file format for vertical grid corrections. Details about all the format can be found in the GDAL documentation. GDAL both reads and writes the format. Using GDAL for construction of new grids is recommended. +Temporal gridshifting +################################################################################ +.. versionadded:: 5.1.0 + +By initializing the vertical gridshift operation with a central epoch, it can be +used as a step function applying the grid offsets only if a coordinate is +transformed from an epoch before grids central epoch to an epoch after. This is +handy in transformations where it is necessary to handle deformations caused by +seismic activity. + +The central epoch of the grid is controlled with :option:`+t_epoch` and the final +epoch of the coordinate is set with :option:`+t_final`. The observation epoch of +the coordinate is part of the coordinate tuple. + +Suppose we want to model the deformation of the 2008 earthquake in Iceland in +a transformation of data from 2005 to 2009:: + + echo 63.992 -21.014 10.0 2005.0 | cct +proj=vgridshift +grids=iceland2008.gtx +t_epoch=2008.4071 +t_final=2009.0 + 63.992 -21.014 10.11 2005.0 + +.. note:: + The timestamp of the resulting coordinate is still 2005.0. The observation + time is always kept unchanged as it would otherwise be impossible to do the + inverse transformation. + +Temporal gridshifting is especially powerful in transformation pipelines where +several gridshifts can be chained together, effectively acting as a series of +step functions that can be applied to a coordinate that is propagated through +time. In the following example we establish a pipeline that allows transformation +of coordinates from any given epoch up until the current date, applying only +those gridshifts that have central epochs between the observation epoch and +the final epoch:: + + +proj=pipeline +t_final=now + +step +proj=vgridshift +grids=earthquake_1.gtx +t_epoch=2010.421 + +step +proj=vgridshift +grids=earthquake_2.gtx +t_epoch=2013.853 + +step +proj=vgridshift +grids=earthquake_3.gtx +t_epoch=2017.713 + +.. note:: + The special epoch *now* is used when specifying the final epoch with + :option:`+t_final`. This results in coordinates being transformed to the + current date. Additionally, :option:`+t_final` is used as a + :ref:`global pipeline parameter<global-pipeline-parameter>`, which means + that it is applied to all the steps in the pipeline. + +In the above transformation, a coordinate with observation epoch 2009.32 would +be subject to all three gridshift steps in the pipeline. A coordinate with +observation epoch 2014.12 would only by offset by the last step in the pipeline. + Parameters ################################################################################ +Required ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + .. option:: +grids=<list> Comma-separated list of grids to load. If a grid is prefixed by an `@` the - grid is consideres optional and PROJ will the not complain if the grid is + grid is considered optional and PROJ will the not complain if the grid is not available. Grids are expected to be in GTX format. + +Optional ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. include:: ../options/t_epoch.rst +.. versionadded:: 5.1.0 +.. include:: ../options/t_final.rst +.. versionadded:: 5.1.0 + diff --git a/docs/source/resource_files.rst b/docs/source/resource_files.rst index b6aa950d..10cbcb7b 100644 --- a/docs/source/resource_files.rst +++ b/docs/source/resource_files.rst @@ -404,7 +404,7 @@ ITRF2005: :: - +init=ITRF2000:ITRF2005 +tobs=2010.5 + +init=ITRF2000:ITRF2005 +t_obs=2010.5 which then expands to @@ -412,8 +412,8 @@ which then expands to +proj=helmert +x=-0.0001 +y=0.0008 +z=0.0058 +s=-0.0004 +dx=0.0002 +dy=-0.0001 +dz=0.0018 +ds=-0.000008 - +epoch=2000.0 +transpose - +tobs=2010.5 + +t_epoch=2000.0 +transpose + +t_obs=2010.5 Below is a list of the init files that are packaged with PROJ. diff --git a/docs/source/usage/differences.rst b/docs/source/usage/differences.rst new file mode 100644 index 00000000..58e92cd8 --- /dev/null +++ b/docs/source/usage/differences.rst @@ -0,0 +1,47 @@ +.. _differences: + +================================================================================ +Known differences between versions +================================================================================ + +Once in a while, a new version of PROJ causes changes in the existing behaviour. +In this section we track deliberate changes to PROJ that break from previous +behaviour. Most times that will be caused by a bug fix. Unfortunately, some bugs +have existed for so long that their faulty behaviour is relied upon by software +that uses PROJ. + +Behavioural changes caused by new bugs are not tracked here, as they should be +fixed in later versions of PROJ. + +Version 5.0.0 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +Longitude wrapping when using custom central meridian +------------------------------------------------------------------------------- + +By default PROJ wraps output longitudes in the range -180 to 180. Previous to +PROJ 5, this was handled incorrectly when a custom central meridian was set with +:option:`+lon`. This caused a change in sign on the resulting easting as seen +below:: + + $ proj +proj=merc +lon_0=110 + -70 0 + -20037508.34 0.00 + 290 0 + 20037508.34 0.00 + +From PROJ 5 on onwards, the same input now results in same coordinates, as seen +from the example below where PROJ 5 is used:: + + $ proj +proj=merc +lon_0=110 + -70 0 + -20037508.34 0.00 + 290 0 + -20037508.34 0.00 + +The change is made on the basis that :math:`\lambda=290^{\circ}` is a full +rotation of the circle larger than :math:`\lambda=-70^{\circ}` and hence +should return the same output for both. + +Adding the :option:`+over` flag to the projection definition provides +the old behaviour. diff --git a/docs/source/usage/index.rst b/docs/source/usage/index.rst index 7de0e147..823e8fe7 100644 --- a/docs/source/usage/index.rst +++ b/docs/source/usage/index.rst @@ -16,4 +16,4 @@ command line applications or the C API that is a part of the software package. projections transformation environmentvars - + differences diff --git a/docs/source/usage/transformation.rst b/docs/source/usage/transformation.rst index bd32bcd9..6a80f9c6 100644 --- a/docs/source/usage/transformation.rst +++ b/docs/source/usage/transformation.rst @@ -130,7 +130,7 @@ coordinate timestamps back to GPS weeks. rx=0.00039 ry=-0.00080 rz=0.00114 dx=-0.0029 dy=-0.0002 dz=-0.0006 ds=0.00001 drx=0.00011 dry=0.00019 drz=-0.00007 - epoch=1988.0 + t_epoch=1988.0 step proj=unitconvert t_in=decimalyear t_out=gps_week diff --git a/src/Makefile.am b/src/Makefile.am index 971ac78c..2fd8aa5c 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -18,11 +18,11 @@ EXTRA_DIST = makefile.vc proj.def bin_cct.cmake bin_gie.cmake bin_cs2cs.cmake \ proj_SOURCES = proj.c gen_cheb.c p_series.c cs2cs_SOURCES = cs2cs.c gen_cheb.c p_series.c -cct_SOURCES = cct.c proj_strtod.c optargpm.h +cct_SOURCES = cct.c proj_strtod.c proj_strtod.h optargpm.h nad2bin_SOURCES = nad2bin.c geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h -gie_SOURCES = gie.c proj_strtod.c optargpm.h +gie_SOURCES = gie.c proj_strtod.c proj_strtod.h optargpm.h multistresstest_SOURCES = multistresstest.c test228_SOURCES = test228.c geodtest_SOURCES = geodtest.c 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 <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; @@ -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 <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"); diff --git a/src/bin_cct.cmake b/src/bin_cct.cmake index a204e7e7..caa261a8 100644 --- a/src/bin_cct.cmake +++ b/src/bin_cct.cmake @@ -1,9 +1,9 @@ -set(CCT_SRC cct.c proj_strtod.c) +set(CCT_SRC cct.c proj_strtod.c proj_strtod.h) set(CCT_INCLUDE optargpm.h) source_group("Source Files\\Bin" FILES ${CCT_SRC}) add_executable(cct ${CCT_SRC} ${CCT_INCLUDE}) target_link_libraries(cct ${PROJ_LIBRARIES}) -install(TARGETS cct - RUNTIME DESTINATION ${BINDIR}) +install(TARGETS cct + RUNTIME DESTINATION ${BINDIR}) diff --git a/src/bin_gie.cmake b/src/bin_gie.cmake index ca6dde0e..b5f8f8ef 100644 --- a/src/bin_gie.cmake +++ b/src/bin_gie.cmake @@ -1,9 +1,9 @@ -set(GIE_SRC gie.c proj_strtod.c) +set(GIE_SRC gie.c proj_strtod.c proj_strtod.h) set(GIE_INCLUDE optargpm.h) source_group("Source Files\\Bin" FILES ${GIE_SRC}) add_executable(gie ${GIE_SRC} ${GIE_INCLUDE}) target_link_libraries(gie ${PROJ_LIBRARIES}) -install(TARGETS gie - RUNTIME DESTINATION ${BINDIR}) +install(TARGETS gie + RUNTIME DESTINATION ${BINDIR}) @@ -80,13 +80,11 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-10-26 #include "proj.h" #include "proj_internal.h" +#include "proj_strtod.h" #include "projects.h" #include "optargpm.h" -/* Prototypes for functions in proj_strtod.c */ -double proj_strtod(const char *str, char **endptr); -double proj_atof(const char *str); static void logger(void *data, int level, const char *msg); static void print(PJ_LOG_LEVEL log_level, const char *fmt, ...); @@ -116,6 +116,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-10-01/2017-10-08 #include "proj.h" #include "proj_internal.h" #include "proj_math.h" +#include "proj_strtod.h" #include "projects.h" #include "optargpm.h" @@ -154,10 +155,6 @@ static const char *gie_tags[] = { static const size_t n_gie_tags = sizeof gie_tags / sizeof gie_tags[0]; -/* from proj_strtod.c */ -double proj_strtod(const char *str, char **endptr); -double proj_atof(const char *str); - int main(int argc, char **argv); static int dispatch (const char *cmnd, const char *args); @@ -705,6 +702,20 @@ static int roundtrip (const char *args) { /***************************************************************************** Check how far we go from the ACCEPTed point when doing successive back/forward transformation pairs. + +Without args, roundtrip defaults to 100 iterations: + + roundtrip + +With one arg, roundtrip will default to a tolerance of T.tolerance: + + roundtrip ntrips + +With two args: + + roundtrip ntrips tolerance + +Always returns 0. ******************************************************************************/ int ntrips; double d, r, ans; @@ -719,7 +730,17 @@ back/forward transformation pairs. } ans = proj_strtod (args, &endp); - ntrips = (int) (endp==args? 100: fabs(ans)); + if (endp==args) { + /* Default to 100 iterations if not args. */ + ntrips = 100; + } else { + if (ans < 1.0 || ans > 1000000.0) { + errmsg (2, "Invalid number of roundtrips: %lf\n", ans); + return another_failing_roundtrip (); + } + ntrips = (int)ans; + } + d = strtod_scaled (endp, 1); d = d==HUGE_VAL? T.tolerance: d; diff --git a/src/proj_strtod.c b/src/proj_strtod.c index fa683465..ad197d2a 100644 --- a/src/proj_strtod.c +++ b/src/proj_strtod.c @@ -85,6 +85,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-01-17/2017-09-18 ***********************************************************************/ +#include "proj_strtod.h" #include <stdlib.h> /* for abs */ #include <string.h> /* for strchr */ @@ -93,9 +94,6 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-01-17/2017-09-18 #include <float.h> /* for HUGE_VAL */ #include <math.h> /* for pow() */ -double proj_strtod(const char *str, char **endptr); -double proj_atof(const char *str); - double proj_strtod(const char *str, char **endptr) { double number = 0, integral_part = 0; @@ -321,8 +319,6 @@ double proj_atof(const char *str) { /* compile/run: gcc -DTEST -o proj_strtod_test proj_strtod.c && proj_strtod_test */ #include <stdio.h> -#include <stdlib.h> -#include <string.h> char *un_underscore (char *s) { static char u[1024]; diff --git a/src/proj_strtod.h b/src/proj_strtod.h new file mode 100644 index 00000000..38c2d1f4 --- /dev/null +++ b/src/proj_strtod.h @@ -0,0 +1,4 @@ +/* Internal header for proj_strtod.c */ + +double proj_strtod(const char *str, char **endptr); +double proj_atof(const char *str); diff --git a/test/gie/deformation.gie b/test/gie/deformation.gie index b6ca3e0f..03071d1e 100644 --- a/test/gie/deformation.gie +++ b/test/gie/deformation.gie @@ -62,4 +62,84 @@ 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 +roundtrip 100 + +accept 12 56 0.0 2011.0 +expect 12 56 0.0 2011.0 +roundtrip 100 + +accept 12 56 0.0 2019.0 +expect 12 56 0.0 2019.0 +roundtrip 100 + +accept 12 56 0.0 +expect 12 56 -36.5966 +roundtrip 100 + + +------------------------------------------------------------------------------- +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 +roundtrip 100 + +accept 12 56 0.0 2011.0 +expect 12 56 0.0 2011.0 +roundtrip 1000 + +accept 12 56 0.0 3011.0 +expect 12 56 0.0 3011.0 +roundtrip 100 + + +------------------------------------------------------------------------------- +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> |