aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/Makefile.am2
-rw-r--r--src/PJ_deformation.c291
-rw-r--r--src/lib_proj.cmake1
-rw-r--r--src/makefile.vc3
-rw-r--r--src/pj_ell_set.c39
-rw-r--r--src/pj_list.h1
-rw-r--r--src/proj_internal.h1
-rw-r--r--test/gie/deformation.gie73
-rwxr-xr-xtravis/install.sh3
9 files changed, 410 insertions, 4 deletions
diff --git a/src/Makefile.am b/src/Makefile.am
index 363f2cce..0bfcfb24 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -87,7 +87,7 @@ libproj_la_SOURCES = \
\
proj_4D_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c \
PJ_vgridshift.c PJ_hgridshift.c PJ_unitconvert.c PJ_molodensky.c \
- pj_internal.c
+ PJ_deformation.c pj_internal.c
install-exec-local:
rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT)
diff --git a/src/PJ_deformation.c b/src/PJ_deformation.c
new file mode 100644
index 00000000..05858900
--- /dev/null
+++ b/src/PJ_deformation.c
@@ -0,0 +1,291 @@
+/***********************************************************************
+
+ Kinematic datum shifting utilizing a deformation model
+
+ Kristian Evers, 2017-10-29
+
+************************************************************************
+
+Perform datum shifts by means of a deformation/velocity model.
+
+ X_out = X_in + (T_ct - T_obs)*DX
+ Y_out = Y_in + (T_ct - T_obs)*DY
+ Z_out = Z_in + (T_ct - T_obs)*DZ
+
+
+Corrections are done in cartesian space. Coordinates of the
+gridded model are in lat/long space and the corrections in the model are
+in cartesian space. Hence the input coordinates needs to be converted
+to lat/long space when searching for corrections in the grid. The
+corrections are then applied to the input coordinates in cartesian
+space.
+
+A full deformation model is described by two grids, one for the horizontal
+components and one for the vertical component. The horizontal grid is
+stored in CTable/CTable2 and the vertical grid is stored in the GTX
+format. Both grids are expected to contain grid-values in units of
+m/year.
+
+************************************************************************
+* Copyright (c) 2017, Kristian Evers
+*
+* Permission is hereby granted, free of charge, to any person obtaining a
+* copy of this software and associated documentation files (the "Software"),
+* to deal in the Software without restriction, including without limitation
+* the rights to use, copy, modify, merge, publish, distribute, sublicense,
+* and/or sell copies of the Software, and to permit persons to whom the
+* Software is furnished to do so, subject to the following conditions:
+*
+* The above copyright notice and this permission notice shall be included
+* in all copies or substantial portions of the Software.
+*
+* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+* DEALINGS IN THE SOFTWARE.
+*
+***********************************************************************/
+#define PJ_LIB__
+#include <errno.h>
+#include <proj.h>
+#include "proj_internal.h"
+#include "projects.h"
+
+PROJ_HEAD(deformation, "Kinematic grid shift");
+
+#define TOL 1e-8
+#define MAX_ITERATIONS 10
+
+struct pj_opaque {
+ double t_obs;
+ double t_epoch;
+ int has_xy_grids, has_z_grids;
+ PJ *cart;
+};
+
+static XYZ get_xygrid_shift(PJ* P, XYZ cartesian) {
+ PJ_COORD geodetic, shift;
+
+ geodetic.lpz = pj_inv3d(cartesian, P->opaque->cart);
+ shift.lp = proj_hgrid_value(P, geodetic.lp);
+
+ return shift.xyz;
+}
+
+static double get_zgrid_shift(PJ* P, XYZ cartesian) {
+ PJ_COORD geodetic;
+
+ geodetic.lpz = pj_inv3d(cartesian, P->opaque->cart);
+
+ return proj_vgrid_value(P, geodetic.lp);
+}
+
+static XYZ reverse_hshift(PJ *P, XYZ input, double dt) {
+ XYZ out, delta, dif;
+ int i = MAX_ITERATIONS;
+
+ delta = get_xygrid_shift(P, input);
+
+ out.x = input.x + dt*delta.x;
+ out.y = input.y - dt*delta.y;
+ out.z = input.z;
+
+ do {
+ delta = get_xygrid_shift(P, out);
+
+ if (delta.x == HUGE_VAL)
+ break;
+
+ dif.x = out.x - dt*delta.x - input.x;
+ dif.y = out.y + dt*delta.y - input.y;
+ out.x += dif.x;
+ out.y += dif.y;
+
+ } while ( --i && hypot(dif.x, dif.y) > TOL );
+
+ return out;
+}
+
+static XYZ forward_3d(LPZ lpz, PJ *P) {
+ struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+ PJ_COORD out, in;
+ XYZ shift;
+ double dt = 0.0;
+ in.lpz = lpz;
+ out = in;
+
+ if (Q->t_obs != HUGE_VAL) {
+ dt = Q->t_obs - Q->t_epoch;
+ } else {
+ out = proj_coord_error(); /* in the 3D case +t_obs must be specified */
+ proj_log_debug(P, "deformation: +t_obs must be specified");
+ return out.xyz;
+ }
+
+ if (Q->has_xy_grids) {
+ shift = get_xygrid_shift(P, in.xyz);
+
+ out.xyz.x += dt * shift.x;
+ out.xyz.y += dt * shift.y;
+ }
+
+ if (Q->has_z_grids)
+ out.xyz.z += dt * get_zgrid_shift(P, in.xyz);
+
+ return out.xyz;
+}
+
+
+static PJ_OBS forward_obs(PJ_OBS in, PJ *P) {
+ struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+ double dt;
+ XYZ shift;
+ PJ_OBS out = in;
+
+ if (Q->t_obs != HUGE_VAL) {
+ dt = Q->t_obs - Q->t_epoch;
+ } else {
+ dt = in.coo.xyzt.t - Q->t_epoch;
+ }
+
+ if (Q->has_xy_grids) {
+ shift = get_xygrid_shift(P, in.coo.xyz);
+ out.coo.xyz.x += dt * shift.x;
+ out.coo.xyz.y += dt * shift.y;
+ }
+
+ if (Q->has_z_grids)
+ out.coo.xyz.z += dt * get_zgrid_shift(P, in.coo.xyz);
+
+ return out;
+}
+
+
+static LPZ reverse_3d(XYZ in, PJ *P) {
+ struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+ PJ_COORD out;
+ double dt = 0.0;
+ out.xyz = in;
+
+ if (Q->t_obs != HUGE_VAL) {
+ dt = Q->t_epoch - Q->t_obs;
+ } else {
+ out = proj_coord_error(); /* in the 3D case +t_obs must be specified */
+ proj_log_debug(P, "deformation: +t_obs must be specified");
+ return out.lpz;
+ }
+
+ if (Q->has_xy_grids) {
+ out.xyz = reverse_hshift(P, in, dt);
+ }
+
+ if (Q->has_z_grids)
+ out.xyz.z = in.z + dt * get_zgrid_shift(P, in);
+
+ return out.lpz;
+}
+
+static PJ_OBS reverse_obs(PJ_OBS in, PJ *P) {
+ struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+ PJ_OBS out = in;
+ double dt;
+
+ if (Q->t_obs != HUGE_VAL) {
+ dt = Q->t_epoch - Q->t_obs;
+ } else {
+ dt = Q->t_epoch - in.coo.xyzt.t;
+ }
+
+ if (Q->has_xy_grids)
+ out.coo.xyz = reverse_hshift(P, in.coo.xyz, dt);
+
+ if (Q->has_z_grids)
+ out.coo.xyz.z = in.coo.xyz.z + dt * get_zgrid_shift(P, in.coo.xyz);
+
+ return out;
+}
+
+static void *destructor(PJ *P, int errlev) {
+ if (0==P)
+ return 0;
+
+ if (0==P->opaque)
+ return pj_default_destructor (P, errlev);
+
+ if (P->opaque->cart)
+ P->opaque->cart->destructor (P->opaque->cart, errlev);
+
+ return pj_default_destructor(P, errlev);
+}
+
+
+PJ *PROJECTION(deformation) {
+ struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+ if (0==Q)
+ return pj_default_destructor(P, ENOMEM);
+ P->opaque = (void *) Q;
+
+ Q->cart = proj_create(P->ctx, "+proj=cart");
+ if (Q->cart == 0)
+ return destructor(P, ENOMEM);
+
+ /* inherit ellipsoid definition from P to Q->cart (simpler than guessing */
+ /* how the ellipsoid was specified in original definition) */
+ pj_inherit_ellipsoid_defs(P, Q->cart);
+
+ Q->has_xy_grids = pj_param(P->ctx, P->params, "txy_grids").i;
+ Q->has_z_grids = pj_param(P->ctx, P->params, "tz_grids").i;
+
+ /* Build gridlists. P->gridlist and P->vgridlist_geoid can be empty if */
+ /* +xy_grids or +z_grids only ask for optional grids */
+ if (!Q->has_xy_grids && !Q->has_z_grids) {
+ proj_log_error(P, "deformation: At least one of either +xy_grids or +z_grids should be specified.");
+ return destructor(P, PJD_ERR_NO_ARGS );
+ }
+
+ if (Q->has_xy_grids) {
+ Q->has_xy_grids = proj_hgrid_init(P, "xy_grids");
+ if (proj_errno(P)) {
+ proj_log_error(P, "deformation: could not find requested grid(s).");
+ return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+ }
+ }
+
+ if (Q->has_z_grids) {
+ Q->has_z_grids = proj_vgrid_init(P, "z_grids");
+ if (proj_errno(P)) {
+ proj_log_error(P, "deformation: could not find requested grid(s).");
+ return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
+ }
+ }
+
+ Q->t_obs = HUGE_VAL;
+ if (pj_param(P->ctx, P->params, "tt_obs").i) {
+ Q->t_obs = pj_param(P->ctx, P->params, "dt_obs").f;
+ }
+
+ if (pj_param(P->ctx, P->params, "tt_epoch").i) {
+ Q->t_epoch = pj_param(P->ctx, P->params, "dt_epoch").f;
+ } else {
+ proj_log_error(P, "deformation: +t_epoch parameter missing.");
+ return destructor(P, PJD_ERR_MISSING_ARGS);
+ }
+
+ P->fwdobs = forward_obs;
+ P->invobs = reverse_obs;
+ P->fwd3d = forward_3d;
+ P->inv3d = reverse_3d;
+ P->fwd = 0;
+ P->inv = 0;
+
+ P->left = PJ_IO_UNITS_METERS;
+ P->right = PJ_IO_UNITS_METERS;
+ P->destructor = destructor;
+
+ return P;
+}
+
+int pj_deformation_selftest (void) {return 10000;}
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 053e9ef6..969ad64b 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -58,6 +58,7 @@ SET(SRC_LIBPROJ_PJ
PJ_collg.c
PJ_comill.c
PJ_crast.c
+ PJ_deformation.c
PJ_denoy.c
PJ_eck1.c
PJ_eck2.c
diff --git a/src/makefile.vc b/src/makefile.vc
index 1330e9bb..15fd6428 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -61,7 +61,8 @@ support = \
pipeline = \
proj_4D_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \
- PJ_vgridshift.obj PJ_hgridshift.obj PJ_unitconvert.obj PJ_molodensky.obj
+ PJ_vgridshift.obj PJ_hgridshift.obj PJ_unitconvert.obj PJ_molodensky.obj \
+ PJ_deformation.obj
geodesic = geodesic.obj
diff --git a/src/pj_ell_set.c b/src/pj_ell_set.c
index 72892ac2..d43feecb 100644
--- a/src/pj_ell_set.c
+++ b/src/pj_ell_set.c
@@ -1,12 +1,49 @@
/* set ellipsoid parameters a and es */
-#include <projects.h>
#include <string.h>
+#include "proj_internal.h"
+#include "projects.h"
#define SIXTH .1666666666666666667 /* 1/6 */
#define RA4 .04722222222222222222 /* 17/360 */
#define RA6 .02215608465608465608 /* 67/3024 */
#define RV4 .06944444444444444444 /* 5/72 */
#define RV6 .04243827160493827160 /* 55/1296 */
+/* copy ellipsoidal parameters from src to dst */
+void pj_inherit_ellipsoid_defs(const PJ *src, PJ *dst) {
+
+ /* The linear parameters */
+ dst->a = src->a;
+ dst->b = src->b;
+ dst->ra = src->ra;
+ dst->rb = src->rb;
+
+ /* The eccentricities */
+ dst->alpha = src->alpha;
+ dst->e = src->e;
+ dst->es = src->es;
+ dst->e2 = src->e2;
+ dst->e2s = src->e2s;
+ dst->e3 = src->e3;
+ dst->e3s = src->e3s;
+ dst->one_es = src->one_es;
+ dst->rone_es = src->rone_es;
+
+ /* The flattenings */
+ dst->f = src->f;
+ dst->f2 = src->f2;
+ dst->n = src->n;
+ dst->rf = src->rf;
+ dst->rf2 = src->rf2;
+ dst->rn = src->rn;
+
+ /* This one's for GRS80 */
+ dst->J = src->J;
+
+ /* es and a before any +proj related adjustment */
+ dst->es_orig = src->es_orig;
+ dst->a_orig = src->a_orig;
+}
+
/* initialize geographic shape parameters */
int pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
int i;
diff --git a/src/pj_list.h b/src/pj_list.h
index bf287219..76c37126 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -25,6 +25,7 @@ PROJ_HEAD(chamb, "Chamberlin Trimetric")
PROJ_HEAD(collg, "Collignon")
PROJ_HEAD(comill, "Compact Miller")
PROJ_HEAD(crast, "Craster Parabolic (Putnins P4)")
+PROJ_HEAD(deformation, "Kinematic grid shift")
PROJ_HEAD(denoy, "Denoyer Semi-Elliptical")
PROJ_HEAD(eck1, "Eckert I")
PROJ_HEAD(eck2, "Eckert II")
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 29e73e70..ba542dea 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -117,6 +117,7 @@ void proj_log_trace (PJ *P, const char *fmt, ...);
/*void proj_log_func (PJ *P, void *app_data, void (*log)(void *, int, const char *));*/
void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION log);
+void pj_inherit_ellipsoid_defs(const PJ *src, PJ *dst);
/* Lowest level: Minimum support for fileapi */
void proj_fileapi_set (PJ *P, void *fileapi);
diff --git a/test/gie/deformation.gie b/test/gie/deformation.gie
new file mode 100644
index 00000000..74a6b25d
--- /dev/null
+++ b/test/gie/deformation.gie
@@ -0,0 +1,73 @@
+BEGIN
+===============================================================================
+Test for the deformation operation - Kinematic Gridshifting
+
+For all the deformation tests the alaska and egm96_15.gtx grids are used even
+though they are not parts of a deformation model, they are in the proper format
+and for testing purposes it doesn't really matter all that much...
+
+The input coordinate is located at lon=60, lam=-160 - somewhere in Alaska.
+
+===============================================================================
+
+
+-------------------------------------------------------------------------------
+Test using only horizontal grid and +tobs parameter
+-------------------------------------------------------------------------------
+OPERATION +proj=deformation +xy_grids=alaska +t_epoch=2016.0 +t_obs=2000.0
+-------------------------------------------------------------------------------
+TOLERANCE 0.000001 m
+ACCEPT -3004295.5882503074 -1093474.1690603832 5500477.1338251457
+EXPECT -3004295.5888766116 -1093474.1688513425 5500477.1338251457
+ROUNDTRIP 1000
+
+-------------------------------------------------------------------------------
+Test using only vertical grid and +tobs parameter
+-------------------------------------------------------------------------------
+OPERATION +proj=deformation +z_grids=egm96_15.gtx +t_epoch=2016.0 +t_obs=2000.0
+-------------------------------------------------------------------------------
+TOLERANCE 0.000001 m
+ACCEPT -3004295.5882503074 -1093474.1690603832 5500477.1338251457
+EXPECT -3004295.5882503074 -1093474.1690603832 5500234.008855661
+ROUNDTRIP 1000
+
+-------------------------------------------------------------------------------
+Test using both horizontal and vertical grids as well as the +tobs parameter
+-------------------------------------------------------------------------------
+OPERATION +proj=deformation +xy_grids=alaska +z_grids=egm96_15.gtx +t_epoch=2016.0 +t_obs=2000.0
+-------------------------------------------------------------------------------
+TOLERANCE 0.000001 m
+ACCEPT -3004295.5882503074 -1093474.1690603832 5500477.1338251457
+EXPECT -3004295.5888766116 -1093474.1688513425 5500234.008855661
+ROUNDTRIP 1000
+
+-------------------------------------------------------------------------------
+Test using only horizontal grid
+-------------------------------------------------------------------------------
+OPERATION +proj=deformation +xy_grids=alaska +t_epoch=2016.0
+-------------------------------------------------------------------------------
+TOLERANCE 0.000001 m
+ACCEPT -3004295.5882503074 -1093474.1690603832 5500477.1338251457 2000.0
+EXPECT -3004295.5888766116 -1093474.1688513425 5500477.1338251457 2000.0
+ROUNDTRIP 1000
+
+-------------------------------------------------------------------------------
+Test using only vertical grid
+-------------------------------------------------------------------------------
+OPERATION +proj=deformation +z_grids=egm96_15.gtx +t_epoch=2016.0
+-------------------------------------------------------------------------------
+TOLERANCE 0.000001 m
+ACCEPT -3004295.5882503074 -1093474.1690603832 5500477.1338251457 2000.0
+EXPECT -3004295.5882503074 -1093474.1690603832 5500234.008855661 2000.0
+ROUNDTRIP 1000
+
+-------------------------------------------------------------------------------
+Test using both horizontal and vertical grids
+-------------------------------------------------------------------------------
+OPERATION +proj=deformation +xy_grids=alaska +z_grids=egm96_15.gtx +t_epoch=2016.0
+-------------------------------------------------------------------------------
+TOLERANCE 0.000001 m
+ACCEPT -3004295.5882503074 -1093474.1690603832 5500477.1338251457 2000.0
+EXPECT -3004295.5888766116 -1093474.1688513425 5500234.008855661 2000.0
+ROUNDTRIP 1000
+END
diff --git a/travis/install.sh b/travis/install.sh
index 1247435b..62c2677f 100755
--- a/travis/install.sh
+++ b/travis/install.sh
@@ -81,7 +81,8 @@ if [ $TRAVIS_OS_NAME == "osx" ]; then
make -j3
make check
PROJ_LIB=$GRIDDIR ./src/proj -VC
-./src/gie ./test/gie/builtins.gie
+PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/builtins.gie
+PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/deformation.gie
# install & run the working GIGS test
# create locations that pyproj understands
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-16 02:35:44 +0000