From 26e65892728d6bfb010677091b0a464d7c7ebc76 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Sat, 12 Nov 2016 16:51:21 +0100 Subject: [PATCH] Pipeline plus api - in continuation of the Genereic Coordinates pull request (#445) This commit reflects continued work on the new rationalized API with geodetic extensions (see rationale in proj.h). It also reflects the parallel work on the transformation pipeline reimplementation, by introducing the PJ_cart cartesian/ellipsoidal converter. See example/pj_obs_api_test.c for demo example code --- examples/pj_obs_api_test.c | 132 +++++++++++++++ src/Makefile.am | 2 +- src/PJ_cart.c | 336 +++++++++++++++++++++++++++++++++++++ src/lib_proj.cmake | 3 +- src/makefile.vc | 2 +- src/pj_ctx.c | 2 - src/pj_list.h | 1 + src/pj_malloc.c | 128 ++++++++++++-- src/pj_obs_api.c | 248 +++++++++++++++++++++++++++ src/pj_observation.c | 114 ------------- src/proj.def | 18 +- src/proj.h | 199 ++++++++++++++++------ src/proj_api.h | 102 +++++++---- src/projects.h | 16 +- 14 files changed, 1072 insertions(+), 231 deletions(-) create mode 100644 examples/pj_obs_api_test.c create mode 100644 src/PJ_cart.c create mode 100644 src/pj_obs_api.c delete mode 100644 src/pj_observation.c diff --git a/examples/pj_obs_api_test.c b/examples/pj_obs_api_test.c new file mode 100644 index 0000000000..49d1224e24 --- /dev/null +++ b/examples/pj_obs_api_test.c @@ -0,0 +1,132 @@ +/******************************************************************************* + 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_apply (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 main (void) { + PJ *p; + PJ_OBS a, b; + int err; + double dist; + + /* Log everything libproj offers to log for you */ + pj_debug_set (0, PJ_LOG_DEBUG_MINOR); + + /* An utm projection on the GRS80 ellipsoid */ + p = pj_create ("+proj=utm +zone=32 +ellps=GRS80"); + 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_apply (p, PJ_FWD, a); + printf ("FWD: %15.9f %15.9f\n", b.coo.enh.e, b.coo.enh.n); + + /* Inverse projection */ + a = pj_apply (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_apply (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_apply (p, PJ_FWD, a); + printf ("FWD: %15.9f %15.9f\n", b.coo.enh.e, b.coo.enh.n); + + dist = pj_obs_dist_2d (a, b); + printf ("Roundtrip deviation, (nm): %15.9f\n", dist*1e9); + + /* Invalid projection */ + a = pj_apply (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_error (p); + printf ("pj_error: %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_apply (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_apply (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_apply (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_apply (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); + return 0; +} diff --git a/src/Makefile.am b/src/Makefile.am index 659441117f..1d53375a90 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -74,7 +74,7 @@ libproj_la_SOURCES = \ jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c geodesic.c \ pj_strtod.c \ \ - pj_observation.c + pj_obs_api.c PJ_cart.c install-exec-local: rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT) diff --git a/src/PJ_cart.c b/src/PJ_cart.c new file mode 100644 index 0000000000..5df95c31d9 --- /dev/null +++ b/src/PJ_cart.c @@ -0,0 +1,336 @@ +/****************************************************************************** + * Project: PROJ.4 + * Purpose: Convert between ellipsoidal, geodetic coordinates and + * cartesian, geocentric coordinates. + * + * Formally, this functionality is also found in the PJ_geocent.c + * code. + * + * Actually, however, the PJ_geocent transformations are carried + * out in concert between 2D stubs in PJ_geocent.c and 3D code + * placed in pj_transform.c. + * + * For pipeline-style datum shifts, we do need direct access + * to the full 3D interface for this functionality. + * + * Hence this code, which may look like "just another PJ_geocent" + * but really is something substantially different. + * + * Author: Thomas Knudsen, thokn@sdfe.dk + * + ****************************************************************************** + * Copyright (c) 2016, 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_LIB__ +#include +#include +#include +#include +#include +#include +PROJ_HEAD(cart, "Geodetic/cartesian conversions"); + + +/************************************************************** + CARTESIAN / GEODETIC CONVERSIONS +*************************************************************** + This material follows: + + Bernhard Hofmann-Wellenhof & Helmut Moritz: + Physical Geodesy, 2nd edition. + Springer, 2005. + + chapter 5.6: Coordinate transformations + (HM, below), + + and + + Wikipedia: Geographic Coordinate Conversion, + https://en.m.wikipedia.org/wiki/Geographic_coordinate_conversion + + (WP, below). + + The cartesian-to-geodetic conversion is based on Bowring's + celebrated method: + + B. R. Bowring: + Transformation from spatial to geographical coordinates + Survey Review 23(181), pp. 323-327, 1976 + + (BB, below), + + but could probably use some TLC from a newer and faster + algorithm: + + Toshio Fukushima: + Transformation from Cartesian to Geodetic Coordinates + Accelerated by Halley’s Method + Journal of Geodesy, February 2006 + + (TF, below). + + + These routines are probably not as robust at those in + geocent.c, at least thay haven't been through as heavy + use as their geocent sisters. Some care has been taken + to avoid singularities, but extreme cases (e.g. setting + es, the squared eccentricity, to 1), will cause havoc. + +**************************************************************/ + + +static void freeup (PJ *P) { + pj_freeup_plain (P); + return; +} + + +/*********************************************************************/ +static double semiminor_axis (double a, double es) { +/*********************************************************************/ + return a * sqrt (1 - es); +} + +/*********************************************************************/ +static double normal_radius_of_curvature (double a, double es, double phi) { +/*********************************************************************/ + double s = sin(phi); + if (es==0) + return a; + /* This is from WP. HM formula 2-149 gives an a,b version */ + return a / sqrt (1 - es*s*s); +} + + +/*********************************************************************/ +static double second_eccentricity_squared (double a, double es) { +/********************************************************************** + Follows the definition, e'2 = (a2-b2)/b2, but uses the identity + (a2-b2) = (a-b)(a+b), for improved numerical precision. +***********************************************************************/ + double b = semiminor_axis (a, es); + return (a - b) * (a + b) / (b*b); +} + + +/*********************************************************************/ +static XYZ cartesian (LPZ geodetic, PJ *P) { +/*********************************************************************/ + double N, h, lam, phi, cosphi = cos(geodetic.phi); + XYZ xyz; + + N = normal_radius_of_curvature(P->a, P->es, geodetic.phi); + phi = geodetic.phi; + lam = geodetic.lam; + h = geodetic.z; + + /* HM formula 5-27 (z formula follows WP) */ + xyz.x = (N + h) * cosphi * cos(lam); + xyz.y = (N + h) * cosphi * sin(lam); + xyz.z = (N * (1 - P->es) + h) * sin(phi); + + /*********************************************************************/ + /* */ + /* For historical reasons, in proj, plane coordinates are measured */ + /* in units of the semimajor axis. Since 3D handling is grafted on */ + /* later, this is not the case for heights. And even though this */ + /* coordinate really is 3D cartesian, the z-part looks like a height */ + /* to proj. Hence, we have the somewhat unusual situation of having */ + /* a point coordinate with differing units between dimensions. */ + /* */ + /* The scaling and descaling is handled by the pj_fwd/inv functions. */ + /* */ + /*********************************************************************/ + xyz.x /= P->a; + xyz.y /= P->a; + + return xyz; +} + + +/*********************************************************************/ +static LPZ geodetic (XYZ cartesian, PJ *P) { +/*********************************************************************/ + double N, b, p, theta, c, s, e2s; + LPZ lpz; + + cartesian.x *= P->a; + cartesian.y *= P->a; + + /* Perpendicular distance from point to Z-axis (HM eq. 5-28) */ + p = hypot (cartesian.x, cartesian.y); + + /* Ancillary ellipsoidal parameters */ + b = semiminor_axis (P->a, P->es); + e2s = second_eccentricity_squared (P->a, P->es); + + /* HM eq. (5-37) */ + theta = atan2 (cartesian.z * P->a, p * b); + + /* HM eq. (5-36) (from BB, 1976) */ + c = cos(theta); + s = sin(theta); + lpz.phi = atan2 (cartesian.z + e2s*b*s*s*s, p - P->es*P->a*c*c*c); + lpz.lam = atan2 (cartesian.y, cartesian.x); + N = normal_radius_of_curvature (P->a, P->es, lpz.phi); + + + c = cos(lpz.phi); + if (fabs(c) < 1e-6) { + /* poleward of 89.99994 deg, we avoid division by zero */ + /* by computing the height as the cartesian z value */ + /* minus the semiminor axis length */ + lpz.z = cartesian.z - (cartesian.z > 0? b: -b); + } + else + lpz.z = p / c - N; + + return lpz; +} + + + +/* In effect, 2 cartesian coordinates of a point on the ellipsoid. Rather pointless, but... */ +static XY cart_forward (LP lp, PJ *P) { + PJ_TRIPLET point; + point.lp = lp; + point.lpz.z = 0; + + point.xyz = cartesian (point.lpz, P); + return point.xy; +} + +/* And the other way round. Still rather pointless, but... */ +static LP cart_reverse (XY xy, PJ *P) { + PJ_TRIPLET point; + point.xy = xy; + point.xyz.z = 0; + + point.lpz = geodetic (point.xyz, P); + return point.lp; +} + + + +/*********************************************************************/ +PJ *PROJECTION(cart) { +/*********************************************************************/ + P->fwd3d = cartesian; + P->inv3d = geodetic; + P->fwd = cart_forward; + P->inv = cart_reverse; + return P; +} + +#ifndef PJ_SELFTEST +/* selftest stub */ +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 *p; + PJ_OBS a, b; + int err; + double dist; + + /* Log everything libproj offers to log for you */ + pj_debug_set (0, PJ_LOG_DEBUG_MINOR); + + /* An utm projection on the GRS80 ellipsoid */ + p = pj_create ("+proj=utm +zone=32 +ellps=GRS80"); + if (0==p) + return 1; + + /* Turn off logging */ + pj_debug_set (0, PJ_LOG_NONE); + + /* 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_apply (p, PJ_FWD, a); + + /* Inverse projection */ + a = pj_apply (p, PJ_INV, b); + + /* Null projection */ + a = pj_apply (p, PJ_IDENT, a); + + /* Forward again, to get two linear items for comparison */ + a = pj_apply (p, PJ_FWD, a); + + dist = pj_obs_dist_2d (a, b); + if (dist > 2e-9) + return 2; + + /* Invalid projection */ + a = pj_apply (p, 42, a); + if (a.coo.lpz.lam!=DBL_MAX) + return 3; + err = pj_error (p); + if (0==err) + return 4; + + /* Clear error */ + pj_error_set (p, 0); + + /* Clean up */ + pj_free (p); + + /* Now do some 3D transformations */ + p = pj_create ("+proj=cart +ellps=GRS80"); + if (0==p) + return 5; + + /* 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_apply (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); + if (dist > 2e-9) + return 6; + + /* Inverse projection: Ellipsoidal-to-3D-Cartesian */ + b = pj_apply (p, PJ_INV, b); + + /* Move p to another context */ + pj_context_renew (p); + b = pj_apply (p, PJ_FWD, b); + + /* Move it back to the default context */ + pj_context_free (p); + b = pj_apply (p, PJ_INV, b); + + pj_free (p); + return 0; +} +#endif diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 54e59ae771..2386a102b6 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -50,6 +50,7 @@ SET(SRC_LIBPROJ_PJ PJ_boggs.c PJ_bonne.c PJ_calcofi.c + PJ_cart.c PJ_cass.c PJ_cc.c PJ_cea.c @@ -193,7 +194,7 @@ SET(SRC_LIBPROJ_CORE pj_mlfn.c pj_msfn.c pj_mutex.c - pj_observation.c + pj_obs_api.c pj_open_lib.c pj_param.c pj_phi2.c diff --git a/src/makefile.vc b/src/makefile.vc index f06a769a3d..51b4320b69 100644 --- a/src/makefile.vc +++ b/src/makefile.vc @@ -59,7 +59,7 @@ support = \ pj_strtod.obj pj_run_selftests.obj pj_generic_selftest.obj pipeline = \ - pj_observation.obj + pj_obs_api.obj PJ_cart.obj geodesic = geodesic.obj diff --git a/src/pj_ctx.c b/src/pj_ctx.c index 43d12255cd..9ae3d4e0e3 100644 --- a/src/pj_ctx.c +++ b/src/pj_ctx.c @@ -192,5 +192,3 @@ projFileAPI *pj_ctx_get_fileapi( projCtx ctx ) { return ctx->fileapi; } - - diff --git a/src/pj_list.h b/src/pj_list.h index a6806556ef..3569e335b3 100644 --- a/src/pj_list.h +++ b/src/pj_list.h @@ -17,6 +17,7 @@ PROJ_HEAD(bipc, "Bipolar conic of western hemisphere") PROJ_HEAD(boggs, "Boggs Eumorphic") PROJ_HEAD(bonne, "Bonne (Werner lat_1=90)") PROJ_HEAD(calcofi, "Cal Coop Ocean Fish Invest Lines/Stations") +PROJ_HEAD(cart, "Geodetic/cartesian conversions") PROJ_HEAD(cass, "Cassini") PROJ_HEAD(cc, "Central Cylindrical") PROJ_HEAD(cea, "Equal Area Cylindrical") diff --git a/src/pj_malloc.c b/src/pj_malloc.c index aab69e9917..b9d8f825be 100644 --- a/src/pj_malloc.c +++ b/src/pj_malloc.c @@ -1,28 +1,68 @@ +/****************************************************************************** + * Project: PROJ.4 + * Purpose: Memory management for proj.4. + * This version includes an implementation of generic destructors, + * for memory deallocation for the large majority of PJ-objects + * that do not allocate anything alse than the PJ-object itself, + * and its associated opaque object - i.e. no additional malloc'ed + * memory inside the opaque object. + * + * Author: Gerald I. Evenden (Original proj.4 author), + * Frank Warmerdam (2000) pj_malloc? + * Thomas Knudsen (2016) - freeup/dealloc parts + * + ****************************************************************************** + * Copyright (c) 2000, Frank Warmerdam + * Copyright (c) 2016, 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. + *****************************************************************************/ + + /* allocate and deallocate memory */ /* These routines are used so that applications can readily replace ** projection system memory allocation/deallocation call with custom ** application procedures. */ + #include #include - void * -pj_malloc(size_t size) { -/* -/ Currently, pj_malloc is a hack to solve an errno problem. -/ The problem is described in more details at -/ https://bugzilla.redhat.com/bugzilla/show_bug.cgi?id=86420. -/ It seems, that pj_init and similar functions incorrectly -/ (under debian/glibs-2.3.2) assume that pj_malloc resets -/ errno after success. pj_malloc tries to mimic this. -*/ +/**********************************************************************/ +void *pj_malloc(size_t size) { +/*********************************************************************** +Currently, pj_malloc is a hack to solve an errno problem. +The problem is described in more details at +https://bugzilla.redhat.com/bugzilla/show_bug.cgi?id=86420. +It seems, that pj_init and similar functions incorrectly +(under debian/glibs-2.3.2) assume that pj_malloc resets +errno after success. pj_malloc tries to mimic this. +***********************************************************************/ int old_errno = errno; void *res = malloc(size); if ( res && !old_errno ) errno = 0; return res; } - void -pj_dalloc(void *ptr) { + +/**********************************************************************/ +void pj_dalloc(void *ptr) { +/**********************************************************************/ free(ptr); } @@ -70,3 +110,67 @@ pointer" to signal an error in a multi level allocation: pj_dalloc (ptr); return 0; } + + + + +/*****************************************************************************/ +void *pj_freeup_msg_plain (PJ *P, int errlev) { /* Destructor */ +/***************************************************************************** + Does memory deallocation for "plain" PJ objects, i.e. that vast majority + of PJs where the opaque object does not contain any additionally + allocated memory. +******************************************************************************/ + if (0==P) + return 0; + + if (0!=errlev) + pj_ctx_set_errno (P->ctx, errlev); + + if (0==P->opaque) + return pj_dealloc (P); + + pj_dealloc (P->opaque); + return pj_dealloc(P); +} + + +/*****************************************************************************/ +void pj_freeup_plain (PJ *P) { +/***************************************************************************** + Adapts pj_freeup_msg_plain to the format defined for the callback in + the PJ object. + + i.e. reduces most instances of projection deallocation code to: + + static void freeup (PJ *P) { + pj_freeup_plain (P); + return; + } + + rather than: + + static void *freeup_msg_add (PJ *P, int errlev) { + if (0==P) + return 0; + pj_ctx_set_errno (P->ctx, errlev); + + if (0==P->opaque) + return pj_dealloc (P); + + (* projection specific deallocation goes here *) + + pj_dealloc (P->opaque); + return pj_dealloc(P); + } + + (* Adapts pipeline_freeup to the format defined for the PJ object *) + static void freeup_msg_add (PJ *P) { + freeup_new_add (P, 0); + return; + } + + ******************************************************************************/ + pj_freeup_msg_plain (P, 0); + return; + } diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c new file mode 100644 index 0000000000..ccec8a64eb --- /dev/null +++ b/src/pj_obs_api.c @@ -0,0 +1,248 @@ +/****************************************************************************** + * Project: PROJ.4 + * Purpose: Implement a (currently minimalistic) proj API based primarily + * on the PJ_OBS generic geodetic data type. + * + * proj thread contexts have not seen widespread use, so one of the + * intentions with this new API is to make them less visible on the + * API surface: Contexts do not have a life by themselves, they are + * visible only through their associated PJs, and the number of + * functions supporting them is limited. + * + * Author: Thomas Knudsen, thokn@sdfe.dk, 2016-06-09/2016-11-06 + * + ****************************************************************************** + * Copyright (c) 2016, Thomas Knudsen/SDFE + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + *****************************************************************************/ +#define PJ_OBS_C +#include +#include +#include +#include + + +/* Used as return value in case of errors */ +const PJ_OBS pj_obs_error = { + /* Cannot use HUGE_VAL here: MSVC misimplements HUGE_VAL as something that is not compile time constant */ + {{DBL_MAX,DBL_MAX,DBL_MAX,DBL_MAX}}, + {{DBL_MAX,DBL_MAX,DBL_MAX}}, + 0, 0 +}; + +/* Used for zero-initializing new objects */ +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; + +/* Euclidean distance between two 2D coordinates stored in PJ_OBSs */ +double pj_obs_dist_2d (PJ_OBS a, PJ_OBS b) { + double *A = a.coo.v, *B = b.coo.v; + return hypot (A[0] - B[0], A[1] - B[1]); +} + +/* Euclidean distance between two 3D coordinates stored in PJ_OBSs */ +double pj_obs_dist_3d (PJ_OBS a, PJ_OBS b) { + double *A = a.coo.v, *B = b.coo.v; + return hypot (hypot (A[0] - B[0], A[1] - B[1]), A[2] - B[2]); +} + + +PJ_OBS pj_fwdobs (PJ_OBS obs, PJ *P) { + 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_error_set (P, EINVAL); + return pj_obs_error; +} + + +PJ_OBS pj_invobs (PJ_OBS obs, PJ *P) { + 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_error_set (P, EINVAL); + return pj_obs_error; +} + + +/* Apply the transformation P to the observation obs */ +PJ_OBS pj_apply (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_error_set (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) { + int i; + PJ_OBS o, u; + + if (0==P) + return HUGE_VAL; + + if (n < 1) { + pj_error_set (P, EINVAL); + return HUGE_VAL; + } + + o.coo = obs.coo; + + for (i = 0; i < n; i++) { + switch (direction) { + case PJ_FWD: + u.coo.xyz = pj_fwd3d (o.coo.lpz, P); + o.coo.lpz = pj_inv3d (u.coo.xyz, P); + break; + case PJ_INV: + u.coo.lpz = pj_inv3d (o.coo.xyz, P); + o.coo.xyz = pj_fwd3d (u.coo.lpz, P); + break; + default: + pj_error_set (P, EINVAL); + return HUGE_VAL; + } + } + + return pj_obs_dist_3d (o, obs); +} + + +PJ *pj_create (const char *definition) { + return pj_init_plus (definition); +} + + +PJ *pj_create_argv (int argc, char **argv) { + return pj_init (argc, argv); +} + + +/* From here: 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 */ + +int pj_error (PJ *P) { + return pj_ctx_get_errno (pj_get_ctx(P)); +} + + +void pj_error_set (PJ *P, int err) { + pj_ctx_set_errno (pj_get_ctx(P), err); +} + + +/* Set debug level 0-3. Higher number means more debug info. 0 turns it off */ +void pj_debug_set (PJ *P, enum pj_debug_level debuglevel) { + PJ_CONTEXT *ctx; + if (0==P) + ctx = pj_get_default_ctx(); + else + ctx = pj_get_ctx (P); + ctx->debug_level = debuglevel; +} + + +/* 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_set (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_error_set (P, ENOMEM); + return 1; + } + + pj_set_ctx (P, ctx); + return 0; +} + +/* Move daughter to mother's context */ +void pj_context_inherit (PJ *mother, PJ *daughter) { + pj_set_ctx (daughter, pj_get_ctx (mother)); +} + + +void pj_context_free (const PJ *P) { + PJ_CONTEXT *ctx; + + if (0==P) + return; + + /* 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 */ + return; + } + + ctx = pj_get_ctx ((PJ *) P); + + /* Otherwise, trying to free the default context is a no-op */ + 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; +} diff --git a/src/pj_observation.c b/src/pj_observation.c deleted file mode 100644 index 0d1f7acf57..0000000000 --- a/src/pj_observation.c +++ /dev/null @@ -1,114 +0,0 @@ -/****************************************************************************** - * Project: PROJ.4 - * Purpose: Implement some service routines for the PJ_OBSERVATION generic - * geodetic data type - * Author: Thomas Knudsen, thokn@sdfe.dk, 2016-06-09/2016-10-11 - * - ****************************************************************************** - * Copyright (c) 2016, Thomas Knudsen/SDFE - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the "Software"), - * to deal in the Software without restriction, including without limitation - * the rights to use, copy, modify, merge, publish, distribute, sublicense, - * and/or sell copies of the Software, and to permit persons to whom the - * Software is furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL - * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING - * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER - * DEALINGS IN THE SOFTWARE. - *****************************************************************************/ -#include -#include - - -/* Constructor for the OBSERVATION object */ -PJ_OBSERVATION pj_observation ( - double x, double y, double z, double t, - double o, double p, double k, - int id, unsigned int flags -) { - PJ_OBSERVATION g; - g.coo.xyzt.x = x; - g.coo.xyzt.y = y; - g.coo.xyzt.z = z; - g.coo.xyzt.t = t; - g.anc.opk.o = o; - g.anc.opk.p = p; - g.anc.opk.k = k; - g.id = id; - g.flags = flags; - return g; -} - - -PJ_OBSERVATION pj_apply (PJ *P, enum pj_direction direction, PJ_OBSERVATION obs) { - - switch (direction) { - case 1: - obs.coo.xyz = pj_fwd3d (obs.coo.lpz, P); - return obs; - case -1: - obs.coo.lpz = pj_inv3d (obs.coo.xyz, P); - return obs; - case 0: - return obs; - default: - pj_ctx_set_errno (P->ctx, EINVAL); - } - - return pj_observation (HUGE_VAL,HUGE_VAL,HUGE_VAL,HUGE_VAL,HUGE_VAL,HUGE_VAL,HUGE_VAL,0,0); -} - - -/* initial attempt, following a suggestion by Kristian Evers */ -double pj_roundtrip(PJ *P, enum pj_direction direction, int n, PJ_OBSERVATION obs) { - double d; - PJ_OBSERVATION o; - - /* multiple roundtrips not supported yet */ - if (n > 1) { - pj_ctx_set_errno (P->ctx, EINVAL); - return HUGE_VAL; - } - - switch (direction) { - case 1: - o.coo.xyz = pj_fwd3d (obs.coo.lpz, P); - break; - case -1: - o.coo.lpz = pj_inv3d (obs.coo.xyz, P); - break; - default: - pj_ctx_set_errno (P->ctx, EINVAL); - return HUGE_VAL; - } - - o.coo.xyz = pj_fwd3d (obs.coo.lpz, P); - o.coo.lpz = pj_inv3d (o.coo.xyz, P); - - /* Should add "both ways" here */ - d = hypot (hypot (o.coo.v[0] - obs.coo.v[0], o.coo.v[1] - obs.coo.v[1]), o.coo.v[2] - obs.coo.v[2]); - return d; -} - - -int pj_show_triplet (FILE *stream, const char *banner, PJ_TRIPLET point) { - int i = 0; - if (banner) - i += fprintf (stream, "%s", banner); - - i += fprintf(stream, "%16.10f ", point.xyz.x); - i += fprintf(stream, "%16.10f ", point.xyz.y); - i += fprintf(stream, "%16.4f", point.xyz.z); - if (banner) - i += fprintf(stream, "\n"); - return i; -} diff --git a/src/proj.def b/src/proj.def index dbf27f96e0..09b32c32ce 100644 --- a/src/proj.def +++ b/src/proj.def @@ -89,7 +89,17 @@ EXPORTS geod_polygon_testpoint @87 geod_polygon_clear @88 pj_run_selftests @89 - pj_observation @90 - pj_show_triplet @91 - pj_apply @92 - + pj_error @90 + pj_create @91 + pj_create_argv @92 + pj_apply @93 + pj_fwdobs @94 + pj_invobs @95 + pj_roundtrip @96 + pj_debug_set @97 + pj_error_set @98 + pj_log_set @99 + pj_context_renew @100 + pj_context_inherit @101 + pj_context_free @102 + pj_fileapi_set @103 diff --git a/src/proj.h b/src/proj.h index 3d8ba35904..4b7cfc8775 100644 --- a/src/proj.h +++ b/src/proj.h @@ -1,23 +1,88 @@ /****************************************************************************** * Project: PROJ.4 - * Purpose: Revised, experimental, "bare essentials" API for PROJ.4. - * Intended as the foundation for added geodetic functionality. + * Purpose: Revised, experimental API for PROJ.4, intended as the foundation + * for added geodetic functionality. * - * Introduces the OBSERVATION data type, for generic coordinate - * and ancillary data handling. + * The original proj API (defined in projects.h) has grown organically + * over the years, but it has also grown somewhat messy. * - * Also introduces the PJ_SPATIOTEMPORAL and PJ_TRIPLET unions - * making it possible to make explicit the previously used - * "implicit type punning", where a XY is turned into a LP by - * re#defining both as UV, behind the back of the user. + * The same has happened with the newer high level API (defined in + * proj_api.h): To support various historical objectives, proj_api.h + * contains a rather complex combination of conditional defines and + * typedefs. Probably for good (historical) reasons, which are not + * always evident from today's perspective. + * + * 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, + * while all data types are being moved to the PJ_ namespace. + * + * Please note that this API is *orthogonal* to the previous APIs: + * Apart from some inclusion guards, projects.h and proj_api.h are not + * touched - if you do not include proj.h, the projects and proj_api + * APIs should work as they always have. + * + * A few implementation details: + * + * Apart from the namespacing efforts, I'm trying to eliminate three + * proj_api elements, which I have found especially confusing. + * + * FIRST and foremost, I try to avoid typedef'ing away pointer + * semantics. I agree that it can be occasionally useful, but I + * prefer having the pointer nature of function arguments being + * explicitly visible. + * + * Hence, projCtx has been replaced by PJ_CONTEXT *. + * and projPJ has been replaced by PJ * + * + * SECOND, I try to eliminate cases of information hiding implemented + * by redefining data types to void pointers. + * + * I prefer using a combination of forward declarations and typedefs. + * Hence: + * typedef void *projCtx; + * Has been replaced by: + * struct projCtx_t; + * typedef struct projCtx_t PJ_CONTEXT; + * This makes it possible for the calling program to know that the + * PJ_CONTEXT data type exists, and handle pointers to that data type + * without having any idea about its internals. + * + * (obviously, in this example, struct projCtx_t should also be + * renamed struct pj_ctx some day...) + * + * THIRD, I try to eliminate implicit type punning. Hence this API + * introduces the PJ_OBS ("observation") data type, for generic + * coordinate and handling of ancillary data. + * + * It includes the PJ_COORD and PJ_TRIPLET unions making it possible + * to make explicit the previously used "implicit type punning", where + * a XY is turned into a LP by re#defining both as UV, behind the back + * of the user. + * + * The PJ_COORD union is used for storing 1D, 2D, 3D and 4D coordinates. + * The PJ_TRIPLET union is used for storing any set of up to 3 related + * observations. At the application code level, the names of these + * unions will usually not be used - they will only be accessed via + * their tag names in the PJ_OBS data type. * * The bare essentials API presented here follows the PROJ.4 * convention of sailing the coordinate to be reprojected, up on * the stack ("call by value"), and symmetrically returning the - * result on the stack. Although the OBSERVATION object is 4 times + * result on the stack. Although the PJ_OBS object is 4 times * as large as the traditional XY and LP objects, timing results * have shown the overhead to be very reasonable. * + * In its current incarnation, the API is focused on compactness: + * Currently it consists of only nine functions. + * + * 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. + * + * See pj_proj_test.c for an example of how to use the API. + * * Author: Thomas Knudsen, * ****************************************************************************** @@ -35,7 +100,7 @@ * * 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 + * 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 @@ -46,11 +111,21 @@ #endif #include +#include +#include + #include #include #include #include +#ifdef PROJECTS_H +#error proj.h must be included before projects.h +#endif +#ifdef PROJ_API_H +#error proj.h must be included before proj_api.h +#endif + #ifndef PROJ_H #define PROJ_H #ifdef __cplusplus @@ -58,33 +133,26 @@ extern "C" { #endif -/****************************************************************************** - proj.h is included by projects.h in order to determine the size of the - PJ_OBSERVATION object. +/* 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 - In order to stomp as little as possible on the traditional proj.4 name - space, proj.h is littered with inclusion guards, leaving only the minimum - possible implementation when included from projects.h - The PJ_OBSERVATION object is fully defined if proj.h is included alone or - in connection with *but before* projects.h (the latter may be needed in - some cases, where it is necessary to access this "bare essentials" API, - while still having direct access to PJ object internals) -******************************************************************************/ +/* first forward declare everything needed */ /* Data type for generic geodetic observations */ -struct PJ_OBSERVATION; -typedef struct PJ_OBSERVATION PJ_OBSERVATION; +struct PJ_OBS; +typedef struct PJ_OBS PJ_OBS; /* Data type for generic geodetic 3D data */ union PJ_TRIPLET; typedef union PJ_TRIPLET PJ_TRIPLET; /* Data type for generic geodetic 3D data plus epoch information */ -union PJ_SPATIOTEMPORAL; -typedef union PJ_SPATIOTEMPORAL PJ_SPATIOTEMPORAL; +union PJ_COORD; +typedef union PJ_COORD PJ_COORD; -#ifndef PROJECTS_H /* Data type for projection/transformation information */ struct PJconsts; typedef struct PJconsts PJ; /* the PJ object herself */ @@ -115,21 +183,18 @@ typedef struct { double lam, phi, z; } LPZ; /* Degrees, minutes, and seconds */ typedef struct { double d, m, s; } PJ_DMS; -/* Geoid undulation (N) and deflections of the vertical (z, e) */ -typedef struct { double N, z, e; } PJ_NZE; +/* Geoid undulation (N) and deflections of the vertical (eta, zeta) */ +typedef struct { double e, z, N; } PJ_EZN; /* Ellipsoidal parameters */ typedef struct { double a, f; } PJ_AF; -#endif /* ndef PROJECTS_H */ -union PJ_SPATIOTEMPORAL { -#ifndef PROJECTS_H +union PJ_COORD { PJ_XYZT xyzt; PJ_UVWT uvwt; PJ_ENHT enht; PJ_LPZT lpzt; PJ_ENH enh; -#endif double v[4]; /* Who cares - it's just a vector! */ XYZ xyz; UVW uvw; @@ -143,12 +208,10 @@ union PJ_SPATIOTEMPORAL { /* Avoid preprocessor renaming and implicit type-punning: Use a union to make it explicit */ union PJ_TRIPLET { -#ifndef PROJECTS_H PJ_OPK opk; PJ_ENH enh; - PJ_NZE nze; + PJ_EZN ezn; PJ_AF af; -#endif double v[3]; /* Who cares - it's just a vector! */ XYZ xyz; LPZ lpz; @@ -158,39 +221,68 @@ union PJ_TRIPLET { UV uv; }; -struct PJ_OBSERVATION { - PJ_SPATIOTEMPORAL coo; /* coordinate data */ +struct PJ_OBS { + PJ_COORD coo; /* coordinate data */ PJ_TRIPLET anc; /* ancillary data */ int id; /* integer ancillary data - e.g. observation number, EPSG code... */ unsigned int flags; /* additional data, intended for flags */ }; +/* The context type - properly namespaced synonym for projCtx */ +struct projCtx_t; +typedef struct projCtx_t PJ_CONTEXT; +typedef int *PJ_FILE; -#ifndef PROJECTS_H +/* 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); -/* Direction: "+" forward, "-" reverse, 0: do nothing */ -enum pj_direction { - PJ_FWD = -1, - PJ_IDENT = 0, - PJ_INV = 1 +/* High level functionality for handling thread contexts */ +enum pj_debug_level { + PJ_LOG_NONE = 0, + PJ_LOG_ERROR = 1, + PJ_LOG_DEBUG_MAJOR = 2, + PJ_LOG_DEBUG_MINOR = 3 }; +void pj_debug_set (PJ *P, enum pj_debug_level debuglevel); +void pj_error_set (PJ *P, int err); +void pj_log_set (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); /* Apply transformation to observation - in forward or inverse direction */ -PJ_OBSERVATION pj_apply (PJ *P, enum pj_direction direction, PJ_OBSERVATION obs); +enum pj_direction { + PJ_FWD = 1, /* Forward */ + PJ_IDENT = 0, /* Do nothing */ + PJ_INV = -1 /* Inverse */ +}; +PJ_OBS pj_apply (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_OBSERVATION obs); +double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs); + +/* Euclidean distance between two 2D coordinates stored in PJ_OBSs */ +double pj_obs_dist_2d (PJ_OBS a, PJ_OBS b); +/* Euclidean distance between two 3D coordinates stored in PJ_OBSs */ +double pj_obs_dist_3d (PJ_OBS a, PJ_OBS b); -int pj_show_triplet (FILE *stream, const char *banner, PJ_TRIPLET point); -/* Constructor for the OBSERVATION object */ -PJ_OBSERVATION pj_observation ( - double x, double y, double z, double t, - double o, double p, double k, - int id, unsigned int flags -); + +#ifndef PJ_OBS_C +extern const PJ_OBS pj_obs_error; +extern const PJ_OBS pj_obs_null; +extern const PJ *pj_shutdown; +#endif #ifndef TODEG #define TODEG(rad) ((rad)*180.0/M_PI) @@ -199,10 +291,9 @@ PJ_OBSERVATION pj_observation ( #define TORAD(deg) ((deg)*M_PI/180.0) #endif -#endif /* ndef PROJECTS_H */ - #ifdef __cplusplus } #endif + #endif /* ndef PROJ_H */ diff --git a/src/proj_api.h b/src/proj_api.h index 24a6f05395..706e35ab11 100644 --- a/src/proj_api.h +++ b/src/proj_api.h @@ -25,45 +25,77 @@ * DEALINGS IN THE SOFTWARE. *****************************************************************************/ + + + /* + * This version number 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. + */ + #ifndef PJ_VERSION + #define PJ_VERSION 493 + #endif + + +/* If we're not asked for PJ_VERSION only, give them everything */ +#ifndef PROJ_API_INCLUDED_FOR_PJ_VERSION_ONLY /* General projections header file */ #ifndef PROJ_API_H #define PROJ_API_H -/* standard inclusions */ -#include -#include - #ifdef __cplusplus extern "C" { #endif -/* - * This version number 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 PJ_VERSION 493 + +/* standard inclusions */ +#include +#include + /* pj_init() and similar functions can be used with a non-C locale */ /* Can be detected too at runtime if the symbol pj_atof exists */ #define PJ_LOCALE_SAFE 1 extern char const pj_release[]; /* global release id string */ +extern int pj_errno; /* global error return code */ #define RAD_TO_DEG 57.295779513082321 #define DEG_TO_RAD .017453292519943296 -extern int pj_errno; /* global error return code */ +#if defined(PROJECTS_H) || defined(PROJ_H) +#define PROJ_API_H_NOT_INVOKED_AS_PRIMARY_API +#endif + + + +#ifndef PROJ_H +/* In proj.h these macros are replaced by the enumeration pj_debug_level */ +#define PJ_LOG_NONE 0 +#define PJ_LOG_ERROR 1 +#define PJ_LOG_DEBUG_MAJOR 2 +#define PJ_LOG_DEBUG_MINOR 3 +#endif + -#if !defined(PROJECTS_H) +#ifdef PROJ_API_H_NOT_INVOKED_AS_PRIMARY_API + /* These make the function declarations below conform with classic proj */ + typedef PJ *projPJ; /* projPJ is a pointer to PJ */ + typedef projCtx_t *projCtx; /* projCtx is a pointer to projCtx_t */ +# define projXY XY +# define projLP LP +# define projXYZ XYZ +# define projLPZ LPZ +#else + /* i.e. proj_api invoked as primary API */ typedef struct { double u, v; } projUV; typedef struct { double u, v, w; } projUVW; typedef void *projPJ; @@ -72,15 +104,12 @@ extern int pj_errno; /* global error return code */ #define projXYZ projUVW #define projLPZ projUVW typedef void *projCtx; -#else - typedef PJ *projPJ; - typedef projCtx_t *projCtx; -# define projXY XY -# define projLP LP -# define projXYZ XYZ -# define projLPZ LPZ #endif + + + +/* If included *after* proj.h finishes, we have alternative names */ /* file reading api, like stdio */ typedef int *PAFile; typedef struct projFileAPI_t { @@ -91,14 +120,20 @@ typedef struct projFileAPI_t { void (*FClose)(PAFile); } projFileAPI; + + /* procedure prototypes */ +projCtx pj_get_default_ctx(void); +projCtx pj_get_ctx( projPJ ); + projXY pj_fwd(projLP, projPJ); projLP pj_inv(projXY, projPJ); projXYZ pj_fwd3d(projLPZ, projPJ); projLPZ pj_inv3d(projXYZ, projPJ); + int pj_transform( projPJ src, projPJ dst, long point_count, int point_offset, double *x, double *y, double *z ); int pj_datum_transform( projPJ src, projPJ dst, long point_count, int point_offset, @@ -128,6 +163,8 @@ projPJ pj_init_ctx( projCtx, int, char ** ); projPJ pj_init_plus_ctx( projCtx, const char * ); char *pj_get_def(projPJ, int); projPJ pj_latlong_from_proj( projPJ ); + + void *pj_malloc(size_t); void pj_dalloc(void *); void *pj_calloc (size_t n, size_t size); @@ -139,8 +176,8 @@ void pj_acquire_lock(void); void pj_release_lock(void); void pj_cleanup_lock(void); -projCtx pj_get_default_ctx(void); -projCtx pj_get_ctx( projPJ ); +int pj_run_selftests (int verbosity); + void pj_set_ctx( projPJ, projCtx ); projCtx pj_ctx_alloc(void); void pj_ctx_free( projCtx ); @@ -168,17 +205,10 @@ char *pj_ctx_fgets(projCtx ctx, char *line, int size, PAFile file); PAFile pj_open_lib(projCtx, const char *, const char *); -int pj_run_selftests (int verbosity); - - -#define PJ_LOG_NONE 0 -#define PJ_LOG_ERROR 1 -#define PJ_LOG_DEBUG_MAJOR 2 -#define PJ_LOG_DEBUG_MINOR 3 #ifdef __cplusplus } #endif #endif /* ndef PROJ_API_H */ - +#endif /* ndef PROJ_API_INCLUDED_FOR_PJ_VERSION_ONLY */ diff --git a/src/projects.h b/src/projects.h index f6ac6acfe1..d07c43c7c3 100644 --- a/src/projects.h +++ b/src/projects.h @@ -168,6 +168,7 @@ typedef struct { double u, v, w; } UVW; /* Forward declarations and typedefs for stuff needed inside the PJ object */ struct PJconsts; +struct PJ_OBS; struct pj_opaque; struct ARG_list; struct FACTORS; @@ -181,6 +182,7 @@ enum pj_io_units { }; #ifndef PROJ_H typedef struct PJconsts PJ; /* the PJ object herself */ +typedef struct PJ_OBS PJ_OBS; #endif struct PJ_REGION_S { @@ -191,11 +193,9 @@ struct PJ_REGION_S { }; -#include "proj.h" /* Need this for sizeof(PJ_OBSERVATION) */ struct projCtx_t; typedef struct projCtx_t projCtx_t; - /* base projection data structure */ struct PJconsts { @@ -230,12 +230,13 @@ struct PJconsts { **************************************************************************************/ + XY (*fwd)(LP, PJ *); LP (*inv)(XY, PJ *); XYZ (*fwd3d)(LPZ, PJ *); LPZ (*inv3d)(XYZ, PJ *); - PJ_OBSERVATION (*fwdobs)(PJ_OBSERVATION, PJ *); - PJ_OBSERVATION (*invobs)(PJ_OBSERVATION, PJ *); + PJ_OBS (*fwdobs)(PJ_OBS, PJ *); + PJ_OBS (*invobs)(PJ_OBS, PJ *); void (*spc)(LP, PJ *, struct FACTORS *); @@ -280,8 +281,8 @@ struct PJconsts { double one_es; /* 1 - e^2 */ double rone_es; /* 1/one_es */ - /* The flattenings */ + /* The flattenings */ double f; /* first flattening */ double f2; /* second flattening */ double n; /* third flattening */ @@ -589,7 +590,6 @@ typedef struct _PJ_GridCatalog { struct _PJ_GridCatalog *next; } PJ_GridCatalog; - /* procedure prototypes */ double dmstor(const char *, char **); double dmstor_ctx(projCtx ctx, const char *, char **); @@ -722,9 +722,13 @@ struct PJ_PRIME_MERIDIANS *pj_get_prime_meridians_ref( void ); double pj_atof( const char* nptr ); double pj_strtod( const char *nptr, char **endptr ); +void pj_freeup_plain (PJ *P); #ifdef __cplusplus } #endif +#ifndef PROJECTS_H_ATEND +#define PROJECTS_H_ATEND +#endif #endif /* end of basic projections header */