aboutsummaryrefslogtreecommitdiff
path: root/src/pj_obs_api.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/pj_obs_api.c')
-rw-r--r--src/pj_obs_api.c485
1 files changed, 325 insertions, 160 deletions
diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c
index 37e55bc3..d8395b0f 100644
--- a/src/pj_obs_api.c
+++ b/src/pj_obs_api.c
@@ -34,118 +34,45 @@
*****************************************************************************/
#define PJ_OBS_API_C
#include <proj.h>
-#include <projects.h>
+#include "proj_internal.h"
+#include "projects.h"
#include <geodesic.h>
-#include <float.h>
-#include <math.h>
+#include <stddef.h>
+#include <errno.h>
-/* 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);
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 23:57:54 +0000