From 11e5226a96f2e77dac210d38b8afc11e34a2c196 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 2ad201bdb0f3408eed0aab07fe255c6ff1cd3249 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 a3fa749bc4f378d005c9e3fd809c0be25de5ffb2 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 a15c64a4f82a9996ee2d6a21de902f6bc80dc4ac 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 076247ec39e8ad1ce9554744b8c65b26f43b2017 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 9649cc099728162c4c8862b7d5a69d0dfee92c1d Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Thu, 26 Oct 2017 23:10:42 +0200 Subject: Repair and improve broken cct output routine (#631) Repair and improve broken cct output routine and do a few minor cleanups --- src/cct.c | 51 ++++++++++++++++++++++++++++----------------------- 1 file changed, 28 insertions(+), 23 deletions(-) (limited to 'src') diff --git a/src/cct.c b/src/cct.c index 72d8c6cc..90b84159 100644 --- a/src/cct.c +++ b/src/cct.c @@ -44,7 +44,7 @@ Hence, in honour of cct (the geodesist) this is cct (the program). ************************************************************************ -Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-09-19 +Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-10-26 ************************************************************************ @@ -73,7 +73,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-09-19 #include "optargpm.h" #include -#include +#include "projects.h" #include #include #include @@ -85,8 +85,6 @@ double proj_atof(const char *str); char *column (char *buf, int n); PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double fixed_time); -int print_output_line (FILE *fout, char *buf, PJ_COORD point); -int main(int argc, char **argv); @@ -237,16 +235,31 @@ int main(int argc, char **argv) { } - /* Loop over all lines of all input files */ + /* Loop over all records of all input files */ while (opt_input_loop (o, optargs_file_format_text)) { void *ret = fgets (buf, 10000, o->input); - int res; opt_eof_handler (o); if (0==ret) { fprintf (stderr, "Read error in record %d\n", (int) o->record_index); continue; } point = parse_input_line (buf, columns_xyzt, fixed_z, fixed_time); + if (HUGE_VAL==point.xyzt.x) { + char *c = column (buf, 1); + + /* if it's a comment or blank line, we reflect it */ + if (c && ((*c=='\0') || (*c=='#'))) { + fprintf (fout, "%s\n", buf); + continue; + } + + /* otherwise, it must be a syntax error */ + fprintf (fout, "# Record %d UNREADABLE: %s", (int) o->record_index, buf); + if (verbose) + fprintf (stderr, "%s: Could not parse file '%s' line %d\n", o->progname, opt_filename (o), opt_record (o)); + continue; + } + if (PJ_IO_UNITS_RADIANS==input_unit) { point.lpzt.lam = proj_torad (point.lpzt.lam); point.lpzt.phi = proj_torad (point.lpzt.phi); @@ -256,13 +269,17 @@ int main(int argc, char **argv) { point.lpzt.lam = proj_todeg (point.lpzt.lam); point.lpzt.phi = proj_todeg (point.lpzt.phi); } - res = print_output_line (fout, buf, point); - if (0==res) { - fprintf (fout, "# UNREADABLE: %s", buf); - if (verbose) - fprintf (stderr, "%s: Could not parse file '%s' line %d\n", o->progname, opt_filename (o), opt_record (o)); + + if (HUGE_VAL==point.xyzt.x) { + /* transformation error (TODO provide existing internal errmsg here) */ + fprintf (fout, "# Record %d TRANSFORMATION ERROR: %s", (int) o->record_index, buf); + continue; } + + /* Time to print the result */ + fprintf (fout, "%20.15f %20.15f %20.15f %20.15f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t); } + if (stdout != fout) fclose (fout); free (o); @@ -322,15 +339,3 @@ PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double errno = prev_errno; return result; } - - -int print_output_line (FILE *fout, char *buf, PJ_COORD point) { - char *c; - if (HUGE_VAL!=point.xyzt.x) - return fprintf (fout, "%20.15f %20.15f %20.15f %20.15f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t); - c = column (buf, 1); - /* reflect comments and blanks */ - if (c && ((*c=='\0') || (*c=='#'))) - return fprintf (fout, "%s\n", buf); - return 0; -} -- cgit v1.2.3 From 5646ff12f32adf78e2bc187e6557ce64e4e04b39 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Fri, 27 Oct 2017 22:50:40 +0200 Subject: Linguistics: Clarify this and that, here and there (#632) * Linguistics: Clarify this and that, here and there * Revert nullification of PJ_cart->fwd, inv --- src/PJ_cart.c | 32 +++++++++++++++----------------- src/PJ_helmert.c | 11 +++++++---- src/PJ_horner.c | 32 +++++++++++++++----------------- src/PJ_molodensky.c | 7 +++++++ src/PJ_pipeline.c | 21 ++++++++++++--------- src/gie.c | 9 +++++---- src/optargpm.h | 16 ++++------------ src/proj.h | 27 +++++++++++---------------- src/proj_4D_api.c | 6 ------ src/proj_strtod.c | 2 +- src/projects.h | 2 +- 11 files changed, 78 insertions(+), 87 deletions(-) (limited to 'src') diff --git a/src/PJ_cart.c b/src/PJ_cart.c index 99e1af7d..67154ab3 100644 --- a/src/PJ_cart.c +++ b/src/PJ_cart.c @@ -120,8 +120,13 @@ static double normal_radius_of_curvature (double a, double es, double phi) { /*********************************************************************/ static double geocentric_radius (double a, double b, double phi) { -/*********************************************************************/ - /* This is from WP2, but uses hypot() for potentially better numerical robustness */ +/********************************************************************* + Return the geocentric radius at latitude phi, of an ellipsoid + with semimajor axis a and semiminor axis b. + + This is from WP2, but uses hypot() for potentially better + numerical robustness +***********************************************************************/ return hypot (a*a*cos (phi), b*b*sin(phi)) / hypot (a*cos(phi), b*sin(phi)); } @@ -129,18 +134,15 @@ static double geocentric_radius (double a, double b, double phi) { /*********************************************************************/ static XYZ cartesian (LPZ geodetic, PJ *P) { /*********************************************************************/ - double N, h, lam, phi, cosphi = cos(geodetic.phi); + double N, cosphi = cos(geodetic.phi); XYZ xyz; N = normal_radius_of_curvature(P->a, P->es, geodetic.phi); - phi = geodetic.phi; - lam = geodetic.lam; - h = geodetic.z; /* HM formula 5-27 (z formula follows WP) */ - xyz.x = (N + h) * cosphi * cos(lam); - xyz.y = (N + h) * cosphi * sin(lam); - xyz.z = (N * (1 - P->es) + h) * sin(phi); + xyz.x = (N + geodetic.z) * cosphi * cos(geodetic.lam); + xyz.y = (N + geodetic.z) * cosphi * sin(geodetic.lam); + xyz.z = (N * (1 - P->es) + geodetic.z) * sin(geodetic.phi); return xyz; } @@ -149,23 +151,19 @@ static XYZ cartesian (LPZ geodetic, PJ *P) { /*********************************************************************/ static LPZ geodetic (XYZ cartesian, PJ *P) { /*********************************************************************/ - double N, b, p, theta, c, s, e2s; + double N, p, theta, c, s; LPZ lpz; /* Perpendicular distance from point to Z-axis (HM eq. 5-28) */ p = hypot (cartesian.x, cartesian.y); - /* Ancillary ellipsoidal parameters */ - b = P->b; - e2s = P->e2s; - /* HM eq. (5-37) */ - theta = atan2 (cartesian.z * P->a, p * b); + theta = atan2 (cartesian.z * P->a, p * P->b); /* HM eq. (5-36) (from BB, 1976) */ c = cos(theta); s = sin(theta); - lpz.phi = atan2 (cartesian.z + e2s*b*s*s*s, p - P->es*P->a*c*c*c); + lpz.phi = atan2 (cartesian.z + P->e2s*P->b*s*s*s, p - P->es*P->a*c*c*c); lpz.lam = atan2 (cartesian.y, cartesian.x); N = normal_radius_of_curvature (P->a, P->es, lpz.phi); @@ -176,7 +174,7 @@ static LPZ geodetic (XYZ cartesian, PJ *P) { /* by computing the height as the cartesian z value */ /* minus the geocentric radius of the Earth at the given */ /* latitude */ - double r = geocentric_radius (P->a, b, lpz.phi); + double r = geocentric_radius (P->a, P->b, lpz.phi); lpz.z = fabs (cartesian.z) - r; } else diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c index 503dc392..c9c30401 100644 --- a/src/PJ_helmert.c +++ b/src/PJ_helmert.c @@ -92,8 +92,9 @@ struct pj_opaque_helmert { #define R21 (Q->R[2][1]) #define R22 (Q->R[2][2]) +/**************************************************************************/ static void update_parameters(PJ *P) { - /************************************************************************** +/*************************************************************************** Update transformation parameters. --------------------------------- @@ -114,7 +115,8 @@ static void update_parameters(PJ *P) { of that parameter. [0] http://itrf.ign.fr/doc_ITRF/Transfo-ITRF2008_ITRFs.txt - **************************************************************************/ + +*******************************************************************************/ struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; double dt = Q->t_obs - Q->epoch; @@ -146,8 +148,9 @@ static void update_parameters(PJ *P) { return; } +/**************************************************************************/ static void build_rot_matrix(PJ *P) { - /************************************************************************** +/*************************************************************************** Build rotation matrix. ---------------------- @@ -208,7 +211,7 @@ static void build_rot_matrix(PJ *P) { as published, the "transpose" option provides the ability to switch between the conventions. - ***************************************************************************/ +***************************************************************************/ struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque; double f, t, p; /* phi/fi , theta, psi */ diff --git a/src/PJ_horner.c b/src/PJ_horner.c index d6d2c51c..a9a2c922 100644 --- a/src/PJ_horner.c +++ b/src/PJ_horner.c @@ -1,19 +1,3 @@ -#define PJ_LIB__ -#include "proj_internal.h" -#include "projects.h" -#include -#include -#include -#include - -PROJ_HEAD(horner, "Horner polynomial evaluation"); - -/* make horner.h interface with proj's memory management */ -#define horner_dealloc(x) pj_dealloc(x) -#define horner_calloc(n,x) pj_calloc(n,x) - -/* The next few hundred lines is mostly cut-and-paste from the horner.h header library */ - /*********************************************************************** Interfacing to a classic piece of geodetic software @@ -91,6 +75,20 @@ PROJ_HEAD(horner, "Horner polynomial evaluation"); * *****************************************************************************/ +#define PJ_LIB__ +#include "proj_internal.h" +#include "projects.h" +#include +#include +#include +#include + +PROJ_HEAD(horner, "Horner polynomial evaluation"); + +/* make horner.h interface with proj's memory management */ +#define horner_dealloc(x) pj_dealloc(x) +#define horner_calloc(n,x) pj_calloc(n,x) + struct horner; typedef struct horner HORNER; @@ -399,7 +397,7 @@ static int parse_coefs (PJ *P, double *coefs, char *param, int ncoefs) { buf = pj_calloc (strlen (param) + 2, sizeof(char)); if (0==buf) { - proj_log_error (P, "Horner: Out of core"); + proj_log_error (P, "Horner: No memory left"); return 0; } diff --git a/src/PJ_molodensky.c b/src/PJ_molodensky.c index f09f07cd..50423fd8 100644 --- a/src/PJ_molodensky.c +++ b/src/PJ_molodensky.c @@ -90,6 +90,13 @@ static double RM (double a, double es, double phi) { E.J Krakiwsky & D.B. Thomson, 1974, GEODETIC POSITION COMPUTATIONS, + + Fredericton NB, Canada: + University of New Brunswick, + Department of Geodesy and Geomatics Engineering, + Lecture Notes No. 39, + 99 pp. + http://www2.unb.ca/gge/Pubs/LN39.pdf **********************************************************/ diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c index 137fdec8..ab2d3420 100644 --- a/src/PJ_pipeline.c +++ b/src/PJ_pipeline.c @@ -45,9 +45,10 @@ It is a very powerful concept, that increases the range of relevance of the proj.4 system substantially. It is, however, not a particularly intrusive - addition to the code base: The implementation is by and large completed by - adding an extra projection called "pipeline" (i.e. this file), which - handles all business. + addition to the PROJ.4 code base: The implementation is by and large completed + by adding an extra projection called "pipeline" (i.e. this file), which + handles all business, and a small amount of added functionality in the + pj_init code, implementing support for multilevel, embedded pipelines. Syntactically, the pipeline system introduces the "+step" keyword (which indicates the start of each transformation step), the "+omit_fwd" and @@ -63,10 +64,11 @@ Where indicate the Helmert arguments: 3 translations (+x=..., +y=..., +z=...), 3 rotations (+rx=..., +ry=..., +rz=...) and a scale factor (+s=...). - Following geodetic conventions, the rotations are given in Milliarcseconds, + Following geodetic conventions, the rotations are given in arcseconds, and the scale factor is given as parts-per-million. - [1] B. W. Kernighan & P. J. Plauger: Software tools. Addison-Wesley, 1976, 338 pp. + [1] B. W. Kernighan & P. J. Plauger: Software tools. + Reading, Massachusetts, Addison-Wesley, 1976, 338 pp. ******************************************************************************** @@ -156,10 +158,11 @@ static LP pipeline_reverse (XY xyz, PJ *P); The P->left and P->right elements indicate isomorphism. - For classic proj style projections, they both have the - value PJ_IO_UNITS_CLASSIC (default initialization), indicating - that the forward driver expects angular input coordinates, and - provides linear output coordinates. + For classic proj style projections, P->left has the value + PJ_IO_UNITS_RADIANS, while P->right has the value + PJ_IO_UNITS_CLASSIC, indicating that the forward driver expects + angular input coordinates, and provides linear output coordinates, + scaled by the P->a semimajor axis length. Newer projections may set P->left and P->right to either PJ_IO_UNITS_METERS, PJ_IO_UNITS_RADIANS or PJ_IO_UNITS_ANY, diff --git a/src/gie.c b/src/gie.c index f0850b61..269a7239 100644 --- a/src/gie.c +++ b/src/gie.c @@ -52,11 +52,12 @@ but all of them *just special functions*, and not particularly more special than the sin(), cos(), tan(), and hypot() already available in the C standard library. -And hence, *they should not be particularly much harder to use*, for a -programmer, than the sin()s, tan()s and hypot()s so readily available. +And hence, according to Gerald, *they should not be particularly much +harder to use*, for a programmer, than the sin()s, tan()s and hypot()s +so readily available. Gerald's ingenuity also showed in the implementation of the vision, -where he devised a highly comprehensible, yet simple, system of key-value +where he devised a comprehensive, yet simple, system of key-value pairs for parameterising a map projection, and the highly flexible PJ struct, storing run-time compiled versions of those key-value pairs, hence making a map projection function call, pj_fwd(PJ, point), as easy @@ -332,7 +333,7 @@ char *column (char *buf, int n) { } - +/* interpret as a numeric followed by a linear decadal prefix - return the properly scaled numeric */ static double strtod_scaled (char *args, double default_scale) { double s; char *endp = args; diff --git a/src/optargpm.h b/src/optargpm.h index c47dfec1..23e375f8 100644 --- a/src/optargpm.h +++ b/src/optargpm.h @@ -10,8 +10,8 @@ For PROJ.4 command line programs, we have a somewhat complex option decoding situation, since we have to navigate in a cocktail of classic single letter style options, prefixed by "-", GNU style long options -prefixwd by "--", transformation specification elements prefixed by "+", -and input file names prefixed by "" nothing. +prefixed by "--", transformation specification elements prefixed by "+", +and input file names prefixed by "" (i.e. nothing). Hence, classic getopt.h style decoding does not cut the mustard, so this is an attempt to catch up and chop the ketchup. @@ -48,8 +48,8 @@ Operator specs are +proj=utm +zone=32 +ellps=GRS80 Operands are bar baz -While claiming neither to save the world, nor to hint at the "shape of -jazz to come", at least optargpm has shown useful in constructing cs2cs +While neither claiming to save the world, nor to hint at the "shape of +jazz to come", at least optargpm has shown useful in constructing cs2cs style transformation filters. Supporting a wide range of option syntax, the getoptpm API is somewhat @@ -424,14 +424,6 @@ OPTARGS *opt_parse (int argc, char **argv, const char *flags, const char *keys, o->argc = argc; o->argv = argv; o->progname = opt_strip_path (argv[0]); -/* o->progname = argv[0]; - last_path_delim = strrchr (argv[0], '\\'); - if (last_path_delim > o->progname) - o->progname = last_path_delim; - last_path_delim = strrchr (argv[0], '/'); - if (last_path_delim > o->progname) - o->progname = last_path_delim; -*/ /* Reset all flags */ for (i = 0; i < (int) strlen (flags); i++) diff --git a/src/proj.h b/src/proj.h index ed885091..dd0672f9 100644 --- a/src/proj.h +++ b/src/proj.h @@ -53,26 +53,21 @@ * it remains as-is for now). * * THIRD, I try to eliminate implicit type punning. Hence this API - * introduces the PJ_OBS ("observation") data type, for generic - * coordinate and handling of ancillary data. + * introduces the PJ_COORD union data type, for generic 4D coordinate + * handling. * - * It includes the PJ_COORD and PJ_TRIPLET unions making it possible - * to make explicit the previously used "implicit type punning", where - * a XY is turned into a LP by re#defining both as UV, behind the back - * of the user. + * PJ_COORD makes it possible to make explicit the previously used + * "implicit type punning", where a XY is turned into a LP by + * re#defining both as UV, behind the back of the user. * * The PJ_COORD union is used for storing 1D, 2D, 3D and 4D coordinates. - * The PJ_TRIPLET union is used for storing any set of up to 3 related - * observations. At the application code level, the names of these - * unions will usually not be used - they will only be accessed via - * their tag names in the PJ_OBS data type. * * The bare essentials API presented here follows the PROJ.4 * convention of sailing the coordinate to be reprojected, up on * the stack ("call by value"), and symmetrically returning the - * result on the stack. Although the PJ_OBS object is 4 times - * as large as the traditional XY and LP objects, timing results - * have shown the overhead to be very reasonable. + * result on the stack. Although the PJ_COORD object is twice as large + * as the traditional XY and LP objects, timing results have shown the + * overhead to be very reasonable. * * Contexts and thread safety * -------------------------- @@ -82,8 +77,8 @@ * context subsystem is unavoidable in a multi-threaded world. * Hence, instead of hiding it away, we move it into the limelight, * highly recommending (but not formally requiring) the bracketing - * with calls to proj_context_create(...)/proj_context_destroy() of - * any code block calling PROJ.4 functions. + * of any code block calling PROJ.4 functions with calls to + * proj_context_create(...)/proj_context_destroy() * * Legacy single threaded code need not do anything, but *may* * implement a bit of future compatibility by using the backward @@ -426,7 +421,7 @@ const PJ_ELLPS *proj_list_ellps(void); const PJ_UNITS *proj_list_units(void); const PJ_PRIME_MERIDIANS *proj_list_prime_meridians(void); -/* These are trivial, and while occasionaly useful in real code, primarily here to */ +/* These are trivial, and while occasionally useful in real code, primarily here to */ /* simplify demo code, and in acknowledgement of the proj-internal discrepancy between */ /* angular units expected by classical proj, and by Charles Karney's geodesics subsystem */ double proj_torad (double angle_in_degrees); diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c index b1aa3883..5484b650 100644 --- a/src/proj_4D_api.c +++ b/src/proj_4D_api.c @@ -3,12 +3,6 @@ * 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 * ****************************************************************************** diff --git a/src/proj_strtod.c b/src/proj_strtod.c index 757dfaf6..fa683465 100644 --- a/src/proj_strtod.c +++ b/src/proj_strtod.c @@ -17,7 +17,7 @@ In the early versions of proj, iirc, a gnu version of strtod was used, mostly to work around cases where the same system library was used for C and Fortran linking, hence making strtod accept "D" and "d" as exponentiation indicators, following Fortran Double Precision constant -syntax. This broke the proj angular syntax accepting a "d" to mean +syntax. This broke the proj angular syntax, accepting a "d" to mean "degree": 12d34'56", meaning 12 degrees 34 minutes and 56 seconds. With an explicit MIT licence, PROJ.4 could not include GPL code any diff --git a/src/projects.h b/src/projects.h index b9d88cf3..1fe595bc 100644 --- a/src/projects.h +++ b/src/projects.h @@ -186,7 +186,7 @@ struct PJ_REGION_S; typedef struct PJ_REGION_S PJ_Region; typedef struct ARG_list paralist; /* parameter list */ enum pj_io_units { - PJ_IO_UNITS_CLASSIC = 0, /* LEFT: Radians RIGHT: Scaled meters */ + PJ_IO_UNITS_CLASSIC = 0, /* Scaled meters (right) */ PJ_IO_UNITS_METERS = 1, /* Meters */ PJ_IO_UNITS_RADIANS = 2 /* Radians */ }; -- cgit v1.2.3 From 234d712d2cf1145ba84c48f51511f92b461e9707 Mon Sep 17 00:00:00 2001 From: Ilya Oshchepkov Date: Sat, 28 Oct 2017 14:13:05 +0300 Subject: Add ellipsoids for the Russian coordinate systems PZ-90 and GSK-2011 --- src/pj_ellps.c | 2 ++ 1 file changed, 2 insertions(+) (limited to 'src') diff --git a/src/pj_ellps.c b/src/pj_ellps.c index 32fbbd42..2a2c5d3b 100644 --- a/src/pj_ellps.c +++ b/src/pj_ellps.c @@ -15,6 +15,7 @@ pj_ellps[] = { {"andrae", "a=6377104.43", "rf=300.0", "Andrae 1876 (Den., Iclnd.)"}, {"aust_SA", "a=6378160.0", "rf=298.25", "Australian Natl & S. Amer. 1969"}, {"GRS67", "a=6378160.0", "rf=298.2471674270", "GRS 67(IUGG 1967)"}, +{"GSK2011", "a=6378136.5", "rf=298.2564151", "GSK-2011"}, {"bessel", "a=6377397.155", "rf=299.1528128", "Bessel 1841"}, {"bess_nam", "a=6377483.865", "rf=299.1528128", "Bessel 1841 (Namibia)"}, {"clrk66", "a=6378206.4", "b=6356583.8", "Clarke 1866"}, @@ -40,6 +41,7 @@ pj_ellps[] = { {"mprts", "a=6397300.", "rf=191.", "Maupertius 1738"}, {"new_intl", "a=6378157.5", "b=6356772.2", "New International 1967"}, {"plessis", "a=6376523.", "b=6355863.", "Plessis 1817 (France)"}, +{"PZ90", "a=6378136.0", "rf=298.25784", "PZ-90"}, {"SEasia", "a=6378155.0", "b=6356773.3205", "Southeast Asia"}, {"walbeck", "a=6376896.0", "b=6355834.8467", "Walbeck"}, {"WGS60", "a=6378165.0", "rf=298.3", "WGS 60"}, -- cgit v1.2.3 From 0d0beff91ddfc2cc4d195a141524a139f3afb756 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Sat, 28 Oct 2017 15:42:00 +0200 Subject: Allow nested pipelines. (#629) Allow nested pipelines when wrapped in +init's. The previous behaviour was to quit pipeline initialization when encountering a nested pipeline definition. With this commit that behaviour is changed so that it is possible to nest pipelines as long as they are defined elsewhere in a init-file. This is useful in init-files where steps in complicated transformations can be grouped in "sub-pipelines". These "sub-pipelines" can then be used as individual steps in a larger and more complicated pipeline. Nested pipelines are governed by the following rules: 1. You can't have more than one literal +proj=pipeline in a proj-string 2. Pipelines can be nested if they are wrapped up in a +init 3. More than one +init is disallowed in non-pipeline proj-strings 4. +inits are expanded as late as possible, that is they will only be expanded in single operations (that can be a part of a pipeline) --- src/PJ_pipeline.c | 2 +- src/pj_init.c | 96 +++++++++++++++++++++++++++++++++---------------------- src/pj_strerrno.c | 3 +- src/projects.h | 1 + 4 files changed, 62 insertions(+), 40 deletions(-) (limited to 'src') diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c index ab2d3420..f1411cb4 100644 --- a/src/PJ_pipeline.c +++ b/src/PJ_pipeline.c @@ -380,7 +380,7 @@ PJ *PROJECTION(pipeline) { if (0==strcmp ("proj=pipeline", argv[i])) { if (-1 != i_pipeline) { - proj_log_error (P, "Pipeline: Nesting invalid"); + proj_log_error (P, "Pipeline: Nesting only allowed when child pipelines are wrapped in +init's"); return destructor (P, PJD_ERR_MALFORMED_PIPELINE); /* ERROR: nested pipelines */ } i_pipeline = i; diff --git a/src/pj_init.c b/src/pj_init.c index 704a8b55..2a945397 100644 --- a/src/pj_init.c +++ b/src/pj_init.c @@ -381,8 +381,9 @@ pj_init_plus_ctx( projCtx ctx, const char *definition ) if( argc+1 == MAX_ARG ) { - pj_ctx_set_errno( ctx, -44 ); - goto bum_call; + pj_dalloc( defn_copy ); + pj_ctx_set_errno( ctx, PJD_ERR_UNPARSEABLE_CS_DEF ); + return 0; } argv[argc++] = defn_copy + i + 1; @@ -410,9 +411,7 @@ pj_init_plus_ctx( projCtx ctx, const char *definition ) /* perform actual initialization */ result = pj_init_ctx( ctx, argc, argv ); -bum_call: pj_dalloc( defn_copy ); - return result; } @@ -441,6 +440,8 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) { int i; int found_def = 0; PJ *PIN = 0; + int n_pipelines = 0; + int n_inits = 0; if (0==ctx) ctx = pj_get_default_ctx (); @@ -452,36 +453,45 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) { pj_ctx_set_errno (ctx, PJD_ERR_NO_ARGS); return 0; } - + + /* count occurrences of pipelines and inits */ + for (i = 0; i < argc; ++i) { + if (!strcmp (argv[i], "+proj=pipeline") || !strcmp(argv[i], "proj=pipeline") ) + n_pipelines++; + if (!strncmp (argv[i], "+init=", 6) || !strncmp(argv[i], "init=", 5)) + n_inits++; + } + + /* can't have nested pipeline directly */ + if (n_pipelines > 1) { + pj_ctx_set_errno (ctx, PJD_ERR_MALFORMED_PIPELINE); + return 0; + } + + /* don't allow more than one +init in non-pipeline operations */ + if (n_pipelines == 0 && n_inits > 1) { + pj_ctx_set_errno (ctx, PJD_ERR_TOO_MANY_INITS); + return 0; + } + /* put arguments into internal linked list */ start = curr = pj_mkparam(argv[0]); if (!curr) return pj_dealloc_params (ctx, start, ENOMEM); - /* build parameter list and expand +init's. Does not take care of a single +init. */ for (i = 1; i < argc; ++i) { curr->next = pj_mkparam(argv[i]); if (!curr->next) return pj_dealloc_params (ctx, start, ENOMEM); - - /* check if +init present */ - if (pj_param(ctx, curr, "tinit").i) { - found_def = 0; - curr = get_init(ctx, &curr, curr->next, pj_param(ctx, curr, "sinit").s, &found_def); - if (!curr) - return pj_dealloc_params (ctx, start, PJD_ERR_NO_ARGS); - - if (!found_def) - return pj_dealloc_params (ctx, start, PJD_ERR_NO_OPTION_IN_INIT_FILE); - - } else { - curr = curr->next; - } + curr = curr->next; } - /* in case the parameter list only consist of a +init parameter - it is expanded here (will not be handled in the above loop). */ - if (pj_param(ctx, start, "tinit").i && argc == 1) { + + /* Only expand +init's in non-pipeline operations. +init's in pipelines are */ + /* expanded in the individual pipeline steps during pipeline initialization. */ + /* Potentially this leads to many nested pipelines, which shouldn't be a */ + /* problem when +inits are expanded as late as possible. */ + if (pj_param(ctx, start, "tinit").i && n_pipelines == 0) { found_def = 0; curr = get_init(ctx, &start, curr, pj_param(ctx, start, "sinit").s, &found_def); if (!curr) @@ -489,10 +499,10 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) { if (!found_def) return pj_dealloc_params (ctx, start, PJD_ERR_NO_OPTION_IN_INIT_FILE); } - + if (ctx->last_errno) return pj_dealloc_params (ctx, start, ctx->last_errno); - + /* find projection selection */ if (!(name = pj_param(ctx, start, "sproj").s)) return pj_default_destructor (PIN, PJD_ERR_PROJ_NOT_NAMED); @@ -500,10 +510,11 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) { if (!s) return pj_dealloc_params (ctx, start, PJD_ERR_UNKNOWN_PROJECTION_ID); - - /* set defaults, unless inhibited */ - if (!(pj_param(ctx, start, "bno_defs").i)) - curr = get_defaults(ctx,&start, curr, name); + + /* set defaults, unless inhibited or we are initializing a pipeline */ + if (!(pj_param(ctx, start, "bno_defs").i) && n_pipelines == 0) + curr = get_defaults(ctx,&start, curr, name); + proj = (PJ *(*)(PJ *)) pj_list[i].proj; /* allocate projection structure */ @@ -511,6 +522,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) { if (0==PIN) return pj_dealloc_params (ctx, start, ENOMEM); + PIN->ctx = ctx; PIN->params = start; PIN->is_latlong = 0; @@ -526,13 +538,21 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) { PIN->vgridlist_geoid_count = 0; /* set datum parameters */ - if (pj_datum_set(ctx, start, PIN)) + if (pj_datum_set(ctx, start, PIN)) return pj_default_destructor (PIN, PJD_ERR_MISSING_ARGS); - /* set ellipsoid/sphere parameters */ - if (pj_ell_set(ctx, start, &PIN->a, &PIN->es)) { - pj_log (ctx, PJ_LOG_DEBUG_MINOR, "pj_init_ctx: Must specify ellipsoid or sphere"); - return pj_default_destructor (PIN, PJD_ERR_MISSING_ARGS); + /* set ellipsoid/sphere parameters. If we are initializing a pipeline we */ + /* override ellipsoid-setter, since it also adds a "+ellps=WGS84" to the */ + /* parameter list which creates problems when the pipeline driver passes */ + /* the global parameters on the it's children. */ + if (n_pipelines > 0) { + PIN->a = 6378137.0; + PIN->es = .00669438002290341575; + } else { + if (pj_ell_set(ctx, start, &PIN->a, &PIN->es)) { + pj_log (ctx, PJ_LOG_DEBUG_MINOR, "pj_init_ctx: Must specify ellipsoid or sphere"); + return pj_default_destructor (PIN, PJD_ERR_MISSING_ARGS); + } } PIN->a_orig = PIN->a; @@ -606,7 +626,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) { if( !(fabs(PIN->long_wrap_center) < 10 * M_TWOPI) ) return pj_default_destructor (PIN, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT); } - + /* axis orientation */ if( (pj_param(ctx, start,"saxis").s) != NULL ) { @@ -619,7 +639,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) { || strchr( axis_legal, axis_arg[1] ) == NULL || strchr( axis_legal, axis_arg[2] ) == NULL) return pj_default_destructor (PIN, PJD_ERR_AXIS); - + /* it would be nice to validate we don't have on axis repeated */ strcpy( PIN->axis, axis_arg ); } @@ -643,7 +663,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) { PIN->k0 = pj_param(ctx, start, "dk").f; else PIN->k0 = 1.; - if (PIN->k0 <= 0.) + if (PIN->k0 <= 0.) return pj_default_destructor (PIN, PJD_ERR_K_LESS_THAN_ZERO); /* set units */ @@ -704,7 +724,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) { && *next_str == '\0' ) value = name; - if (!value) + if (!value) return pj_default_destructor (PIN, PJD_ERR_UNKNOWN_PRIME_MERIDIAN); PIN->from_greenwich = dmstor_ctx(ctx,value,NULL); } diff --git a/src/pj_strerrno.c b/src/pj_strerrno.c index 89a5a025..fe137b87 100644 --- a/src/pj_strerrno.c +++ b/src/pj_strerrno.c @@ -62,7 +62,8 @@ pj_err_list[] = { "missing required arguments", /* -54 */ "lat_0 = 0", /* -55 */ "ellipsoidal usage unsupported", /* -56 */ - + "only one +init allowed for non-pipeline operations", /* -57 */ + /* When adding error messages, remember to update ID defines in projects.h, and transient_error array in pj_transform */ }; diff --git a/src/projects.h b/src/projects.h index 1fe595bc..6ea0c914 100644 --- a/src/projects.h +++ b/src/projects.h @@ -532,6 +532,7 @@ struct FACTORS { #define PJD_ERR_MISSING_ARGS -54 #define PJD_ERR_LAT_0_IS_ZERO -55 #define PJD_ERR_ELLIPSOIDAL_UNSUPPORTED -56 +#define PJD_ERR_TOO_MANY_INITS -57 struct projFileAPI_t; -- cgit v1.2.3