diff --git a/examples/pj_obs_api_mini_demo.c b/examples/pj_obs_api_mini_demo.c new file mode 100644 index 0000000000..855c37f803 --- /dev/null +++ b/examples/pj_obs_api_mini_demo.c @@ -0,0 +1,66 @@ +/******************************************************************************* + Tiny test of an evolving new API, demonstrating simple examples of + 2D and 3D transformations. + + The main transformation setup object is PJ, well known from the two + former proj APIs (projects.h and proj_api.h) + + The main data objects PJ_COORD and PJ_OBS are new, but encapsulate + the older LP, XY etc. objects in a framework for storing a 3D + observation taken at a 4D point in space-time, and including some + additional metadata (e.g. a serial number or an epsg code). + + PJ_OBS and PJ_COORD use unions to enforce explicit statement of what + kind of coordinates are expected at a given spot in the code, where + the old API uses type punning, implicitly assuming that "you know what + you do". For backward compatibility, the new API is not really type + safe in the sense that you cannot use a cartesian coordinate where a + geographic is expected - but it makes it possible to explicitly state + in the code whet the programmer expected and intended. + + The 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. + + A series of experiments have, however, shown that they, for (mostly) + historical reasons, are very hard to eliminate totally. But we have + reduced their API surface presence to a constructor and a destructor, + plus an extra argument to the PJ constructor, pj_create(). + + For single threaded programs, the calls to the context constructor + and destructor may be left out, and the default context selected + by passing a null-pointer to pj_create. + + Thomas Knudsen, 2016-10-30/2017-07-06 +*******************************************************************************/ +#include +#include + +int main (void) { + PJ_CONTEXT *C; + PJ *P; + PJ_COORD a, b; + + /* or you may set C=0 if you are sure you will use PJ objects from only one thread */ + C = proj_context_create(); + + P = proj_create (C, "+proj=utm +zone=32 +ellps=GRS80"); + if (0==P) + return puts ("Oops"), 0; + + /* a coordinate union representing Copenhagen: 55d N, 12d E */ + /* note: PROJ.4 works in radians, hence the proj_torad() calls */ + a = proj_coord (proj_torad(12), proj_torad(55), 0, 0); + + /* transform to UTM zone 32, then back to geographical */ + /* note the use of union selectors to indicate what kind of coordinates are expected */ + b = proj_trans_coord (P, PJ_FWD, a); + printf ("easting: %g, northing: %g\n", b.en.e, b.en.n); + b = proj_trans_coord (P, PJ_INV, b); + printf ("longitude: %g, latitude: %g\n", b.lp.lam, b.lp.phi); + + /* Clean up */ + proj_destroy (P); + proj_context_destroy (C); /* may be omitted in the single threaded case */ + return 0; +} diff --git a/examples/pj_obs_api_test.c b/examples/pj_obs_api_test.c deleted file mode 100644 index bb792c62e6..0000000000 --- a/examples/pj_obs_api_test.c +++ /dev/null @@ -1,243 +0,0 @@ -/******************************************************************************* - Tiny test of an evolving new API, demonstrating simple examples of - 2D and 3D transformations. - - The main transformation setup object is PJ, well known from the two - former proj APIs (projects.h and proj_api.h) - - The main data object PJ_OBS is new, but encapsulates the older - LP, XY etc. objects in a framework for storing a 3D observation taken at - a 4D point in space-time, and including some additional metadata (e.g. - a serial number or an epsg code). - - PJ_OBS uses unions to enforce explicit statement of what kind of - coordinates are expected at a given spot in the code. - - The primary elements of the API are: - - pj_create (char *desc): - Create a new PJ object from a text description. - pj_trans (PJ *P, int direction, PJ_OBS obs): - Forward or inverse transformation of obs. - pj_error (PJ *P): - Check error status of P. - pj_free (PJ *P): - Clean up when done with the transformation object. - - These 4 functions should cover most of the common use cases. Additional - functionality for handling thread contexts are provided by the functions - pj_debug_set, pj_error_set, pj_log_set, pj_context_renew, - pj_context_inherit, pj_context_free, pj_fileapi_set. - - The 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 only - visible through their associated PJs, and the number of functions - supporting them is limited. - - Compilation example: - gcc -pedantic -W -Wall -Wextra -I../src -o pj_proj_test pj_proj_test.c -L../../../build/proj.4/lib -lproj_4_9 - - Thomas Knudsen, 2016-10-30/11-03 -*******************************************************************************/ -#include - -int pj_pipeline_selftest (void); - -int main (void) { - PJ *P; - PJ_OBS a, b; - char *args[3] = {"proj=utm", "zone=32", "ellps=GRS80"}; - int err; - double dist; - XY cph_utm32; - - /* Log everything libproj offers to log for you */ - pj_log_level (0, PJ_LOG_TRACE); - - /* An utm projection on the GRS80 ellipsoid */ - P = pj_create ("+proj=utm +zone=32 +ellps=GRS80"); - if (0==P) - return puts ("Oops"), 0; - - /* Clean up */ - pj_free (P); - - /* Same projection, now using argc/argv style initialization */ - P = pj_create_argv (3, args); - if (0==P) - return puts ("Oops"), 0; - - /* zero initialize everything, then set (longitude, latitude) to (12, 55) */ - a = pj_obs_null (); - /* a.coo.lp: The coordinate part of a, interpreted as a classic LP pair */ - a.coo.lp.lam = TORAD(12); - a.coo.lp.phi = TORAD(55); - - /* Forward projection */ - b = pj_trans (P, PJ_FWD, a); - printf ("FWD: %15.9f %15.9f\n", b.coo.enh.e, b.coo.enh.n); - cph_utm32 = b.coo.xy; - - /* Inverse projection */ - a = pj_trans (P, PJ_INV, b); - printf ("INV: %15.9f %15.9f\n", TODEG(a.coo.lpz.lam), TODEG(a.coo.lpz.phi)); - - /* Null projection */ - a = pj_trans (P, PJ_IDENT, a); - printf ("IDENT: %15.9f %15.9f\n", TODEG(a.coo.lpz.lam), TODEG(a.coo.lpz.phi)); - - /* Forward again, to get two linear items for comparison */ - a = pj_trans (P, PJ_FWD, a); - printf ("FWD: %15.9f %15.9f\n", b.coo.enh.e, b.coo.enh.n); - - dist = pj_xy_dist (a.coo.xy, b.coo.xy); - printf ("Roundtrip deviation, (nm): %15.9f\n", dist*1e9); - - /* should be identical - checking whether null-init is done */ - dist = pj_xyz_dist (a.coo.xyz, b.coo.xyz); - printf ("Roundtrip deviation, (nm): %15.9f\n", dist*1e9); - - /* Invalid projection */ - a = pj_trans (P, 42, a); - if (a.coo.lpz.lam!=DBL_MAX) - printf ("%15.9f %15.9f\n", a.coo.lpz.lam, a.coo.lpz.phi); - err = pj_err_level (P, PJ_ERR_TELL); - printf ("pj_err_level: %d\n", err); - - /* Clean up */ - pj_free (P); - - /* Now do some 3D transformations */ - P = pj_create ("+proj=cart +ellps=GRS80"); - if (0==P) - return puts ("Oops"), 0; - - /* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */ - a = b = pj_obs_null (); - a.coo.lpz.lam = TORAD(12); - a.coo.lpz.phi = TORAD(55); - a.coo.lpz.z = 100; - - /* Forward projection: 3D-Cartesian-to-Ellipsoidal */ - b = pj_trans (P, PJ_FWD, a); - printf ("FWD: %15.9f %15.9f %15.9f\n", b.coo.xyz.x, b.coo.xyz.y, b.coo.xyz.z); - - /* Check roundtrip precision for 10000 iterations each way */ - dist = pj_roundtrip (P, PJ_FWD, 10000, a); - printf ("Roundtrip deviation, fwd (nm): %15.9f\n", dist*1e9*1e5); - dist = pj_roundtrip (P, PJ_INV, 10000, b); - printf ("Roundtrip deviation, inv (nm): %15.9f\n", dist*1e9); - - /* Inverse projection: Ellipsoidal-to-3D-Cartesian */ - b = pj_trans (P, PJ_INV, b); - printf ("INV: %15.9f %15.9f %15.9f\n", TODEG(b.coo.lpz.lam), TODEG(b.coo.lpz.phi), b.coo.lpz.z); - - /* Move p to another context */ - pj_context_renew (P); - b = pj_trans (P, PJ_FWD, b); - printf ("CTX1: %15.9f %15.9f %15.9f\n", b.coo.xyz.x, b.coo.xyz.y, b.coo.xyz.z); - - /* Move it back to the default context */ - pj_context_free (P); - b = pj_trans (P, PJ_INV, b); - printf ("CTX0: %15.9f %15.9f %15.9f\n", TODEG(b.coo.lpz.lam), TODEG(b.coo.lpz.phi), b.coo.lpz.z); - - pj_free (P); - - /*************************************************************************** - - P I P E L I N E T E S T S - - ***************************************************************************/ - - /* forward-reverse geo->utm->geo */ - P = pj_create ( - "+proj=pipeline +ellps=GRS80 +zone=32 +step " - "+proj=utm +step " - "+proj=utm +inv" - ); - if (0==P) - return puts ("Oops"), 0; - /* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */ - a = b = pj_obs_null (); - a.coo.lpz.lam = TORAD(12); - a.coo.lpz.phi = TORAD(55); - a.coo.lpz.z = 100; - printf ("PRE: %15.9f %15.9f\n", a.coo.lpz.lam, a.coo.lpz.phi); - - /* Forward projection */ - b = pj_trans (P, PJ_FWD, a); - printf ("FWD: %15.9f %15.9f\n", TODEG(b.coo.lpz.lam), TODEG(b.coo.lpz.phi)); - - /* Inverse projection */ - a = pj_trans (P, PJ_INV, b); - printf ("INV: %15.9f %15.9f\n", TODEG(a.coo.lpz.lam), TODEG(a.coo.lpz.phi)); - - pj_free (P); - - - /* And now the back-to-back situation utm->geo->utm */ - P = pj_create ( - "+proj=pipeline +ellps=GRS80 +zone=32 +step " - "+proj=utm +inv +step " - "+proj=utm"); - if (0==P) - return puts ("Oops"), 0; - - /* zero initialize everything, then set (easting, northing) to utm(12, 55) */ - a = b = pj_obs_null (); - a.coo.xy = cph_utm32; - printf ("PRE: %15.9f %15.9f\n", a.coo.xy.x, a.coo.xy.y); - - /* Forward projection */ - b = pj_trans (P, PJ_FWD, a); - printf ("FWD: %15.9f %15.9f\n", b.coo.xy.x, b.coo.xy.y); - - /* Inverse projection */ - a = pj_trans (P, PJ_INV, b); - printf ("INV: %15.9f %15.9f\n", a.coo.xy.x, a.coo.xy.y); - - pj_free (P); - - - - /* Finally testing a corner case: A rather pointless one-step pipeline geo->utm */ - P = pj_create ("+proj=pipeline +ellps=GRS80 +zone=32 +step +proj=utm"); - if (0==P) - return puts ("Oops"), 0; - - /* zero initialize everything, then set (easting, northing) to utm(12, 55) */ - a = b = pj_obs_null (); - a.coo.lpz.lam = TORAD(12); - a.coo.lpz.phi = TORAD(55); - printf ("PRE: %15.9f %15.9f\n", TODEG(a.coo.lp.lam), TODEG(a.coo.lp.phi)); - printf ("EXP: %15.9f %15.9f\n", cph_utm32.x, cph_utm32.y); - - /* Forward projection */ - b = pj_trans (P, PJ_FWD, a); - printf ("FWD: %15.9f %15.9f\n", b.coo.xy.x, b.coo.xy.y); - - /* Inverse projection */ - a = pj_trans (P, PJ_INV, b); - printf ("INV: %15.9f %15.9f\n", TODEG(a.coo.lp.lam), TODEG(a.coo.lp.phi)); - - /* Geodesic distance between two points with angular 2D coordinates */ - a.coo.lp.lam = TORAD(12); - a.coo.lp.phi = TORAD(60); - b.coo.lp.lam = TORAD(12); - b.coo.lp.phi = TORAD(61); - dist = pj_lp_dist (P, a.coo.lp, b.coo.lp); - printf ("1 deg at 60N: %15.9f\n", dist); - - a.coo.lp.lam = TORAD(12); - a.coo.lp.phi = TORAD(0.); - b.coo.lp.lam = TORAD(12); - b.coo.lp.phi = TORAD(1.); - dist = pj_lp_dist (P, a.coo.lp, b.coo.lp); - printf ("1 deg at Equator: %15.9f\n", dist); - - pj_free (P); - pj_pipeline_selftest (); - return 0; -} diff --git a/src/Makefile.am b/src/Makefile.am index 25a7eaf402..f261541298 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -32,7 +32,7 @@ lib_LTLIBRARIES = libproj.la libproj_la_LDFLAGS = -no-undefined -version-info 12:0:0 libproj_la_SOURCES = \ - pj_list.h \ + pj_list.h proj_internal.h\ PJ_aeqd.c PJ_gnom.c PJ_laea.c PJ_mod_ster.c \ PJ_nsper.c PJ_nzmg.c PJ_ortho.c PJ_stere.c PJ_sterea.c \ PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c PJ_isea.c \ @@ -75,7 +75,7 @@ libproj_la_SOURCES = \ pj_strtod.c \ \ pj_obs_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c \ - PJ_vgridshift.c PJ_hgridshift.c PJ_unitconvert.c + PJ_vgridshift.c PJ_hgridshift.c PJ_unitconvert.c pj_internal.c install-exec-local: rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT) diff --git a/src/PJ_cart.c b/src/PJ_cart.c index f3379a1b86..26acac2dce 100644 --- a/src/PJ_cart.c +++ b/src/PJ_cart.c @@ -41,7 +41,7 @@ *****************************************************************************/ #define PJ_LIB__ -#include +#include "proj_internal.h" #include #include #include @@ -232,117 +232,196 @@ int pj_cart_selftest (void) {return 0;} #else /* Testing quite a bit of the pj_obs_api as a side effect (inspired by pj_obs_api_test.c) */ int pj_cart_selftest (void) { + PJ_CONTEXT *ctx; PJ *P; - PJ_OBS a, b; + PJ_OBS a, b, obs[2]; int err; - double dist; + size_t n, sz; + double dist, h, t; char *args[3] = {"proj=utm", "zone=32", "ellps=GRS80"}; - + char *arg = {" +proj=utm +zone=32 +ellps=GRS80"}; + char *arg_def; /* An utm projection on the GRS80 ellipsoid */ - P = pj_create ("+proj=utm +zone=32 +ellps=GRS80"); + P = proj_create (0, arg); if (0==P) return 1; + /* note arg is handcrafted to go undisturbed through get def reconstruction */ + arg_def = proj_definition_retrieve (P); + if (0!=strcmp(arg, arg_def)) + return 44; + proj_release (arg_def); + /* Clean up */ - pj_free (P); + proj_destroy (P); /* Same projection, now using argc/argv style initialization */ - P = pj_create_argv (3, args); + P = proj_create_argv (0, 3, args); if (0==P) return 2; /* zero initialize everything, then set (longitude, latitude) to (12, 55) */ - a = pj_obs_null; + a = proj_obs_null; /* a.coo.lp: The coordinate part of a, interpreted as a classic LP pair */ - a.coo.lp.lam = TORAD(12); - a.coo.lp.phi = TORAD(55); + a.coo.lp.lam = PJ_TORAD(12); + a.coo.lp.phi = PJ_TORAD(55); /* Forward projection */ - b = pj_trans (P, PJ_FWD, a); + b = proj_trans_obs (P, PJ_FWD, a); /* Inverse projection */ - a = pj_trans (P, PJ_INV, b); + a = proj_trans_obs (P, PJ_INV, b); /* Null projection */ - a = pj_trans (P, PJ_IDENT, a); + a = proj_trans_obs (P, PJ_IDENT, a); /* Forward again, to get two linear items for comparison */ - a = pj_trans (P, PJ_FWD, a); + a = proj_trans_obs (P, PJ_FWD, a); - dist = pj_xy_dist (a.coo.xy, b.coo.xy); + dist = proj_xy_dist (a.coo.xy, b.coo.xy); if (dist > 2e-9) return 3; + /* Clear any previous error */ + proj_errno_set (P, 0); + /* Invalid projection */ - a = pj_trans (P, 42, a); + a = proj_trans_obs (P, 42, a); if (a.coo.lpz.lam!=HUGE_VAL) return 4; - err = pj_err_level (P, PJ_ERR_TELL); + err = proj_errno (P); if (0==err) return 5; - /* Clear error */ - pj_err_level (P, 0); + /* Clear error again */ + proj_errno_set (P, 0); /* Clean up */ - pj_free (P); + proj_destroy (P); /* Now do some 3D transformations */ - P = pj_create ("+proj=cart +ellps=GRS80"); + P = proj_create (0, "+proj=cart +ellps=GRS80"); if (0==P) return 6; /* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */ - a = b = pj_obs_null; - a.coo.lpz.lam = TORAD(12); - a.coo.lpz.phi = TORAD(55); + a = b = proj_obs_null; + a.coo.lpz.lam = PJ_TORAD(12); + a.coo.lpz.phi = PJ_TORAD(55); a.coo.lpz.z = 100; /* Forward projection: 3D-Cartesian-to-Ellipsoidal */ - b = pj_trans (P, PJ_FWD, a); + b = proj_trans_obs (P, PJ_FWD, a); /* Check roundtrip precision for 10000 iterations each way */ - dist = pj_roundtrip (P, PJ_FWD, 10000, a); - dist = pj_roundtrip (P, PJ_INV, 10000, b); + dist = proj_roundtrip (P, PJ_FWD, 10000, a); + dist = proj_roundtrip (P, PJ_INV, 10000, b); if (dist > 2e-9) return 7; /* Test at the North Pole */ - a = b = pj_obs_null; - a.coo.lpz.lam = TORAD(0); - a.coo.lpz.phi = TORAD(90); + a = b = proj_obs_null; + a.coo.lpz.lam = PJ_TORAD(0); + a.coo.lpz.phi = PJ_TORAD(90); a.coo.lpz.z = 100; /* Forward projection: Ellipsoidal-to-3D-Cartesian */ - dist = pj_roundtrip (P, PJ_FWD, 1, a); + dist = proj_roundtrip (P, PJ_FWD, 1, a); if (dist > 1e-12) return 8; /* Test at the South Pole */ - a = b = pj_obs_null; - a.coo.lpz.lam = TORAD(0); - a.coo.lpz.phi = TORAD(-90); + a = b = proj_obs_null; + a.coo.lpz.lam = PJ_TORAD(0); + a.coo.lpz.phi = PJ_TORAD(-90); a.coo.lpz.z = 100; /* Forward projection: Ellipsoidal-to-3D-Cartesian */ - dist = pj_roundtrip (P, PJ_FWD, 1, a); + dist = proj_roundtrip (P, PJ_FWD, 1, a); if (dist > 1e-12) return 9; /* Inverse projection: 3D-Cartesian-to-Ellipsoidal */ - b = pj_trans (P, PJ_INV, b); + b = proj_trans_obs (P, PJ_INV, b); /* Move p to another context */ - pj_context_renew (P); - b = pj_trans (P, PJ_FWD, b); + ctx = proj_context_create (); + if (ctx==pj_get_default_ctx()) + return 10; + proj_context_set (P, ctx); + if (ctx != P->ctx) + return 11; + b = proj_trans_obs (P, PJ_FWD, b); /* Move it back to the default context */ - pj_context_free (P); - b = pj_trans (P, PJ_INV, b); + proj_context_set (P, 0); + if (pj_get_default_ctx() != P->ctx) + return 12; + proj_context_destroy (ctx); + + /* We go on with the work - now back on the default context */ + b = proj_trans_obs (P, PJ_INV, b); + proj_destroy (P); + + + /* Testing the proj_transform nightmare */ + + /* An utm projection on the GRS80 ellipsoid */ + P = proj_create (0, "+proj=utm +zone=32 +ellps=GRS80"); + if (0==P) + return 13; + + 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); + sz = sizeof (PJ_OBS); + + /* Forward projection */ + a = proj_trans_obs (P, PJ_FWD, obs[0]); + b = proj_trans_obs (P, PJ_FWD, obs[1]); + + n = proj_transform ( + P, PJ_FWD, + &(obs[0].coo.lpz.lam), sz, 2, + &(obs[0].coo.lpz.phi), sz, 2, + &(obs[0].coo.lpz.z), sz, 2, + 0, sz, 0 + ); + if (2!=n) + return 14; + if (a.coo.lpz.lam != obs[0].coo.lpz.lam) return 15; + if (a.coo.lpz.phi != obs[0].coo.lpz.phi) return 16; + if (a.coo.lpz.z != obs[0].coo.lpz.z) return 17; + if (b.coo.lpz.lam != obs[1].coo.lpz.lam) return 18; + if (b.coo.lpz.phi != obs[1].coo.lpz.phi) return 19; + if (b.coo.lpz.z != obs[1].coo.lpz.z) return 20; + + /* now test the case of constant z */ + 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); + h = 27; + t = 33; + n = proj_transform ( + P, PJ_FWD, + &(obs[0].coo.lpz.lam), sz, 2, + &(obs[0].coo.lpz.phi), sz, 2, + &h, 0, 1, + &t, 0, 1 + ); + if (2!=n) + return 21; + if (a.coo.lpz.lam != obs[0].coo.lpz.lam) return 22; + if (a.coo.lpz.phi != obs[0].coo.lpz.phi) return 23; + if (45 != obs[0].coo.lpz.z) return 24; + if (b.coo.lpz.lam != obs[1].coo.lpz.lam) return 25; + if (b.coo.lpz.phi != obs[1].coo.lpz.phi) return 26; + if (50 != obs[1].coo.lpz.z) return 27; /* NOTE: unchanged */ + if (50==h) return 28; + + /* Clean up */ + proj_destroy (P); - pj_free (P); return 0; } #endif diff --git a/src/PJ_chamb.c b/src/PJ_chamb.c index eb3ed479c0..e9e9848bf7 100644 --- a/src/PJ_chamb.c +++ b/src/PJ_chamb.c @@ -156,7 +156,7 @@ int pj_chamb_selftest (void) { double tolerance_lp = 1e-10; double tolerance_xy = 1e-7; - char s_args[] = {"+proj=chamb +a=6400000 +lat_1=0.5 +lat_2=2"}; + char s_args[] = {"+proj=chamb +R=6400000 +lat_1=0.5 +lat_2=2"}; LP fwd_in[] = { { 2, 1}, diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c index d5d5a6834e..e291c189de 100644 --- a/src/PJ_helmert.c +++ b/src/PJ_helmert.c @@ -43,7 +43,7 @@ Last update: 2017-05-15 ***********************************************************************/ #define PJ_LIB__ -#include +#include "proj_internal.h" #include #include #include @@ -154,16 +154,16 @@ static void update_parameters(PJ *P) { Q->theta = Q->theta_0 + Q->dtheta * dt; /* debugging output */ - if (pj_err_level(P, PJ_ERR_TELL) >= PJ_LOG_TRACE) { - pj_log_trace(P, "Transformation parameters for observation epoch %g:", Q->t_obs); - pj_log_trace(P, "x: %g", Q->xyz.x); - pj_log_trace(P, "y: %g", Q->xyz.y); - pj_log_trace(P, "z: %g", Q->xyz.z); - pj_log_trace(P, "s: %g", Q->scale*1e-6); - pj_log_trace(P, "rx: %g", Q->opk.o); - pj_log_trace(P, "ry: %g", Q->opk.p); - pj_log_trace(P, "rz: %g", Q->opk.k); - pj_log_trace(P, "theta: %g", Q->theta); + if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_TRACE) { + proj_log_trace(P, "Transformation parameters for observation epoch %g:", Q->t_obs); + proj_log_trace(P, "x: %g", Q->xyz.x); + proj_log_trace(P, "y: %g", Q->xyz.y); + proj_log_trace(P, "z: %g", Q->xyz.z); + proj_log_trace(P, "s: %g", Q->scale*1e-6); + proj_log_trace(P, "rx: %g", Q->opk.o); + proj_log_trace(P, "ry: %g", Q->opk.p); + proj_log_trace(P, "rz: %g", Q->opk.k); + proj_log_trace(P, "theta: %g", Q->theta); } return; } @@ -312,11 +312,11 @@ static void build_rot_matrix(PJ *P) { } /* some debugging output */ - if (pj_err_level(P, PJ_ERR_TELL) >= PJ_LOG_TRACE) { - pj_log_trace(P, "Rotation Matrix:"); - pj_log_trace(P, " | % 6.6g % 6.6g % 6.6g |", R00, R01, R02); - pj_log_trace(P, " | % 6.6g % 6.6g % 6.6g |", R10, R11, R12); - pj_log_trace(P, " | % 6.6g % 6.6g % 6.6g |", R20, R21, R22); + if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_TRACE) { + proj_log_trace(P, "Rotation Matrix:"); + proj_log_trace(P, " | % 6.6g % 6.6g % 6.6g |", R00, R01, R02); + proj_log_trace(P, " | % 6.6g % 6.6g % 6.6g |", R10, R11, R12); + proj_log_trace(P, " | % 6.6g % 6.6g % 6.6g |", R20, R21, R22); } return; @@ -579,16 +579,16 @@ PJ *PROJECTION(helmert) { Q->theta = Q->theta_0; /* Let's help with debugging */ - if (pj_err_level(P, PJ_ERR_TELL) >= PJ_LOG_DEBUG) { - pj_log_debug(P, "Helmert parameters:"); - pj_log_debug(P, "x= % 3.5f y= % 3.5f z= % 3.5f", Q->xyz.x, Q->xyz.y, Q->xyz.z); - pj_log_debug(P, "rx= % 3.5f ry= % 3.5f rz= % 3.5f", + if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) { + proj_log_debug(P, "Helmert parameters:"); + proj_log_debug(P, "x= % 3.5f y= % 3.5f z= % 3.5f", Q->xyz.x, Q->xyz.y, Q->xyz.z); + proj_log_debug(P, "rx= % 3.5f ry= % 3.5f rz= % 3.5f", Q->opk.o / ARCSEC_TO_RAD, Q->opk.p / ARCSEC_TO_RAD, Q->opk.k / ARCSEC_TO_RAD); - pj_log_debug(P, "s=% 3.5f approximate=% d transpose=% d", + proj_log_debug(P, "s=% 3.5f approximate=% d transpose=% d", Q->scale, Q->approximate, Q->transpose); - pj_log_debug(P, "dx= % 3.5f dy= % 3.5f dz= % 3.5f", Q->dxyz.x, Q->dxyz.y, Q->dxyz.z); - pj_log_debug(P, "drx=% 3.5f dry=% 3.5f drz=% 3.5f", Q->dopk.o, Q->dopk.p, Q->dopk.k); - pj_log_debug(P, "ds=% 3.5f epoch=% 5.5f tobs=% 5.5f", Q->dscale, Q->epoch, Q->t_obs); + proj_log_debug(P, "dx= % 3.5f dy= % 3.5f dz= % 3.5f", Q->dxyz.x, Q->dxyz.y, Q->dxyz.z); + proj_log_debug(P, "drx=% 3.5f dry=% 3.5f drz=% 3.5f", Q->dopk.o, Q->dopk.p, Q->dopk.k); + proj_log_debug(P, "ds=% 3.5f epoch=% 5.5f tobs=% 5.5f", Q->dscale, Q->epoch, Q->t_obs); } if ((Q->opk.o==0) && (Q->opk.p==0) && (Q->opk.k==0) && (Q->scale==0)) { @@ -619,23 +619,26 @@ static int test (char *args, PJ_TRIPLET in, PJ_TRIPLET expect, double tol) { return 5; out.xyz = pj_fwd3d (in.lpz, P); - if (pj_xyz_dist (out.xyz, expect.xyz) > tol) { - pj_log_error(P, "Tolerance of forward calculation not met"); - pj_log_error(P, " Expect: %10.10f, %10.10f, %10.10f", expect.xyz.x, expect.xyz.y, expect.xyz.z); - pj_log_error(P, " Out: %10.10f, %10.10f, %10.10f", out.xyz.x, out.xyz.y, out.xyz.z); - pj_log_level(NULL, 0); + if (proj_xyz_dist (out.xyz, expect.xyz) > tol) { + proj_log_error(P, "Tolerance of forward calculation not met"); + proj_log_error(P, " Expect: %10.10f, %10.10f, %10.10f", expect.xyz.x, expect.xyz.y, expect.xyz.z); + proj_log_error(P, " Out: %10.10f, %10.10f, %10.10f", out.xyz.x, out.xyz.y, out.xyz.z); + proj_log_level(NULL, 0); + proj_destroy (P); return 1; } out.lpz = pj_inv3d (out.xyz, P); - if (pj_xyz_dist (out.xyz, in.xyz) > tol) { - pj_log_error(P, "Tolerance of inverse calculation not met"); - pj_log_error(P, " In: %10.10f, %10.10f, %10.10f", in.xyz.x, in.xyz.y, in.xyz.z); - pj_log_error(P, " Out: %10.10f, %10.10f, %10.10f", out.xyz.x, out.xyz.y, out.xyz.z); - pj_log_level(NULL, 0); + if (proj_xyz_dist (out.xyz, in.xyz) > tol) { + proj_log_error(P, "Tolerance of inverse calculation not met"); + proj_log_error(P, " In: %10.10f, %10.10f, %10.10f", in.xyz.x, in.xyz.y, in.xyz.z); + proj_log_error(P, " Out: %10.10f, %10.10f, %10.10f", out.xyz.x, out.xyz.y, out.xyz.z); + proj_log_level(NULL, 0); + proj_destroy (P); return 2; } + proj_destroy (P); return 0; } @@ -675,7 +678,7 @@ int pj_helmert_selftest (void) { PJ_TRIPLET in3 = {{3370658.378, 711877.314, 5349787.086}}; /* ITRF2000@2017.0 */ PJ_TRIPLET expect3 = {{3370658.18890, 711877.42370, 5349787.12430}}; /* ITRF93@2017.0 */ char args3[] = { - " +proj=helmert" + " +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" " +dx=-0.0029 +dy=-0.0002 +dz=-0.0006 +ds=0.00001" @@ -690,8 +693,9 @@ int pj_helmert_selftest (void) { PJ_OBS in4 = {{{3370658.378, 711877.314, 5349787.086, 2017.0}}, {{ 0, 0, 0}}, 0, 0}; PJ_OBS out; - PJ *helmert = pj_create( - " +proj=helmert" + PJ *helmert = proj_create( + 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" " +dx=-0.0029 +dy=-0.0002 +dz=-0.0006 +ds=0.00001" @@ -709,7 +713,7 @@ int pj_helmert_selftest (void) { to computed coordinates, hence the test tolerance is set accordingly. */ PJ_TRIPLET in5 = {{2546506.957, 542256.609, 0}}; PJ_TRIPLET expect5 = {{766563.675, 165282.277, 0}}; - char args5[] = " +proj=helmert +x=-9597.3572 +y=.6112 +s=0.304794780637 +theta=-1.244048"; + char args5[] = " +proj=helmert +ellps=GRS80 +x=-9597.3572 +y=.6112 +s=0.304794780637 +theta=-1.244048"; /* Run tests 1-3 & 5 */ ret = test (args1, in1, expect1, 1e-4); if (ret) return ret; @@ -718,24 +722,24 @@ int pj_helmert_selftest (void) { ret = test (args5, in5, expect5, 0.001); if (ret) return ret + 40; /* Run test 4 */ - out = pj_trans(helmert, PJ_FWD, in4); - if (pj_xyz_dist (out.coo.xyz, expect4a) > 1e-4) { - pj_log_error(helmert, "Tolerance of test 4a not met!"); - pj_log_error(helmert, " In: %10.10f, %10.10f, %10.10f", in4.coo.xyz.x, in4.coo.xyz.y, in4.coo.xyz.z); - pj_log_error(helmert, " Out: %10.10f, %10.10f, %10.10f", out.coo.xyz.x, out.coo.xyz.y, out.coo.xyz.z); + out = proj_trans_obs (helmert, PJ_FWD, in4); + if (proj_xyz_dist (out.coo.xyz, expect4a) > 1e-4) { + proj_log_error(helmert, "Tolerance of test 4a not met!"); + proj_log_error(helmert, " In: %10.10f, %10.10f, %10.10f", in4.coo.xyz.x, in4.coo.xyz.y, in4.coo.xyz.z); + proj_log_error(helmert, " Out: %10.10f, %10.10f, %10.10f", out.coo.xyz.x, out.coo.xyz.y, out.coo.xyz.z); return 31; } in4.coo.xyzt.t = 2018.0; - out = pj_trans(helmert, PJ_FWD, in4); - if (pj_xyz_dist (out.coo.xyz, expect4b) > 1e-4) { - pj_log_error(helmert, "Tolerance of test 4b not met!"); - pj_log_error(helmert, " In: %10.10f, %10.10f, %10.10f", in4.coo.xyz.x, in4.coo.xyz.y, in4.coo.xyz.z); - pj_log_error(helmert, " Out: %10.10f, %10.10f, %10.10f", out.coo.xyz.x, out.coo.xyz.y, out.coo.xyz.z); + out = proj_trans_obs (helmert, PJ_FWD, in4); + if (proj_xyz_dist (out.coo.xyz, expect4b) > 1e-4) { + proj_log_error(helmert, "Tolerance of test 4b not met!"); + proj_log_error(helmert, " In: %10.10f, %10.10f, %10.10f", in4.coo.xyz.x, in4.coo.xyz.y, in4.coo.xyz.z); + proj_log_error(helmert, " Out: %10.10f, %10.10f, %10.10f", out.coo.xyz.x, out.coo.xyz.y, out.coo.xyz.z); return 32; } - pj_free(helmert); + proj_destroy(helmert); return 0; } diff --git a/src/PJ_hgridshift.c b/src/PJ_hgridshift.c index 1b9d272701..41c2b62949 100644 --- a/src/PJ_hgridshift.c +++ b/src/PJ_hgridshift.c @@ -1,5 +1,5 @@ #define PJ_LIB__ -#include +#include "proj_internal.h" #include PROJ_HEAD(hgridshift, "Horizontal grid shift"); @@ -70,7 +70,7 @@ static PJ_OBS reverse_obs(PJ_OBS obs, PJ *P) { PJ *PROJECTION(hgridshift) { if (!pj_param(P->ctx, P->params, "tgrids").i) { - pj_log_error(P, "hgridshift: +grids parameter missing."); + proj_log_error(P, "hgridshift: +grids parameter missing."); return freeup_msg(P, -1); } @@ -80,7 +80,7 @@ PJ *PROJECTION(hgridshift) { /* Was gridlist compiled properly? */ if ( pj_ctx_get_errno(P->ctx) ) { - pj_log_error(P, "hgridshift: could not find required grid(s)."); + proj_log_error(P, "hgridshift: could not find required grid(s)."); return freeup_msg(P, -38); } @@ -108,38 +108,36 @@ int pj_hgridshift_selftest (void) { double dist; /* fail on purpose: +grids parameter is mandatory*/ - P = pj_create("+proj=hgridshift"); + P = proj_create(0, "+proj=hgridshift"); if (0!=P) return 99; /* fail on purpose: open non-existing grid */ - P = pj_create("+proj=hgridshift +grids=nonexistinggrid.gsb"); + P = proj_create(0, "+proj=hgridshift +grids=nonexistinggrid.gsb"); if (0!=P) return 999; /* Failure most likely means the grid is missing */ - P = pj_create ("+proj=hgridshift +grids=nzgd2kgrid0005.gsb +ellps=GRS80"); + P = proj_create (0, "+proj=hgridshift +grids=nzgd2kgrid0005.gsb +ellps=GRS80"); if (0==P) return 10; - a = pj_obs_null; - a.coo.lpz.lam = TORAD(173); - a.coo.lpz.phi = TORAD(-45); + a = proj_obs_null; + a.coo.lpz.lam = PJ_TORAD(173); + a.coo.lpz.phi = PJ_TORAD(-45); - dist = pj_roundtrip (P, PJ_FWD, 1, a); + dist = proj_roundtrip (P, PJ_FWD, 1, a); if (dist > 0.00000001) return 1; - expect.coo.lpz.lam = TORAD(172.999892181021551); - expect.coo.lpz.phi = TORAD(-45.001620431954613); - b = pj_trans(P, PJ_FWD, a); - if (pj_xy_dist(expect.coo.xy, b.coo.xy) > 1e-4) + expect.coo.lpz.lam = PJ_TORAD(172.999892181021551); + expect.coo.lpz.phi = PJ_TORAD(-45.001620431954613); + b = proj_trans_obs(P, PJ_FWD, a); + if (proj_xy_dist(expect.coo.xy, b.coo.xy) > 1e-4) return 2; - - pj_free(P); - + proj_destroy(P); return 0; } #endif diff --git a/src/PJ_horner.c b/src/PJ_horner.c index 1b055a4e79..56980f724c 100644 --- a/src/PJ_horner.c +++ b/src/PJ_horner.c @@ -1,5 +1,5 @@ #define PJ_LIB__ -#include +#include "proj_internal.h" #include #include #include @@ -94,7 +94,7 @@ PROJ_HEAD(horner, "Horner polynomial evaluation"); struct horner; typedef struct horner HORNER; -static UV horner (const HORNER *transformation, enum pj_direction, UV position); +static UV horner (const HORNER *transformation, enum proj_direction, UV position); static HORNER *horner_alloc (size_t order, int complex_polynomia); static void horner_free (HORNER *h); @@ -176,7 +176,7 @@ static HORNER *horner_alloc (size_t order, int complex_polynomia) { /**********************************************************************/ -static UV horner (const HORNER *transformation, enum pj_direction direction, UV position) { +static UV horner (const HORNER *transformation, enum proj_direction direction, UV position) { /*********************************************************************** A reimplementation of the classic Engsager/Poder 2D Horner polynomial @@ -301,7 +301,7 @@ static PJ_OBS horner_reverse_obs (PJ_OBS point, PJ *P) { /**********************************************************************/ -static UV complex_horner (const HORNER *transformation, enum pj_direction direction, UV position) { +static UV complex_horner (const HORNER *transformation, enum proj_direction direction, UV position) { /*********************************************************************** A reimplementation of a classic Engsager/Poder Horner complex @@ -401,7 +401,7 @@ static int parse_coefs (PJ *P, double *coefs, char *param, int ncoefs) { buf = pj_calloc (strlen (param) + 2, sizeof(char)); if (0==buf) { - pj_log_error (P, "Horner: Out of core"); + proj_log_error (P, "Horner: Out of core"); return 0; } @@ -415,7 +415,7 @@ static int parse_coefs (PJ *P, double *coefs, char *param, int ncoefs) { for (i = 0; i < ncoefs; i++) { if (i > 0) { if ( next == 0 || ','!=*next) { - pj_log_error (P, "Horner: Malformed polynomium set %s. need %d coefs", param, ncoefs); + proj_log_error (P, "Horner: Malformed polynomium set %s. need %d coefs", param, ncoefs); return 0; } init = ++next; @@ -443,7 +443,7 @@ PJ *PROJECTION(horner) { if (pj_param (P->ctx, P->params, "tdeg").i) /* degree specified? */ degree = pj_param(P->ctx, P->params, "ideg").i; else { - pj_log_debug (P, "Horner: Must specify polynomial degree, (+deg=n)"); + proj_log_debug (P, "Horner: Must specify polynomial degree, (+deg=n)"); return horner_freeup (P); } @@ -525,48 +525,49 @@ int pj_horner_selftest (void) { double dist; /* Real polynonia relating the technical coordinate system TC32 to "System 45 Bornholm" */ - P = pj_create (tc32_utm32); + P = proj_create (0, tc32_utm32); if (0==P) return 10; - a = b = pj_obs_null; + a = b = proj_obs_null; a.coo.uv.v = 6125305.4245; a.coo.uv.u = 878354.8539; /* Check roundtrip precision for 1 iteration each way, starting in forward direction */ - dist = pj_roundtrip (P, PJ_FWD, 1, a); + dist = proj_roundtrip (P, PJ_FWD, 1, a); if (dist > 0.01) return 1; /* The complex polynomial transformation between the "System Storebaelt" and utm32/ed50 */ - P = pj_create (sb_utm32); + P = proj_create (0, sb_utm32); if (0==P) return 11; /* Test value: utm32_ed50(620000, 6130000) = sb_ed50(495136.8544, 6130821.2945) */ - a = b = c = pj_obs_null; + a = b = c = proj_obs_null; a.coo.uv.v = 6130821.2945; a.coo.uv.u = 495136.8544; c.coo.uv.v = 6130000.0000; c.coo.uv.u = 620000.0000; /* Forward projection */ - b = pj_trans (P, PJ_FWD, a); - dist = pj_xy_dist (b.coo.xy, c.coo.xy); + b = proj_trans_obs (P, PJ_FWD, a); + dist = proj_xy_dist (b.coo.xy, c.coo.xy); if (dist > 0.001) return 2; /* Inverse projection */ - b = pj_trans (P, PJ_INV, c); - dist = pj_xy_dist (b.coo.xy, a.coo.xy); + b = proj_trans_obs (P, PJ_INV, c); + dist = proj_xy_dist (b.coo.xy, a.coo.xy); if (dist > 0.001) return 3; /* Check roundtrip precision for 1 iteration each way */ - dist = pj_roundtrip (P, PJ_FWD, 1, a); + dist = proj_roundtrip (P, PJ_FWD, 1, a); if (dist > 0.01) return 4; + proj_destroy(P); return 0; } #endif diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c index 387c0ca862..6e42c0cba0 100644 --- a/src/PJ_pipeline.c +++ b/src/PJ_pipeline.c @@ -96,7 +96,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-20 ********************************************************************************/ #define PJ_LIB__ -#include +#include "proj_internal.h" #include #include @@ -183,7 +183,7 @@ static PJ_OBS pipeline_forward_obs (PJ_OBS point, PJ *P) { last_step = P->opaque->steps + 1; for (i = first_step; i != last_step; i++) { - pj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %20.20s", + proj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %20.20s", i-first_step, point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z, P->opaque->pipeline[i]->descr ); @@ -191,13 +191,13 @@ static PJ_OBS pipeline_forward_obs (PJ_OBS point, PJ *P) { if (P->opaque->omit_forward[i]) continue; if (P->opaque->reverse_step[i]) - point = pj_trans (P->opaque->pipeline[i], PJ_INV, point); + point = proj_trans_obs (P->opaque->pipeline[i], PJ_INV, point); else - point = pj_trans (P->opaque->pipeline[i], PJ_FWD, point); + point = proj_trans_obs (P->opaque->pipeline[i], PJ_FWD, point); if (P->opaque->depth < PIPELINE_STACK_SIZE) P->opaque->stack[P->opaque->depth++] = point; } - pj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z); + proj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z); P->opaque->depth = 0; /* Clear the stack */ return point; @@ -210,20 +210,20 @@ static PJ_OBS pipeline_reverse_obs (PJ_OBS point, PJ *P) { first_step = P->opaque->steps; last_step = 0; for (i = first_step; i != last_step; i--) { - pj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %.4f %.4f", + proj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %.4f %.4f", i - 1, point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z, P->opaque->pipeline[i]->a, P->opaque->pipeline[i]->rf ); if (P->opaque->omit_inverse[i]) continue; if (P->opaque->reverse_step[i]) - point = pj_trans (P->opaque->pipeline[i], PJ_FWD, point); + point = proj_trans_obs (P->opaque->pipeline[i], PJ_FWD, point); else - point = pj_trans (P->opaque->pipeline[i], PJ_INV, point); + point = proj_trans_obs (P->opaque->pipeline[i], PJ_INV, point); if (P->opaque->depth < PIPELINE_STACK_SIZE) P->opaque->stack[P->opaque->depth++] = point; } - pj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z); + proj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z); P->opaque->depth = 0; /* Clear the stack */ return point; @@ -232,7 +232,7 @@ static PJ_OBS pipeline_reverse_obs (PJ_OBS point, PJ *P) { /* Delegate the work to pipeline_forward_obs() */ static XYZ pipeline_forward_3d (LPZ lpz, PJ *P) { - PJ_OBS point = pj_obs_null; + PJ_OBS point = proj_obs_null; point.coo.lpz = lpz; point = pipeline_forward_obs (point, P); return point.coo.xyz; @@ -240,21 +240,21 @@ static XYZ pipeline_forward_3d (LPZ lpz, PJ *P) { /* Delegate the work to pipeline_reverse_obs() */ static LPZ pipeline_reverse_3d (XYZ xyz, PJ *P) { - PJ_OBS point = pj_obs_null; + PJ_OBS point = proj_obs_null; point.coo.xyz = xyz; point = pipeline_reverse_obs (point, P); return point.coo.lpz; } static XY pipeline_forward (LP lp, PJ *P) { - PJ_OBS point = pj_obs_null; + PJ_OBS point = proj_obs_null; point.coo.lp = lp; point = pipeline_forward_obs (point, P); return point.coo.xy; } static LP pipeline_reverse (XY xy, PJ *P) { - PJ_OBS point = pj_obs_null; + PJ_OBS point = proj_obs_null; point.coo.xy = xy; point = pipeline_reverse_obs (point, P); return point.coo.lp; @@ -275,7 +275,8 @@ static void *pipeline_freeup (PJ *P, int errlev) { /* Destructor */ if (0==P) return 0; - pj_err_level (P, errlev); + if (errlev) + proj_errno_set (P, errlev); if (0==P->opaque) return pj_dealloc (P); @@ -382,7 +383,7 @@ PJ *PROJECTION(pipeline) { for (i = 0; i < argc; i++) { if (0==strcmp ("step", argv[i])) { if (-1==i_pipeline) { - pj_log_error (P, "Pipeline: +step before +proj=pipeline"); + proj_log_error (P, "Pipeline: +step before +proj=pipeline"); return pipeline_freeup (P, -50); } if (0==nsteps) @@ -393,7 +394,7 @@ PJ *PROJECTION(pipeline) { if (0==strcmp ("proj=pipeline", argv[i])) { if (-1 != i_pipeline) { - pj_log_error (P, "Pipeline: Nesting invalid"); + proj_log_error (P, "Pipeline: Nesting invalid"); return pipeline_freeup (P, -50); /* ERROR: nested pipelines */ } i_pipeline = i; @@ -419,7 +420,7 @@ PJ *PROJECTION(pipeline) { PJ *next_step = 0; /* Build a set of setup args for the current step */ - pj_log_trace (P, "Pipeline: Building arg list for step no. %d", i); + proj_log_trace (P, "Pipeline: Building arg list for step no. %d", i); /* First add the step specific args */ for (j = i_current_step + 1; 0 != strcmp ("step", argv[j]); j++) @@ -445,21 +446,21 @@ PJ *PROJECTION(pipeline) { P->opaque->reverse_step[i+1] = 1; } - pj_log_trace (P, "Pipeline: init - %s, %d", current_argv[0], current_argc); + proj_log_trace (P, "Pipeline: init - %s, %d", current_argv[0], current_argc); for (j = 1; j < current_argc; j++) - pj_log_trace (P, " %s", current_argv[j]); + proj_log_trace (P, " %s", current_argv[j]); next_step = pj_init (current_argc, current_argv); - pj_log_trace (P, "Pipeline: Step %d at %p", i, next_step); + proj_log_trace (P, "Pipeline: Step %d at %p", i, next_step); if (0==next_step) { - pj_log_error (P, "Pipeline: Bad step definition: %s", current_argv[0]); + proj_log_error (P, "Pipeline: Bad step definition: %s", current_argv[0]); return pipeline_freeup (P, -50); /* ERROR: bad pipeline def */ } P->opaque->pipeline[i+1] = next_step; - pj_log_trace (P, "Pipeline: step done"); + proj_log_trace (P, "Pipeline: step done"); } - pj_log_trace (P, "Pipeline: %d steps built. Determining i/o characteristics", nsteps); + proj_log_trace (P, "Pipeline: %d steps built. Determining i/o characteristics", nsteps); /* Determine forward input (= reverse output) data type */ @@ -468,7 +469,7 @@ PJ *PROJECTION(pipeline) { if (0==P->opaque->omit_forward[i+1]) break; if (i==nsteps) { - pj_log_error (P, "Pipeline: No forward steps"); + proj_log_error (P, "Pipeline: No forward steps"); return pipeline_freeup (P, -50); } @@ -491,7 +492,7 @@ PJ *PROJECTION(pipeline) { if (0==P->opaque->omit_inverse[i+1]) break; if (i==-1) { - pj_log (pj_get_ctx (P), PJ_LOG_ERROR, "Pipeline: No reverse steps"); + proj_log_error (P, "Pipeline: No reverse steps"); return pipeline_freeup (P, -50); } @@ -506,7 +507,7 @@ PJ *PROJECTION(pipeline) { else P->right = PJ_IO_UNITS_METERS; } - pj_log_trace (P, "Pipeline: Units - left: [%s], right: [%s]\n", + proj_log_trace (P, "Pipeline: Units - left: [%s], right: [%s]\n", P->left ==PJ_IO_UNITS_RADIANS? "angular": "linear", P->right==PJ_IO_UNITS_RADIANS? "angular": "linear"); @@ -525,94 +526,96 @@ int pj_pipeline_selftest (void) { double dist; /* forward-reverse geo->utm->geo */ - P = pj_create ("+proj=pipeline +zone=32 +step +proj=utm +ellps=GRS80 +step +proj=utm +ellps=GRS80 +inv"); + P = proj_create (0, "+proj=pipeline +zone=32 +step +proj=utm +ellps=GRS80 +step +proj=utm +ellps=GRS80 +inv"); if (0==P) return 1000; /* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 0) */ - a = b = pj_obs_null; - a.coo.lpz.lam = TORAD(12); - a.coo.lpz.phi = TORAD(55); + a = b = proj_obs_null; + a.coo.lpz.lam = PJ_TORAD(12); + a.coo.lpz.phi = PJ_TORAD(55); a.coo.lpz.z = 0; /* Forward projection */ - b = pj_trans (P, PJ_FWD, a); - if (pj_lp_dist (P, a.coo.lp, b.coo.lp) > 1e-4) + b = proj_trans_obs (P, PJ_FWD, a); + if (proj_lp_dist (P, a.coo.lp, b.coo.lp) > 1e-4) return 1001; /* Inverse projection (still same result: pipeline is symmetrical) */ - a = pj_trans (P, PJ_INV, b); - if (pj_lp_dist (P, a.coo.lp, b.coo.lp) > 1e-4) + a = proj_trans_obs (P, PJ_INV, b); + if (proj_lp_dist (P, a.coo.lp, b.coo.lp) > 1e-4) return 1002; - pj_free (P); + proj_destroy (P); /* And now the back-to-back situation utm->geo->utm */ - P = pj_create ("+proj=pipeline +zone=32 +step +proj=utm +ellps=GRS80 +inv +step +proj=utm +ellps=GRS80"); + P = proj_create (0, "+proj=pipeline +zone=32 +step +proj=utm +ellps=GRS80 +inv +step +proj=utm +ellps=GRS80"); if (0==P) return 2000; /* zero initialize everything, then set (easting, northing) to utm(12, 55) */ - a = b = pj_obs_null; + a = b = proj_obs_null; a.coo.xy = cph_utm32; /* Forward projection */ - b = pj_trans (P, PJ_FWD, a); - if (pj_xy_dist (a.coo.xy, b.coo.xy) > 1e-4) + b = proj_trans_obs (P, PJ_FWD, a); + if (proj_xy_dist (a.coo.xy, b.coo.xy) > 1e-4) return 2001; /* Inverse projection */ - a = pj_trans (P, PJ_INV, b); - if (pj_xy_dist (a.coo.xy, b.coo.xy) > 1e-4) + a = proj_trans_obs (P, PJ_INV, b); + if (proj_xy_dist (a.coo.xy, b.coo.xy) > 1e-4) return 2001; - if (pj_xyz_dist (a.coo.xyz, b.coo.xyz) > 1e-4) + if (proj_xyz_dist (a.coo.xyz, b.coo.xyz) > 1e-4) return 2002; - pj_free (P); + proj_destroy (P); /* Finally testing a corner case: A rather pointless one-step pipeline geo->utm */ - P = pj_create ("+proj=pipeline +zone=32 +step +proj=utm +ellps=GRS80 "); + P = proj_create (0, "+proj=pipeline +zone=32 +step +proj=utm +ellps=GRS80 "); if (0==P) return 3000; - a = b = pj_obs_null; - a.coo.lpz.lam = TORAD(12); - a.coo.lpz.phi = TORAD(55); + a = b = proj_obs_null; + a.coo.lpz.lam = PJ_TORAD(12); + a.coo.lpz.phi = PJ_TORAD(55); /* Forward projection */ - b = pj_trans (P, PJ_FWD, a); - if (pj_xy_dist (cph_utm32, b.coo.xy) > 1e-4) + b = proj_trans_obs (P, PJ_FWD, a); + if (proj_xy_dist (cph_utm32, b.coo.xy) > 1e-4) return 3001; /* Inverse projection */ - b = pj_trans (P, PJ_INV, b); - if (pj_lp_dist (P, a.coo.lp, b.coo.lp) > 1e-4) + b = proj_trans_obs (P, PJ_INV, b); + if (proj_lp_dist (P, a.coo.lp, b.coo.lp) > 1e-4) return 3002; /* Since we use pj_lp_dist to determine success above, we should also test that it works */ /* Geodesic distance between two points with angular 2D coordinates */ - a.coo.lp.lam = TORAD(12); - a.coo.lp.phi = TORAD(60); - b.coo.lp.lam = TORAD(12); - b.coo.lp.phi = TORAD(61); - dist = pj_lp_dist (P, a.coo.lp, b.coo.lp); + a.coo.lp.lam = PJ_TORAD(12); + a.coo.lp.phi = PJ_TORAD(60); + b.coo.lp.lam = PJ_TORAD(12); + b.coo.lp.phi = PJ_TORAD(61); + dist = proj_lp_dist (P, a.coo.lp, b.coo.lp); if (fabs (111420.727870234 - dist) > 1e-4) return 4001; - a.coo.lp.lam = TORAD(12); - a.coo.lp.phi = TORAD(0.); - b.coo.lp.lam = TORAD(12); - b.coo.lp.phi = TORAD(1.); - dist = pj_lp_dist (P, a.coo.lp, b.coo.lp); + a.coo.lp.lam = PJ_TORAD(12); + a.coo.lp.phi = PJ_TORAD(0.); + b.coo.lp.lam = PJ_TORAD(12); + b.coo.lp.phi = PJ_TORAD(1.); + dist = proj_lp_dist (P, a.coo.lp, b.coo.lp); if (fabs (110574.388554153 - dist) > 1e-4) return 4002; + proj_destroy (P); /* test a pipeline with several +init steps */ - P = pj_create( + P = proj_create( + 0, "+proj=pipeline " "+step +init=epsg:25832 +inv " "+step +init=epsg:25833 " @@ -625,13 +628,13 @@ int pj_pipeline_selftest (void) { a.coo.xy.x = 700000.0; a.coo.xy.y = 6000000.0; - b = pj_trans(P, PJ_FWD, a); - dist = pj_xy_dist(a.coo.xy, b.coo.xy); + b = proj_trans_obs(P, PJ_FWD, a); + dist = proj_xy_dist(a.coo.xy, b.coo.xy); if (dist > 1e-7) return 5001; - pj_free (P); + proj_destroy (P); return 0; } #endif diff --git a/src/PJ_unitconvert.c b/src/PJ_unitconvert.c index be9d2f705e..27672c8d07 100644 --- a/src/PJ_unitconvert.c +++ b/src/PJ_unitconvert.c @@ -65,8 +65,9 @@ Last update: 2017-05-16 #define PJ_LIB__ #include -#include +#include "proj_internal.h" #include +#include PROJ_HEAD(unitconvert, "Unit conversion"); @@ -348,7 +349,7 @@ PJ *PROJECTION(unitconvert) { if (!s) return freeup_msg(P, -8); /* unknown unit conversion id */ Q->xy_in_id = i; - pj_log_debug(P, "xy_in unit: %s", pj_units[i].name); + proj_log_debug(P, "xy_in unit: %s", pj_units[i].name); } if ((name = pj_param (P->ctx, P->params, "sxy_out").s) != NULL) { @@ -357,7 +358,7 @@ PJ *PROJECTION(unitconvert) { if (!s) return freeup_msg(P, -8); /* unknown unit conversion id */ Q->xy_out_id = i; - pj_log_debug(P, "xy_out unit: %s", pj_units[i].name); + proj_log_debug(P, "xy_out unit: %s", pj_units[i].name); } if ((name = pj_param (P->ctx, P->params, "sz_in").s) != NULL) { @@ -366,7 +367,7 @@ PJ *PROJECTION(unitconvert) { if (!s) return freeup_msg(P, -8); /* unknown unit conversion id */ Q->z_in_id = i; - pj_log_debug(P, "z_in unit: %s", pj_units[i].name); + proj_log_debug(P, "z_in unit: %s", pj_units[i].name); } if ((name = pj_param (P->ctx, P->params, "sz_out").s) != NULL) { @@ -375,7 +376,7 @@ PJ *PROJECTION(unitconvert) { if (!s) return freeup_msg(P, -8); /* unknown unit conversion id */ Q->z_out_id = i; - pj_log_debug(P, "z_out unit: %s", pj_units[i].name); + proj_log_debug(P, "z_out unit: %s", pj_units[i].name); } @@ -385,7 +386,7 @@ PJ *PROJECTION(unitconvert) { if (!s) return freeup_msg(P, -8); /* unknown unit conversion id */ Q->t_in_id = i; - pj_log_debug(P, "t_in unit: %s", time_units[i].name); + proj_log_debug(P, "t_in unit: %s", time_units[i].name); } s = 0; @@ -395,7 +396,7 @@ PJ *PROJECTION(unitconvert) { return freeup_msg(P, -8); /* unknown unit conversion id */ } Q->t_out_id = i; - pj_log_debug(P, "t_out unit: %s", time_units[i].name); + proj_log_debug(P, "t_out unit: %s", time_units[i].name); } return P; @@ -409,7 +410,7 @@ int pj_unitconvert_selftest (void) {return 0;} static int test_time(char* args, double tol, double t_in, double t_exp) { PJ_OBS in, out; - PJ *P = pj_create(args); + PJ *P = proj_create(0, args); int ret = 0; if (P == 0) @@ -417,46 +418,46 @@ static int test_time(char* args, double tol, double t_in, double t_exp) { in.coo.xyzt.t = t_in; - out = pj_trans(P, PJ_FWD, in); + out = proj_trans_obs(P, PJ_FWD, in); if (fabs(out.coo.xyzt.t - t_exp) > tol) { - pj_log_error(P, "out: %10.10g, expect: %10.10g", out.coo.xyzt.t, t_exp); + proj_log_error(P, "out: %10.10g, expect: %10.10g", out.coo.xyzt.t, t_exp); ret = 1; } - out = pj_trans(P, PJ_INV, out); + out = proj_trans_obs(P, PJ_INV, out); if (fabs(out.coo.xyzt.t - t_in) > tol) { - pj_log_error(P, "out: %10.10g, expect: %10.10g", out.coo.xyzt.t, t_in); + proj_log_error(P, "out: %10.10g, expect: %10.10g", out.coo.xyzt.t, t_in); ret = 2; } pj_free(P); - pj_log_level(NULL, 0); + proj_log_level(NULL, 0); return ret; } static int test_xyz(char* args, double tol, PJ_TRIPLET in, PJ_TRIPLET exp) { PJ_OBS out, obs_in; - PJ *P = pj_create(args); + PJ *P = proj_create(0, args); int ret = 0; if (P == 0) return 5; obs_in.coo.xyz = in.xyz; - out = pj_trans(P, PJ_FWD, obs_in); - if (pj_xyz_dist(out.coo.xyz, exp.xyz) > tol) { + out = proj_trans_obs(P, PJ_FWD, obs_in); + if (proj_xyz_dist(out.coo.xyz, exp.xyz) > tol) { printf("exp: %10.10g, %10.10g, %10.10g\n", exp.xyz.x, exp.xyz.y, exp.xyz.z); printf("out: %10.10g, %10.10g, %10.10g\n", out.coo.xyz.x, out.coo.xyz.y, out.coo.xyz.z); ret = 1; } - out = pj_trans(P, PJ_INV, out); - if (pj_xyz_dist(out.coo.xyz, in.xyz) > tol) { + out = proj_trans_obs(P, PJ_INV, out); + if (proj_xyz_dist(out.coo.xyz, in.xyz) > tol) { printf("exp: %g, %g, %g\n", in.xyz.x, in.xyz.y, in.xyz.z); printf("out: %g, %g, %g\n", out.coo.xyz.x, out.coo.xyz.y, out.coo.xyz.z); ret += 2; } - pj_free(P); - pj_log_level(NULL, 0); + proj_destroy(P); + proj_log_level(NULL, 0); return ret; } diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c index 72d975172d..e790cd4d39 100644 --- a/src/PJ_vgridshift.c +++ b/src/PJ_vgridshift.c @@ -1,5 +1,5 @@ #define PJ_LIB__ -#include +#include "proj_internal.h" #include PROJ_HEAD(vgridshift, "Vertical grid shift"); @@ -73,7 +73,7 @@ static PJ_OBS reverse_obs(PJ_OBS obs, PJ *P) { PJ *PROJECTION(vgridshift) { if (!pj_param(P->ctx, P->params, "tgrids").i) { - pj_log_error(P, "vgridshift: +grids parameter missing."); + proj_log_error(P, "vgridshift: +grids parameter missing."); return freeup_msg(P, -1); } @@ -83,7 +83,7 @@ PJ *PROJECTION(vgridshift) { /* Was gridlist compiled properly? */ if ( pj_ctx_get_errno(P->ctx) ) { - pj_log_error(P, "vgridshift: could not find required grid(s)."); + proj_log_error(P, "vgridshift: could not find required grid(s)."); pj_dalloc(P->gridlist); P->gridlist = NULL; return freeup_msg(P, -38); @@ -113,32 +113,32 @@ int pj_vgridshift_selftest (void) { double dist; /* fail on purpose: +grids parameter it mandatory*/ - P = pj_create("+proj=vgridshift"); + P = proj_create(0, "+proj=vgridshift"); if (0!=P) return 99; /* fail on purpose: open non-existing grid */ - P = pj_create("+proj=vgridshift +grids=nonexistinggrid.gtx"); + P = proj_create(0, "+proj=vgridshift +grids=nonexistinggrid.gtx"); if (0!=P) return 999; /* Failure most likely means the grid is missing */ - P = pj_create ("+proj=vgridshift +grids=egm96_15.gtx +ellps=GRS80"); + P = proj_create (0, "+proj=vgridshift +grids=egm96_15.gtx +ellps=GRS80"); if (0==P) return 10; - a = pj_obs_null; - a.coo.lpz.lam = TORAD(12.5); - a.coo.lpz.phi = TORAD(55.5); + a = proj_obs_null; + a.coo.lpz.lam = PJ_TORAD(12.5); + a.coo.lpz.phi = PJ_TORAD(55.5); - dist = pj_roundtrip (P, PJ_FWD, 1, a); + dist = proj_roundtrip (P, PJ_FWD, 1, a); if (dist > 0.00000001) return 1; expect = a; expect.coo.lpz.z = -36.021305084228515625; - b = pj_trans(P, PJ_FWD, a); - if (pj_xyz_dist(expect.coo.xyz, b.coo.xyz) > 1e-10) + b = proj_trans_obs(P, PJ_FWD, a); + if (proj_xyz_dist(expect.coo.xyz, b.coo.xyz) > 1e-10) return 2; diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 9ed1623c63..39dab3a4fc 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -201,6 +201,8 @@ SET(SRC_LIBPROJ_CORE pj_msfn.c pj_mutex.c pj_obs_api.c + pj_internal.c + proj_internal.h pj_open_lib.c pj_param.c pj_phi2.c diff --git a/src/makefile.vc b/src/makefile.vc index 8584d41980..b37d3b7b29 100644 --- a/src/makefile.vc +++ b/src/makefile.vc @@ -56,7 +56,8 @@ support = \ pj_utils.obj pj_gridlist.obj pj_gridinfo.obj \ proj_mdist.obj pj_mutex.obj pj_initcache.obj \ pj_ctx.obj pj_fileapi.obj pj_log.obj pj_apply_vgridshift.obj \ - pj_strtod.obj pj_run_selftests.obj pj_generic_selftest.obj + pj_strtod.obj pj_run_selftests.obj pj_generic_selftest.obj \ + pj_internal.obj pipeline = \ pj_obs_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \ diff --git a/src/pj_ctx.c b/src/pj_ctx.c index 9ae3d4e0e3..89b2816f67 100644 --- a/src/pj_ctx.c +++ b/src/pj_ctx.c @@ -60,8 +60,13 @@ void pj_set_ctx( projPJ pj, projCtx ctx ) projCtx pj_get_default_ctx() { + /* If already initialized, don't bother locking */ + if( default_context_initialized ) + return &default_context; + pj_acquire_lock(); + /* Ask again, since it may have been initialized in another thread */ if( !default_context_initialized ) { default_context.last_errno = 0; diff --git a/src/pj_fwd.c b/src/pj_fwd.c index 598c6dd2b4..02f18d1c77 100644 --- a/src/pj_fwd.c +++ b/src/pj_fwd.c @@ -1,21 +1,21 @@ /* general forward projection */ #define PJ_LIB__ +#include #include -#include # define EPS 1.0e-12 XY /* forward projection entry */ pj_fwd(LP lp, PJ *P) { XY xy; XY err; double t; + int last_errno; /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */ err.x = err.y = HUGE_VAL; if (0==P->fwd) return err; - - P->ctx->last_errno = pj_errno = errno = 0; + last_errno = proj_errno_reset (P); /* Check input coordinates if angular */ if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) { @@ -39,7 +39,7 @@ pj_fwd(LP lp, PJ *P) { /* Do the transformation */ xy = (*P->fwd)(lp, P); - if ( P->ctx->last_errno ) + if ( proj_errno (P) ) return err; /* Classic proj.4 functions return plane coordinates in units of the semimajor axis */ @@ -53,5 +53,6 @@ pj_fwd(LP lp, PJ *P) { xy.y = P->fr_meter * (xy.y + P->y0); /* z is not scaled since this is handled by vto_meter outside */ + proj_errno_restore (P, last_errno); return xy; } diff --git a/src/pj_fwd3d.c b/src/pj_fwd3d.c index 0fb2f63313..be683745b7 100644 --- a/src/pj_fwd3d.c +++ b/src/pj_fwd3d.c @@ -1,4 +1,5 @@ #define PJ_LIB__ +#include #include #include # define EPS 1.0e-12 @@ -9,6 +10,7 @@ XYZ pj_fwd3d(LPZ lpz, PJ *P) { XYZ xyz; XYZ err; double t; + int last_errno; /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */ err.x = err.y = err.z = HUGE_VAL; @@ -16,8 +18,7 @@ XYZ pj_fwd3d(LPZ lpz, PJ *P) { if (0==P->fwd3d) return err; - - P->ctx->last_errno = pj_errno = errno = 0; + last_errno = proj_errno_reset(P); /* Check input coordinates if angular */ if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) { @@ -55,5 +56,6 @@ XYZ pj_fwd3d(LPZ lpz, PJ *P) { xyz.y = P->fr_meter * (xyz.y + P->y0); /* z is not scaled since this is handled by vto_meter outside */ + proj_errno_restore (P, last_errno); return xyz; } diff --git a/src/pj_internal.c b/src/pj_internal.c new file mode 100644 index 0000000000..f3ca35665d --- /dev/null +++ b/src/pj_internal.c @@ -0,0 +1,236 @@ +/****************************************************************************** + * 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. + * + * Author: Thomas Knudsen, thokn@sdfe.dk, 2017-07-05 + * + ****************************************************************************** + * 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_INTERNAL_C +#include "proj_internal.h" +#include +#include + +#include +#include +#include + + + + +/* 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 +}; + + +/* Work around non-constness of MSVC HUGE_VAL by providing functions rather than constants */ +PJ_COORD proj_coord_error (void) { + PJ_COORD c; + c.v[0] = c.v[1] = c.v[2] = c.v[3] = HUGE_VAL; + return c; +} + +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; +} + + + +PJ_OBS pj_fwdobs (PJ_OBS obs, PJ *P) { + if (0!=P->fwdobs) { + obs = P->fwdobs (obs, P); + return obs; + } + if (0!=P->fwd3d) { + obs.coo.xyz = pj_fwd3d (obs.coo.lpz, P); + return obs; + } + if (0!=P->fwd) { + obs.coo.xy = pj_fwd (obs.coo.lp, P); + return obs; + } + proj_errno_set (P, EINVAL); + return proj_obs_error (); +} + + +PJ_OBS pj_invobs (PJ_OBS obs, PJ *P) { + if (0!=P->invobs) { + obs = P->invobs (obs, P); + return obs; + } + if (0!=P->inv3d) { + obs.coo.lpz = pj_inv3d (obs.coo.xyz, P); + return obs; + } + if (0!=P->inv) { + obs.coo.lp = pj_inv (obs.coo.xy, P); + return obs; + } + proj_errno_set (P, EINVAL); + return proj_obs_error (); +} + + + +PJ_COORD pj_fwdcoord (PJ_COORD coo, PJ *P) { + if (0!=P->fwdcoord) + return P->fwdcoord (coo, P); + if (0!=P->fwdobs) { + PJ_OBS obs = proj_obs_null; + obs.coo = coo; + obs = P->fwdobs (obs, P); + return obs.coo; + } + if (0!=P->fwd3d) { + coo.xyz = pj_fwd3d (coo.lpz, P); + return coo; + } + if (0!=P->fwd) { + coo.xy = pj_fwd (coo.lp, P); + return coo; + } + proj_errno_set (P, EINVAL); + return proj_coord_error (); +} + + +PJ_COORD pj_invcoord (PJ_COORD coo, PJ *P) { + if (0!=P->invcoord) + return P->invcoord (coo, P); + if (0!=P->invobs) { + PJ_OBS obs = proj_obs_null; + obs.coo = coo; + obs = P->invobs (obs, P); + return obs.coo; + } + if (0!=P->inv3d) { + coo.lpz = pj_inv3d (coo.xyz, P); + return coo; + } + if (0!=P->inv) { + coo.lp = pj_inv (coo.xy, P); + return coo; + } + proj_errno_set (P, EINVAL); + return proj_coord_error (); +} + + + + + +/* Move P to a new context - or to the default context if 0 is specified */ +void proj_context_set (PJ *P, PJ_CONTEXT *ctx) { + if (0==ctx) + ctx = pj_get_default_ctx (); + pj_set_ctx (P, ctx); + return; +} + +void proj_context_inherit (PJ *parent, PJ *child) { + if (0==parent) + pj_set_ctx (child, pj_get_default_ctx()); + else + pj_set_ctx (child, pj_get_ctx(parent)); + return; +} + + + + + +/* stuff below is *not* considered API, and will be moved to an "internal plumbing toolset" */ + + + +void proj_context_errno_set (PJ_CONTEXT *ctx, int err) { + if (0==ctx) + ctx = pj_get_default_ctx(); + pj_ctx_set_errno (ctx, err); + return; +} + + + +/* Set logging level 0-3. Higher number means more debug info. 0 turns it off */ +enum proj_log_level proj_log_level (PJ_CONTEXT *ctx, enum proj_log_level log_level) { + enum proj_log_level previous; + if (0==ctx) + ctx = pj_get_default_ctx(); + if (0==ctx) + return PJ_LOG_TELL; + previous = ctx->debug_level; + if (PJ_LOG_TELL==log_level) + return previous; + ctx->debug_level = log_level; + return previous; +} + + + +/* logging */ + +/* pj_vlog resides in pj_log.c and relates to pj_log as vsprintf relates to sprintf */ +void pj_vlog( projCtx ctx, int level, const char *fmt, va_list args ); + +void proj_log_error (PJ *P, const char *fmt, ...) { + va_list args; + va_start( args, fmt ); + pj_vlog (pj_get_ctx (P), PJ_LOG_ERROR , fmt, args); + va_end( args ); +} + +void proj_log_debug (PJ *P, const char *fmt, ...) { + va_list args; + va_start( args, fmt ); + pj_vlog (pj_get_ctx (P), PJ_LOG_DEBUG_MAJOR , fmt, args); + va_end( args ); +} + +void proj_log_trace (PJ *P, const char *fmt, ...) { + va_list args; + va_start( args, fmt ); + pj_vlog (pj_get_ctx (P), PJ_LOG_DEBUG_MINOR , fmt, args); + va_end( args ); +} + +/* Put a new logging function into P's context. The opaque object app_data is passed as first arg at each call to the logger */ +void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION log) { + if (0==ctx) + pj_get_default_ctx (); + if (0==ctx) + return; + ctx->app_data = app_data; + if (0!=log) + ctx->logger = log; +} diff --git a/src/pj_inv.c b/src/pj_inv.c index d5393261a4..55fc917f19 100644 --- a/src/pj_inv.c +++ b/src/pj_inv.c @@ -1,5 +1,6 @@ /* general inverse projection */ #define PJ_LIB__ +#include #include #include # define EPS 1.0e-12 @@ -8,6 +9,7 @@ LP pj_inv(XY xy, PJ *P) { LP lp; LP err; + int last_errno; /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */ err.lam = err.phi = HUGE_VAL; @@ -20,7 +22,8 @@ LP pj_inv(XY xy, PJ *P) { pj_ctx_set_errno( P->ctx, -15); return err; } - P->ctx->last_errno = errno = pj_errno = 0; + + last_errno = proj_errno_reset (P); /* de-scale and de-offset */ xy.x = (xy.x * P->to_meter - P->x0); @@ -52,5 +55,6 @@ LP pj_inv(XY xy, PJ *P) { lp.phi = atan(P->one_es * tan(lp.phi)); } + proj_errno_restore (P, last_errno); return lp; -} +} \ No newline at end of file diff --git a/src/pj_inv3d.c b/src/pj_inv3d.c index c6052ad4b3..a01cfa7e4a 100644 --- a/src/pj_inv3d.c +++ b/src/pj_inv3d.c @@ -1,4 +1,5 @@ #define PJ_LIB__ +#include #include #include # define EPS 1.0e-12 @@ -8,6 +9,7 @@ LPZ pj_inv3d (XYZ xyz, PJ *P) { LPZ lpz; LPZ err; + int last_errno; /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */ err.lam = err.phi = err.z = HUGE_VAL; @@ -21,7 +23,7 @@ LPZ pj_inv3d (XYZ xyz, PJ *P) { return err; } - P->ctx->last_errno = errno = pj_errno = 0; + last_errno = proj_errno_reset (P); /* de-scale and de-offset */ /* z is not de-scaled since that is handled by vto_meter before we get here */ @@ -53,5 +55,6 @@ LPZ pj_inv3d (XYZ xyz, PJ *P) { lpz.phi = atan(P->one_es * tan(lpz.phi)); } + proj_errno_restore (P, last_errno); return lpz; } diff --git a/src/pj_log.c b/src/pj_log.c index 5b59e7a50d..2525d050b1 100644 --- a/src/pj_log.c +++ b/src/pj_log.c @@ -45,9 +45,9 @@ void pj_stderr_logger( void *app_data, int level, const char *msg ) /************************************************************************/ /* pj_vlog() */ /************************************************************************/ - +void pj_vlog( projCtx ctx, int level, const char *fmt, va_list args ); /* Workhorse for the log functions - relates to pj_log as vsprintf relates to sprintf */ -static void pj_vlog( projCtx ctx, int level, const char *fmt, va_list args ) +void pj_vlog( projCtx ctx, int level, const char *fmt, va_list args ) { char *msg_buf; @@ -84,38 +84,3 @@ void pj_log( projCtx ctx, int level, const char *fmt, ... ) pj_vlog( ctx, level, fmt, args ); va_end( args ); } - - - -/************************************************************************/ -/* logging functions for the proj.h/obs_api */ -/************************************************************************/ - -void pj_log_error (PJ *P, const char *fmt, ...) - -{ - va_list args; - va_start( args, fmt ); - pj_vlog (pj_get_ctx (P), PJ_LOG_ERROR , fmt, args); - va_end( args ); -} - - -void pj_log_debug (PJ *P, const char *fmt, ...) - -{ - va_list args; - va_start( args, fmt ); - pj_vlog (pj_get_ctx (P), PJ_LOG_DEBUG_MAJOR , fmt, args); - va_end( args ); -} - - -void pj_log_trace (PJ *P, const char *fmt, ...) - -{ - va_list args; - va_start( args, fmt ); - pj_vlog (pj_get_ctx (P), PJ_LOG_DEBUG_MINOR , fmt, args); - va_end( args ); -} diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c index 37e55bc341..d8395b0f44 100644 --- a/src/pj_obs_api.c +++ b/src/pj_obs_api.c @@ -34,118 +34,45 @@ *****************************************************************************/ #define PJ_OBS_API_C #include -#include +#include "proj_internal.h" +#include "projects.h" #include -#include -#include +#include +#include -/* Used for zero-initializing new objects */ -const PJ_COORD pj_coo_null = {{0, 0, 0, 0}}; -const PJ_OBS pj_obs_null = { - {{0, 0, 0, 0}}, - {{0, 0, 0}}, - 0, 0 -}; -/* Magic object signaling proj system shutdown mode to routines taking a PJ * arg */ -const PJ *pj_shutdown = (PJ *) &pj_shutdown; +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 between two points with angular 2D coordinates */ -double pj_lp_dist (PJ *P, LP a, LP b) { +double proj_lp_dist (PJ *P, LP a, LP b) { double s12, azi1, azi2; /* Note: the geodesic code takes arguments in degrees */ - geod_inverse (P->geod, TODEG(a.phi), TODEG(a.lam), TODEG(b.phi), TODEG(b.lam), &s12, &azi1, &azi2); + 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 pj_xy_dist (XY a, XY b) { +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 pj_xyz_dist (XYZ a, XYZ b) { +double proj_xyz_dist (XYZ a, XYZ b) { return hypot (hypot (a.x - b.x, a.y - b.y), a.z - b.z); } -/* Work around non-constness of MSVC HUGE_VAL by providing functions rather than constants */ -static PJ_COORD pj_coo_error (void) { - PJ_COORD c; - c.v[0] = c.v[1] = c.v[2] = c.v[3] = HUGE_VAL; - return c; -} - -static PJ_OBS pj_obs_error (void) { - PJ_OBS obs; - obs.coo = pj_coo_error (); - obs.anc.v[0] = obs.anc.v[1] = obs.anc.v[2] = HUGE_VAL; - obs.id = obs.flags = 0; - return obs; -} - - - -static PJ_OBS pj_fwdobs (PJ_OBS obs, PJ *P) { - if (0!=P->fwdobs) { - obs = P->fwdobs (obs, P); - return obs; - } - if (0!=P->fwd3d) { - obs.coo.xyz = pj_fwd3d (obs.coo.lpz, P); - return obs; - } - if (0!=P->fwd) { - obs.coo.xy = pj_fwd (obs.coo.lp, P); - return obs; - } - pj_err_level (P, EINVAL); - return pj_obs_error (); -} - - -static PJ_OBS pj_invobs (PJ_OBS obs, PJ *P) { - if (0!=P->invobs) { - obs = P->invobs (obs, P); - return obs; - } - if (0!=P->inv3d) { - obs.coo.lpz = pj_inv3d (obs.coo.xyz, P); - return obs; - } - if (0!=P->inv) { - obs.coo.lp = pj_inv (obs.coo.xy, P); - return obs; - } - pj_err_level (P, EINVAL); - return pj_obs_error (); -} - - -/* Apply the transformation P to the observation obs */ -PJ_OBS pj_trans (PJ *P, enum 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; - } - - pj_err_level (P, EINVAL); - return pj_obs_error (); -} - /* Measure numerical deviation after n roundtrips fwd-inv (or inv-fwd) */ -double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs) { +double proj_roundtrip (PJ *P, enum proj_direction direction, int n, PJ_OBS obs) { int i; PJ_OBS o, u; @@ -153,7 +80,7 @@ double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs) { return HUGE_VAL; if (n < 1) { - pj_err_level (P, EINVAL); + proj_errno_set (P, EINVAL); return HUGE_VAL; } @@ -170,116 +97,354 @@ double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs) { o = pj_fwdobs (u, P); break; default: - pj_err_level (P, EINVAL); + proj_errno_set (P, EINVAL); return HUGE_VAL; } } - return pj_xyz_dist (o.coo.xyz, obs.coo.xyz); + return proj_xyz_dist (o.coo.xyz, obs.coo.xyz); } -PJ *pj_create (const char *definition) { - return pj_init_plus (definition); -} +/* Apply the transformation P to the coordinate coo */ +PJ_OBS proj_trans_obs (PJ *P, enum proj_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; + } -PJ *pj_create_argv (int argc, char **argv) { - return pj_init (argc, argv); + proj_errno_set (P, EINVAL); + return proj_obs_error (); } -/* Below: Minimum viable support for contexts. The first four functions */ -/* relate to error reporting, debugging, and logging, hence being generically */ -/* useful. The remaining is a compact implementation of the more low level */ -/* proj_api.h thread contexts, which may or may not be useful */ -/* Set error level 0-3, or query current level */ -int pj_err_level (PJ *P, int err_level) { - int previous; - PJ_CONTEXT *ctx; +/* Apply the transformation P to the coordinate coo */ +PJ_COORD proj_trans_coord (PJ *P, enum proj_direction direction, PJ_COORD coo) { if (0==P) - ctx = pj_get_default_ctx(); - else - ctx = pj_get_ctx (P); - previous = pj_ctx_get_errno (ctx); - if (PJ_ERR_TELL==err_level) - return previous; - pj_ctx_set_errno (ctx, err_level); - return previous; -} + 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; + } -/* Set logging level 0-3. Higher number means more debug info. 0 turns it off */ -enum pj_log_level pj_log_level (PJ *P, enum pj_log_level log_level) { - enum pj_log_level previous; - PJ_CONTEXT *ctx; - if (0==P) - ctx = pj_get_default_ctx(); - else - ctx = pj_get_ctx (P); - previous = ctx->debug_level; - if (PJ_LOG_TELL==log_level) - return previous; - ctx->debug_level = log_level; - return previous; + proj_errno_set (P, EINVAL); + return proj_coord_error (); } -/* Put a new logging function into P's context. The opaque object app_data is passed as first arg at each call to the logger */ -void pj_log_func (PJ *P, void *app_data, void (*log)(void *, int, const char *)) { - PJ_CONTEXT *ctx = pj_get_ctx (P); - ctx->app_data = app_data; - if (0!=log) - ctx->logger = log; -} - -/* Move P to a new context - initially a copy of the default context */ -int pj_context_renew (PJ *P) { - PJ_CONTEXT *ctx = pj_ctx_alloc (); - if (0==ctx) { - pj_err_level (P, ENOMEM); - return 1; +/*************************************************************************************/ +size_t proj_transform ( + PJ *P, + enum proj_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; } - pj_set_ctx (P, ctx); - 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; } -/* Move daughter to mother's context */ -void pj_context_inherit (PJ *mother, PJ *daughter) { - pj_set_ctx (daughter, pj_get_ctx (mother)); + +PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) { + if (0==ctx) + ctx = pj_get_default_ctx (); + return pj_init_plus_ctx (ctx, definition); } +PJ *proj_create_argv (PJ_CONTEXT *ctx, int argc, char **argv) { + if (0==ctx) + ctx = pj_get_default_ctx (); + return pj_init_ctx (ctx, argc, argv); +} -void pj_context_free (const PJ *P) { - PJ_CONTEXT *ctx; +PJ *proj_destroy (PJ *P) { + pj_free (P); + return 0; +} +/* For now, if PJ itself is clean, we return the thread local error level. */ +/* This may change as OBS_API error reporting matures */ +int proj_errno (PJ *P) { if (0==P) + return pj_ctx_get_errno (pj_get_default_ctx ()); + if (0 != P->last_errno) + return P->last_errno; + return pj_ctx_get_errno (pj_get_ctx (P)); +} + +/*****************************************************************************/ +void proj_errno_set (PJ *P, int err) { +/****************************************************************************** + Sets errno in the PJ, and bubbles it up to the context and pj_errno levels + through the low level pj_ctx interface. +******************************************************************************/ + if (0==P) { + errno = EINVAL; return; + } + + /* Use proj_errno_reset to explicitly clear the error status */ + if (0==err) + return; + + /* set local error level */ + P->last_errno = err; + /* and let it bubble up */ + proj_context_errno_set (pj_get_ctx (P), err); + errno = err; + return; +} + +/*****************************************************************************/ +void proj_errno_restore (PJ *P, int err) { +/****************************************************************************** + Reduce some mental impedance in the canonical reset/restore use case: + Basically, proj_errno_restore() is a synonym for proj_errno_set(), + but the use cases are very different (_set: indicate an error to higher + level user code, _restore: pass previously set error indicators in case + of no errors at this level). + + Hence, although the inner working is identical, we provide both options, + to avoid some rather confusing real world code. + + See usage example under proj_errno_reset () +******************************************************************************/ + proj_errno_set (P, err); +} - /* During shutdown it should be OK to free the default context - but this needs more work */ - if (pj_shutdown==P) { - /* The obvious solution "pj_ctx_free (pj_get_default_ctx ())" fails */ +/*****************************************************************************/ +int proj_errno_reset (PJ *P) { +/****************************************************************************** + Clears errno in the PJ, and bubbles it up to the context and + pj_errno levels through the low level pj_ctx interface. + + Returns the previous value of the errno, for convenient reset/restore + operations: + + void foo (PJ *P) { + int last_errno = proj_errno_reset (P); + + do_something_with_P (P); + + (* failure - keep latest error status *) + if (proj_errno(P)) + return; + (* success - restore previous error status *) + proj_errno_restore (P, last_errno); return; } +******************************************************************************/ + int last_errno; + if (0==P) { + errno = EINVAL; + return EINVAL; + } + last_errno = proj_errno (P); + + /* set local error level */ + P->last_errno = 0; + /* and let it bubble up */ + pj_ctx_set_errno (pj_get_ctx (P), 0); + errno = 0; + return last_errno; +} - ctx = pj_get_ctx ((PJ *) P); - /* Otherwise, trying to free the default context is a no-op */ +/* Create a new context - or provide a pointer to the default context */ +PJ_CONTEXT *proj_context_create (void) { + return pj_ctx_alloc (); +} + + +void proj_context_destroy (PJ_CONTEXT *ctx) { + if (0==ctx) + return; + + /* Trying to free the default context is a no-op (since it is statically allocated) */ if (pj_get_default_ctx ()==ctx) return; - /* The common (?) case: free the context and move the PJ to the default context */ pj_ctx_free (ctx); - pj_set_ctx ((PJ *) P, pj_get_default_ctx ()); } -/* Minimum support for fileapi - which may never have been used... */ -void pj_fileapi_set (PJ *P, void *fileapi) { - PJ_CONTEXT *ctx = pj_get_ctx (P); - ctx->fileapi = (projFileAPI *) fileapi; +/* Build a fully expanded proj_create() compatible representation of P */ +char *proj_definition_retrieve (PJ *P) { + if (0==P) + return 0; + return pj_get_def(P, 0); } + +/* ...and get rid of it safely */ +void *proj_release (void *buffer) { + return pj_dealloc (buffer); +} + + +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);} + + +/* The shape of jazz to come! */ + +int proj_transform_obs (PJ *P, enum proj_direction direction, size_t n, PJ_OBS *obs); +int proj_transform_coord (PJ *P, enum proj_direction direction, size_t n, PJ_COORD *coord); +PJ_OBS proj_obs (double x, double y, double z, double t, double o, double p, double k, int id, unsigned int flags); +PJ *proj_create_crs_to_crs (PJ_CONTEXT *ctx, const char *def_from, const char *def_to); diff --git a/src/proj.def b/src/proj.def index 3792d3c032..314c438536 100644 --- a/src/proj.def +++ b/src/proj.def @@ -89,21 +89,45 @@ EXPORTS geod_polygon_testpoint @87 geod_polygon_clear @88 pj_run_selftests @89 - pj_create @90 - pj_create_argv @91 - pj_trans @92 - pj_roundtrip @95 - pj_log_level @96 - pj_err_level @97 - pj_log_func @98 - pj_context_renew @99 - pj_context_inherit @100 - pj_context_free @101 - pj_fileapi_set @102 - pj_log_error @103 - pj_log_debug @104 - pj_log_trace @105 - pj_lp_dist @106 - pj_xy_dist @107 - pj_xyz_dist @108 - pj_find_file @109 + + proj_create @90 + proj_create_argv @91 + proj_destroy @92 + + proj_trans_obs @93 + proj_trans_coord @94 + proj_transform @95 + proj_roundtrip @96 + + proj_coord @97 + proj_coord_error @98 + proj_obs_error @99 + + proj_errno @100 + proj_errno_set @101 + proj_errno_reset @102 + proj_errno_restore @103 + proj_context_errno_set @104 + + proj_context_create @105 + proj_context_set @106 + proj_context_inherit @107 + proj_context_destroy @108 + + proj_lp_dist @109 + proj_xy_dist @110 + proj_xyz_dist @111 + + proj_log_level @112 + proj_log_func @113 + proj_log_error @114 + proj_log_debug @115 + proj_log_trace @116 + + proj_definition_retrieve @117 + proj_release @118 + proj_torad @119 + proj_todeg @120 + + + pj_find_file @121 diff --git a/src/proj.h b/src/proj.h index f6b3bd2a35..1977acde36 100644 --- a/src/proj.h +++ b/src/proj.h @@ -14,7 +14,7 @@ * * This is an evolving attempt at creating a re-rationalized API * with primary design goals focused on sanitizing the namespaces. - * Hence, all symbols exposed are being moved to the pj_ namespace, + * Hence, all symbols exposed are being moved to the proj_ namespace, * while all data types are being moved to the PJ_ namespace. * * Please note that this API is *orthogonal* to the previous APIs: @@ -49,7 +49,8 @@ * without having any idea about its internals. * * (obviously, in this example, struct projCtx_t should also be - * renamed struct pj_ctx some day...) + * renamed struct pj_ctx some day, but for backwards compatibility + * 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 @@ -73,20 +74,30 @@ * as large as the traditional XY and LP objects, timing results * have shown the overhead to be very reasonable. * - * In its current incarnation, the API is focused on compactness: - * Currently it consists of only nine functions. + * Contexts and thread safety + * -------------------------- * - * Hence, due to the proj_ctx subsystem being little used, it has - * been replaced by a tiny set of 3 functions, making it possible, - * but not very convenient, to do almost everything possible with - * the original ctx API. + * After a year of experiments (and previous experience from the + * trlib transformation library) it has become clear that the + * 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. * - * See pj_proj_test.c for an example of how to use the API. + * Legacy single threaded code need not do anything, but *may* + * implement a bit of future compatibility by using the backward + * 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. * * Author: Thomas Knudsen, + * Benefitting from a large number of comments and suggestions + * by (primarily) Kristian Evers and Even Rouault. * ****************************************************************************** - * Copyright (c) 2016, Thomas Knudsen / SDFE + * 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"), @@ -106,18 +117,15 @@ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER * DEALINGS IN THE SOFTWARE. *****************************************************************************/ +/* #ifdef _MSC_VER +#ifndef _USE_MATH_DEFINES #define _USE_MATH_DEFINES #endif -#include - -#include -#include +#endif +#include For M_PI */ +#include /* For size_t */ -#include -#include -#include -#include #ifdef PROJECTS_H #error proj.h must be included before projects.h @@ -132,14 +140,33 @@ extern "C" { #endif +/************************************************************************ + These version numbers should be updated with every release! + The format of PJ_VERSION is + + * Before version 4.10.0: PJ_VERSION=MNP + where M, N, and P are the major, minor, and patch numbers; + e.g., PJ_VERSION=493 for version 4.9.3. + + * Version 4.10.0 and later: PJ_VERSION=MMMNNNPP later + where MMM, NNN, PP are the major, minor, and patch numbers. + The minor and patch numbers are padded with leading zeros if + necessary); + e.g., PJ_VERSION=401000 for version 4.10.0. +************************************************************************/ +#define PROJ_VERSION_MAJOR 10 +#define PROJ_VERSION_MINOR 0 +#define PROJ_VERSION_MICRO 0 + +#ifndef PJ_VERSION +#define PJ_VERSION 10##000##00 +#endif +#ifndef PROJ_VERSION +#define PROJ_VERSION PJ_VERSION +#endif -/* We need to access the PJ_VERSION symbol from proj_api.h */ -#define PROJ_API_INCLUDED_FOR_PJ_VERSION_ONLY -#include -#undef PROJ_API_INCLUDED_FOR_PJ_VERSION_ONLY extern char const pj_release[]; /* global release id string */ -extern int pj_errno; /* global error return code */ /* first forward declare everything needed */ @@ -182,6 +209,9 @@ typedef struct { double lam, phi, z; } LPZ; /* Ancillary pairs and triplets for geodetic computations */ +/* easting and northing */ +typedef struct { double e, n; } PJ_EN; + /* Degrees, minutes, and seconds */ typedef struct { double d, m, s; } PJ_DMS; @@ -198,6 +228,7 @@ union PJ_COORD { PJ_ENHT enht; PJ_LPZT lpzt; PJ_ENH enh; + PJ_EN en; double v[4]; /* It's just a vector */ XYZ xyz; UVW uvw; @@ -227,6 +258,7 @@ union PJ_PAIR { LP lp; UV uv; PJ_AF af; + PJ_EN en; double v[2]; /* Yes - It's really just a vector! */ }; @@ -240,91 +272,79 @@ struct PJ_OBS { /* The context type - properly namespaced synonym for projCtx */ struct projCtx_t; typedef struct projCtx_t PJ_CONTEXT; -typedef int *PJ_FILE; + + + +/* A P I */ + + +/* Functionality for handling thread contexts */ +PJ_CONTEXT *proj_context_create (void); +void proj_context_destroy (PJ_CONTEXT *ctx); /* Manage the transformation definition object PJ */ -PJ *pj_create (const char *definition); -PJ *pj_create_argv (int argc, char **argv); -void pj_free (PJ *P); -int pj_error (PJ *P); +PJ *proj_create (PJ_CONTEXT *ctx, const char *definition); +PJ *proj_create_argv (PJ_CONTEXT *ctx, int argc, char **argv); +PJ *proj_destroy (PJ *P); /* Apply transformation to observation - in forward or inverse direction */ -enum pj_direction { +enum proj_direction { PJ_FWD = 1, /* Forward */ PJ_IDENT = 0, /* Do nothing */ PJ_INV = -1 /* Inverse */ }; -PJ_OBS pj_trans (PJ *P, enum pj_direction direction, PJ_OBS obs); -/* Measure internal consistency - in forward or inverse direction */ -double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs); -/* Geodesic distance between two points with angular 2D coordinates */ -double pj_lp_dist (PJ *P, LP a, LP b); +PJ_OBS proj_trans_obs (PJ *P, enum proj_direction direction, PJ_OBS obs); +PJ_COORD proj_trans_coord (PJ *P, enum proj_direction direction, PJ_COORD coord); -/* Euclidean distance between two points with linear 2D coordinates */ -double pj_xy_dist (XY a, XY b); -/* Euclidean distance between two points with linear 3D coordinates */ -double pj_xyz_dist (XYZ a, XYZ b); +size_t proj_transform ( + PJ *P, + enum proj_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 +); -#ifndef PJ_OBS_API_C -/* Part of MSVC workaround: Make pj_*_null look function-like for symmetry with pj_*_error */ -#define pj_coo_null(x) pj_coo_null -#define pj_obs_null(x) pj_obs_null -extern const PJ_COORD pj_coo_null; -extern const PJ_OBS pj_obs_null; -extern const PJ *pj_shutdown; -#endif - -#ifndef TODEG -#define TODEG(rad) ((rad)*180.0/M_PI) -#endif -#ifndef TORAD -#define TORAD(deg) ((deg)*M_PI/180.0) -#endif +/* not a constructor, but an initializer */ +PJ_COORD proj_coord (double x, double y, double z, double t); +/* Measure internal consistency - in forward or inverse direction */ +double proj_roundtrip (PJ *P, enum proj_direction direction, int n, PJ_OBS obs); +/* Geodesic distance between two points with angular 2D coordinates */ +double proj_lp_dist (PJ *P, LP a, LP b); +/* Euclidean distance between two points with linear 2D coordinates */ +double proj_xy_dist (XY a, XY b); +/* Euclidean distance between two points with linear 3D coordinates */ +double proj_xyz_dist (XYZ a, XYZ b); -/* High level functionality for handling thread contexts */ -enum pj_log_level { - PJ_LOG_NONE = 0, - PJ_LOG_ERROR = 1, - PJ_LOG_DEBUG = 2, - PJ_LOG_TRACE = 3, - PJ_LOG_TELL = 4, - PJ_LOG_DEBUG_MAJOR = 2, /* for proj_api.h compatibility */ - PJ_LOG_DEBUG_MINOR = 3 /* for proj_api.h compatibility */ -}; - /* Set or read error level */ -#define PJ_ERR_TELL -56789 -int pj_err_level (PJ *P, int err_level); - -/* Set logging level 0-3. Higher number means more debug info. 0 turns it off */ -enum pj_log_level pj_log_level (PJ *P, enum pj_log_level log_level); - -void pj_log_error (PJ *P, const char *fmt, ...); -void pj_log_debug (PJ *P, const char *fmt, ...); -void pj_log_trace (PJ *P, const char *fmt, ...); -void pj_log_func (PJ *P, void *app_data, void (*log)(void *, int, const char *)); - - -/* Lower level functionality for handling thread contexts */ -int pj_context_renew (PJ *P); -void pj_context_inherit (PJ *mother, PJ *daughter); -void pj_context_free (const PJ *P); - -/* Lowest level: Minimum support for fileapi */ -void pj_fileapi_set (PJ *P, void *fileapi); - +int proj_errno (PJ *P); +void proj_errno_set (PJ *P, int err); +int proj_errno_reset (PJ *P); +void proj_errno_restore (PJ *P, int err); + +/* Build a fully expanded proj_create() compatible representation of P */ +char *proj_definition_retrieve (PJ *P); +/* ...and get rid of it safely */ +void *proj_release (void *buffer); + + +/* These are trivial, and while occasionaly 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); +double proj_todeg (double angle_in_radians); #ifdef __cplusplus } diff --git a/src/proj_internal.h b/src/proj_internal.h new file mode 100644 index 0000000000..da457077aa --- /dev/null +++ b/src/proj_internal.h @@ -0,0 +1,115 @@ +/****************************************************************************** + * Project: PROJ.4 + * Purpose: Internal plumbing for the PROJ.4 library. + * + * Author: Thomas Knudsen, + * + ****************************************************************************** + * 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 COORD 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. + *****************************************************************************/ + +#ifdef _MSC_VER +#ifndef _USE_MATH_DEFINES +#define _USE_MATH_DEFINES +#endif +#endif +#include /* For M_PI */ + +#include + +#ifdef PROJECTS_H +#error proj_internal.h must be included before projects.h +#endif +#ifdef PROJ_API_H +#error proj_internal.h must be included before proj_api.h +#endif + +#ifndef PROJ_INTERNAL_H +#define PROJ_INTERNAL_H +#ifdef __cplusplus +extern "C" { +#endif + + + +#ifndef PJ_TODEG +#define PJ_TODEG(rad) ((rad)*180.0/M_PI) +#endif +#ifndef PJ_TORAD +#define PJ_TORAD(deg) ((deg)*M_PI/180.0) +#endif + + + +PJ_COORD proj_coord_error (void); +PJ_OBS proj_obs_error (void); +#ifndef PJ_INTERNAL_C +extern const PJ_COORD proj_coord_null; +extern const PJ_OBS proj_obs_null; +#endif +/* Part of MSVC workaround: Make proj_*_null look function-like for symmetry with proj_*_error */ +#define proj_coord_null(x) proj_coord_null +#define proj_obs_null(x) proj_obs_null + + +void proj_context_errno_set (PJ_CONTEXT *ctx, int err); +void proj_context_set (PJ *P, PJ_CONTEXT *ctx); +void proj_context_inherit (PJ *parent, PJ *child); + + +PJ_OBS pj_fwdobs (PJ_OBS obs, PJ *P); +PJ_OBS pj_invobs (PJ_OBS obs, PJ *P); +PJ_COORD pj_fwdcoord (PJ_COORD coo, PJ *P); +PJ_COORD pj_invcoord (PJ_COORD coo, PJ *P); + + + +/* High level functionality for handling thread contexts */ +enum proj_log_level { + PJ_LOG_NONE = 0, + PJ_LOG_ERROR = 1, + PJ_LOG_DEBUG = 2, + PJ_LOG_TRACE = 3, + PJ_LOG_TELL = 4, + PJ_LOG_DEBUG_MAJOR = 2, /* for proj_api.h compatibility */ + PJ_LOG_DEBUG_MINOR = 3 /* for proj_api.h compatibility */ +}; + +/* Set logging level 0-3. Higher number means more debug info. 0 turns it off */ +enum proj_log_level proj_log_level (PJ_CONTEXT *ctx, enum proj_log_level log_level); +typedef void (*PJ_LOG_FUNCTION)(void *, int, const char *); + +void proj_log_error (PJ *P, const char *fmt, ...); +void proj_log_debug (PJ *P, const char *fmt, ...); +void proj_log_trace (PJ *P, const char *fmt, ...); +/*void proj_log_func (PJ *P, void *app_data, void (*log)(void *, int, const char *));*/ +void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION log); + + +/* Lowest level: Minimum support for fileapi */ +void proj_fileapi_set (PJ *P, void *fileapi); + + +#ifdef __cplusplus +} +#endif + +#endif /* ndef PROJ_INTERNAL_H */ diff --git a/src/projects.h b/src/projects.h index c1bf1b3a25..1c99506f21 100644 --- a/src/projects.h +++ b/src/projects.h @@ -175,6 +175,7 @@ typedef struct { double u, v, w; } UVW; /* Forward declarations and typedefs for stuff needed inside the PJ object */ struct PJconsts; struct PJ_OBS; +union PJ_COORD; struct geod_geodesic; struct pj_opaque; struct ARG_list; @@ -190,6 +191,7 @@ 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 struct PJ_REGION_S { @@ -208,11 +210,12 @@ struct PJconsts { /************************************************************************************* - G E N E R A L C O N T E X T + G E N E R A L P A R A M E T E R S T R U C T ************************************************************************************** TODO: Need some description here - especially about the thread context... + This is the struct behind the PJ typedef **************************************************************************************/ @@ -245,11 +248,67 @@ struct PJconsts { LPZ (*inv3d)(XYZ, PJ *); PJ_OBS (*fwdobs)(PJ_OBS, PJ *); PJ_OBS (*invobs)(PJ_OBS, PJ *); + PJ_COORD (*fwdcoord)(PJ_COORD, PJ *); + PJ_COORD (*invcoord)(PJ_COORD, PJ *); void (*spc)(LP, PJ *, struct FACTORS *); void (*pfree)(PJ *); + /************************************************************************************* + + E R R O R R E P O R T I N G + + ************************************************************************************** + + Currently, we're doing error reporting through the context->last_errno indicator. + + It is, however, not entirely sure this will be the right way to do it in all + cases: During allocation/initialization, it is certainly nice to have a higher + level error indicator, since we primarily signal "something went wrong", by + returning 0 from the pj_init family of functions - and with a null return we + cannot pass messages through internal state in the PJ object. + + Historically, the errno variable has been used for that kind of messages, but + apparently, thread safety was added to PROJ.4 at a time where it was not clear + that errno is actually thread local. + + Additionally, errno semantics has historically been misinterpreted in parts of + pj_init.c, a misinterpretation, that was mitigated by a hack in pj_malloc.c + some 15 years ago. + + This PJ-local errno is a first step towards a more structured approach to + error reporting, being implemented in the OBS_API (cf. pj_obs_api.[ch]), and + related plumbing efforts. + + In due course this will let us get rid of the pj_malloc.c hack, and allow us + to introduce a more layered error reporting structure, where errors are + reported where they occur, and bubble up to the higher levels (context->errno, + then thread local errno), so the highest level indicate "something went wrong + somewhere", then the more localized ones can be used for pinpointing: + + errno: "something went wrong somewhere on this thread", + context->last_errno: "something went wrong in some PROJ.4 related code", + PJ->last_errno: "It was in THIS PJ something went wrong", + pj_strerrno: "This was what went wrong". + + Which will be quite helpful, once fully implemented, especially for + debugging complex transformation pipelines, while still maintaining backward + compatibility in the messaging system. + + Note that there is even a global pj_errno, which is here and there accessed + without acquiring lock. This, and the practise of resetting the thread local + errno, should be given some consideration during the cleanup of the error + reporting system. + + The name "last_errno", rather than "errno" is used partially for alignment + with the context->last_errno, partially because in a multithreaded environment, + errno is a macro, and spurious spaces turning "errno" into a separate token + will expose it to macro expansion, to the tune of much confusion and agony. + + **************************************************************************************/ + int last_errno; + /************************************************************************************* @@ -345,6 +404,11 @@ struct PJconsts { D A T U M S A N D H E I G H T S Y S T E M S + ************************************************************************************** + + It may be possible, and meaningful, to move the list parts of this up to the + PJ_CONTEXT level. + **************************************************************************************/ int datum_type; /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */ @@ -374,9 +438,15 @@ struct PJconsts { struct _pj_gi *last_after_grid; /* TODO: Description needed */ PJ_Region last_after_region; /* TODO: Description needed */ double last_after_date; /* TODO: Description needed */ -}; + /************************************************************************************* + + E N D O F G E N E R A L P A R A M E T E R S T R U C T + + **************************************************************************************/ + +}; @@ -389,7 +459,6 @@ struct ARG_list { }; - typedef union { double f; int i; char *s; } PROJVALUE; @@ -404,6 +473,7 @@ struct PJ_ELLPS { char *ell; /* elliptical parameter */ char *name; /* comments */ }; + struct PJ_UNITS { char *id; /* units keyword */ char *to_meter; /* multiply by value to get meters */ @@ -490,8 +560,6 @@ extern struct PJ_LIST pj_list[]; extern struct PJ_SELFTEST_LIST pj_selftest_list[]; #endif - - #ifndef PJ_ELLPS__ extern struct PJ_ELLPS pj_ellps[]; #endif