From b8f765def1c54ddd0e8c61fd4619f58aa35165ff Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Tue, 24 Oct 2017 18:34:40 +0200 Subject: Extend proj_strtod test case collection and improve its strtod-replication --- src/proj_strtod.c | 66 ++++++++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 53 insertions(+), 13 deletions(-) (limited to 'src') diff --git a/src/proj_strtod.c b/src/proj_strtod.c index d4063b0b..757dfaf6 100644 --- a/src/proj_strtod.c +++ b/src/proj_strtod.c @@ -121,10 +121,16 @@ double proj_strtod(const char *str, char **endptr) { /* Empty string? */ if (0==*p) { - errno = EINVAL; if (endptr) - *endptr = p; - return HUGE_VAL; + *endptr = (char *) str; + return 0; + } + + /* non-numeric? */ + if (0==strchr("0123456789+-._", *p)) { + if (endptr) + *endptr = (char *) str; + return 0; } /* Then handle optional prefixed sign and skip prefix zeros */ @@ -137,9 +143,15 @@ double proj_strtod(const char *str, char **endptr) { if (isdigit(*p) || '_'==*p || '.'==*p) break; if (endptr) - *endptr = p; - errno = EINVAL; - return HUGE_VAL; + *endptr = (char *) str; + return 0; + } + + /* stray sign, as in "+/-"? */ + if (0!=sign && (0==strchr ("0123456789._", *p) || 0==*p)) { + if (endptr) + *endptr = (char *) str; + return 0; } /* skip prefixed zeros before '.' */ @@ -147,8 +159,11 @@ double proj_strtod(const char *str, char **endptr) { p++; /* zero? */ - if (0==*p || 0==strchr ("0123456789eE.", *p)) - return 0; + if ((0==*p) || 0==strchr ("0123456789eE.", *p) || isspace(*p)) { + if (endptr) + *endptr = p; + return sign==-1? -0: 0; + } /* Now expect a (potentially zero-length) string of digits */ while (isdigit(*p) || ('_'==*p)) { @@ -228,8 +243,15 @@ double proj_strtod(const char *str, char **endptr) { number = -number; /* Do we have an exponent part? */ - if (*p == 'e' || *p == 'E') { + while (*p == 'e' || *p == 'E') { p++; + + /* Just a stray "e", as in 100elephants? */ + if (0==*p || 0==strchr ("0123456789+-_", *p)) { + p--; + break; + } + while ('_'==*p) p++; /* Does it have a sign? */ @@ -263,6 +285,7 @@ double proj_strtod(const char *str, char **endptr) { if (-1==sign) n = -n; exponent += n; + break; } if (endptr) @@ -351,14 +374,31 @@ int main (int argc, char **argv) { errno = 0; - test ("1"); - test ("0"); + test (""); + test (" "); + test (" abcde"); + test (" edcba"); + test ("abcde"); + test ("edcba"); + test ("+"); + test ("-"); + test ("+ "); + test ("- "); + test (" + "); + test (" - "); + test ("e 1"); + test ("e1"); + test ("0 66"); test ("1."); test ("0."); test ("1.0"); test ("0.0"); - test ("1.0"); - test ("0.0"); + test ("1 "); + test ("0 "); + test ("-0 "); + test ("0_ "); + test ("0_"); + test ("1e"); test ("_1.0"); test ("_0.0"); test ("1_.0"); -- cgit v1.2.3 From 395ec20b5a661e77d048212182841e5684dd6757 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Tue, 24 Oct 2017 23:36:13 +0200 Subject: Make gie roundtrips compatible with updated proj_strtod In order to mimic strtod, proj_strtod now returns 0 and not HUGE_VAL on nonnumeric input. Hence, we must check the return pointers to identify an error. --- src/gie.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/gie.c b/src/gie.c index 6ff157ac..33fc4d82 100644 --- a/src/gie.c +++ b/src/gie.c @@ -457,11 +457,11 @@ static int accept (char *args) { static int roundtrip (char *args) { int ntrips; double d, r, ans; - char *endp; + char *endp, *endq; ans = proj_strtod (args, &endp); ntrips = (int) (ans==HUGE_VAL? 100: fabs(ans)); - d = proj_strtod (endp, &endp); - d = d==HUGE_VAL? T.tolerance: d / 1000; + d = proj_strtod (endp, &endq); + d = (endp==endq)? T.tolerance: d / 1000; r = proj_roundtrip (T.P, PJ_FWD, ntrips, T.a); if (r > d) { if (T.verbosity > -1) { -- cgit v1.2.3 From fd95842b9fafb0e140bc867af71f358c92258ff8 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Wed, 25 Oct 2017 10:39:56 +0200 Subject: Remove PJ_OBS from the API surface, rename pj_obs_api.c to pj_4D_api.c (#625) * Remove PJ_OBS from the API surface, rename pj_obs_api.c to pj_4D_api.c * Repair proj.def --- src/Makefile.am | 2 +- src/PJ_cart.c | 10 +- src/PJ_helmert.c | 4 +- src/PJ_latlong.c | 3 +- src/lib_proj.cmake | 2 +- src/makefile.vc | 2 +- src/pj_internal.c | 55 +++- src/pj_obs_api.c | 818 ---------------------------------------------------- src/proj.def | 75 +++-- src/proj.h | 19 +- src/proj_4D_api.c | 752 +++++++++++++++++++++++++++++++++++++++++++++++ src/proj_internal.h | 9 + src/projects.h | 6 +- 13 files changed, 863 insertions(+), 894 deletions(-) delete mode 100644 src/pj_obs_api.c create mode 100644 src/proj_4D_api.c (limited to 'src') diff --git a/src/Makefile.am b/src/Makefile.am index a039151d..363f2cce 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -85,7 +85,7 @@ 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 \ + 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 diff --git a/src/PJ_cart.c b/src/PJ_cart.c index 37aa3b97..99e1af7d 100644 --- a/src/PJ_cart.c +++ b/src/PJ_cart.c @@ -387,10 +387,10 @@ int pj_cart_selftest (void) { b = proj_trans_obs (P, PJ_FWD, obs[1]); n = proj_transform ( - P, PJ_FWD, + P, PJ_FWD, &(obs[0].coo.lpz.lam), sz, 2, &(obs[0].coo.lpz.phi), sz, 2, - &(obs[0].coo.lpz.z), sz, 2, + &(obs[0].coo.lpz.z), sz, 2, 0, sz, 0 ); if (2!=n) @@ -408,10 +408,10 @@ int pj_cart_selftest (void) { h = 27; t = 33; n = proj_transform ( - P, PJ_FWD, + P, PJ_FWD, &(obs[0].coo.lpz.lam), sz, 2, &(obs[0].coo.lpz.phi), sz, 2, - &h, 0, 1, + &h, 0, 1, &t, 0, 1 ); if (2!=n) @@ -429,7 +429,7 @@ int pj_cart_selftest (void) { obs[0].coo = proj_coord (PJ_TORAD(12), PJ_TORAD(55), 45, 0); obs[1].coo = proj_coord (PJ_TORAD(12), PJ_TORAD(56), 50, 0); - if (proj_transform_obs(P, PJ_FWD, 2, obs)) + if (proj_transform_coord(P, PJ_FWD, 2, (PJ_COORD *) obs)) return 30; if (a.coo.lpz.lam != obs[0].coo.lpz.lam) return 31; diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c index ffbdd01a..503dc392 100644 --- a/src/PJ_helmert.c +++ b/src/PJ_helmert.c @@ -668,11 +668,11 @@ int pj_helmert_selftest (void) { matrix is updated when necessary. Test coordinates from GNSStrans. */ XYZ expect4a = {3370658.18890, 711877.42370, 5349787.12430}; XYZ expect4b = {3370658.18087, 711877.42750, 5349787.12648}; - PJ_OBS in4 = {{{3370658.378, 711877.314, 5349787.086, 2017.0}}, {{ 0, 0, 0}}, 0, 0}; + PJ_OBS in4 = {{{3370658.378, 711877.314, 5349787.086, 2017.0}}}; PJ_OBS out; PJ *helmert = proj_create( - 0, + 0, " +proj=helmert +ellps=GRS80" " +x=0.0127 +y=0.0065 +z=-0.0209 +s=0.00195" " +rx=-0.00039 +ry=0.00080 +rz=-0.00114" diff --git a/src/PJ_latlong.c b/src/PJ_latlong.c index 1677142a..7ee41e2a 100644 --- a/src/PJ_latlong.c +++ b/src/PJ_latlong.c @@ -1,6 +1,6 @@ /****************************************************************************** * Project: PROJ.4 - * Purpose: Stub projection implementation for lat/long coordinates. We + * Purpose: Stub projection implementation for lat/long coordinates. We * don't actually change the coordinates, but we want proj=latlong * to act sort of like a projection. * Author: Frank Warmerdam, warmerdam@pobox.com @@ -29,6 +29,7 @@ /* very loosely based upon DMA code by Bradford W. Drew */ #define PJ_LIB__ +#include "proj_internal.h" #include #include "projects.h" diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 1be10362..053e9ef6 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -201,7 +201,7 @@ SET(SRC_LIBPROJ_CORE pj_mlfn.c pj_msfn.c pj_mutex.c - pj_obs_api.c + proj_4D_api.c pj_internal.c proj_internal.h pj_open_lib.c diff --git a/src/makefile.vc b/src/makefile.vc index fdf03bd3..1330e9bb 100644 --- a/src/makefile.vc +++ b/src/makefile.vc @@ -60,7 +60,7 @@ support = \ pj_internal.obj pipeline = \ - pj_obs_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \ + 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 geodesic = geodesic.obj diff --git a/src/pj_internal.c b/src/pj_internal.c index 31c299ac..5eb98afb 100644 --- a/src/pj_internal.c +++ b/src/pj_internal.c @@ -1,8 +1,9 @@ /****************************************************************************** * Project: PROJ.4 - * Purpose: This is primarily material originating from pj_obs_api.c, - * that does not fit into the API category. Hence this pile of - * tubings and fittings for PROJ.4 internal plumbing. + * Purpose: This is primarily material originating from pj_obs_api.c + * (now proj_4D_api.c), that does not fit into the API + * category. Hence this pile of tubings and fittings for + * PROJ.4 internal plumbing. * * Author: Thomas Knudsen, thokn@sdfe.dk, 2017-07-05 * @@ -42,12 +43,52 @@ /* Used for zero-initializing new objects */ const PJ_COORD proj_coord_null = {{0, 0, 0, 0}}; const PJ_OBS proj_obs_null = { - {{0, 0, 0, 0}}, - {{0, 0, 0}}, - 0, 0 + {{0, 0, 0, 0}} }; + +/* Initialize PJ_OBS struct */ +PJ_OBS proj_obs (double x, double y, double z, double t) { + PJ_OBS res; + res.coo = proj_coord (x, y, z, t); + return res; +} + + + + + + + + + + + + + + +/* Apply the transformation P to the coordinate coo */ +PJ_OBS proj_trans_obs (PJ *P, PJ_DIRECTION direction, PJ_OBS obs) { + if (0==P) + return obs; + + switch (direction) { + case PJ_FWD: + return pj_fwdobs (obs, P); + case PJ_INV: + return pj_invobs (obs, P); + case PJ_IDENT: + return obs; + default: + break; + } + + proj_errno_set (P, EINVAL); + return proj_obs_error (); +} + + /* Work around non-constness of MSVC HUGE_VAL by providing functions rather than constants */ PJ_COORD proj_coord_error (void) { PJ_COORD c; @@ -58,8 +99,6 @@ PJ_COORD proj_coord_error (void) { PJ_OBS proj_obs_error (void) { PJ_OBS obs; obs.coo = proj_coord_error (); - obs.anc.v[0] = obs.anc.v[1] = obs.anc.v[2] = HUGE_VAL; - obs.id = obs.flags = 0; return obs; } diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c deleted file mode 100644 index 55da9fa2..00000000 --- a/src/pj_obs_api.c +++ /dev/null @@ -1,818 +0,0 @@ -/****************************************************************************** - * Project: PROJ.4 - * Purpose: Implement a (currently minimalistic) proj API based primarily - * on the PJ_OBS generic geodetic data type. - * - * proj thread contexts have not seen widespread use, so one of the - * intentions with this new API is to make them less visible on the - * API surface: Contexts do not have a life by themselves, they are - * visible only through their associated PJs, and the number of - * functions supporting them is limited. - * - * Author: Thomas Knudsen, thokn@sdfe.dk, 2016-06-09/2016-11-06 - * - ****************************************************************************** - * Copyright (c) 2016, 2017 Thomas Knudsen/SDFE - * - * 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_OBS_API_C -#include -#include "proj_internal.h" -#include "projects.h" -#include "geodesic.h" -#include -#include - - -/* Initialize PJ_COORD struct */ -PJ_COORD proj_coord (double x, double y, double z, double t) { - PJ_COORD res; - res.v[0] = x; - res.v[1] = y; - res.v[2] = z; - res.v[3] = t; - return res; -} - -/* Initialize PJ_OBS struct */ -PJ_OBS proj_obs (double x, double y, double z, double t, double o, double p, double k, int id, unsigned int flags) { - PJ_OBS res; - res.coo.v[0] = x; - res.coo.v[1] = y; - res.coo.v[2] = z; - res.coo.v[3] = t; - res.anc.v[0] = o; - res.anc.v[1] = p; - res.anc.v[2] = k; - res.id = id; - res.flags = flags; - - return res; -} - - - -/* Geodesic distance (in meter) between two points with angular 2D coordinates */ -double proj_lp_dist (const PJ *P, LP a, LP b) { - double s12, azi1, azi2; - /* Note: the geodesic code takes arguments in degrees */ - geod_inverse (P->geod, PJ_TODEG(a.phi), PJ_TODEG(a.lam), PJ_TODEG(b.phi), PJ_TODEG(b.lam), &s12, &azi1, &azi2); - return s12; -} - -/* Euclidean distance between two points with linear 2D coordinates */ -double proj_xy_dist (XY a, XY b) { - return hypot (a.x - b.x, a.y - b.y); -} - -/* Euclidean distance between two points with linear 3D coordinates */ -double proj_xyz_dist (XYZ a, XYZ b) { - return hypot (hypot (a.x - b.x, a.y - b.y), a.z - b.z); -} - - - -/* Measure numerical deviation after n roundtrips fwd-inv (or inv-fwd) */ -double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD coo) { - int i; - PJ_COORD o, u; - enum pj_io_units unit; - - if (0==P) - return HUGE_VAL; - - if (n < 1) { - proj_errno_set (P, EINVAL); - return HUGE_VAL; - } - - o = coo; - - switch (direction) { - case PJ_FWD: - for (i = 0; i < n; i++) { - u = pj_fwdcoord (o, P); - o = pj_invcoord (u, P); - } - break; - case PJ_INV: - for (i = 0; i < n; i++) { - u = pj_invcoord (o, P); - o = pj_fwdcoord (u, P); - } - break; - default: - proj_errno_set (P, EINVAL); - return HUGE_VAL; - } - - /* left when forward, because we do a roundtrip, and end where we begin */ - unit = direction==PJ_FWD? P->left: P->right; - if (unit==PJ_IO_UNITS_RADIANS) - return hypot (proj_lp_dist (P, coo.lp, o.lp), coo.lpz.z - o.lpz.z); - - return proj_xyz_dist (coo.xyz, coo.xyz); -} - - - - - - - - - - - - - -/* Apply the transformation P to the coordinate coo */ -PJ_OBS proj_trans_obs (PJ *P, PJ_DIRECTION direction, PJ_OBS obs) { - if (0==P) - return obs; - - switch (direction) { - case PJ_FWD: - return pj_fwdobs (obs, P); - case PJ_INV: - return pj_invobs (obs, P); - case PJ_IDENT: - return obs; - default: - break; - } - - proj_errno_set (P, EINVAL); - return proj_obs_error (); -} - - - -/* Apply the transformation P to the coordinate coo */ -PJ_COORD proj_trans_coord (PJ *P, PJ_DIRECTION direction, PJ_COORD coo) { - if (0==P) - return coo; - - switch (direction) { - case PJ_FWD: - return pj_fwdcoord (coo, P); - case PJ_INV: - return pj_invcoord (coo, P); - case PJ_IDENT: - return coo; - default: - break; - } - - proj_errno_set (P, EINVAL); - return proj_coord_error (); -} - - - -/*************************************************************************************/ -size_t proj_transform ( - PJ *P, - PJ_DIRECTION direction, - double *x, size_t sx, size_t nx, - double *y, size_t sy, size_t ny, - double *z, size_t sz, size_t nz, - double *t, size_t st, size_t nt -) { -/************************************************************************************** - - Transform a series of coordinates, where the individual coordinate dimension - may be represented by an array that is either - - 1. fully populated - 2. a null pointer and/or a length of zero, which will be treated as a - fully populated array of zeroes - 3. of length one, i.e. a constant, which will be treated as a fully - populated array of that constant value - - The strides, sx, sy, sz, st, represent the step length, in bytes, between - consecutive elements of the corresponding array. This makes it possible for - proj_transform to handle transformation of a large class of application - specific data structures, without necessarily understanding the data structure - format, as in: - - typedef struct {double x, y; int quality_level; char surveyor_name[134];} XYQS; - XYQS survey[345]; - double height = 23.45; - PJ *P = {...}; - size_t stride = sizeof (XYQS); - ... - proj_transform ( - P, PJ_INV, sizeof(XYQS), - &(survey[0].x), stride, 345, (* We have 345 eastings *) - &(survey[0].y), stride, 345, (* ...and 345 northings. *) - &height, 1, (* The height is the constant 23.45 m *) - 0, 0 (* and the time is the constant 0.00 s *) - ); - - This is similar to the inner workings of the pj_transform function, but the - stride functionality has been generalized to work for any size of basic unit, - not just a fixed number of doubles. - - In most cases, the stride will be identical for x, y,z, and t, since they will - typically be either individual arrays (stride = sizeof(double)), or strided - views into an array of application specific data structures (stride = sizeof (...)). - - But in order to support cases where x, y, z, and t come from heterogeneous - sources, individual strides, sx, sy, sz, st, are used. - - Caveat: Since proj_transform does its work *in place*, this means that even the - supposedly constants (i.e. length 1 arrays) will return from the call in altered - state. Hence, remember to reinitialize between repeated calls. - - Return value: Number of transformations completed. - -**************************************************************************************/ - PJ_COORD coord = proj_coord_null; - size_t i, nmin; - double null_broadcast = 0; - if (0==P) - return 0; - - /* ignore lengths of null arrays */ - if (0==x) nx = 0; - if (0==y) ny = 0; - if (0==z) nz = 0; - if (0==t) nt = 0; - - /* and make the nullities point to some real world memory for broadcasting nulls */ - if (0==nx) x = &null_broadcast; - if (0==ny) y = &null_broadcast; - if (0==nz) z = &null_broadcast; - if (0==nt) t = &null_broadcast; - - /* nothing to do? */ - if (0==nx+ny+nz+nt) - return 0; - - /* arrays of length 1 are constants, which we broadcast along the longer arrays */ - /* so we need to find the length of the shortest non-unity array to figure out */ - /* how many coordinate pairs we must transform */ - nmin = (nx > 1)? nx: (ny > 1)? ny: (nz > 1)? nz: (nt > 1)? nt: 1; - if ((nx > 1) && (nx < nmin)) nmin = nx; - if ((ny > 1) && (ny < nmin)) nmin = ny; - if ((nz > 1) && (nz < nmin)) nmin = nz; - if ((nt > 1) && (nt < nmin)) nmin = nt; - - /* Check validity of direction flag */ - switch (direction) { - case PJ_FWD: - case PJ_INV: - break; - case PJ_IDENT: - return nmin; - default: - proj_errno_set (P, EINVAL); - return 0; - } - - /* Arrays of length==0 are broadcast as the constant 0 */ - /* Arrays of length==1 are broadcast as their single value */ - /* Arrays of length >1 are iterated over (for the first nmin values) */ - /* The slightly convolved incremental indexing is used due */ - /* to the stride, which may be any size supported by the platform */ - for (i = 0; i < nmin; i++) { - coord.xyzt.x = *x; - coord.xyzt.y = *y; - coord.xyzt.z = *z; - coord.xyzt.t = *t; - - if (PJ_FWD==direction) - coord = pj_fwdcoord (coord, P); - else - coord = pj_invcoord (coord, P); - - /* in all full length cases, we overwrite the input with the output */ - if (nx > 1) { - *x = coord.xyzt.x; - x = (double *) ( ((char *) x) + sx); - } - if (ny > 1) { - *y = coord.xyzt.y; - y = (double *) ( ((char *) y) + sy); - } - if (nz > 1) { - *z = coord.xyzt.z; - z = (double *) ( ((char *) z) + sz); - } - if (nt > 1) { - *t = coord.xyzt.t; - t = (double *) ( ((char *) t) + st); - } - } - /* Last time around, we update the length 1 cases with their transformed alter egos */ - /* ... or would we rather not? Then what about the nmin==1 case? */ - /* perhaps signalling the non-array case by setting all strides to 0? */ - if (nx==1) - *x = coord.xyzt.x; - if (ny==1) - *y = coord.xyzt.y; - if (nz==1) - *z = coord.xyzt.z; - if (nt==1) - *t = coord.xyzt.t; - - return i; -} - -/*****************************************************************************/ -int proj_transform_obs (PJ *P, PJ_DIRECTION direction, size_t n, PJ_OBS *obs) { -/****************************************************************************** - Batch transform an array of PJ_OBS. - - Returns 0 if all observations are transformed without error, otherwise - returns error number. -******************************************************************************/ - size_t i; - for (i=0; i= 511) - continue; - - if (strlen(info.searchpath) != 0) { - strcat(info.searchpath, ";"); - len = strlen(paths[i]); - strncat(info.searchpath, paths[i], sizeof(info.searchpath)-len-1); - } else { - pj_strlcpy(info.searchpath, paths[i], sizeof(info.searchpath)); - } - } - - return info; -} - - -/*****************************************************************************/ -PJ_PROJ_INFO proj_pj_info(PJ *P) { -/****************************************************************************** - Basic info about a particular instance of a projection object. - - Returns PJ_PROJ_INFO struct. - -******************************************************************************/ - PJ_PROJ_INFO info; - char *def; - - memset(&info, 0, sizeof(PJ_PROJ_INFO)); - - /* Expected accuracy of the transformation. Hardcoded for now, will be improved */ - /* later. Most likely to be used when a transformation is set up with */ - /* proj_create_crs_to_crs in a future version that leverages the EPSG database. */ - info.accuracy = -1.0; - - if (!P) { - return info; - } - - /* projection id */ - if (pj_param(P->ctx, P->params, "tproj").i) - pj_strlcpy(info.id, pj_param(P->ctx, P->params, "sproj").s, sizeof(info.id)); - - /* projection description */ - pj_strlcpy(info.description, P->descr, sizeof(info.description)); - - /* projection definition */ - def = pj_get_def(P, 0); /* pj_get_def takes a non-const PJ pointer */ - pj_strlcpy(info.definition, &def[1], sizeof(info.definition)); /* def includes a leading space */ - pj_dealloc(def); - - /* this does not take into account that a pipeline potentially does not */ - /* have an inverse. */ - info.has_inverse = (P->inv != 0 || P->inv3d != 0 || P->invobs != 0); - - return info; -} - - -/*****************************************************************************/ -PJ_GRID_INFO proj_grid_info(const char *gridname) { -/****************************************************************************** - Information about a named datum grid. - - Returns PJ_GRID_INFO struct. - -******************************************************************************/ - PJ_GRID_INFO info; - - /*PJ_CONTEXT *ctx = proj_context_create(); */ - PJ_CONTEXT *ctx = pj_get_default_ctx(); - PJ_GRIDINFO *gridinfo = pj_gridinfo_init(ctx, gridname); - memset(&info, 0, sizeof(PJ_GRID_INFO)); - - /* in case the grid wasn't found */ - if (gridinfo->filename == NULL) { - pj_gridinfo_free(ctx, gridinfo); - strcpy(info.format, "missing"); - return info; - } - - /* name of grid */ - pj_strlcpy(info.gridname, gridname, sizeof(info.gridname)); - - /* full path of grid */ - pj_find_file(ctx, gridname, info.filename, sizeof(info.filename)); - - /* grid format */ - pj_strlcpy(info.format, gridinfo->format, sizeof(info.format)); - - /* grid size */ - info.n_lon = gridinfo->ct->lim.lam; - info.n_lat = gridinfo->ct->lim.phi; - - /* cell size */ - info.cs_lon = gridinfo->ct->del.lam; - info.cs_lat = gridinfo->ct->del.phi; - - /* bounds of grid */ - info.lowerleft = gridinfo->ct->ll; - info.upperright.lam = info.lowerleft.lam + info.n_lon*info.cs_lon; - info.upperright.phi = info.lowerleft.phi + info.n_lat*info.cs_lat; - - pj_gridinfo_free(ctx, gridinfo); - - return info; -} - -/*****************************************************************************/ -PJ_INIT_INFO proj_init_info(const char *initname){ -/****************************************************************************** - Information about a named init file. - - Maximum length of initname is 64. - - Returns PJ_INIT_INFO struct. - - If the init file is not found all members of - the return struct are set to 0. If the init file is found, but it the - metadata is missing, the value is set to "Unknown". - -******************************************************************************/ - int file_found, def_found=0; - char param[80], key[74]; - paralist *start, *next; - PJ_INIT_INFO info; - PJ_CONTEXT *ctx = pj_get_default_ctx(); - - memset(&info, 0, sizeof(PJ_INIT_INFO)); - - file_found = pj_find_file(ctx, initname, info.filename, sizeof(info.filename)); - if (!file_found || strlen(initname) > 64) { - return info; - } - - pj_strlcpy(info.name, initname, sizeof(info.name)); - strcpy(info.origin, "Unknown"); - strcpy(info.version, "Unknown"); - strcpy(info.lastupdate, "Unknown"); - - pj_strlcpy(key, initname, 64); /* make room for ":metadata\0" at the end */ - strncat(key, ":metadata", 9); - strcpy(param, "+init="); - strncat(param, key, 73); - - start = pj_mkparam(param); - next = pj_get_init(ctx, &start, start, key, &def_found); - - if (pj_param(ctx, start, "tversion").i) { - pj_strlcpy(info.version, pj_param(ctx, start, "sversion").s, sizeof(info.version)); - } - - if (pj_param(ctx, start, "torigin").i) { - pj_strlcpy(info.origin, pj_param(ctx, start, "sorigin").s, sizeof(info.origin)); - } - - if (pj_param(ctx, start, "tlastupdate").i) { - pj_strlcpy(info.lastupdate, pj_param(ctx, start, "slastupdate").s, sizeof(info.lastupdate)); - } - - for ( ; start; start = next) { - next = start->next; - pj_dalloc(start); - } - - return info; -} - - -/*****************************************************************************/ -PJ_DERIVS proj_derivatives(PJ *P, const LP lp) { -/****************************************************************************** - Derivatives of coordinates. - - returns PJ_DERIVS. If unsuccessfull error number is set and the returned - struct contains NULL data. - -******************************************************************************/ - PJ_DERIVS derivs; - - if (pj_deriv(lp, 1e-5, P, &derivs)) { - /* errno set in pj_derivs */ - memset(&derivs, 0, sizeof(PJ_DERIVS)); - } - - return derivs; -} - - -/*****************************************************************************/ -PJ_FACTORS proj_factors(PJ *P, const LP lp) { -/****************************************************************************** - Cartographic characteristics at point lp. - - Characteristics include meridian, parallel and areal scales, angular - distortion, meridian/parallel, meridian convergence and scale error. - - returns PJ_FACTORS. If unsuccessfull error number is set and the returned - struct contains NULL data. - -******************************************************************************/ - PJ_FACTORS factors; - - /* pj_factors rely code being zero */ - factors.code = 0; - - if (pj_factors(lp, P, 0.0, &factors)) { - /* errno set in pj_factors */ - memset(&factors, 0, sizeof(PJ_FACTORS)); - } - - return factors; -} - - -const PJ_ELLPS *proj_list_ellps(void) { - return pj_get_ellps_ref(); -} - -const PJ_UNITS *proj_list_units(void) { - return pj_get_units_ref(); -} - -const PJ_OPERATIONS *proj_list_operations(void) { - return pj_get_list_ref(); -} - -const PJ_PRIME_MERIDIANS *proj_list_prime_meridians(void) { - return pj_get_prime_meridians_ref(); -} - -double proj_torad (double angle_in_degrees) { return PJ_TORAD (angle_in_degrees);} -double proj_todeg (double angle_in_radians) { return PJ_TODEG (angle_in_radians);} - - -double proj_dmstor(const char *is, char **rs) { - return dmstor(is, rs); -} - -char* proj_rtodms(char *s, double r, int pos, int neg) { - return rtodms(s, r, pos, neg); -} diff --git a/src/proj.def b/src/proj.def index 598f2824..86deb24d 100644 --- a/src/proj.def +++ b/src/proj.def @@ -100,50 +100,49 @@ EXPORTS proj_trans_obs @95 proj_trans_coord @96 proj_transform @97 - proj_transform_obs @98 - proj_transform_coord @99 - proj_roundtrip @100 + proj_transform_coord @98 + proj_roundtrip @99 - proj_coord @101 - proj_obs @102 - proj_coord_error @103 - proj_obs_error @104 + proj_coord @100 + proj_obs @101 + proj_coord_error @102 + proj_obs_error @103 - proj_errno @105 - proj_errno_set @106 - proj_errno_reset @107 - proj_errno_restore @108 - proj_context_errno_set @109 + proj_errno @104 + proj_errno_set @105 + proj_errno_reset @106 + proj_errno_restore @107 + proj_context_errno_set @108 - proj_context_create @110 - proj_context_set @111 - proj_context_inherit @112 - proj_context_destroy @113 + proj_context_create @109 + proj_context_set @110 + proj_context_inherit @111 + proj_context_destroy @112 - proj_lp_dist @114 - proj_xy_dist @115 - proj_xyz_dist @116 + proj_lp_dist @113 + proj_xy_dist @114 + proj_xyz_dist @115 - proj_log_level @117 - proj_log_func @118 - proj_log_error @119 - proj_log_debug @120 - proj_log_trace @121 + proj_log_level @116 + proj_log_func @117 + proj_log_error @118 + proj_log_debug @119 + proj_log_trace @120 - proj_info @122 - proj_pj_info @123 - proj_grid_info @124 - proj_init_info @125 + proj_info @121 + proj_pj_info @122 + proj_grid_info @123 + proj_init_info @124 - proj_torad @126 - proj_todeg @127 - proj_rtodms @128 - proj_dmstor @129 + proj_torad @125 + proj_todeg @126 + proj_rtodms @127 + proj_dmstor @128 - proj_derivatives @130 - proj_factors @131 + proj_derivatives @129 + proj_factors @130 - proj_list_operations @132 - proj_list_ellps @133 - proj_list_units @134 - proj_list_prime_meridians @135 + proj_list_operations @131 + proj_list_ellps @132 + proj_list_units @133 + proj_list_prime_meridians @134 diff --git a/src/proj.h b/src/proj.h index aa9f3b17..ed885091 100644 --- a/src/proj.h +++ b/src/proj.h @@ -90,7 +90,7 @@ * compatible call proj_context_create(0), which will not create * a new context, but simply provide a pointer to the default one. * - * See pj_obs_api_test.c for an example of how to use the API. + * See proj_4D_api_test.c for examples of how to use the API. * * Author: Thomas Knudsen, * Benefitting from a large number of comments and suggestions @@ -170,10 +170,6 @@ extern char const pj_release[]; /* global release id string */ /* first forward declare everything needed */ -/* Data type for generic geodetic observations */ -struct PJ_OBS; -typedef struct PJ_OBS PJ_OBS; - /* Data type for generic geodetic 3D data */ union PJ_TRIPLET; typedef union PJ_TRIPLET PJ_TRIPLET; @@ -299,14 +295,6 @@ union PJ_PAIR { }; -struct PJ_OBS { - PJ_COORD coo; /* coordinate data */ - PJ_TRIPLET anc; /* ancillary data */ - int id; /* integer ancillary data - e.g. observation number, EPSG code... */ - unsigned int flags; /* additional data, intended for flags */ -}; - - #define PJ_IS_ANAL_XL_YL 01 /* derivatives of lon analytic */ #define PJ_IS_ANAL_XP_YP 02 /* derivatives of lat analytic */ #define PJ_IS_ANAL_HK 04 /* h and k analytic */ @@ -385,7 +373,6 @@ enum PJ_DIRECTION { }; typedef enum PJ_DIRECTION PJ_DIRECTION; -PJ_OBS proj_trans_obs (PJ *P, PJ_DIRECTION direction, PJ_OBS obs); PJ_COORD proj_trans_coord (PJ *P, PJ_DIRECTION direction, PJ_COORD coord); @@ -398,16 +385,14 @@ size_t proj_transform ( double *t, size_t st, size_t nt ); -int proj_transform_obs (PJ *P, PJ_DIRECTION direction, size_t n, PJ_OBS *obs); int proj_transform_coord (PJ *P, PJ_DIRECTION direction, size_t n, PJ_COORD *coord); /* Initializers */ PJ_COORD proj_coord (double x, double y, double z, double t); -PJ_OBS proj_obs (double x, double y, double z, double t, double o, double p, double k, int id, unsigned int flags); /* Measure internal consistency - in forward or inverse direction */ double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD coo); - + /* Geodesic distance between two points with angular 2D coordinates */ double proj_lp_dist (const PJ *P, LP a, LP b); diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c new file mode 100644 index 00000000..b1aa3883 --- /dev/null +++ b/src/proj_4D_api.c @@ -0,0 +1,752 @@ +/****************************************************************************** + * Project: PROJ.4 + * Purpose: Implement a (currently minimalistic) proj API based primarily + * on the PJ_COORD 4D geodetic spatiotemporal data type. + * + * proj thread contexts have not seen widespread use, so one of the + * intentions with this new API is to make them less visible on the + * API surface: Contexts do not have a life by themselves, they are + * visible only through their associated PJs, and the number of + * functions supporting them is limited. + * + * Author: Thomas Knudsen, thokn@sdfe.dk, 2016-06-09/2016-11-06 + * + ****************************************************************************** + * Copyright (c) 2016, 2017 Thomas Knudsen/SDFE + * + * 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. + *****************************************************************************/ +#include +#include "proj_internal.h" +#include "projects.h" +#include "geodesic.h" +#include +#include + + +/* Initialize PJ_COORD struct */ +PJ_COORD proj_coord (double x, double y, double z, double t) { + PJ_COORD res; + res.v[0] = x; + res.v[1] = y; + res.v[2] = z; + res.v[3] = t; + return res; +} + + + +/* Geodesic distance (in meter) between two points with angular 2D coordinates */ +double proj_lp_dist (const PJ *P, LP a, LP b) { + double s12, azi1, azi2; + /* Note: the geodesic code takes arguments in degrees */ + geod_inverse (P->geod, PJ_TODEG(a.phi), PJ_TODEG(a.lam), PJ_TODEG(b.phi), PJ_TODEG(b.lam), &s12, &azi1, &azi2); + return s12; +} + +/* Euclidean distance between two points with linear 2D coordinates */ +double proj_xy_dist (XY a, XY b) { + return hypot (a.x - b.x, a.y - b.y); +} + +/* Euclidean distance between two points with linear 3D coordinates */ +double proj_xyz_dist (XYZ a, XYZ b) { + return hypot (hypot (a.x - b.x, a.y - b.y), a.z - b.z); +} + + + +/* Measure numerical deviation after n roundtrips fwd-inv (or inv-fwd) */ +double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD coo) { + int i; + PJ_COORD o, u; + enum pj_io_units unit; + + if (0==P) + return HUGE_VAL; + + if (n < 1) { + proj_errno_set (P, EINVAL); + return HUGE_VAL; + } + + o = coo; + + switch (direction) { + case PJ_FWD: + for (i = 0; i < n; i++) { + u = pj_fwdcoord (o, P); + o = pj_invcoord (u, P); + } + break; + case PJ_INV: + for (i = 0; i < n; i++) { + u = pj_invcoord (o, P); + o = pj_fwdcoord (u, P); + } + break; + default: + proj_errno_set (P, EINVAL); + return HUGE_VAL; + } + + /* left when forward, because we do a roundtrip, and end where we begin */ + unit = direction==PJ_FWD? P->left: P->right; + if (unit==PJ_IO_UNITS_RADIANS) + return hypot (proj_lp_dist (P, coo.lp, o.lp), coo.lpz.z - o.lpz.z); + + return proj_xyz_dist (coo.xyz, coo.xyz); +} + + + +/* Apply the transformation P to the coordinate coo */ +PJ_COORD proj_trans_coord (PJ *P, PJ_DIRECTION direction, PJ_COORD coo) { + if (0==P) + return coo; + + switch (direction) { + case PJ_FWD: + return pj_fwdcoord (coo, P); + case PJ_INV: + return pj_invcoord (coo, P); + case PJ_IDENT: + return coo; + default: + break; + } + + proj_errno_set (P, EINVAL); + return proj_coord_error (); +} + + + +/*************************************************************************************/ +size_t proj_transform ( + PJ *P, + PJ_DIRECTION direction, + double *x, size_t sx, size_t nx, + double *y, size_t sy, size_t ny, + double *z, size_t sz, size_t nz, + double *t, size_t st, size_t nt +) { +/************************************************************************************** + + Transform a series of coordinates, where the individual coordinate dimension + may be represented by an array that is either + + 1. fully populated + 2. a null pointer and/or a length of zero, which will be treated as a + fully populated array of zeroes + 3. of length one, i.e. a constant, which will be treated as a fully + populated array of that constant value + + The strides, sx, sy, sz, st, represent the step length, in bytes, between + consecutive elements of the corresponding array. This makes it possible for + proj_transform to handle transformation of a large class of application + specific data structures, without necessarily understanding the data structure + format, as in: + + typedef struct {double x, y; int quality_level; char surveyor_name[134];} XYQS; + XYQS survey[345]; + double height = 23.45; + PJ *P = {...}; + size_t stride = sizeof (XYQS); + ... + proj_transform ( + P, PJ_INV, sizeof(XYQS), + &(survey[0].x), stride, 345, (* We have 345 eastings *) + &(survey[0].y), stride, 345, (* ...and 345 northings. *) + &height, 1, (* The height is the constant 23.45 m *) + 0, 0 (* and the time is the constant 0.00 s *) + ); + + This is similar to the inner workings of the pj_transform function, but the + stride functionality has been generalized to work for any size of basic unit, + not just a fixed number of doubles. + + In most cases, the stride will be identical for x, y,z, and t, since they will + typically be either individual arrays (stride = sizeof(double)), or strided + views into an array of application specific data structures (stride = sizeof (...)). + + But in order to support cases where x, y, z, and t come from heterogeneous + sources, individual strides, sx, sy, sz, st, are used. + + Caveat: Since proj_transform does its work *in place*, this means that even the + supposedly constants (i.e. length 1 arrays) will return from the call in altered + state. Hence, remember to reinitialize between repeated calls. + + Return value: Number of transformations completed. + +**************************************************************************************/ + PJ_COORD coord = proj_coord_null; + size_t i, nmin; + double null_broadcast = 0; + if (0==P) + return 0; + + /* ignore lengths of null arrays */ + if (0==x) nx = 0; + if (0==y) ny = 0; + if (0==z) nz = 0; + if (0==t) nt = 0; + + /* and make the nullities point to some real world memory for broadcasting nulls */ + if (0==nx) x = &null_broadcast; + if (0==ny) y = &null_broadcast; + if (0==nz) z = &null_broadcast; + if (0==nt) t = &null_broadcast; + + /* nothing to do? */ + if (0==nx+ny+nz+nt) + return 0; + + /* arrays of length 1 are constants, which we broadcast along the longer arrays */ + /* so we need to find the length of the shortest non-unity array to figure out */ + /* how many coordinate pairs we must transform */ + nmin = (nx > 1)? nx: (ny > 1)? ny: (nz > 1)? nz: (nt > 1)? nt: 1; + if ((nx > 1) && (nx < nmin)) nmin = nx; + if ((ny > 1) && (ny < nmin)) nmin = ny; + if ((nz > 1) && (nz < nmin)) nmin = nz; + if ((nt > 1) && (nt < nmin)) nmin = nt; + + /* Check validity of direction flag */ + switch (direction) { + case PJ_FWD: + case PJ_INV: + break; + case PJ_IDENT: + return nmin; + default: + proj_errno_set (P, EINVAL); + return 0; + } + + /* Arrays of length==0 are broadcast as the constant 0 */ + /* Arrays of length==1 are broadcast as their single value */ + /* Arrays of length >1 are iterated over (for the first nmin values) */ + /* The slightly convolved incremental indexing is used due */ + /* to the stride, which may be any size supported by the platform */ + for (i = 0; i < nmin; i++) { + coord.xyzt.x = *x; + coord.xyzt.y = *y; + coord.xyzt.z = *z; + coord.xyzt.t = *t; + + if (PJ_FWD==direction) + coord = pj_fwdcoord (coord, P); + else + coord = pj_invcoord (coord, P); + + /* in all full length cases, we overwrite the input with the output */ + if (nx > 1) { + *x = coord.xyzt.x; + x = (double *) ( ((char *) x) + sx); + } + if (ny > 1) { + *y = coord.xyzt.y; + y = (double *) ( ((char *) y) + sy); + } + if (nz > 1) { + *z = coord.xyzt.z; + z = (double *) ( ((char *) z) + sz); + } + if (nt > 1) { + *t = coord.xyzt.t; + t = (double *) ( ((char *) t) + st); + } + } + /* Last time around, we update the length 1 cases with their transformed alter egos */ + /* ... or would we rather not? Then what about the nmin==1 case? */ + /* perhaps signalling the non-array case by setting all strides to 0? */ + if (nx==1) + *x = coord.xyzt.x; + if (ny==1) + *y = coord.xyzt.y; + if (nz==1) + *z = coord.xyzt.z; + if (nt==1) + *t = coord.xyzt.t; + + return i; +} + +/*****************************************************************************/ +int proj_transform_coord (PJ *P, PJ_DIRECTION direction, size_t n, PJ_COORD *coord) { +/****************************************************************************** + Batch transform an array of PJ_COORD. + + Returns 0 if all coordinates are transformed without error, otherwise + returns error number. +******************************************************************************/ + size_t i; + for (i=0; i= 511) + continue; + + if (strlen(info.searchpath) != 0) { + strcat(info.searchpath, ";"); + len = strlen(paths[i]); + strncat(info.searchpath, paths[i], sizeof(info.searchpath)-len-1); + } else { + pj_strlcpy(info.searchpath, paths[i], sizeof(info.searchpath)); + } + } + + return info; +} + + +/*****************************************************************************/ +PJ_PROJ_INFO proj_pj_info(PJ *P) { +/****************************************************************************** + Basic info about a particular instance of a projection object. + + Returns PJ_PROJ_INFO struct. + +******************************************************************************/ + PJ_PROJ_INFO info; + char *def; + + memset(&info, 0, sizeof(PJ_PROJ_INFO)); + + /* Expected accuracy of the transformation. Hardcoded for now, will be improved */ + /* later. Most likely to be used when a transformation is set up with */ + /* proj_create_crs_to_crs in a future version that leverages the EPSG database. */ + info.accuracy = -1.0; + + if (!P) { + return info; + } + + /* projection id */ + if (pj_param(P->ctx, P->params, "tproj").i) + pj_strlcpy(info.id, pj_param(P->ctx, P->params, "sproj").s, sizeof(info.id)); + + /* projection description */ + pj_strlcpy(info.description, P->descr, sizeof(info.description)); + + /* projection definition */ + def = pj_get_def(P, 0); /* pj_get_def takes a non-const PJ pointer */ + pj_strlcpy(info.definition, &def[1], sizeof(info.definition)); /* def includes a leading space */ + pj_dealloc(def); + + /* this does not take into account that a pipeline potentially does not */ + /* have an inverse. */ + info.has_inverse = (P->inv != 0 || P->inv3d != 0 || P->invobs != 0); + + return info; +} + + +/*****************************************************************************/ +PJ_GRID_INFO proj_grid_info(const char *gridname) { +/****************************************************************************** + Information about a named datum grid. + + Returns PJ_GRID_INFO struct. + +******************************************************************************/ + PJ_GRID_INFO info; + + /*PJ_CONTEXT *ctx = proj_context_create(); */ + PJ_CONTEXT *ctx = pj_get_default_ctx(); + PJ_GRIDINFO *gridinfo = pj_gridinfo_init(ctx, gridname); + memset(&info, 0, sizeof(PJ_GRID_INFO)); + + /* in case the grid wasn't found */ + if (gridinfo->filename == NULL) { + pj_gridinfo_free(ctx, gridinfo); + strcpy(info.format, "missing"); + return info; + } + + /* name of grid */ + pj_strlcpy(info.gridname, gridname, sizeof(info.gridname)); + + /* full path of grid */ + pj_find_file(ctx, gridname, info.filename, sizeof(info.filename)); + + /* grid format */ + pj_strlcpy(info.format, gridinfo->format, sizeof(info.format)); + + /* grid size */ + info.n_lon = gridinfo->ct->lim.lam; + info.n_lat = gridinfo->ct->lim.phi; + + /* cell size */ + info.cs_lon = gridinfo->ct->del.lam; + info.cs_lat = gridinfo->ct->del.phi; + + /* bounds of grid */ + info.lowerleft = gridinfo->ct->ll; + info.upperright.lam = info.lowerleft.lam + info.n_lon*info.cs_lon; + info.upperright.phi = info.lowerleft.phi + info.n_lat*info.cs_lat; + + pj_gridinfo_free(ctx, gridinfo); + + return info; +} + +/*****************************************************************************/ +PJ_INIT_INFO proj_init_info(const char *initname){ +/****************************************************************************** + Information about a named init file. + + Maximum length of initname is 64. + + Returns PJ_INIT_INFO struct. + + If the init file is not found all members of + the return struct are set to 0. If the init file is found, but it the + metadata is missing, the value is set to "Unknown". + +******************************************************************************/ + int file_found, def_found=0; + char param[80], key[74]; + paralist *start, *next; + PJ_INIT_INFO info; + PJ_CONTEXT *ctx = pj_get_default_ctx(); + + memset(&info, 0, sizeof(PJ_INIT_INFO)); + + file_found = pj_find_file(ctx, initname, info.filename, sizeof(info.filename)); + if (!file_found || strlen(initname) > 64) { + return info; + } + + pj_strlcpy(info.name, initname, sizeof(info.name)); + strcpy(info.origin, "Unknown"); + strcpy(info.version, "Unknown"); + strcpy(info.lastupdate, "Unknown"); + + pj_strlcpy(key, initname, 64); /* make room for ":metadata\0" at the end */ + strncat(key, ":metadata", 9); + strcpy(param, "+init="); + strncat(param, key, 73); + + start = pj_mkparam(param); + next = pj_get_init(ctx, &start, start, key, &def_found); + + if (pj_param(ctx, start, "tversion").i) { + pj_strlcpy(info.version, pj_param(ctx, start, "sversion").s, sizeof(info.version)); + } + + if (pj_param(ctx, start, "torigin").i) { + pj_strlcpy(info.origin, pj_param(ctx, start, "sorigin").s, sizeof(info.origin)); + } + + if (pj_param(ctx, start, "tlastupdate").i) { + pj_strlcpy(info.lastupdate, pj_param(ctx, start, "slastupdate").s, sizeof(info.lastupdate)); + } + + for ( ; start; start = next) { + next = start->next; + pj_dalloc(start); + } + + return info; +} + + +/*****************************************************************************/ +PJ_DERIVS proj_derivatives(PJ *P, const LP lp) { +/****************************************************************************** + Derivatives of coordinates. + + returns PJ_DERIVS. If unsuccessfull error number is set and the returned + struct contains NULL data. + +******************************************************************************/ + PJ_DERIVS derivs; + + if (pj_deriv(lp, 1e-5, P, &derivs)) { + /* errno set in pj_derivs */ + memset(&derivs, 0, sizeof(PJ_DERIVS)); + } + + return derivs; +} + + +/*****************************************************************************/ +PJ_FACTORS proj_factors(PJ *P, const LP lp) { +/****************************************************************************** + Cartographic characteristics at point lp. + + Characteristics include meridian, parallel and areal scales, angular + distortion, meridian/parallel, meridian convergence and scale error. + + returns PJ_FACTORS. If unsuccessfull error number is set and the returned + struct contains NULL data. + +******************************************************************************/ + PJ_FACTORS factors; + + /* pj_factors rely code being zero */ + factors.code = 0; + + if (pj_factors(lp, P, 0.0, &factors)) { + /* errno set in pj_factors */ + memset(&factors, 0, sizeof(PJ_FACTORS)); + } + + return factors; +} + + +const PJ_ELLPS *proj_list_ellps(void) { + return pj_get_ellps_ref(); +} + +const PJ_UNITS *proj_list_units(void) { + return pj_get_units_ref(); +} + +const PJ_OPERATIONS *proj_list_operations(void) { + return pj_get_list_ref(); +} + +const PJ_PRIME_MERIDIANS *proj_list_prime_meridians(void) { + return pj_get_prime_meridians_ref(); +} + +double proj_torad (double angle_in_degrees) { return PJ_TORAD (angle_in_degrees);} +double proj_todeg (double angle_in_radians) { return PJ_TODEG (angle_in_radians);} + + +double proj_dmstor(const char *is, char **rs) { + return dmstor(is, rs); +} + +char* proj_rtodms(char *s, double r, int pos, int neg) { + return rtodms(s, r, pos, neg); +} diff --git a/src/proj_internal.h b/src/proj_internal.h index 0b74b563..5773637c 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -50,6 +50,12 @@ extern "C" { +struct PJ_OBS { + PJ_COORD coo; /* coordinate data */ +}; +typedef struct PJ_OBS PJ_OBS; + + #ifndef PJ_TODEG #define PJ_TODEG(rad) ((rad)*180.0/M_PI) #endif @@ -59,6 +65,9 @@ extern "C" { +PJ_OBS proj_obs (double x, double y, double z, double t); +PJ_OBS proj_trans_obs (PJ *P, PJ_DIRECTION direction, PJ_OBS obs); + PJ_COORD proj_coord_error (void); PJ_OBS proj_obs_error (void); #ifndef PJ_INTERNAL_C diff --git a/src/projects.h b/src/projects.h index cbb980ca..b9d88cf3 100644 --- a/src/projects.h +++ b/src/projects.h @@ -174,6 +174,9 @@ typedef struct { double u, v, w; } UVW; /* Forward declarations and typedefs for stuff needed inside the PJ object */ struct PJconsts; struct PJ_OBS; +#ifndef PROJ_INTERNAL_H +typedef struct PJ_OBS PJ_OBS; +#endif union PJ_COORD; struct geod_geodesic; struct pj_opaque; @@ -189,7 +192,6 @@ enum pj_io_units { }; #ifndef PROJ_H typedef struct PJconsts PJ; /* the PJ object herself */ -typedef struct PJ_OBS PJ_OBS; typedef union PJ_COORD PJ_COORD; #endif @@ -259,7 +261,7 @@ struct PJconsts { void (*spc)(LP, PJ *, struct FACTORS *); void *(*destructor)(PJ *, int); - + /************************************************************************************* -- cgit v1.2.3 From 064efb77179850a552c56831047e8e739aed4067 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Wed, 25 Oct 2017 13:11:48 +0200 Subject: Repair gie and cct after breakage due to proj_strtod update (#628) * Repair gie and cct after breakage due to proj_strtod update * Remove unused variables --- src/cct.c | 20 +++++++++---- src/gie.c | 99 +++++++++++++++++++++++++++++++++++++++++---------------------- 2 files changed, 80 insertions(+), 39 deletions(-) (limited to 'src') diff --git a/src/cct.c b/src/cct.c index 35d4ec3f..72d8c6cc 100644 --- a/src/cct.c +++ b/src/cct.c @@ -289,22 +289,32 @@ char *column (char *buf, int n) { return buf; } +/* column to double */ +static double cold (char *args, int col) { + char *endp; + char *target; + double d; + target = column (args, col); + d = proj_strtod (target, &endp); + if (endp==target) + return HUGE_VAL; + return d; +} PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double fixed_time) { PJ_COORD err = proj_coord (HUGE_VAL, HUGE_VAL, HUGE_VAL, HUGE_VAL); PJ_COORD result = err; int prev_errno = errno; - char *endptr = 0; errno = 0; result.xyzt.z = fixed_height; result.xyzt.t = fixed_time; - result.xyzt.x = proj_strtod (column (buf, columns[0]), &endptr); - result.xyzt.y = proj_strtod (column (buf, columns[1]), &endptr); + result.xyzt.x = cold (buf, columns[0]); + result.xyzt.y = cold (buf, columns[1]); if (result.xyzt.z==HUGE_VAL) - result.xyzt.z = proj_strtod (column (buf, columns[2]), &endptr); + result.xyzt.z = cold (buf, columns[2]); if (result.xyzt.t==HUGE_VAL) - result.xyzt.t = proj_strtod (column (buf, columns[3]), &endptr); + result.xyzt.t = cold (buf, columns[3]); if (0!=errno) return err; diff --git a/src/gie.c b/src/gie.c index 33fc4d82..f0850b61 100644 --- a/src/gie.c +++ b/src/gie.c @@ -132,6 +132,7 @@ static int get_inp (FILE *f, char *inp, int size); static int get_cmnd (char *inp, char *cmnd, int len); static char *get_args (char *inp); static int dispatch (char *cmnd, char *args); +static char *column (char *buf, int n); @@ -314,6 +315,52 @@ static int process_file (char *fname) { } +/* return a pointer to the n'th column of buf */ +char *column (char *buf, int n) { + int i; + if (n <= 0) + return buf; + for (i = 0; i < n; i++) { + while (isspace(*buf)) + buf++; + if (i == n - 1) + break; + while ((0 != *buf) && !isspace(*buf)) + buf++; + } + return buf; +} + + + +static double strtod_scaled (char *args, double default_scale) { + double s; + char *endp = args; + s = proj_strtod (args, &endp); + if (args==endp) + return HUGE_VAL; + + endp = column (args, 2); + + if (0==strcmp(endp, "km")) + s *= 1000; + else if (0==strcmp(endp, "m")) + s *= 1; + else if (0==strcmp(endp, "dm")) + s /= 10; + else if (0==strcmp(endp, "cm")) + s /= 100; + else if (0==strcmp(endp, "mm")) + s /= 1000; + else if (0==strcmp(endp, "um")) + s /= 1e6; + else if (0==strcmp(endp, "nm")) + s /= 1e9; + else + s *= default_scale; + return s; +} + @@ -334,31 +381,11 @@ static int banner (char *args) { static int tolerance (char *args) { - char *endp = args; - T.tolerance = proj_strtod (endp, &endp); + T.tolerance = strtod_scaled (args, 1); if (HUGE_VAL==T.tolerance) { T.tolerance = 0.0005; return 1; } - while (isspace (*endp)) - endp++; - - if (0==strcmp(endp, "km")) - T.tolerance *= 1000; - else if (0==strcmp(endp, "m")) - T.tolerance *= 1; - else if (0==strcmp(endp, "dm")) - T.tolerance /= 10; - else if (0==strcmp(endp, "cm")) - T.tolerance /= 100; - else if (0==strcmp(endp, "mm")) - T.tolerance /= 1000; - else if (0==strcmp(endp, "um")) - T.tolerance /= 1e6; - else if (0==strcmp(endp, "nm")) - T.tolerance /= 1e9; - else - T.tolerance /= 1000; /* If no unit, assume mm */ return 0; } @@ -436,15 +463,17 @@ static PJ_COORD torad_if_needed (PJ *P, PJ_DIRECTION dir, PJ_COORD a) { static int accept (char *args) { int n, i; - char *endp = args; + char *endp, *prev = args; T.a = proj_coord (0,0,0,0); n = 4; for (i = 0; i < 4; i++) { - T.a.v[i] = proj_strtod (endp, &endp); - if (HUGE_VAL==T.a.v[i]) { + T.a.v[i] = proj_strtod (prev, &endp); + if (prev==endp) { n--; T.a.v[i] = 0; + break; } + prev = endp; } T.a = torad_if_needed (T.P, T.dir, T.a); if (T.verbosity > 3) @@ -457,11 +486,11 @@ static int accept (char *args) { static int roundtrip (char *args) { int ntrips; double d, r, ans; - char *endp, *endq; + char *endp; ans = proj_strtod (args, &endp); - ntrips = (int) (ans==HUGE_VAL? 100: fabs(ans)); - d = proj_strtod (endp, &endq); - d = (endp==endq)? T.tolerance: d / 1000; + ntrips = (int) (endp==args? 100: fabs(ans)); + d = strtod_scaled (endp, 1); + d = d==HUGE_VAL? T.tolerance: d; r = proj_roundtrip (T.P, PJ_FWD, ntrips, T.a); if (r > d) { if (T.verbosity > -1) { @@ -469,7 +498,7 @@ static int roundtrip (char *args) { banner (T.operation); fprintf (T.fout, "%s", T.op_ko? " -----\n": delim); fprintf (T.fout, " FAILURE in %s(%d):\n", opt_strip_path (T.curr_file), (int) lineno); - fprintf (T.fout, " roundtrip deivation: %.3f mm, expected: %.3f mm\n", 1000*r, 1000*d); + fprintf (T.fout, " roundtrip deviation: %.3f mm, expected: %.3f mm\n", 1000*r, 1000*d); } T.op_ko++; T.total_ko++; @@ -485,18 +514,20 @@ static int roundtrip (char *args) { static int expect (char *args) { double d; enum pj_io_units unit; - char *endp = args; + char *endp, *prev = args; int i; T.e = proj_coord (0,0,0,0); T.b = proj_coord (0,0,0,0); T.nargs = 4; for (i = 0; i < 4; i++) { - T.e.v[i] = proj_strtod (endp, &endp); - if (HUGE_VAL==T.e.v[i]) { + T.e.v[i] = proj_strtod (prev, &endp); + if (prev==endp) { T.nargs--; T.e.v[i] = 0; + break; } + prev = endp; } T.e = torad_if_needed (T.P, T.dir==PJ_FWD? PJ_INV:PJ_FWD, T.e); @@ -525,8 +556,8 @@ static int expect (char *args) { d = proj_xyz_dist (T.b.xyz, T.e.xyz); if (d > T.tolerance) { - if (d > 10) - d = 9.999999; + if (d > 1e6) + d = 999999.999999; if (T.verbosity > -1) { if (0==T.op_ko && T.verbosity < 2) banner (T.operation); -- cgit v1.2.3 From 47dd489e1def98f583a3288e4e3f101ae741b4e2 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Wed, 25 Oct 2017 20:46:56 +0200 Subject: Remove superfluous internal selftests from projection files matching patterns PJ_f....c and PJ_e....c --- src/PJ_fahey.c | 44 +------------ src/PJ_fouc_s.c | 44 +------------ src/PJ_gall.c | 44 +------------ src/PJ_geos.c | 59 +---------------- src/PJ_gins8.c | 30 +-------- src/PJ_gn_sinu.c | 197 ++----------------------------------------------------- src/PJ_gnom.c | 44 +------------ src/PJ_goode.c | 46 +------------ src/PJ_gstmerc.c | 45 +------------ 9 files changed, 14 insertions(+), 539 deletions(-) (limited to 'src') diff --git a/src/PJ_fahey.c b/src/PJ_fahey.c index 42318f8f..4fcb7849 100644 --- a/src/PJ_fahey.c +++ b/src/PJ_fahey.c @@ -37,46 +37,4 @@ PJ *PROJECTION(fahey) { return P; } -#ifndef PJ_SELFTEST -int pj_fahey_selftest (void) {return 0;} -#else - -int pj_fahey_selftest (void) { - double tolerance_lp = 1e-10; - double tolerance_xy = 1e-7; - - char s_args[] = {"+proj=fahey +a=6400000 +lat_1=0.5 +lat_2=2"}; - - LP fwd_in[] = { - { 2, 1}, - { 2,-1}, - {-2, 1}, - {-2,-1} - }; - - XY s_fwd_expect[] = { - { 182993.34464912376, 101603.19356988439}, - { 182993.34464912376, -101603.19356988439}, - {-182993.34464912376, 101603.19356988439}, - {-182993.34464912376, -101603.19356988439}, - }; - - XY inv_in[] = { - { 200, 100}, - { 200,-100}, - {-200, 100}, - {-200,-100} - }; - - LP s_inv_expect[] = { - {0.0021857886080359551, 0.00098424601668238403}, - {0.0021857886080359551, -0.00098424601668238403}, - {-0.0021857886080359551, 0.00098424601668238403}, - {-0.0021857886080359551, -0.00098424601668238403}, - }; - - return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect); -} - - -#endif +int pj_fahey_selftest (void) {return 10000;} diff --git a/src/PJ_fouc_s.c b/src/PJ_fouc_s.c index 343e5878..c581e690 100644 --- a/src/PJ_fouc_s.c +++ b/src/PJ_fouc_s.c @@ -67,46 +67,4 @@ PJ *PROJECTION(fouc_s) { } -#ifndef PJ_SELFTEST -int pj_fouc_s_selftest (void) {return 0;} -#else - -int pj_fouc_s_selftest (void) { - double tolerance_lp = 1e-10; - double tolerance_xy = 1e-7; - - char s_args[] = {"+proj=fouc_s +a=6400000 +lat_1=0.5 +lat_2=2"}; - - LP fwd_in[] = { - { 2, 1}, - { 2,-1}, - {-2, 1}, - {-2,-1} - }; - - XY s_fwd_expect[] = { - { 223402.14425527424, 111695.40119861449}, - { 223402.14425527424, -111695.40119861449}, - {-223402.14425527424, 111695.40119861449}, - {-223402.14425527424, -111695.40119861449}, - }; - - XY inv_in[] = { - { 200, 100}, - { 200,-100}, - {-200, 100}, - {-200,-100} - }; - - LP s_inv_expect[] = { - { 0.0017904931097838226, 0.000895246554928339}, - { 0.0017904931097838226, -0.000895246554928339}, - {-0.0017904931097838226, 0.000895246554928339}, - {-0.0017904931097838226, -0.000895246554928339}, - }; - - return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect); -} - - -#endif +int pj_fouc_s_selftest (void) {return 10000;} diff --git a/src/PJ_gall.c b/src/PJ_gall.c index 01a56e33..8a003073 100644 --- a/src/PJ_gall.c +++ b/src/PJ_gall.c @@ -41,46 +41,4 @@ PJ *PROJECTION(gall) { } -#ifndef PJ_SELFTEST -int pj_gall_selftest (void) {return 0;} -#else - -int pj_gall_selftest (void) { - double tolerance_lp = 1e-10; - double tolerance_xy = 1e-7; - - char s_args[] = {"+proj=gall +a=6400000 +lat_1=0.5 +lat_2=2"}; - - LP fwd_in[] = { - { 2, 1}, - { 2,-1}, - {-2, 1}, - {-2,-1} - }; - - XY s_fwd_expect[] = { - { 157969.17113451968, 95345.249178385886}, - { 157969.17113451968, -95345.249178385886}, - {-157969.17113451968, 95345.249178385886}, - {-157969.17113451968, -95345.249178385886}, - }; - - XY inv_in[] = { - { 200, 100}, - { 200,-100}, - {-200, 100}, - {-200,-100} - }; - - LP s_inv_expect[] = { - { 0.0025321396391918614, 0.001048846580346495}, - { 0.0025321396391918614, -0.001048846580346495}, - {-0.0025321396391918614, 0.001048846580346495}, - {-0.0025321396391918614, -0.001048846580346495}, - }; - - return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect); -} - - -#endif +int pj_gall_selftest (void) {return 10000;} diff --git a/src/PJ_geos.c b/src/PJ_geos.c index 5fd3e56b..7d527409 100644 --- a/src/PJ_geos.c +++ b/src/PJ_geos.c @@ -234,61 +234,4 @@ PJ *PROJECTION(geos) { } -#ifndef PJ_SELFTEST -int pj_geos_selftest (void) {return 0;} -#else - -int pj_geos_selftest (void) { - double tolerance_lp = 1e-10; - double tolerance_xy = 1e-7; - - char e_args[] = {"+proj=geos +ellps=GRS80 +lat_1=0.5 +lat_2=2 +h=35785831"}; - char s_args[] = {"+proj=geos +R=6400000 +lat_1=0.5 +lat_2=2 +h=35785831"}; - - LP fwd_in[] = { - { 2, 1}, - { 2,-1}, - {-2, 1}, - {-2,-1} - }; - - XY e_fwd_expect[] = { - { 222527.07036580026, 110551.30341332949}, - { 222527.07036580026, -110551.30341332949}, - {-222527.07036580026, 110551.30341332949}, - {-222527.07036580026, -110551.30341332949}, - }; - - XY s_fwd_expect[] = { - { 223289.45763579503, 111677.65745653701}, - { 223289.45763579503, -111677.65745653701}, - {-223289.45763579503, 111677.65745653701}, - {-223289.45763579503, -111677.65745653701}, - }; - - XY inv_in[] = { - { 200, 100}, - { 200,-100}, - {-200, 100}, - {-200,-100} - }; - - LP e_inv_expect[] = { - { 0.0017966305689715385, 0.00090436947723267452}, - { 0.0017966305689715385, -0.00090436947723267452}, - {-0.0017966305689715385, 0.00090436947723267452}, - {-0.0017966305689715385, -0.00090436947723267452}, - }; - - LP s_inv_expect[] = { - { 0.0017904931105078943, 0.00089524655504237148}, - { 0.0017904931105078943, -0.00089524655504237148}, - {-0.0017904931105078943, 0.00089524655504237148}, - {-0.0017904931105078943, -0.00089524655504237148}, - }; - - return pj_generic_selftest (e_args, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, s_fwd_expect, inv_in, e_inv_expect, s_inv_expect); -} - - -#endif +int pj_geos_selftest (void) {return 10000;} diff --git a/src/PJ_gins8.c b/src/PJ_gins8.c index b27ec092..26db4f31 100644 --- a/src/PJ_gins8.c +++ b/src/PJ_gins8.c @@ -31,32 +31,4 @@ PJ *PROJECTION(gins8) { } -#ifndef PJ_SELFTEST -int pj_gins8_selftest (void) {return 0;} -#else - -int pj_gins8_selftest (void) { - double tolerance_lp = 1e-10; - double tolerance_xy = 1e-7; - - char s_args[] = {"+proj=gins8 +a=6400000 +lat_1=0.5 +lat_2=2"}; - - LP fwd_in[] = { - { 2, 1}, - { 2,-1}, - {-2, 1}, - {-2,-1} - }; - - XY s_fwd_expect[] = { - { 194350.25093959007, 111703.90763533533}, - { 194350.25093959007, -111703.90763533533}, - {-194350.25093959007, 111703.90763533533}, - {-194350.25093959007, -111703.90763533533}, - }; - - return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, 0, 0, 0); -} - - -#endif +int pj_gins8_selftest (void) {return 10000;} diff --git a/src/PJ_gn_sinu.c b/src/PJ_gn_sinu.c index d13f2834..bfe26b53 100644 --- a/src/PJ_gn_sinu.c +++ b/src/PJ_gn_sinu.c @@ -120,7 +120,7 @@ PJ *PROJECTION(sinu) { if (!(Q->en = pj_enfn(P->es))) return pj_default_destructor (P, ENOMEM); - + if (P->es != 0.0) { P->inv = e_inverse; P->fwd = e_forward; @@ -184,194 +184,7 @@ PJ *PROJECTION(gn_sinu) { } -#ifndef PJ_SELFTEST -int pj_sinu_selftest (void) {return 0;} -#else - -int pj_sinu_selftest (void) { - double tolerance_lp = 1e-10; - double tolerance_xy = 1e-7; - - char e_args[] = {"+proj=sinu +ellps=GRS80 +lat_1=0.5 +lat_2=2"}; - char s_args[] = {"+proj=sinu +R=6400000 +lat_1=0.5 +lat_2=2"}; - - LP fwd_in[] = { - { 2, 1}, - { 2,-1}, - {-2, 1}, - {-2,-1} - }; - - XY e_fwd_expect[] = { - { 222605.29953946592, 110574.38855415257}, - { 222605.29953946592, -110574.38855415257}, - {-222605.29953946592, 110574.38855415257}, - {-222605.29953946592, -110574.38855415257}, - }; - - XY s_fwd_expect[] = { - { 223368.11902663155, 111701.07212763709}, - { 223368.11902663155, -111701.07212763709}, - {-223368.11902663155, 111701.07212763709}, - {-223368.11902663155, -111701.07212763709}, - }; - - XY inv_in[] = { - { 200, 100}, - { 200,-100}, - {-200, 100}, - {-200,-100} - }; - - LP e_inv_expect[] = { - { 0.0017966305684613522, 0.00090436947707945409}, - { 0.0017966305684613522, -0.00090436947707945409}, - {-0.0017966305684613522, 0.00090436947707945409}, - {-0.0017966305684613522, -0.00090436947707945409}, - }; - - LP s_inv_expect[] = { - { 0.0017904931100023887, 0.00089524655489191132}, - { 0.0017904931100023887, -0.00089524655489191132}, - {-0.0017904931100023887, 0.00089524655489191132}, - {-0.0017904931100023887, -0.00089524655489191132}, - }; - - return pj_generic_selftest (e_args, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, s_fwd_expect, inv_in, e_inv_expect, s_inv_expect); -} - - -#endif - -#ifndef PJ_SELFTEST -int pj_eck6_selftest (void) {return 0;} -#else - -int pj_eck6_selftest (void) { - double tolerance_lp = 1e-10; - double tolerance_xy = 1e-7; - - char s_args[] = {"+proj=eck6 +a=6400000 +lat_1=0.5 +lat_2=2"}; - - LP fwd_in[] = { - { 2, 1}, - { 2,-1}, - {-2, 1}, - {-2,-1} - }; - - XY s_fwd_expect[] = { - { 197021.60562899226, 126640.42073317352}, - { 197021.60562899226, -126640.42073317352}, - {-197021.60562899226, 126640.42073317352}, - {-197021.60562899226, -126640.42073317352}, - }; - - XY inv_in[] = { - { 200, 100}, - { 200,-100}, - {-200, 100}, - {-200,-100} - }; - - LP s_inv_expect[] = { - { 0.002029978749734037, 0.00078963032910382171}, - { 0.002029978749734037, -0.00078963032910382171}, - {-0.002029978749734037, 0.00078963032910382171}, - {-0.002029978749734037, -0.00078963032910382171}, - }; - - return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect); -} - - -#endif - -#ifndef PJ_SELFTEST -int pj_mbtfps_selftest (void) {return 0;} -#else - -int pj_mbtfps_selftest (void) { - double tolerance_lp = 1e-10; - double tolerance_xy = 1e-7; - - char s_args[] = {"+proj=mbtfps +a=6400000 +lat_1=0.5 +lat_2=2"}; - - LP fwd_in[] = { - { 2, 1}, - { 2,-1}, - {-2, 1}, - {-2,-1} - }; - - XY s_fwd_expect[] = { - { 204740.11747857218, 121864.72971934026}, - { 204740.11747857218, -121864.72971934026}, - {-204740.11747857218, 121864.72971934026}, - {-204740.11747857218, -121864.72971934026}, - }; - - XY inv_in[] = { - { 200, 100}, - { 200,-100}, - {-200, 100}, - {-200,-100} - }; - - LP s_inv_expect[] = { - { 0.0019534152166442065, 0.00082057965689633387}, - { 0.0019534152166442065, -0.00082057965689633387}, - {-0.0019534152166442065, 0.00082057965689633387}, - {-0.0019534152166442065, -0.00082057965689633387}, - }; - - return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect); -} - - -#endif - - -#ifndef PJ_SELFTEST -int pj_gn_sinu_selftest (void) {return 0;} -#else - -int pj_gn_sinu_selftest (void) { - double tolerance_lp = 1e-10; - double tolerance_xy = 1e-7; - - char s_args[] = {"+proj=gn_sinu +a=6400000 +lat_1=0.5 +lat_2=2 +m=1 +n=2"}; - - LP fwd_in[] = { - { 2, 1}, - { 2,-1}, - {-2, 1}, - {-2,-1} - }; - - XY s_fwd_expect[] = { - { 223385.13250469571, 111698.23644718733}, - { 223385.13250469571, -111698.23644718733}, - {-223385.13250469571, 111698.23644718733}, - {-223385.13250469571, -111698.23644718733}, - }; - - XY inv_in[] = { - { 200, 100}, - { 200,-100}, - {-200, 100}, - {-200,-100} - }; - - LP s_inv_expect[] = { - { 0.0017904931098931057, 0.00089524655491012516}, - { 0.0017904931098931057, -0.00089524655491012516}, - {-0.0017904931098931057, 0.00089524655491012516}, - {-0.0017904931098931057, -0.00089524655491012516}, - }; - - return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect); -} - - -#endif +int pj_sinu_selftest (void) {return 10000;} +int pj_eck6_selftest (void) {return 10000;} +int pj_mbtfps_selftest (void) {return 10000;} +int pj_gn_sinu_selftest (void) {return 10000;} diff --git a/src/PJ_gnom.c b/src/PJ_gnom.c index 1d3f3386..2aedefec 100644 --- a/src/PJ_gnom.c +++ b/src/PJ_gnom.c @@ -136,46 +136,4 @@ PJ *PROJECTION(gnom) { } -#ifndef PJ_SELFTEST -int pj_gnom_selftest (void) {return 0;} -#else - -int pj_gnom_selftest (void) { - double tolerance_lp = 1e-10; - double tolerance_xy = 1e-7; - - char s_args[] = {"+proj=gnom +a=6400000 +lat_1=0.5 +lat_2=2"}; - - LP fwd_in[] = { - { 2, 1}, - { 2,-1}, - {-2, 1}, - {-2,-1} - }; - - XY s_fwd_expect[] = { - { 223492.92474718543, 111780.50920659291}, - { 223492.92474718543, -111780.50920659291}, - {-223492.92474718543, 111780.50920659291}, - {-223492.92474718543, -111780.50920659291}, - }; - - XY inv_in[] = { - { 200, 100}, - { 200,-100}, - {-200, 100}, - {-200,-100} - }; - - LP s_inv_expect[] = { - { 0.0017904931092009798, 0.00089524655438192376}, - { 0.0017904931092009798, -0.00089524655438192376}, - {-0.0017904931092009798, 0.00089524655438192376}, - {-0.0017904931092009798, -0.00089524655438192376}, - }; - - return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect); -} - - -#endif +int pj_gnom_selftest (void) {return 10000;} diff --git a/src/PJ_goode.c b/src/PJ_goode.c index 3bfeb21f..3c5d172d 100644 --- a/src/PJ_goode.c +++ b/src/PJ_goode.c @@ -71,7 +71,7 @@ PJ *PROJECTION(goode) { Q->moll->ctx = P->ctx; if (!(Q->sinu = pj_sinu(Q->sinu)) || !(Q->moll = pj_moll(Q->moll))) return destructor (P, ENOMEM); - + P->fwd = s_forward; P->inv = s_inverse; @@ -79,46 +79,4 @@ PJ *PROJECTION(goode) { } -#ifndef PJ_SELFTEST -int pj_goode_selftest (void) {return 0;} -#else - -int pj_goode_selftest (void) { - double tolerance_lp = 1e-10; - double tolerance_xy = 1e-7; - - char s_args[] = {"+proj=goode +a=6400000 +lat_1=0.5 +lat_2=2"}; - - LP fwd_in[] = { - { 2, 1}, - { 2,-1}, - {-2, 1}, - {-2,-1} - }; - - XY s_fwd_expect[] = { - { 223368.11902663155, 111701.07212763709}, - { 223368.11902663155, -111701.07212763709}, - {-223368.11902663155, 111701.07212763709}, - {-223368.11902663155, -111701.07212763709}, - }; - - XY inv_in[] = { - { 200, 100}, - { 200,-100}, - {-200, 100}, - {-200,-100} - }; - - LP s_inv_expect[] = { - { 0.0017904931100023887, 0.00089524655489191132}, - { 0.0017904931100023887, -0.00089524655489191132}, - {-0.0017904931100023887, 0.00089524655489191132}, - {-0.0017904931100023887, -0.00089524655489191132}, - }; - - return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect); -} - - -#endif +int pj_goode_selftest (void) {return 10000;} diff --git a/src/PJ_gstmerc.c b/src/PJ_gstmerc.c index c2846761..efb44b88 100644 --- a/src/PJ_gstmerc.c +++ b/src/PJ_gstmerc.c @@ -69,47 +69,4 @@ PJ *PROJECTION(gstmerc) { } -#ifndef PJ_SELFTEST -int pj_gstmerc_selftest (void) {return 0;} -#else - -int pj_gstmerc_selftest (void) { - double tolerance_lp = 1e-10; - double tolerance_xy = 1e-7; - - char s_args[] = {"+proj=gstmerc +R=6400000 +lat_1=0.5 +lat_2=2"}; - - LP fwd_in[] = { - { 2, 1}, - { 2,-1}, - {-2, 1}, - {-2,-1} - }; - - - XY s_fwd_expect[] = { - { 223413.46640632182, 111769.14504058557}, - { 223413.46640632182, -111769.14504058668}, - {-223413.46640632302, 111769.14504058557}, - {-223413.46640632302, -111769.14504058668}, - }; - - XY inv_in[] = { - { 200, 100}, - { 200,-100}, - {-200, 100}, - {-200,-100} - }; - - LP s_inv_expect[] = { - { 0.0017904931097109673, 0.0008952465544509083}, - { 0.0017904931097109673, -0.0008952465544509083}, - {-0.0017904931097109673, 0.0008952465544509083}, - {-0.0017904931097109673, -0.0008952465544509083}, - }; - - return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect); -} - - -#endif +int pj_gstmerc_selftest (void) {return 10000;} -- cgit v1.2.3 From 3a2bd267a67d41a461946a6f7b0a99262f47d8a7 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Mon, 23 Oct 2017 11:38:56 +0200 Subject: Refactor grid shift functions. This refactoring of the grid shift functions allows for easier access to the actual grid values, as well as making it possible to perform horizontal grid shift on a single coordinate without using the pj_apply_gridshift* functions. The latter simplifies the execution path of the forward and inverse functions in PJ_hgridshift.c. This commit introduces proj_*grid_init, proj_*grid_value and proj_hgrid_apply. The init functions initialises horizontal and vertical grids respectivelive (wrappers for pj_gridlist_from_nadgrids with simpler parameters). The proj_*grid_value functions returns the specific grid value at coordinate lp. The proj_hgrid_apply function applies the grid offset to the input coordinate and outputs the adjusted coordinate. --- src/PJ_hgridshift.c | 52 ++------ src/PJ_vgridshift.c | 25 ++-- src/pj_apply_gridshift.c | 275 +++++++++++++++++++++++++++-------------- src/pj_apply_vgridshift.c | 305 ++++++++++++++++++++++++++++------------------ src/proj_internal.h | 7 +- 5 files changed, 398 insertions(+), 266 deletions(-) (limited to 'src') diff --git a/src/PJ_hgridshift.c b/src/PJ_hgridshift.c index 0adc9e00..659039ab 100644 --- a/src/PJ_hgridshift.c +++ b/src/PJ_hgridshift.c @@ -4,7 +4,6 @@ PROJ_HEAD(hgridshift, "Horizontal grid shift"); - static XYZ forward_3d(LPZ lpz, PJ *P) { PJ_TRIPLET point; point.lpz = lpz; @@ -12,9 +11,7 @@ static XYZ forward_3d(LPZ lpz, PJ *P) { 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 ); + point.lp = proj_hgrid_apply(P, point.lp, PJ_FWD); } return point.xyz; @@ -28,9 +25,7 @@ static LPZ reverse_3d(XYZ xyz, PJ *P) { 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 ); + point.lp = proj_hgrid_apply(P, point.lp, PJ_INV); } return point.lpz; @@ -52,29 +47,8 @@ static PJ_OBS reverse_obs(PJ_OBS obs, PJ *P) { -#if 0 -static XY forward_xy(LP lp, PJ *P) { - PJ_TRIPLET point; - point.lp = lp; - point.lpz.z = 0; - point.xyz = forward_3d (point.lpz, P); - return point.xy; -} - - -static LP reverse_lp(XY xy, PJ *P) { - PJ_TRIPLET point; - point.xy = xy; - point.xyz.z = 0; - point.lpz = reverse_3d (point.xyz, P); - return point.lp; -} -#endif - - - PJ *PROJECTION(hgridshift) { - + P->fwdobs = forward_obs; P->invobs = reverse_obs; P->fwd3d = forward_3d; @@ -90,14 +64,12 @@ PJ *PROJECTION(hgridshift) { return pj_default_destructor (P, PJD_ERR_NO_ARGS); } - /* 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) ); + proj_hgrid_init(P, "grids"); /* Was gridlist compiled properly? */ - if ( pj_ctx_get_errno(pj_get_ctx(P)) ) { + if ( proj_errno(P) ) { proj_log_error(P, "hgridshift: could not find required grid(s)."); - return pj_default_destructor (P, PJD_ERR_FAILED_TO_LOAD_GRID); + return pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); } return P; @@ -119,26 +91,28 @@ int pj_hgridshift_selftest (void) { proj_destroy (P); return 99; } - + /* fail on purpose: open non-existing grid */ P = proj_create(PJ_DEFAULT_CTX, "+proj=hgridshift +grids=@nonexistinggrid.gsb,anothernonexistinggrid.gsb"); if (0!=P) { proj_destroy (P); return 999; } - + /* Failure most likely means the grid is missing */ P = proj_create(PJ_DEFAULT_CTX, "+proj=hgridshift +grids=nzgd2kgrid0005.gsb +ellps=GRS80"); if (0==P) return 10; - + a = proj_obs_null; a.coo.lpz.lam = PJ_TORAD(173); a.coo.lpz.phi = PJ_TORAD(-45); - + dist = proj_roundtrip (P, PJ_FWD, 1, a.coo); - if (dist > 0.00000001) + if (dist > 0.00000001) { + printf("dist: %f\n",dist); return 1; + } expect.coo.lpz.lam = PJ_TORAD(172.999892181021551); expect.coo.lpz.phi = PJ_TORAD(-45.001620431954613); diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c index ededd544..691f791b 100644 --- a/src/PJ_vgridshift.c +++ b/src/PJ_vgridshift.c @@ -9,14 +9,10 @@ static XYZ forward_3d(LPZ lpz, PJ *P) { PJ_TRIPLET point; point.lpz = lpz; - if (P->gridlist != NULL) { + if (P->vgridlist_geoid != 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->gridlist), - &(P->gridlist_count), - 1, 1, 0, - &point.xyz.x, &point.xyz.y, &point.xyz.z ); + point.xyz.z -= proj_vgrid_value(P, point.lp); } return point.xyz; @@ -27,14 +23,10 @@ static LPZ reverse_3d(XYZ xyz, PJ *P) { PJ_TRIPLET point; point.xyz = xyz; - if (P->gridlist != NULL) { + if (P->vgridlist_geoid != 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->gridlist), - &(P->gridlist_count), - 0, 1, 0, - &point.xyz.x, &point.xyz.y, &point.xyz.z ); + point.xyz.z += proj_vgrid_value(P, point.lp); } return point.lpz; @@ -61,14 +53,13 @@ PJ *PROJECTION(vgridshift) { return pj_default_destructor(P, PJD_ERR_NO_ARGS); } - /* 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) ); + /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */ + proj_vgrid_init(P, "grids"); /* Was gridlist compiled properly? */ - if ( pj_ctx_get_errno(P->ctx) ) { + if ( proj_errno(P) ) { proj_log_error(P, "vgridshift: could not find required grid(s)."); - return pj_default_destructor(P, -38); + return pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID); } P->fwdobs = forward_obs; diff --git a/src/pj_apply_gridshift.c b/src/pj_apply_gridshift.c index 91e2de26..45887abd 100644 --- a/src/pj_apply_gridshift.c +++ b/src/pj_apply_gridshift.c @@ -2,7 +2,7 @@ * Project: PROJ.4 * Purpose: Apply datum shifts based on grid shift files (normally NAD27 to * NAD83 or the reverse). This module is responsible for keeping - * a list of loaded grids, and calling with each one that is + * a list of loaded grids, and calling with each one that is * allowed for a given datum (expressed as the nadgrids= parameter). * Author: Frank Warmerdam, warmerdam@pobox.com * @@ -30,9 +30,10 @@ #define PJ_LIB__ -#include #include #include +#include "proj_internal.h" +#include "projects.h" /************************************************************************/ /* pj_apply_gridshift() */ @@ -43,7 +44,7 @@ /* it to honour our public api. */ /************************************************************************/ -int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse, +int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse, long point_count, int point_offset, double *x, double *y, double *z ) @@ -51,16 +52,16 @@ int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse, PJ_GRIDINFO **gridlist; int grid_count; int ret; - + gridlist = pj_gridlist_from_nadgrids( ctx, nadgrids, &grid_count ); if( gridlist == NULL || grid_count == 0 ) return ctx->last_errno; - ret = pj_apply_gridshift_3( ctx, gridlist, grid_count, inverse, + ret = pj_apply_gridshift_3( ctx, gridlist, grid_count, inverse, point_count, point_offset, x, y, z ); - /* + /* ** Note this frees the array of grid list pointers, but not the grids ** which is as intended. The grids themselves live on. */ @@ -72,13 +73,13 @@ int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse, /************************************************************************/ /* pj_apply_gridshift_2() */ /* */ -/* This implementation takes uses the gridlist from a coordinate */ +/* This implementation uses the gridlist from a coordinate */ /* system definition. If the gridlist has not yet been */ /* populated in the coordinate system definition we set it up */ /* now. */ /************************************************************************/ -int pj_apply_gridshift_2( PJ *defn, int inverse, +int pj_apply_gridshift_2( PJ *defn, int inverse, long point_count, int point_offset, double *x, double *y, double *z ) @@ -86,10 +87,10 @@ int pj_apply_gridshift_2( PJ *defn, int inverse, if( defn->catalog_name != NULL ) return pj_gc_apply_gridshift( defn, inverse, point_count, point_offset, x, y, z ); - + if( defn->gridlist == NULL ) { - defn->gridlist = + defn->gridlist = pj_gridlist_from_nadgrids( pj_get_ctx( defn ), pj_param(defn->ctx, defn->params,"snadgrids").s, &(defn->gridlist_count) ); @@ -97,12 +98,76 @@ int pj_apply_gridshift_2( PJ *defn, int inverse, if( defn->gridlist == NULL || defn->gridlist_count == 0 ) return defn->ctx->last_errno; } - + return pj_apply_gridshift_3( pj_get_ctx( defn ), - defn->gridlist, defn->gridlist_count, inverse, + defn->gridlist, defn->gridlist_count, inverse, point_count, point_offset, x, y, z ); } +/************************************************************************/ +/* find_ctable() */ +/* */ +/* Determine which grid is the correct given an input coordinate. */ +/************************************************************************/ + +static struct CTABLE* find_ctable(projCtx ctx, LP input, int grid_count, PJ_GRIDINFO **tables) { + int itable; + double epsilon; + struct CTABLE *ct = NULL; + + /* keep trying till we find a table that works */ + for( itable = 0; itable < grid_count; itable++ ) + { + + PJ_GRIDINFO *gi = tables[itable]; + ct = gi->ct; + epsilon = (fabs(ct->del.phi)+fabs(ct->del.lam))/10000.0; + /* skip tables that don't match our point at all. */ + if ( ct->ll.phi - epsilon > input.phi + || ct->ll.lam - epsilon > input.lam + || (ct->ll.phi + (ct->lim.phi-1) * ct->del.phi + epsilon < input.phi) + || (ct->ll.lam + (ct->lim.lam-1) * ct->del.lam + epsilon < input.lam) ) { + continue; + } + + /* If we have child nodes, check to see if any of them apply. */ + while( gi->child ) + { + PJ_GRIDINFO *child; + + for( child = gi->child; child != NULL; child = child->next ) + { + struct CTABLE *ct1 = child->ct; + epsilon = (fabs(ct1->del.phi)+fabs(ct1->del.lam))/10000.0; + + if( ct1->ll.phi - epsilon > input.phi + || ct1->ll.lam - epsilon > input.lam + || (ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi + epsilon < input.phi) + || (ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam + epsilon < input.lam) ) { + continue; + } + break; + } + + /* If we didn't find a child then nothing more to do */ + if( child == NULL ) break; + + /* Otherwise use the child, first checking it's children */ + gi = child; + ct = child->ct; + } + /* load the grid shift info if we don't have it. */ + if( ct->cvs == NULL) { + if (!pj_gridinfo_load( ctx, gi ) ) { + pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); + return NULL; + } + } + /* if we get this far we have found a suitable grid */ + break; + } + return ct; +} /************************************************************************/ /* pj_apply_gridshift_3() */ @@ -113,16 +178,16 @@ int pj_apply_gridshift_2( PJ *defn, int inverse, int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count, int inverse, long point_count, int point_offset, double *x, double *y, double *z ) - { int i; + struct CTABLE *ct; static int debug_count = 0; (void) z; if( tables == NULL || grid_count == 0 ) { - pj_ctx_set_errno( ctx, -38); - return -38; + pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return PJD_ERR_FAILED_TO_LOAD_GRID; } ctx->last_errno = 0; @@ -138,100 +203,39 @@ int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count, output.phi = HUGE_VAL; output.lam = HUGE_VAL; - /* keep trying till we find a table that works */ - for( itable = 0; itable < grid_count; itable++ ) - { - PJ_GRIDINFO *gi = tables[itable]; - struct CTABLE *ct = gi->ct; - double epsilon = (fabs(ct->del.phi)+fabs(ct->del.lam))/10000.0; - - /* skip tables that don't match our point at all. */ - if( ct->ll.phi - epsilon > input.phi - || ct->ll.lam - epsilon > input.lam - || (ct->ll.phi + (ct->lim.phi-1) * ct->del.phi + epsilon - < input.phi) - || (ct->ll.lam + (ct->lim.lam-1) * ct->del.lam + epsilon - < input.lam) ) - continue; - - /* If we have child nodes, check to see if any of them apply. */ - while( gi->child ) - { - PJ_GRIDINFO *child; - - for( child = gi->child; child != NULL; child = child->next ) - { - struct CTABLE *ct1 = child->ct; - epsilon = - (fabs(ct1->del.phi)+fabs(ct1->del.lam))/10000.0; - - if( ct1->ll.phi - epsilon > input.phi - || ct1->ll.lam - epsilon > input.lam - || (ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi + epsilon - < input.phi) - || (ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam + epsilon - < input.lam) ) - continue; - - break; - } - - /* If we didn't find a child then nothing more to do */ - - if( child == NULL ) break; - - /* Otherwise use the child, first checking it's children */ - - gi = child; - ct = child->ct; - } + ct = find_ctable(ctx, input, grid_count, tables); + output = nad_cvt( input, inverse, ct ); - /* load the grid shift info if we don't have it. */ - if( ct->cvs == NULL && !pj_gridinfo_load( ctx, gi ) ) - { - pj_ctx_set_errno( ctx, -38 ); - return -38; - } - - output = nad_cvt( input, inverse, ct ); - if( output.lam != HUGE_VAL ) - { - if( debug_count++ < 20 ) - pj_log( ctx, PJ_LOG_DEBUG_MINOR, - "pj_apply_gridshift(): used %s", ct->id ); - break; - } - } + if ( output.lam != HUGE_VAL && debug_count++ < 20 ) + pj_log( ctx, PJ_LOG_DEBUG_MINOR, "pj_apply_gridshift(): used %s", ct->id ); - if( output.lam == HUGE_VAL ) + if ( output.lam == HUGE_VAL ) { if( ctx->debug_level >= PJ_LOG_DEBUG_MAJOR ) { pj_log( ctx, PJ_LOG_DEBUG_MAJOR, "pj_apply_gridshift(): failed to find a grid shift table for\n" " location (%.7fdW,%.7fdN)", - x[io] * RAD_TO_DEG, + x[io] * RAD_TO_DEG, y[io] * RAD_TO_DEG ); for( itable = 0; itable < grid_count; itable++ ) { PJ_GRIDINFO *gi = tables[itable]; if( itable == 0 ) - pj_log( ctx, PJ_LOG_DEBUG_MAJOR, - " tried: %s", gi->gridname ); + pj_log( ctx, PJ_LOG_DEBUG_MAJOR, " tried: %s", gi->gridname ); else - pj_log( ctx, PJ_LOG_DEBUG_MAJOR, - ",%s", gi->gridname ); + pj_log( ctx, PJ_LOG_DEBUG_MAJOR, ",%s", gi->gridname ); } } - /* - * We don't actually have any machinery currently to set the - * following macro, so this is mostly kept here to make it clear - * how we ought to operate if we wanted to make it super clear + /* + * We don't actually have any machinery currently to set the + * following macro, so this is mostly kept here to make it clear + * how we ought to operate if we wanted to make it super clear * that an error has occurred when points are outside our available - * datum shift areas. But if this is on, we will find that "low - * value" points on the fringes of some datasets will completely - * fail causing lots of problems when it is more or less ok to + * datum shift areas. But if this is on, we will find that "low + * value" points on the fringes of some datasets will completely + * fail causing lots of problems when it is more or less ok to * just not apply a datum shift. So rather than deal with * that we just fallback to no shift. (see also bug #45). */ @@ -252,3 +256,92 @@ int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **tables, int grid_count, return 0; } +/**********************************************/ +int proj_hgrid_init(PJ* P, const char *grids) { +/********************************************** + + Initizalize and populate list of horizontal + grids. + + Takes a PJ-object and the plus-parameter + name that is used in the proj-string to + specify the grids to load, e.g. "+grids". + The + should be left out here. + + Returns the number of loaded grids. + +***********************************************/ + + /* prepend "s" to the "grids" string to allow usage with pj_param */ + char *sgrids = (char *) pj_malloc( (strlen(grids)+1) *sizeof(char) ); + sprintf(sgrids, "%s%s", "s", grids); + + if (P->gridlist == NULL) { + P->gridlist = pj_gridlist_from_nadgrids( + P->ctx, + pj_param(P->ctx, P->params, sgrids).s, + &(P->gridlist_count) + ); + + if( P->gridlist == NULL || P->gridlist_count == 0 ) { + pj_dealloc(sgrids); + return 0; + } + } + + if (P->gridlist_count == 0) { + proj_errno_set(P, PJD_ERR_FAILED_TO_LOAD_GRID); + } + + pj_dealloc(sgrids); + return P->gridlist_count; +} + +/********************************************/ +/* proj_hgrid_value() */ +/* */ +/* Return coordinate offset in grid */ +/********************************************/ +LP proj_hgrid_value(PJ *P, LP lp) { + struct CTABLE *ct; + LP out; + + ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist); + + /* normalize input to ll origin */ + lp.lam -= ct->ll.lam; + lp.phi -= ct->ll.phi; + lp.lam = adjlon(lp.lam - M_PI) + M_PI; + + out = nad_intr(lp, ct); + + if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) { + pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA); + } + + return out; +} + +LP proj_hgrid_apply(PJ *P, LP lp, PJ_DIRECTION direction) { + struct CTABLE *ct; + int inverse; + LP out; + + out.lam = HUGE_VAL; out.phi = HUGE_VAL; + + ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist); + + if (ct == NULL || ct->cvs == NULL) { + pj_ctx_set_errno( P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); + return out; + } + + inverse = direction == PJ_FWD ? 0 : 1; + out = nad_cvt(lp, inverse, ct); + + if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) + pj_ctx_set_errno(P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + + return out; + +} diff --git a/src/pj_apply_vgridshift.c b/src/pj_apply_vgridshift.c index 35047a19..3c7cc210 100644 --- a/src/pj_apply_vgridshift.c +++ b/src/pj_apply_vgridshift.c @@ -28,23 +28,120 @@ #define PJ_LIB__ -#include #include #include +#include "proj_internal.h" +#include "projects.h" + +static double pj_read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GRIDINFO **tables, struct CTABLE *ct) { + int itable = 0; + double value = HUGE_VAL; + double grid_x, grid_y; + int grid_ix, grid_iy; + int grid_ix2, grid_iy2; + float *cvs; + /* do not deal with NaN coordinates */ + if( input.phi != input.phi || input.lam != input.lam ) + itable = *gridlist_count_p; + + /* keep trying till we find a table that works */ + for ( ; itable < *gridlist_count_p; itable++ ) + { + PJ_GRIDINFO *gi = tables[itable]; + + ct = gi->ct; + + /* skip tables that don't match our point at all. */ + if( ct->ll.phi > input.phi || ct->ll.lam > input.lam + || ct->ll.phi + (ct->lim.phi-1) * ct->del.phi < input.phi + || ct->ll.lam + (ct->lim.lam-1) * ct->del.lam < input.lam ) + continue; + + /* If we have child nodes, check to see if any of them apply. */ + while( gi->child != NULL ) + { + PJ_GRIDINFO *child; + + for( child = gi->child; child != NULL; child = child->next ) + { + struct CTABLE *ct1 = child->ct; + + if( ct1->ll.phi > input.phi || ct1->ll.lam > input.lam + || ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi < input.phi + || ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam < input.lam) + continue; + + break; + } + + /* we didn't find a more refined child node to use, so go with current grid */ + if( child == NULL ) + { + break; + } + + /* Otherwise let's try for childrens children .. */ + gi = child; + ct = child->ct; + } + + /* load the grid shift info if we don't have it. */ + if( ct->cvs == NULL && !pj_gridinfo_load( pj_get_ctx(defn), gi ) ) + { + pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID ); + return PJD_ERR_FAILED_TO_LOAD_GRID; + } + + } + + /* Interpolation a location within the grid */ + grid_x = (input.lam - ct->ll.lam) / ct->del.lam; + grid_y = (input.phi - ct->ll.phi) / ct->del.phi; + grid_ix = (int) floor(grid_x); + grid_iy = (int) floor(grid_y); + grid_x -= grid_ix; + grid_y -= grid_iy; + + grid_ix2 = grid_ix + 1; + if( grid_ix2 >= ct->lim.lam ) + grid_ix2 = ct->lim.lam - 1; + grid_iy2 = grid_iy + 1; + if( grid_iy2 >= ct->lim.phi ) + grid_iy2 = ct->lim.phi - 1; + + cvs = (float *) ct->cvs; + value = cvs[grid_ix + grid_iy * ct->lim.lam] + * (1.0-grid_x) * (1.0-grid_y) + + cvs[grid_ix2 + grid_iy * ct->lim.lam] + * (grid_x) * (1.0-grid_y) + + cvs[grid_ix + grid_iy2 * ct->lim.lam] + * (1.0-grid_x) * (grid_y) + + cvs[grid_ix2 + grid_iy2 * ct->lim.lam] + * (grid_x) * (grid_y); + + /* nodata? */ + /* GTX official nodata value if -88.88880f, but some grids also */ + /* use other big values for nodata (e.g naptrans2008.gtx has */ + /* nodata values like -2147479936), so test them too */ + if( value > 1000 || value < -1000 || value == -88.88880f ) + value = HUGE_VAL; + + + return value; +} /************************************************************************/ /* pj_apply_vgridshift() */ /* */ -/* This implementation takes uses the gridlist from a coordinate */ +/* This implementation takes uses the gridlist from a coordinate */ /* system definition. If the gridlist has not yet been */ /* populated in the coordinate system definition we set it up */ /* now. */ /************************************************************************/ - int pj_apply_vgridshift( PJ *defn, const char *listname, - PJ_GRIDINFO ***gridlist_p, + PJ_GRIDINFO ***gridlist_p, int *gridlist_count_p, - int inverse, + int inverse, long point_count, int point_offset, double *x, double *y, double *z ) @@ -52,22 +149,23 @@ int pj_apply_vgridshift( PJ *defn, const char *listname, int i; static int debug_count = 0; PJ_GRIDINFO **tables; + struct CTABLE ct; if( *gridlist_p == NULL ) { - *gridlist_p = - pj_gridlist_from_nadgrids( pj_get_ctx(defn), + *gridlist_p = + pj_gridlist_from_nadgrids( pj_get_ctx(defn), pj_param(defn->ctx,defn->params,listname).s, gridlist_count_p ); if( *gridlist_p == NULL || *gridlist_count_p == 0 ) return defn->ctx->last_errno; } - + if( *gridlist_count_p == 0 ) { - pj_ctx_set_errno( defn->ctx, -38); - return -38; + pj_ctx_set_errno( defn->ctx, PJD_ERR_FAILED_TO_LOAD_GRID); + return PJD_ERR_FAILED_TO_LOAD_GRID; } tables = *gridlist_p; @@ -75,127 +173,37 @@ int pj_apply_vgridshift( PJ *defn, const char *listname, for( i = 0; i < point_count; i++ ) { + double value; long io = i * point_offset; LP input; - int itable = 0; - double value = HUGE_VAL; input.phi = y[io]; input.lam = x[io]; - /* do not deal with NaN coordinates */ - if( input.phi != input.phi || input.lam != input.lam ) - itable = *gridlist_count_p; + value = pj_read_vgrid_value(defn, input, gridlist_count_p, tables, &ct); - /* keep trying till we find a table that works */ - for( ; itable < *gridlist_count_p; itable++ ) + if( inverse ) + z[io] -= value; + else + z[io] += value; + if( value != HUGE_VAL ) { - PJ_GRIDINFO *gi = tables[itable]; - struct CTABLE *ct = gi->ct; - double grid_x, grid_y; - int grid_ix, grid_iy; - int grid_ix2, grid_iy2; - float *cvs; - - /* skip tables that don't match our point at all. */ - if( ct->ll.phi > input.phi || ct->ll.lam > input.lam - || ct->ll.phi + (ct->lim.phi-1) * ct->del.phi < input.phi - || ct->ll.lam + (ct->lim.lam-1) * ct->del.lam < input.lam ) - continue; - - /* If we have child nodes, check to see if any of them apply. */ - while( gi->child != NULL ) - { - PJ_GRIDINFO *child; - - for( child = gi->child; child != NULL; child = child->next ) - { - struct CTABLE *ct1 = child->ct; - - if( ct1->ll.phi > input.phi || ct1->ll.lam > input.lam - || ct1->ll.phi+(ct1->lim.phi-1)*ct1->del.phi < input.phi - || ct1->ll.lam+(ct1->lim.lam-1)*ct1->del.lam < input.lam) - continue; - - break; - } - - /* we didn't find a more refined child node to use, so go with current grid */ - if( child == NULL ) - { - break; - } - - /* Otherwise let's try for childrens children .. */ - gi = child; - ct = child->ct; - } - - /* load the grid shift info if we don't have it. */ - if( ct->cvs == NULL && !pj_gridinfo_load( pj_get_ctx(defn), gi ) ) - { - pj_ctx_set_errno( defn->ctx, -38 ); - return -38; - } - - /* Interpolation a location within the grid */ - grid_x = (input.lam - ct->ll.lam) / ct->del.lam; - grid_y = (input.phi - ct->ll.phi) / ct->del.phi; - grid_ix = (int) floor(grid_x); - grid_iy = (int) floor(grid_y); - grid_x -= grid_ix; - grid_y -= grid_iy; - - grid_ix2 = grid_ix + 1; - if( grid_ix2 >= ct->lim.lam ) - grid_ix2 = ct->lim.lam - 1; - grid_iy2 = grid_iy + 1; - if( grid_iy2 >= ct->lim.phi ) - grid_iy2 = ct->lim.phi - 1; - - cvs = (float *) ct->cvs; - value = cvs[grid_ix + grid_iy * ct->lim.lam] - * (1.0-grid_x) * (1.0-grid_y) - + cvs[grid_ix2 + grid_iy * ct->lim.lam] - * (grid_x) * (1.0-grid_y) - + cvs[grid_ix + grid_iy2 * ct->lim.lam] - * (1.0-grid_x) * (grid_y) - + cvs[grid_ix2 + grid_iy2 * ct->lim.lam] - * (grid_x) * (grid_y); - - /* nodata? */ - /* GTX official nodata value if -88.88880f, but some grids also */ - /* use other big values for nodata (e.g naptrans2008.gtx has */ - /* nodata values like -2147479936), so test them too */ - if( value > 1000 || value < -1000 || value == -88.88880f ) - value = HUGE_VAL; - else - { - if( inverse ) - z[io] -= value; - else - z[io] += value; - } - - if( value != HUGE_VAL ) - { - if( debug_count++ < 20 ) - pj_log( defn->ctx, PJ_LOG_DEBUG_MINOR, - "pj_apply_gridshift(): used %s", - ct->id ); + if( debug_count++ < 20 ) { + proj_log_trace(defn, "pj_apply_gridshift(): used %s", ct.id); break; } } if( value == HUGE_VAL ) { + int itable; char gridlist[3000]; - pj_log( defn->ctx, PJ_LOG_DEBUG_MAJOR, - "pj_apply_vgridshift(): failed to find a grid shift table for\n" - " location (%.7fdW,%.7fdN)", - x[io] * RAD_TO_DEG, - y[io] * RAD_TO_DEG ); + proj_log_debug(defn, + "pj_apply_vgridshift(): failed to find a grid shift table for\n" + " location (%.7fdW,%.7fdN)", + x[io] * RAD_TO_DEG, + y[io] * RAD_TO_DEG ); gridlist[0] = '\0'; for( itable = 0; itable < *gridlist_count_p; itable++ ) @@ -212,10 +220,10 @@ int pj_apply_vgridshift( PJ *defn, const char *listname, else sprintf( gridlist+strlen(gridlist), ",%s", gi->gridname ); } - pj_log( defn->ctx, PJ_LOG_DEBUG_MAJOR, - "%s", gridlist ); - + + proj_log_debug(defn, "%s", gridlist); pj_ctx_set_errno( defn->ctx, PJD_ERR_GRID_AREA ); + return PJD_ERR_GRID_AREA; } } @@ -223,3 +231,64 @@ int pj_apply_vgridshift( PJ *defn, const char *listname, return 0; } +/**********************************************/ +int proj_vgrid_init(PJ* P, const char *grids) { +/********************************************** + + Initizalize and populate gridlist. + + Takes a PJ-object and the plus-parameter + name that is used in the proj-string to + specify the grids to load, e.g. "+grids". + The + should be left out here. + + Returns the number of loaded grids. + +***********************************************/ + + /* prepend "s" to the "grids" string to allow usage with pj_param */ + char *sgrids = (char *) pj_malloc( (strlen(grids)+1+1) *sizeof(char) ); + sprintf(sgrids, "%s%s", "s", grids); + + if (P->vgridlist_geoid == NULL) { + P->vgridlist_geoid = pj_gridlist_from_nadgrids( + P->ctx, + pj_param(P->ctx, P->params, sgrids).s, + &(P->vgridlist_geoid_count) + ); + + if( P->vgridlist_geoid == NULL || P->vgridlist_geoid_count == 0 ) { + pj_dealloc(sgrids); + return 0; + } + } + + if (P->vgridlist_geoid_count == 0) { + proj_errno_set(P, PJD_ERR_FAILED_TO_LOAD_GRID); + } + + pj_dealloc(sgrids); + return P->vgridlist_geoid_count; +} + +/***********************************************/ +double proj_vgrid_value(PJ *P, LP lp){ +/*********************************************** + + Read grid value at position lp in grids loaded + with proj_grid_init. + + Returns the grid value of the given coordinate. + +************************************************/ + + struct CTABLE used_grid; + double value; + memset(&used_grid, 0, sizeof(struct CTABLE)); + + value = pj_read_vgrid_value(P, lp, &(P->vgridlist_geoid_count), P->vgridlist_geoid, &used_grid); + proj_log_trace(P, "proj_vgrid_value: (%f, %f) = %f", lp.lam*RAD_TO_DEG, lp.phi*RAD_TO_DEG, value); + + return value; +} + diff --git a/src/proj_internal.h b/src/proj_internal.h index 5773637c..29e73e70 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -89,7 +89,12 @@ PJ_OBS pj_invobs (PJ_OBS obs, PJ *P); PJ_COORD pj_fwdcoord (PJ_COORD coo, PJ *P); PJ_COORD pj_invcoord (PJ_COORD coo, PJ *P); - +/* Grid functionality */ +int proj_vgrid_init(PJ *P, const char *grids); +int proj_hgrid_init(PJ *P, const char *grids); +double proj_vgrid_value(PJ *P, LP lp); +LP proj_hgrid_value(PJ *P, LP lp); +LP proj_hgrid_apply(PJ *P, LP lp, PJ_DIRECTION direction); /* High level functionality for handling thread contexts */ enum proj_log_level { -- cgit v1.2.3 From 5877ffac2fd4c2f89deb5b9d962a1969192856c4 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Mon, 23 Oct 2017 11:52:59 +0200 Subject: Addition of 'deformation': Kinematic grid shifting. Kinematic deformation models are used in some geodetic transformations. This commit introduces the ability to do transformations involving a gridded deformation/velocity model. For practical reasons a gridded deformation model needs to be split into two seperate files, one for the horizontal components and one for the vertical component. For this we use formats already known to PROJ.4, namely the CTable/CTable2 and the GTX formats. Grids are specified in the proj-string with +xy_grids and +z_grids. Grid values are expected to be in m/year. The kinematic part of the operation is controlled by the +t_epoch parameter, which is the central epoch of the transformation. An observation epoch is also needed. It can be specified either in the PJ_OBS input as the fourth element in the coordinate, or in the proj-string with +t_obs. If +t_obs is present in the proj-string it takes presedence over the value in the PJ_OBS coordinate. --- src/Makefile.am | 2 +- src/PJ_deformation.c | 291 +++++++++++++++++++++++++++++++++++++++++++++++++++ src/lib_proj.cmake | 1 + src/makefile.vc | 3 +- src/pj_ell_set.c | 39 ++++++- src/pj_list.h | 1 + src/proj_internal.h | 1 + 7 files changed, 335 insertions(+), 3 deletions(-) create mode 100644 src/PJ_deformation.c (limited to 'src') 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 +#include +#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 #include +#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); -- cgit v1.2.3