From 56a754ccd83ff41d1c8edee5fe3548d96f4daaf6 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Mon, 27 Mar 2017 11:39:38 +0200 Subject: [PATCH] Horisontal and vertical gridshift drivers Until now gridshifts has not been working with the new API in proj.h since parsing of +nadgrids and +geoidgrids is build into pj_transform(). This commit introduces the possibility to do both horizontal and vertical gridshift with the pipeline API. The vgridshift and hgridshift kernels are simple wrappers for pj_apply_gridshift3() and pj_apply_vgridshift() that are also used by pj_transform(). Introduced in PR #492. --- src/Makefile.am | 3 +- src/PJ_hgridshift.c | 145 +++++++++++++++++++++++++++++++++++++++++++ src/PJ_vgridshift.c | 147 ++++++++++++++++++++++++++++++++++++++++++++ src/lib_proj.cmake | 2 + src/makefile.vc | 3 +- src/pj_list.h | 2 + 6 files changed, 300 insertions(+), 2 deletions(-) create mode 100644 src/PJ_hgridshift.c create mode 100644 src/PJ_vgridshift.c diff --git a/src/Makefile.am b/src/Makefile.am index a2508fe8fc..f9dee2a9ef 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -74,7 +74,8 @@ libproj_la_SOURCES = \ jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c geodesic.c \ pj_strtod.c \ \ - pj_obs_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c + pj_obs_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c \ + PJ_vgridshift.c PJ_hgridshift.c install-exec-local: rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT) diff --git a/src/PJ_hgridshift.c b/src/PJ_hgridshift.c new file mode 100644 index 0000000000..1b9d272701 --- /dev/null +++ b/src/PJ_hgridshift.c @@ -0,0 +1,145 @@ +#define PJ_LIB__ +#include +#include + +PROJ_HEAD(hgridshift, "Horizontal grid shift"); + +static void *freeup_msg (PJ *P, int errlev) { + if (0==P) + return 0; + + if (0!=P->ctx) + pj_ctx_set_errno (P->ctx, errlev); + + return pj_dealloc(P); +} + + +static void freeup (PJ *P) { + freeup_msg (P, 0); + return; +} + + +static XYZ forward_3d(LPZ lpz, PJ *P) { + PJ_TRIPLET point; + point.lpz = lpz; + + if (P->gridlist != NULL) { + /* Only try the gridshift if at least one grid is loaded, + * otherwise just pass the coordinate through unchanged. */ + pj_apply_gridshift_3( P->ctx, P->gridlist, + P->gridlist_count, 1, 1, 0, + &point.xyz.x, &point.xyz.y, &point.xyz.z ); + } + + return point.xyz; +} + + +static LPZ reverse_3d(XYZ xyz, PJ *P) { + PJ_TRIPLET point; + point.xyz = xyz; + + if (P->gridlist != NULL) { + /* Only try the gridshift if at least one grid is loaded, + * otherwise just pass the coordinate through unchanged. */ + pj_apply_gridshift_3( P->ctx, P->gridlist, + P->gridlist_count, 0, 1, 0, + &point.xyz.x, &point.xyz.y, &point.xyz.z ); + } + + return point.lpz; +} + + +static PJ_OBS forward_obs(PJ_OBS obs, PJ *P) { + PJ_OBS point; + point.coo.xyz = forward_3d (obs.coo.lpz, P); + return point; +} + + +static PJ_OBS reverse_obs(PJ_OBS obs, PJ *P) { + PJ_OBS point; + point.coo.lpz = reverse_3d (obs.coo.xyz, P); + return point; +} + + +PJ *PROJECTION(hgridshift) { + + if (!pj_param(P->ctx, P->params, "tgrids").i) { + pj_log_error(P, "hgridshift: +grids parameter missing."); + return freeup_msg(P, -1); + } + + /* Build gridlist. P->gridlist can be empty if +grids only ask for optional grids. */ + P->gridlist = pj_gridlist_from_nadgrids( P->ctx, pj_param(P->ctx, P->params, "sgrids").s, + &(P->gridlist_count) ); + + /* Was gridlist compiled properly? */ + if ( pj_ctx_get_errno(P->ctx) ) { + pj_log_error(P, "hgridshift: could not find required grid(s)."); + return freeup_msg(P, -38); + } + + 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_RADIANS; + P->right = PJ_IO_UNITS_RADIANS; + + return P; +} + + +#ifndef PJ_SELFTEST +/* selftest stub */ +int pj_hgridshift_selftest (void) {return 0;} +#else +int pj_hgridshift_selftest (void) { + PJ *P; + PJ_OBS expect, a, b; + double dist; + + /* fail on purpose: +grids parameter is mandatory*/ + P = pj_create("+proj=hgridshift"); + if (0!=P) + return 99; + + /* fail on purpose: open non-existing grid */ + P = pj_create("+proj=hgridshift +grids=nonexistinggrid.gsb"); + if (0!=P) + return 999; + + + /* Failure most likely means the grid is missing */ + P = pj_create ("+proj=hgridshift +grids=nzgd2kgrid0005.gsb +ellps=GRS80"); + if (0==P) + return 10; + + a = pj_obs_null; + a.coo.lpz.lam = TORAD(173); + a.coo.lpz.phi = TORAD(-45); + + dist = pj_roundtrip (P, PJ_FWD, 1, a); + if (dist > 0.00000001) + return 1; + + expect.coo.lpz.lam = TORAD(172.999892181021551); + expect.coo.lpz.phi = TORAD(-45.001620431954613); + b = pj_trans(P, PJ_FWD, a); + if (pj_xy_dist(expect.coo.xy, b.coo.xy) > 1e-4) + return 2; + + + pj_free(P); + + return 0; +} +#endif diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c new file mode 100644 index 0000000000..c91db968a7 --- /dev/null +++ b/src/PJ_vgridshift.c @@ -0,0 +1,147 @@ +#define PJ_LIB__ +#include +#include + +PROJ_HEAD(vgridshift, "Vertical grid shift"); + +static void *freeup_msg (PJ *P, int errlev) { + if (0==P) + return 0; + + if (0!=P->ctx) + pj_ctx_set_errno (P->ctx, errlev); + + return pj_dealloc(P); +} + + +static void freeup (PJ *P) { + freeup_msg (P, 0); + return; +} + + +static XYZ forward_3d(LPZ lpz, PJ *P) { + PJ_TRIPLET point; + point.lpz = lpz; + + if (P->gridlist != NULL) { + /* Only try the gridshift if at least one grid is loaded, + * otherwise just pass the coordinate through unchanged. */ + pj_apply_vgridshift( P, "sgrids", + &(P->vgridlist_geoid), + &(P->vgridlist_geoid_count), + 1, 1, 0, + &point.xyz.x, &point.xyz.y, &point.xyz.z ); + } + + return point.xyz; +} + + +static LPZ reverse_3d(XYZ xyz, PJ *P) { + PJ_TRIPLET point; + point.xyz = xyz; + + if (P->gridlist != NULL) { + /* Only try the gridshift if at least one grid is loaded, + * otherwise just pass the coordinate through unchanged. */ + pj_apply_vgridshift( P, "sgrids", + &(P->vgridlist_geoid), + &(P->vgridlist_geoid_count), + 0, 1, 0, + &point.xyz.x, &point.xyz.y, &point.xyz.z ); + } + + return point.lpz; +} + + +static PJ_OBS forward_obs(PJ_OBS obs, PJ *P) { + PJ_OBS point; + point.coo.xyz = forward_3d (obs.coo.lpz, P); + return point; +} + +static PJ_OBS reverse_obs(PJ_OBS obs, PJ *P) { + PJ_OBS point; + point.coo.lpz = reverse_3d (obs.coo.xyz, P); + return point; +} + + +PJ *PROJECTION(vgridshift) { + + if (!pj_param(P->ctx, P->params, "tgrids").i) { + pj_log_error(P, "vgridshift: +grids parameter missing."); + return freeup_msg(P, -1); + } + + /* Build gridlist. P->gridlist can be empty if +grids only ask for optional grids. */ + P->gridlist = pj_gridlist_from_nadgrids( P->ctx, pj_param(P->ctx, P->params, "sgrids").s, + &(P->gridlist_count) ); + + /* Was gridlist compiled properly? */ + if ( pj_ctx_get_errno(P->ctx) ) { + pj_log_error(P, "vgridshift: could not find required grid(s)."); + return freeup_msg(P, -38); + } + + 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_RADIANS; + P->right = PJ_IO_UNITS_RADIANS; + + return P; +} + + +#ifndef PJ_SELFTEST +/* selftest stub */ +int pj_vgridshift_selftest (void) {return 0;} +#else +int pj_vgridshift_selftest (void) { + PJ *P; + PJ_OBS expect, a, b; + double dist; + + /* fail on purpose: +grids parameter it mandatory*/ + P = pj_create("+proj=vgridshift"); + if (0!=P) + return 99; + + /* fail on purpose: open non-existing grid */ + P = pj_create("+proj=vgridshift +grids=nonexistinggrid.gtx"); + if (0!=P) + return 999; + + /* Failure most likely means the grid is missing */ + P = pj_create ("+proj=vgridshift +grids=egm96_15.gtx +ellps=GRS80"); + if (0==P) + return 10; + + a = pj_obs_null; + a.coo.lpz.lam = TORAD(12.5); + a.coo.lpz.phi = TORAD(55.5); + + dist = pj_roundtrip (P, PJ_FWD, 1, a); + if (dist > 0.00000001) + return 1; + + expect = a; + expect.coo.lpz.z = -36.021305084228515625; + b = pj_trans(P, PJ_FWD, a); + if (pj_xyz_dist(expect.coo.xyz, b.coo.xyz) > 1e-10) + return 2; + + + pj_free(P); + + return 0; +} +#endif diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 2fcc4e3162..d1f75cb9cf 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -78,6 +78,7 @@ SET(SRC_LIBPROJ_PJ PJ_hammer.c PJ_hatano.c PJ_helmert.c + PJ_hgridshift.c PJ_horner.c PJ_igh.c PJ_isea.c @@ -139,6 +140,7 @@ SET(SRC_LIBPROJ_PJ PJ_vandg.c PJ_vandg2.c PJ_vandg4.c + PJ_vgridshift.c PJ_wag2.c PJ_wag3.c PJ_wag7.c diff --git a/src/makefile.vc b/src/makefile.vc index e1eacf044f..b0d8291d0f 100644 --- a/src/makefile.vc +++ b/src/makefile.vc @@ -59,7 +59,8 @@ support = \ pj_strtod.obj pj_run_selftests.obj pj_generic_selftest.obj pipeline = \ - pj_obs_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj + pj_obs_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \ + PJ_vgridshift.obj PJ_hgridshift.obj geodesic = geodesic.obj diff --git a/src/pj_list.h b/src/pj_list.h index 96db13d71a..7ed6bdacef 100644 --- a/src/pj_list.h +++ b/src/pj_list.h @@ -53,6 +53,7 @@ PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area") PROJ_HEAD(healpix, "HEALPix") PROJ_HEAD(rhealpix, "rHEALPix") PROJ_HEAD(helmert, "3- and 7-parameter Helmert shift") +PROJ_HEAD(hgridshift, "Horizontal grid shift") PROJ_HEAD(horner, "Horner polynomial evaluation") PROJ_HEAD(igh, "Interrupted Goode Homolosine") PROJ_HEAD(imw_p, "International Map of the World Polyconic") @@ -141,6 +142,7 @@ PROJ_HEAD(vandg2, "van der Grinten II") PROJ_HEAD(vandg3, "van der Grinten III") PROJ_HEAD(vandg4, "van der Grinten IV") PROJ_HEAD(vitk1, "Vitkovsky I") +PROJ_HEAD(vgridshift, "Vertical grid shift") PROJ_HEAD(wag1, "Wagner I (Kavraisky VI)") PROJ_HEAD(wag2, "Wagner II") PROJ_HEAD(wag3, "Wagner III")