aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Makefile.am2
-rw-r--r--src/PJ_cart.c65
-rw-r--r--src/PJ_pipeline.c603
-rw-r--r--src/lib_proj.cmake1
-rw-r--r--src/makefile.vc2
-rw-r--r--src/pj_ell_set.c188
-rw-r--r--src/pj_init.c41
-rw-r--r--src/pj_list.h1
-rw-r--r--src/pj_log.c65
-rw-r--r--src/pj_obs_api.c46
-rw-r--r--src/pj_strerrno.c126
-rw-r--r--src/proj.def10
-rw-r--r--src/proj.h95
-rw-r--r--src/proj_api.h16
-rw-r--r--src/projects.h3
15 files changed, 1014 insertions, 250 deletions
diff --git a/src/Makefile.am b/src/Makefile.am
index 1d53375a..69e0e44e 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_obs_api.c PJ_cart.c
+ pj_obs_api.c PJ_cart.c PJ_pipeline.c
install-exec-local:
rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT)
diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index 5df95c31..4ec20ced 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -197,8 +197,9 @@ static LPZ geodetic (XYZ cartesian, PJ *P) {
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 */
+ /* by computing the height as the cartesian z value */
+ /* minus the semiminor axis length */
+ /* (potential improvement: approx. by 2nd order poly) */
lpz.z = cartesian.z - (cartesian.z > 0? b: -b);
}
else
@@ -247,21 +248,31 @@ 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 *P;
PJ_OBS a, b;
int err;
double dist;
+ char *args[3] = {"proj=utm", "zone=32", "ellps=GRS80"};
+
/* Log everything libproj offers to log for you */
- pj_debug_set (0, PJ_LOG_DEBUG_MINOR);
+ 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)
+ P = pj_create ("+proj=utm +zone=32 +ellps=GRS80");
+ if (0==P)
return 1;
+ /* 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;
+
/* Turn off logging */
- pj_debug_set (0, PJ_LOG_NONE);
+ pj_log_level (0, PJ_LOG_NONE);
/* zero initialize everything, then set (longitude, latitude) to (12, 55) */
a = pj_obs_null;
@@ -270,38 +281,38 @@ int pj_cart_selftest (void) {
a.coo.lp.phi = TORAD(55);
/* Forward projection */
- b = pj_apply (p, PJ_FWD, a);
+ b = pj_trans (P, PJ_FWD, a);
/* Inverse projection */
- a = pj_apply (p, PJ_INV, b);
+ a = pj_trans (P, PJ_INV, b);
/* Null projection */
- a = pj_apply (p, PJ_IDENT, a);
+ a = pj_trans (P, PJ_IDENT, a);
/* Forward again, to get two linear items for comparison */
- a = pj_apply (p, PJ_FWD, a);
+ a = pj_trans (P, PJ_FWD, a);
- dist = pj_obs_dist_2d (a, b);
+ dist = pj_xy_dist (a.coo.xy, b.coo.xy);
if (dist > 2e-9)
return 2;
/* Invalid projection */
- a = pj_apply (p, 42, a);
+ a = pj_trans (P, 42, a);
if (a.coo.lpz.lam!=DBL_MAX)
return 3;
- err = pj_error (p);
+ err = pj_error (P);
if (0==err)
return 4;
/* Clear error */
- pj_error_set (p, 0);
+ pj_error_set (P, 0);
/* Clean up */
- pj_free (p);
+ pj_free (P);
/* Now do some 3D transformations */
- p = pj_create ("+proj=cart +ellps=GRS80");
- if (0==p)
+ P = pj_create ("+proj=cart +ellps=GRS80");
+ if (0==P)
return 5;
/* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
@@ -311,26 +322,26 @@ int pj_cart_selftest (void) {
a.coo.lpz.z = 100;
/* Forward projection: 3D-Cartesian-to-Ellipsoidal */
- b = pj_apply (p, PJ_FWD, a);
+ b = pj_trans (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 = 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);
+ b = pj_trans (P, PJ_INV, b);
/* Move p to another context */
- pj_context_renew (p);
- b = pj_apply (p, PJ_FWD, b);
+ pj_context_renew (P);
+ b = pj_trans (P, PJ_FWD, b);
/* Move it back to the default context */
- pj_context_free (p);
- b = pj_apply (p, PJ_INV, b);
+ pj_context_free (P);
+ b = pj_trans (P, PJ_INV, b);
- pj_free (p);
+ pj_free (P);
return 0;
}
#endif
diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c
new file mode 100644
index 00000000..ef86ccd3
--- /dev/null
+++ b/src/PJ_pipeline.c
@@ -0,0 +1,603 @@
+/*******************************************************************************
+
+ Transformation pipeline manager
+
+ Thomas Knudsen, 2016-05-20/2016-11-20
+
+********************************************************************************
+
+ Geodetic transformations are typically organized in a number of
+ steps. For example, a datum shift could be carried out through
+ these steps:
+
+ 1. Convert (latitude, longitude, ellipsoidal height) to
+ 3D geocentric cartesian coordinates (X, Y, Z)
+ 2. Transform the (X, Y, Z) coordinates to the new datum, using a
+ 7 parameter Helmert transformation.
+ 3. Convert (X, Y, Z) back to (latitude, longitude, ellipsoidal height)
+
+ If the height system used is orthometric, rather than ellipsoidal,
+ another step is needed at each end of the process:
+
+ 1. Add the local geoid undulation (N) to the orthometric height
+ to obtain the ellipsoidal (i.e. geometric) height.
+ 2. Convert (latitude, longitude, ellipsoidal height) to
+ 3D geocentric cartesian coordinates (X, Y, Z)
+ 3. Transform the (X, Y, Z) coordinates to the new datum, using a
+ 7 parameter Helmert transformation.
+ 4. Convert (X, Y, Z) back to (latitude, longitude, ellipsoidal height)
+ 5. Subtract the local geoid undulation (N) from the ellipsoidal height
+ to obtain the orthometric height.
+
+ Additional steps can be added for e.g. change of vertical datum, so the
+ list can grow fairly long. None of the steps are, however, particularly
+ complex, and data flow is strictly from top to bottom.
+
+ Hence, in principle, the first example above could be implemented using
+ Unix pipelines:
+
+ cat my_coordinates | geographic_to_xyz | helmert | xyz_to_geographic > my_transformed_coordinates
+
+ in the grand tradition of Software Tools [1].
+
+ The proj pipeline driver implements a similar concept: Stringing together
+ a number of steps, feeding the output of one step to the input of the next.
+
+ It is a very powerful concept, that increases the range of relevance of the
+ proj.4 system substantially. It is, however, not a particularly intrusive
+ addition to the code base: The implementation is by and large completed by
+ adding an extra projection called "pipeline" (i.e. this file), which
+ handles all business.
+
+ Syntactically, the pipeline system introduces the "+step" keyword (which
+ indicates the start of each transformation step), the "+omit_fwd" and
+ "+omit_inv" keywords (which indicate that a given transformation step
+ should be omitted when the pipeline is executed in forward, resp. inverse
+ direction), and reintroduces the +inv keyword (indicating that a given
+ transformation step should run in reverse, i.e. forward, when the pipeline
+ is executed in inverse direction, and vice versa).
+
+ Hence, the first transformation example above, can be implemented as:
+
+ +proj=pipeline +step proj=cart +step proj=helmert <ARGS> +step proj=cart +inv
+
+ Where <ARGS> indicate the Helmert arguments: 3 translations (+x=..., +y=...,
+ +z=...), 3 rotations (+rx=..., +ry=..., +rz=...) and a scale factor (+s=...).
+ Following geodetic conventions, the rotations are given in Milliarcseconds,
+ and the scale factor is given as parts-per-million.
+
+ [1] B. W. Kernighan & P. J. Plauger: Software tools. Addison-Wesley, 1976, 338 pp.
+
+********************************************************************************
+
+Thomas Knudsen, thokn@sdfe.dk, 2016-05-20
+
+********************************************************************************
+* 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 <proj.h>
+#include <projects.h>
+
+#include <assert.h>
+#include <stddef.h>
+#include <errno.h>
+PROJ_HEAD(pipeline, "Transformation pipeline manager");
+
+/* Projection specific elements for the PJ object */
+#define PIPELINE_STACK_SIZE 100
+struct pj_opaque {
+ int reversible;
+ int steps;
+ int depth;
+ int verbose;
+ int *reverse_step;
+ int *omit_forward;
+ int *omit_inverse;
+ PJ_OBS stack[PIPELINE_STACK_SIZE];
+ PJ **pipeline;
+};
+
+
+static PJ_OBS pipeline_forward_obs (PJ_OBS, PJ *P);
+static PJ_OBS pipeline_reverse_obs (PJ_OBS, PJ *P);
+static XYZ pipeline_forward_3d (LPZ lpz, PJ *P);
+static LPZ pipeline_reverse_3d (XYZ xyz, PJ *P);
+static XY pipeline_forward (LP lpz, PJ *P);
+static LP pipeline_reverse (XY xyz, PJ *P);
+
+
+
+/********************************************************************
+
+ ISOMORPHIC TRANSFORMATIONS
+
+*********************************************************************
+
+ In this context, an isomorphic transformation is a proj PJ
+ object returning the same kind of coordinates that it
+ receives, i.e. a transformation from angular to angular or
+ linear to linear coordinates.
+
+ The degrees-to-radians operation is an example of the former,
+ while the latter is typical for most of the datum shift
+ operations used in geodesy, e.g. the Helmert 7-parameter shift.
+
+ Isomorphic transformations trips up the pj_inv/pj_fwd
+ functions, which try to check input sanity and scale output to
+ internal proj units, under the assumption that if input is of
+ angular type, then output must be of linear (or vice versa).
+
+ Hence, to avoid having pj_inv/pj_fwd stomping on output (or
+ choking on input), we need a way to tell them that they should
+ accept whatever is being handed to them.
+
+ The P->left and P->right elements indicate isomorphism.
+
+ For classic proj style projections, they both have the
+ value PJ_IO_UNITS_CLASSIC (default initialization), indicating
+ that the forward driver expects angular input coordinates, and
+ provides linear output coordinates.
+
+ Newer projections may set P->left and P->right to either
+ PJ_IO_UNITS_METERS, PJ_IO_UNITS_RADIANS or PJ_IO_UNITS_ANY,
+ to indicate their I/O style.
+
+ For the forward driver, left indicates input coordinate
+ type, while right indicates output coordinate type.
+
+ For the inverse driver, left indicates output coordinate
+ type, while right indicates input coordinate type.
+
+*********************************************************************/
+
+
+
+
+static PJ_OBS pipeline_forward_obs (PJ_OBS point, PJ *P) {
+ int i, first_step, last_step, incr;
+
+ first_step = 1;
+ last_step = P->opaque->steps + 1;
+ incr = 1;
+
+ for (i = first_step; i != last_step; i += incr) {
+ if (P->opaque->omit_forward[i])
+ continue;
+ if (P->opaque->reverse_step[i])
+ point = pj_trans (P->opaque->pipeline[i], PJ_INV, point);
+ else
+ point = pj_trans (P->opaque->pipeline[i], PJ_FWD, point);
+ if (P->opaque->depth < PIPELINE_STACK_SIZE)
+ P->opaque->stack[P->opaque->depth++] = point;
+ }
+
+ P->opaque->depth = 0; /* Clear the stack */
+ return point;
+}
+
+static PJ_OBS pipeline_reverse_obs (PJ_OBS point, PJ *P) {
+ int i, first_step, last_step, incr;
+
+ first_step = P->opaque->steps;
+ last_step = 0;
+ incr = -1;
+
+ for (i = first_step; i != last_step; i += incr) {
+ if (P->opaque->omit_inverse[i])
+ continue;
+ if (P->opaque->reverse_step[i])
+ point = pj_trans (P->opaque->pipeline[i], PJ_FWD, point);
+ else
+ point = pj_trans (P->opaque->pipeline[i], PJ_INV, point);
+ if (P->opaque->depth < PIPELINE_STACK_SIZE)
+ P->opaque->stack[P->opaque->depth++] = point;
+ }
+
+ P->opaque->depth = 0; /* Clear the stack */
+ return point;
+}
+
+
+/* Delegate the work to pipeline_forward_obs() */
+static XYZ pipeline_forward_3d (LPZ lpz, PJ *P) {
+ PJ_OBS point = pj_obs_null;
+ point.coo.lpz = lpz;
+ point = pipeline_forward_obs (point, P);
+ return point.coo.xyz;
+}
+
+/* Delegate the work to pipeline_reverse_obs() */
+static LPZ pipeline_reverse_3d (XYZ xyz, PJ *P) {
+ PJ_OBS point = pj_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;
+ 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;
+ point.coo.xy = xy;
+ point = pipeline_reverse_obs (point, P);
+ return point.coo.lp;
+}
+
+static void freeup(PJ *P) { /* Destructor */
+ if (P==0)
+ return;
+ /* Projection specific deallocation goes here */
+ pj_dealloc (P->opaque);
+ pj_dealloc (P);
+ return;
+}
+
+
+static void *pipeline_freeup (PJ *P, int errlev) { /* Destructor */
+ int i;
+ if (0==P)
+ return 0;
+
+ pj_error_set (P, errlev);
+
+ if (0==P->opaque)
+ return pj_dealloc (P);
+
+ for (i = 0; i < P->opaque->steps; i++)
+ pj_free (P->opaque->pipeline[i+1]);
+
+ pj_dealloc (P->opaque->reverse_step);
+ pj_dealloc (P->opaque->omit_forward);
+ pj_dealloc (P->opaque->omit_inverse);
+ pj_dealloc (P->opaque->pipeline);
+
+ pj_dealloc (P->opaque);
+ return pj_dealloc(P);
+}
+
+
+/* Adapts pipeline_freeup to the format defined for the PJ object */
+static void pipeline_freeup_wrapper (PJ *P) {
+ pipeline_freeup (P, 0);
+ return;
+}
+
+
+static PJ *pj_create_pipeline (PJ *P, size_t steps) {
+
+ /* Room for the pipeline: An array of PJ * with room for sentinels at both ends */
+ P->opaque->pipeline = pj_calloc (steps + 2, sizeof(PJ *));
+ if (0==P->opaque->pipeline)
+ return 0;
+
+ P->opaque->steps = steps;
+
+ P->opaque->reverse_step = pj_calloc (steps + 2, sizeof(int));
+ if (0==P->opaque->reverse_step)
+ return 0;
+
+ P->opaque->omit_forward = pj_calloc (steps + 2, sizeof(int));
+ if (0==P->opaque->omit_forward)
+ return 0;
+
+ P->opaque->omit_inverse = pj_calloc (steps + 2, sizeof(int));
+ if (0==P->opaque->omit_inverse)
+ return 0;
+
+ return P;
+}
+
+/* count the number of args in pipeline definition */
+size_t argc_params (paralist *params) {
+ size_t argc = 0;
+ for (; params != 0; params = params->next)
+ argc++;
+ return ++argc; /* one extra for the sentinel */
+}
+
+/* Sentinel for argument list */
+static char argv_sentinel[5] = "step";
+
+/* turn paralist int argc/argv style argument list */
+char **argv_params (paralist *params) {
+ char **argv;
+ size_t argc = 0;
+ argv = pj_calloc (argc_params (params), sizeof (char *));
+ if (0==argv)
+ return 0;
+ for (; params != 0; params = params->next)
+ argv[argc++] = params->param;
+ argv[argc++] = argv_sentinel;
+ return argv;
+}
+
+
+PJ *PROJECTION(pipeline) {
+ int i, nsteps = 0, argc;
+ int i_pipeline = -1, i_first_step = -1, i_current_step;
+ char **argv, **current_argv;
+
+ P->fwdobs = pipeline_forward_obs;
+ P->invobs = pipeline_reverse_obs;
+ P->fwd3d = pipeline_forward_3d;
+ P->inv3d = pipeline_reverse_3d;
+ P->fwd = pipeline_forward;
+ P->inv = pipeline_reverse;
+ P->pfree = pipeline_freeup_wrapper;
+
+ P->opaque = pj_calloc (1, sizeof(struct pj_opaque));
+ if (0==P->opaque)
+ return 0;
+
+ argc = argc_params (P->params);
+ argv = argv_params (P->params);
+ if (0==argv)
+ return pipeline_freeup (P, ENOMEM);
+
+ /* The elements of current_argv are not used - we just use argv_params */
+ /* as allocator for a "large enough" container needed later */
+ current_argv = argv_params (P->params);
+ if (0==current_argv)
+ return pipeline_freeup (P, ENOMEM);
+
+ /* Do some syntactic sanity checking */
+ 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");
+ return pipeline_freeup (P, -50);
+ }
+ if (0==nsteps)
+ i_first_step = i;
+ nsteps++;
+ continue;
+ }
+
+ if (0==strcmp ("proj=pipeline", argv[i])) {
+ if (-1 != i_pipeline) {
+ pj_log_error (P, "Pipeline: Nesting invalid");
+ return pipeline_freeup (P, -50); /* ERROR: nested pipelines */
+ }
+ i_pipeline = i;
+ }
+ }
+ nsteps--; /* Last instance of +step is just a sentinel */
+ P->opaque->steps = nsteps;
+
+ if (-1==i_pipeline)
+ return pipeline_freeup (P, -50); /* ERROR: no pipeline def */
+
+ if (0==nsteps)
+ return pipeline_freeup (P, -50); /* ERROR: no pipeline def */
+
+ /* Make room for the pipeline and execution indicators */
+ if (0==pj_create_pipeline (P, nsteps))
+ return pipeline_freeup (P, ENOMEM);
+
+ /* Now loop over all steps, building a new set of arguments for each init */
+ for (i_current_step = i_first_step, i = 0; i < nsteps; i++) {
+ int j;
+ int current_argc = 0;
+ 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);
+
+ /* First add the step specific args */
+ for (j = i_current_step + 1; 0 != strcmp ("step", argv[j]); j++)
+ current_argv[current_argc++] = argv[j];
+
+ i_current_step = j;
+
+ /* Then add the global args */
+ for (j = i_pipeline + 1; 0 != strcmp ("step", argv[j]); j++)
+ current_argv[current_argc++] = argv[j];
+
+ /* Finally handle non-symmetric steps and inverted steps */
+ for (j = 0; j < current_argc; j++) {
+ if (0==strcmp("omit_inv", current_argv[j])) {
+ P->opaque->omit_inverse[i+1] = 1;
+ P->opaque->omit_forward[i+1] = 0;
+ }
+ if (0==strcmp("omit_fwd", current_argv[j])) {
+ P->opaque->omit_inverse[i+1] = 0;
+ P->opaque->omit_forward[i+1] = 1;
+ }
+ if (0==strcmp("inv", current_argv[j]))
+ P->opaque->reverse_step[i+1] = 1;
+ }
+
+ pj_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]);
+
+ next_step = pj_init (current_argc, current_argv);
+ pj_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]);
+ return pipeline_freeup (P, -50); /* ERROR: bad pipeline def */
+ }
+ P->opaque->pipeline[i+1] = next_step;
+ pj_log_trace (P, "Pipeline: step done");
+ }
+
+ pj_log_trace (P, "Pipeline: %d steps built. Determining i/o characteristics", nsteps);
+
+ /* Determine forward input (= reverse output) data type */
+
+ /* First locate the first forward-active pipeline step */
+ for (i = 0; i < nsteps; i++)
+ if (0==P->opaque->omit_forward[i+1])
+ break;
+ if (i==nsteps) {
+ pj_log_error (P, "Pipeline: No forward steps");
+ return pipeline_freeup (P, -50);
+ }
+
+ if (P->opaque->reverse_step[i + 1])
+ P->left = P->opaque->pipeline[i + 1]->right;
+ else
+ P->left = P->opaque->pipeline[i + 1]->left;
+
+ if (P->left==PJ_IO_UNITS_CLASSIC) {
+ if (P->opaque->reverse_step[i + 1])
+ P->left = PJ_IO_UNITS_METERS;
+ else
+ P->left = PJ_IO_UNITS_RADIANS;
+ }
+
+ /* Now, correspondingly determine forward output (= reverse input) data type */
+
+ /* First locate the last reverse-active pipeline step */
+ for (i = nsteps - 1; i >= 0; i--)
+ 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");
+ return pipeline_freeup (P, -50);
+ }
+
+ if (P->opaque->reverse_step[i + 1])
+ P->right = P->opaque->pipeline[i + 1]->left;
+ else
+ P->right = P->opaque->pipeline[i + 1]->right;
+
+ if (P->right==PJ_IO_UNITS_CLASSIC) {
+ if (P->opaque->reverse_step[i + 1])
+ P->right = PJ_IO_UNITS_RADIANS;
+ else
+ P->right = PJ_IO_UNITS_METERS;
+ }
+ pj_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");
+
+ return P;
+}
+
+#ifndef PJ_SELFTEST
+/* selftest stub */
+int pj_pipeline_selftest (void) {return 0;}
+#else
+
+int pj_pipeline_selftest (void) {
+ PJ *P;
+ PJ_OBS a, b;
+ XY cph_utm32 = {691875.63214, 6098907.82501};
+ double dist;
+
+ /* 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 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.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)
+ 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)
+ return 1002;
+
+ 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 2000;
+
+ /* zero initialize everything, then set (easting, northing) to utm(12, 55) */
+ a = b = pj_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)
+ return 2001;
+
+ /* Inverse projection */
+ a = pj_trans (P, PJ_INV, b);
+ if (pj_xy_dist (a.coo.xy, b.coo.xy) > 1e-4)
+ return 2001;
+ if (pj_xyz_dist (a.coo.xyz, b.coo.xyz) > 1e-4)
+ return 2002;
+
+ 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 3000;
+
+ a = b = pj_obs_null;
+ a.coo.lpz.lam = TORAD(12);
+ a.coo.lpz.phi = TORAD(55);
+
+ /* Forward projection */
+ b = pj_trans (P, PJ_FWD, a);
+ if (pj_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)
+ 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);
+ 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);
+ if (fabs (110574.388554153 - dist) > 1e-4)
+ return 4002;
+
+ pj_free (P);
+ return 0;
+}
+#endif
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 2386a102..9e1ac4f7 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -111,6 +111,7 @@ SET(SRC_LIBPROJ_PJ
PJ_omerc.c
PJ_ortho.c
PJ_patterson.c
+ PJ_pipeline.c
PJ_poly.c
PJ_putp2.c
PJ_putp3.c
diff --git a/src/makefile.vc b/src/makefile.vc
index 51b4320b..a2869217 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_obs_api.obj PJ_cart.obj
+ pj_obs_api.obj PJ_cart.obj PJ_pipeline.obj
geodesic = geodesic.obj
diff --git a/src/pj_ell_set.c b/src/pj_ell_set.c
index 7ca88ed1..87fb6e41 100644
--- a/src/pj_ell_set.c
+++ b/src/pj_ell_set.c
@@ -1,104 +1,108 @@
/* set ellipsoid parameters a and es */
#include <projects.h>
#include <string.h>
-#define SIXTH .1666666666666666667 /* 1/6 */
-#define RA4 .04722222222222222222 /* 17/360 */
-#define RA6 .02215608465608465608 /* 67/3024 */
-#define RV4 .06944444444444444444 /* 5/72 */
-#define RV6 .04243827160493827160 /* 55/1296 */
- int /* initialize geographic shape parameters */
-pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
- int i;
- double b=0.0, e;
- char *name;
- paralist *start = 0;
+#define SIXTH .1666666666666666667 /* 1/6 */
+#define RA4 .04722222222222222222 /* 17/360 */
+#define RA6 .02215608465608465608 /* 67/3024 */
+#define RV4 .06944444444444444444 /* 5/72 */
+#define RV6 .04243827160493827160 /* 55/1296 */
- /* clear any previous error */
- pj_ctx_set_errno(ctx,0);
+/* initialize geographic shape parameters */
+int pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
+ int i;
+ double b=0.0, e;
+ char *name;
+ paralist *start = 0;
- /* check for varying forms of ellipsoid input */
- *a = *es = 0.;
- /* R takes precedence */
- if (pj_param(ctx, pl, "tR").i)
- *a = pj_param(ctx,pl, "dR").f;
- else { /* probable elliptical figure */
+ /* clear any previous error */
+ pj_ctx_set_errno(ctx,0);
- /* check if ellps present and temporarily append its values to pl */
+ /* check for varying forms of ellipsoid input */
+ *a = *es = 0.;
+
+ /* R takes precedence */
+ if (pj_param(ctx, pl, "tR").i)
+ *a = pj_param(ctx,pl, "dR").f;
+
+ /* probable elliptical figure */
+ else {
+ /* check if ellps present and temporarily append its values to pl */
if ((name = pj_param(ctx,pl, "sellps").s) != NULL) {
- char *s;
+ char *s;
+
+ for (start = pl; start && start->next ; start = start->next) ;
+ for (i = 0; (s = pj_ellps[i].id) && strcmp(name, s) ; ++i) ;
+ if (!s) { pj_ctx_set_errno( ctx, -9); return 1; }
+ start->next = pj_mkparam(pj_ellps[i].major);
+ start->next->next = pj_mkparam(pj_ellps[i].ell);
+ }
- for (start = pl; start && start->next ; start = start->next) ;
- for (i = 0; (s = pj_ellps[i].id) && strcmp(name, s) ; ++i) ;
- if (!s) { pj_ctx_set_errno( ctx, -9); return 1; }
- start->next = pj_mkparam(pj_ellps[i].major);
- start->next->next = pj_mkparam(pj_ellps[i].ell);
- }
- *a = pj_param(ctx,pl, "da").f;
- if (pj_param(ctx,pl, "tes").i) /* eccentricity squared */
- *es = pj_param(ctx,pl, "des").f;
- else if (pj_param(ctx,pl, "te").i) { /* eccentricity */
- e = pj_param(ctx,pl, "de").f;
- *es = e * e;
- } else if (pj_param(ctx,pl, "trf").i) { /* recip flattening */
- *es = pj_param(ctx,pl, "drf").f;
- if (!*es) {
- pj_ctx_set_errno( ctx, -10);
- goto bomb;
- }
- *es = 1./ *es;
- *es = *es * (2. - *es);
- } else if (pj_param(ctx,pl, "tf").i) { /* flattening */
- *es = pj_param(ctx,pl, "df").f;
- *es = *es * (2. - *es);
- } else if (pj_param(ctx,pl, "tb").i) { /* minor axis */
- b = pj_param(ctx,pl, "db").f;
- *es = 1. - (b * b) / (*a * *a);
- } /* else *es == 0. and sphere of radius *a */
- if (!b)
- b = *a * sqrt(1. - *es);
- /* following options turn ellipsoid into equivalent sphere */
- if (pj_param(ctx,pl, "bR_A").i) { /* sphere--area of ellipsoid */
- *a *= 1. - *es * (SIXTH + *es * (RA4 + *es * RA6));
- *es = 0.;
- } else if (pj_param(ctx,pl, "bR_V").i) { /* sphere--vol. of ellipsoid */
- *a *= 1. - *es * (SIXTH + *es * (RV4 + *es * RV6));
- *es = 0.;
- } else if (pj_param(ctx,pl, "bR_a").i) { /* sphere--arithmetic mean */
- *a = .5 * (*a + b);
- *es = 0.;
- } else if (pj_param(ctx,pl, "bR_g").i) { /* sphere--geometric mean */
- *a = sqrt(*a * b);
- *es = 0.;
- } else if (pj_param(ctx,pl, "bR_h").i) { /* sphere--harmonic mean */
- *a = 2. * *a * b / (*a + b);
- *es = 0.;
- } else if ((i = pj_param(ctx,pl, "tR_lat_a").i) || /* sphere--arith. */
- pj_param(ctx,pl, "tR_lat_g").i) { /* or geom. mean at latitude */
- double tmp;
+ *a = pj_param(ctx,pl, "da").f;
+ if (pj_param(ctx,pl, "tes").i) /* eccentricity squared */
+ *es = pj_param(ctx,pl, "des").f;
+ else if (pj_param(ctx,pl, "te").i) { /* eccentricity */
+ e = pj_param(ctx,pl, "de").f;
+ *es = e * e;
+ } else if (pj_param(ctx,pl, "trf").i) { /* recip flattening */
+ *es = pj_param(ctx,pl, "drf").f;
+ if (!*es) {
+ pj_ctx_set_errno( ctx, -10);
+ goto bomb;
+ }
+ *es = 1./ *es;
+ *es = *es * (2. - *es);
+ } else if (pj_param(ctx,pl, "tf").i) { /* flattening */
+ *es = pj_param(ctx,pl, "df").f;
+ *es = *es * (2. - *es);
+ } else if (pj_param(ctx,pl, "tb").i) { /* minor axis */
+ b = pj_param(ctx,pl, "db").f;
+ *es = 1. - (b * b) / (*a * *a);
+ } /* else *es == 0. and sphere of radius *a */
+ if (!b)
+ b = *a * sqrt(1. - *es);
+ /* following options turn ellipsoid into equivalent sphere */
+ if (pj_param(ctx,pl, "bR_A").i) { /* sphere--area of ellipsoid */
+ *a *= 1. - *es * (SIXTH + *es * (RA4 + *es * RA6));
+ *es = 0.;
+ } else if (pj_param(ctx,pl, "bR_V").i) { /* sphere--vol. of ellipsoid */
+ *a *= 1. - *es * (SIXTH + *es * (RV4 + *es * RV6));
+ *es = 0.;
+ } else if (pj_param(ctx,pl, "bR_a").i) { /* sphere--arithmetic mean */
+ *a = .5 * (*a + b);
+ *es = 0.;
+ } else if (pj_param(ctx,pl, "bR_g").i) { /* sphere--geometric mean */
+ *a = sqrt(*a * b);
+ *es = 0.;
+ } else if (pj_param(ctx,pl, "bR_h").i) { /* sphere--harmonic mean */
+ *a = 2. * *a * b / (*a + b);
+ *es = 0.;
+ } else if ((i = pj_param(ctx,pl, "tR_lat_a").i) || /* sphere--arith. */
+ pj_param(ctx,pl, "tR_lat_g").i) { /* or geom. mean at latitude */
+ double tmp;
- tmp = sin(pj_param(ctx,pl, i ? "rR_lat_a" : "rR_lat_g").f);
- if (fabs(tmp) > M_HALFPI) {
+ tmp = sin(pj_param(ctx,pl, i ? "rR_lat_a" : "rR_lat_g").f);
+ if (fabs(tmp) > M_HALFPI) {
pj_ctx_set_errno(ctx,-11);
- goto bomb;
- }
- tmp = 1. - *es * tmp * tmp;
- *a *= i ? .5 * (1. - *es + tmp) / ( tmp * sqrt(tmp)) :
- sqrt(1. - *es) / tmp;
- *es = 0.;
- }
+ goto bomb;
+ }
+ tmp = 1. - *es * tmp * tmp;
+ *a *= i ? .5 * (1. - *es + tmp) / ( tmp * sqrt(tmp)) :
+ sqrt(1. - *es) / tmp;
+ *es = 0.;
+ }
bomb:
- if (start) { /* clean up temporary extension of list */
- pj_dalloc(start->next->next);
- pj_dalloc(start->next);
- start->next = 0;
- }
- if (ctx->last_errno)
- return 1;
- }
- /* some remaining checks */
- if (*es < 0.)
- { pj_ctx_set_errno( ctx, -12); return 1; }
- if (*a <= 0.)
- { pj_ctx_set_errno( ctx, -13); return 1; }
- return 0;
+ if (start) { /* clean up temporary extension of list */
+ pj_dalloc(start->next->next);
+ pj_dalloc(start->next);
+ start->next = 0;
+ }
+ if (ctx->last_errno)
+ return 1;
+ }
+ /* some remaining checks */
+ if (*es < 0.)
+ { pj_ctx_set_errno( ctx, -12); return 1; }
+ if (*a <= 0.)
+ { pj_ctx_set_errno( ctx, -13); return 1; }
+ return 0;
}
diff --git a/src/pj_init.c b/src/pj_init.c
index 79b64bb3..d5fe2433 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -29,6 +29,7 @@
#define PJ_LIB__
#include <projects.h>
+#include <geodesic.h>
#include <stdio.h>
#include <string.h>
#include <errno.h>
@@ -448,12 +449,42 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
PIN->a_orig = PIN->a;
PIN->es_orig = PIN->es;
- PIN->e = sqrt(PIN->es);
+ /* Compute some ancillary ellipsoidal parameters */
+
+ PIN->e = sqrt(PIN->es); /* eccentricity */
+ PIN->alpha = asin (PIN->e); /* angular eccentricity */
+
+ /* second eccentricity */
+ PIN->e2 = tan (PIN->alpha);
+ PIN->e2s = PIN->e2 * PIN->e2;
+
+ /* third eccentricity */
+ PIN->e3 = sin (PIN->alpha) / sqrt(2 - sin (PIN->alpha)*sin (PIN->alpha));
+ PIN->e3s = PIN->e3 * PIN->e3;
+
+ /* flattening */
+ PIN->f = 1 - cos (PIN->alpha); /* = 1 - sqrt (1 - PIN->es); */
+ PIN->rf = PIN->f? 1/PIN->f: HUGE_VAL;
+
+ /* second flattening */
+ PIN->f2 = 1/cos (PIN->alpha) - 1;
+ PIN->rf = PIN->f2? 1/PIN->f2: HUGE_VAL;
+
+ /* third flattening */
+ PIN->n = pow (tan (PIN->alpha/2), 2);
+ PIN->rn = PIN->n? 1/PIN->n: HUGE_VAL;
+
+ /* ...and a few more */
+ PIN->b = (1 - PIN->f)*PIN->a;
+ PIN->rb = 1. / PIN->b;
PIN->ra = 1. / PIN->a;
+
PIN->one_es = 1. - PIN->es;
if (PIN->one_es == 0.) { pj_ctx_set_errno( ctx, -6 ); goto bum_call; }
PIN->rone_es = 1./PIN->one_es;
+
+
/* Now that we have ellipse information check for WGS84 datum */
if( PIN->datum_type == PJD_3PARAM
&& PIN->datum_params[0] == 0.0
@@ -588,6 +619,11 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
else
PIN->from_greenwich = 0.0;
+ /* Private object for the geodesic functions */
+ PIN->geod = pj_calloc (1, sizeof (struct geod_geodesic));
+ if (0!=PIN->geod)
+ geod_init(PIN->geod, PIN->a, (1 - sqrt (1 - PIN->es)));
+
/* projection specific initialization */
if (!(PIN = (*proj)(PIN)) || ctx->last_errno) {
bum_call: /* cleanup error return */
@@ -635,6 +671,9 @@ pj_free(PJ *P) {
if( P->catalog != NULL )
pj_dalloc( P->catalog );
+ if( P->geod != NULL )
+ pj_dalloc( P->geod );
+
/* free projection parameters */
P->pfree(P);
}
diff --git a/src/pj_list.h b/src/pj_list.h
index 3569e335..f49fd574 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -101,6 +101,7 @@ PROJ_HEAD(ortel, "Ortelius Oval")
PROJ_HEAD(ortho, "Orthographic")
PROJ_HEAD(pconic, "Perspective Conic")
PROJ_HEAD(patterson, "Patterson Cylindrical")
+PROJ_HEAD(pipeline, "Transformation pipeline manager")
PROJ_HEAD(poly, "Polyconic (American)")
PROJ_HEAD(putp1, "Putnins P1")
PROJ_HEAD(putp2, "Putnins P2")
diff --git a/src/pj_log.c b/src/pj_log.c
index 54b83a05..0c711e52 100644
--- a/src/pj_log.c
+++ b/src/pj_log.c
@@ -42,13 +42,13 @@ void pj_stderr_logger( void *app_data, int level, const char *msg )
}
/************************************************************************/
-/* pj_log() */
+/* pj_vlog() */
/************************************************************************/
-void pj_log( projCtx ctx, int level, const char *fmt, ... )
+/* 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 )
{
- va_list args;
char *msg_buf;
if( level > ctx->debug_level )
@@ -58,14 +58,63 @@ void pj_log( projCtx ctx, int level, const char *fmt, ... )
if( msg_buf == NULL )
return;
- va_start( args, fmt );
-
/* we should use vsnprintf where available once we add configure detect.*/
vsprintf( msg_buf, fmt, args );
- va_end( args );
-
ctx->logger( ctx->app_data, level, msg_buf );
-
+
free( msg_buf );
}
+
+
+/************************************************************************/
+/* pj_log() */
+/************************************************************************/
+
+void pj_log( projCtx ctx, int level, const char *fmt, ... )
+
+{
+ va_list args;
+
+ if( level > ctx->debug_level )
+ return;
+
+ va_start( args, 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 ccec8a64..716c6095 100644
--- a/src/pj_obs_api.c
+++ b/src/pj_obs_api.c
@@ -35,6 +35,7 @@
#define PJ_OBS_C
#include <proj.h>
#include <projects.h>
+#include <geodesic.h>
#include <float.h>
#include <math.h>
@@ -57,20 +58,30 @@ const PJ_OBS pj_obs_null = {
/* 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]);
+/* Geodesic distance between two points with angular 2D coordinates */
+double pj_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);
+ return s12;
}
-/* 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]);
+/* Euclidean distance between two points with linear 2D coordinates */
+double pj_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) {
+ return hypot (hypot (a.x - b.x, a.y - b.y), a.z - b.z);
}
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;
@@ -85,6 +96,10 @@ PJ_OBS pj_fwdobs (PJ_OBS obs, PJ *P) {
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;
@@ -99,7 +114,7 @@ PJ_OBS pj_invobs (PJ_OBS obs, PJ *P) {
/* Apply the transformation P to the observation obs */
-PJ_OBS pj_apply (PJ *P, enum pj_direction direction, PJ_OBS obs) {
+PJ_OBS pj_trans (PJ *P, enum pj_direction direction, PJ_OBS obs) {
if (0==P)
return obs;
@@ -150,7 +165,7 @@ double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs) {
}
}
- return pj_obs_dist_3d (o, obs);
+ return pj_xyz_dist (o.coo.xyz, obs.coo.xyz);
}
@@ -179,14 +194,19 @@ void pj_error_set (PJ *P, int 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) {
+/* 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);
- ctx->debug_level = debuglevel;
+ previous = ctx->debug_level;
+ if (PJ_LOG_TELL==log_level)
+ return previous;
+ ctx->debug_level = log_level;
+ return previous;
}
diff --git a/src/pj_strerrno.c b/src/pj_strerrno.c
index b46f143c..9dc180b5 100644
--- a/src/pj_strerrno.c
+++ b/src/pj_strerrno.c
@@ -4,82 +4,82 @@
#include <errno.h>
#include <string.h>
- static char *
+ static char *
pj_err_list[] = {
- "no arguments in initialization list", /* -1 */
- "no options found in 'init' file", /* -2 */
- "no colon in init= string", /* -3 */
- "projection not named", /* -4 */
- "unknown projection id", /* -5 */
- "effective eccentricity = 1.", /* -6 */
- "unknown unit conversion id", /* -7 */
- "invalid boolean param argument", /* -8 */
- "unknown elliptical parameter name", /* -9 */
- "reciprocal flattening (1/f) = 0", /* -10 */
- "|radius reference latitude| > 90", /* -11 */
- "squared eccentricity < 0", /* -12 */
- "major axis or radius = 0 or not given", /* -13 */
- "latitude or longitude exceeded limits", /* -14 */
- "invalid x or y", /* -15 */
- "improperly formed DMS value", /* -16 */
- "non-convergent inverse meridional dist", /* -17 */
- "non-convergent inverse phi2", /* -18 */
- "acos/asin: |arg| >1.+1e-14", /* -19 */
- "tolerance condition error", /* -20 */
- "conic lat_1 = -lat_2", /* -21 */
- "lat_1 >= 90", /* -22 */
- "lat_1 = 0", /* -23 */
- "lat_ts >= 90", /* -24 */
- "no distance between control points", /* -25 */
- "projection not selected to be rotated", /* -26 */
- "W <= 0 or M <= 0", /* -27 */
- "lsat not in 1-5 range", /* -28 */
- "path not in range", /* -29 */
- "h <= 0", /* -30 */
- "k <= 0", /* -31 */
- "lat_0 = 0 or 90 or alpha = 90", /* -32 */
- "lat_1=lat_2 or lat_1=0 or lat_2=90", /* -33 */
- "elliptical usage required", /* -34 */
- "invalid UTM zone number", /* -35 */
- "arg(s) out of range for Tcheby eval", /* -36 */
- "failed to find projection to be rotated", /* -37 */
- "failed to load datum shift file", /* -38 */
- "both n & m must be spec'd and > 0", /* -39 */
- "n <= 0, n > 1 or not specified", /* -40 */
- "lat_1 or lat_2 not specified", /* -41 */
- "|lat_1| == |lat_2|", /* -42 */
- "lat_0 is pi/2 from mean lat", /* -43 */
- "unparseable coordinate system definition", /* -44 */
- "geocentric transformation missing z or ellps", /* -45 */
- "unknown prime meridian conversion id", /* -46 */
- "illegal axis orientation combination", /* -47 */
- "point not within available datum shift grids", /* -48 */
- "invalid sweep axis, choose x or y", /* -49 */
+ "no arguments in initialization list", /* -1 */
+ "no options found in 'init' file", /* -2 */
+ "no colon in init= string", /* -3 */
+ "projection not named", /* -4 */
+ "unknown projection id", /* -5 */
+ "effective eccentricity = 1.", /* -6 */
+ "unknown unit conversion id", /* -7 */
+ "invalid boolean param argument", /* -8 */
+ "unknown elliptical parameter name", /* -9 */
+ "reciprocal flattening (1/f) = 0", /* -10 */
+ "|radius reference latitude| > 90", /* -11 */
+ "squared eccentricity < 0", /* -12 */
+ "major axis or radius = 0 or not given", /* -13 */
+ "latitude or longitude exceeded limits", /* -14 */
+ "invalid x or y", /* -15 */
+ "improperly formed DMS value", /* -16 */
+ "non-convergent inverse meridional dist", /* -17 */
+ "non-convergent inverse phi2", /* -18 */
+ "acos/asin: |arg| >1.+1e-14", /* -19 */
+ "tolerance condition error", /* -20 */
+ "conic lat_1 = -lat_2", /* -21 */
+ "lat_1 >= 90", /* -22 */
+ "lat_1 = 0", /* -23 */
+ "lat_ts >= 90", /* -24 */
+ "no distance between control points", /* -25 */
+ "projection not selected to be rotated", /* -26 */
+ "W <= 0 or M <= 0", /* -27 */
+ "lsat not in 1-5 range", /* -28 */
+ "path not in range", /* -29 */
+ "h <= 0", /* -30 */
+ "k <= 0", /* -31 */
+ "lat_0 = 0 or 90 or alpha = 90", /* -32 */
+ "lat_1=lat_2 or lat_1=0 or lat_2=90", /* -33 */
+ "elliptical usage required", /* -34 */
+ "invalid UTM zone number", /* -35 */
+ "arg(s) out of range for Tcheby eval", /* -36 */
+ "failed to find projection to be rotated", /* -37 */
+ "failed to load datum shift file", /* -38 */
+ "both n & m must be spec'd and > 0", /* -39 */
+ "n <= 0, n > 1 or not specified", /* -40 */
+ "lat_1 or lat_2 not specified", /* -41 */
+ "|lat_1| == |lat_2|", /* -42 */
+ "lat_0 is pi/2 from mean lat", /* -43 */
+ "unparseable coordinate system definition", /* -44 */
+ "geocentric transformation missing z or ellps", /* -45 */
+ "unknown prime meridian conversion id", /* -46 */
+ "illegal axis orientation combination", /* -47 */
+ "point not within available datum shift grids", /* -48 */
+ "invalid sweep axis, choose x or y", /* -49 */
+ "malformed pipeline", /* -50 */
};
- char *
-pj_strerrno(int err)
-{
+
+char *pj_strerrno(int err) {
static char note[50];
- if (err > 0)
+ if (0==err)
+ return 0;
+
+ if (err > 0) {
#ifdef HAVE_STRERROR
return strerror(err);
#else
- {
sprintf(note,"no system list, errno: %d\n", err);
return note;
- }
#endif
- else if (err < 0) {
+ }
+
+ if (err < 0) {
size_t adjusted_err = - err - 1;
if (adjusted_err < (sizeof(pj_err_list) / sizeof(char *)))
return(pj_err_list[adjusted_err]);
- else
- {
- sprintf( note, "invalid projection system error (%d)",
- err );
+ else {
+ sprintf( note, "invalid projection system error (%d)", err );
return note;
}
- } else
- return NULL;
+ }
}
diff --git a/src/proj.def b/src/proj.def
index 09b32c32..3a5e3de9 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -92,14 +92,20 @@ EXPORTS
pj_error @90
pj_create @91
pj_create_argv @92
- pj_apply @93
+ pj_trans @93
pj_fwdobs @94
pj_invobs @95
pj_roundtrip @96
- pj_debug_set @97
+ pj_log_level @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
+ pj_log_error @104
+ pj_log_debug @105
+ pj_log_trace @106
+ pj_lp_dist @107
+ pj_xy_dist @108
+ pj_xyz_dist @109
diff --git a/src/proj.h b/src/proj.h
index 4b7cfc87..6c33a916 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -138,6 +138,8 @@ extern "C" {
#include <proj_api.h>
#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 */
@@ -189,13 +191,14 @@ typedef struct { double e, z, N; } PJ_EZN;
/* Ellipsoidal parameters */
typedef struct { double a, f; } PJ_AF;
+/* Avoid preprocessor renaming and implicit type-punning: Use unions to make it explicit */
union PJ_COORD {
PJ_XYZT xyzt;
PJ_UVWT uvwt;
PJ_ENHT enht;
PJ_LPZT lpzt;
PJ_ENH enh;
- double v[4]; /* Who cares - it's just a vector! */
+ double v[4]; /* It's just a vector */
XYZ xyz;
UVW uvw;
LPZ lpz;
@@ -204,28 +207,34 @@ union PJ_COORD {
LP lp;
};
-
-
-/* Avoid preprocessor renaming and implicit type-punning: Use a union to make it explicit */
union PJ_TRIPLET {
PJ_OPK opk;
PJ_ENH enh;
PJ_EZN ezn;
- PJ_AF af;
- double v[3]; /* Who cares - it's just a vector! */
+ PJ_DMS dms;
+ double v[3]; /* It's just a vector */
XYZ xyz;
LPZ lpz;
UVW uvw;
XY xy;
LP lp;
UV uv;
+ PJ_AF af;
+};
+
+union PJ_PAIR {
+ XY xy;
+ LP lp;
+ UV uv;
+ PJ_AF af;
+ double v[2]; /* Yes - It's really just a vector! */
};
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 */
+ 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 */
@@ -240,24 +249,6 @@ PJ *pj_create_argv (int argc, char **argv);
void pj_free (PJ *P);
int pj_error (PJ *P);
-/* 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 */
enum pj_direction {
@@ -265,17 +256,19 @@ enum pj_direction {
PJ_IDENT = 0, /* Do nothing */
PJ_INV = -1 /* Inverse */
};
-PJ_OBS pj_apply (PJ *P, enum pj_direction direction, PJ_OBS obs);
+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);
-/* Euclidean distance between two 2D coordinates stored in PJ_OBSs */
-double pj_obs_dist_2d (PJ_OBS a, PJ_OBS b);
+/* Geodesic distance between two points with angular 2D coordinates */
+double pj_lp_dist (PJ *P, LP a, LP b);
-/* Euclidean distance between two 3D coordinates stored in PJ_OBSs */
-double pj_obs_dist_3d (PJ_OBS a, PJ_OBS b);
+/* 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);
#ifndef PJ_OBS_C
@@ -292,6 +285,42 @@ extern const PJ *pj_shutdown;
#endif
+
+
+
+
+
+/* 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 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_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);
+
+
#ifdef __cplusplus
}
#endif
diff --git a/src/proj_api.h b/src/proj_api.h
index 706e35ab..863795f0 100644
--- a/src/proj_api.h
+++ b/src/proj_api.h
@@ -64,11 +64,8 @@ extern "C" {
/* 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
+#define RAD_TO_DEG 57.295779513082321
+#define DEG_TO_RAD .017453292519943296
#if defined(PROJECTS_H) || defined(PROJ_H)
@@ -78,7 +75,10 @@ extern int pj_errno; /* global error return code */
#ifndef PROJ_H
-/* In proj.h these macros are replaced by the enumeration pj_debug_level */
+extern char const pj_release[]; /* global release id string */
+extern int pj_errno; /* global error return code */
+
+/* In proj.h these macros are replaced by the enumeration pj_log_level */
#define PJ_LOG_NONE 0
#define PJ_LOG_ERROR 1
#define PJ_LOG_DEBUG_MAJOR 2
@@ -90,7 +90,7 @@ extern int pj_errno; /* global error return code */
/* 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 projXY XY
# define projLP LP
# define projXYZ XYZ
# define projLPZ LPZ
@@ -107,8 +107,6 @@ extern int pj_errno; /* global error return code */
#endif
-
-
/* If included *after* proj.h finishes, we have alternative names */
/* file reading api, like stdio */
typedef int *PAFile;
diff --git a/src/projects.h b/src/projects.h
index d07c43c7..47a3364a 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -169,6 +169,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 geod_geodesic;
struct pj_opaque;
struct ARG_list;
struct FACTORS;
@@ -212,6 +213,7 @@ struct PJconsts {
projCtx_t *ctx;
const char *descr; /* From pj_list.h or individual PJ_*.c file */
paralist *params; /* Parameter list */
+ struct geod_geodesic *geod; /* For geodesic computations */
struct pj_opaque *opaque; /* Projection specific parameters, Defined in PJ_*.c */
@@ -272,6 +274,7 @@ struct PJconsts {
/* The eccentricities */
+ double alpha; /* angular eccentricity */
double e; /* first eccentricity */
double es; /* first eccentricity squared */
double e2; /* second eccentricity */
generated by cgit v1.2.3 (git 2.25.1) at 2026-03-01 02:27:39 +0000