aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--appveyor.yml1
-rw-r--r--src/PJ_healpix.c3
-rw-r--r--src/PJ_pipeline.c2
-rw-r--r--src/pj_ell_set.c490
-rw-r--r--src/pj_init.c2
-rw-r--r--src/pj_malloc.c10
-rw-r--r--src/pj_param.c6
-rw-r--r--src/proj.def3
-rw-r--r--src/proj.h4
-rw-r--r--src/proj_4D_api.c43
-rw-r--r--src/projects.h12
-rw-r--r--test/gie/ellipsoid.gie121
-rwxr-xr-xtravis/install.sh1
13 files changed, 657 insertions, 41 deletions
diff --git a/appveyor.yml b/appveyor.yml
index 9f66b22e..43669a40 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -62,6 +62,7 @@ test_script:
- gie.exe ..\test\gie\more_builtins.gie
- gie.exe ..\test\gie\deformation.gie
- gie.exe ..\test\gie\axisswap.gie
+ - gie.exe ..\test\gie\ellipsoid.gie
deploy: off
diff --git a/src/PJ_healpix.c b/src/PJ_healpix.c
index d8e1a921..5becfd74 100644
--- a/src/PJ_healpix.c
+++ b/src/PJ_healpix.c
@@ -30,6 +30,7 @@
*****************************************************************************/
# define PJ_LIB__
# include <errno.h>
+# include "proj_internal.h"
# include <proj.h>
# include "projects.h"
@@ -624,7 +625,7 @@ PJ *PROJECTION(healpix) {
return destructor(P, ENOMEM);
Q->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */
P->a = P->a*sqrt(0.5*Q->qp); /* Set P->a to authalic radius. */
- P->ra = 1.0/P->a;
+ pj_calc_ellps_params (P, P->a, P->es); /* Ensure we have a consistent parameter set */
P->fwd = e_healpix_forward;
P->inv = e_healpix_inverse;
} else {
diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c
index e1568423..19b8e18d 100644
--- a/src/PJ_pipeline.c
+++ b/src/PJ_pipeline.c
@@ -317,7 +317,7 @@ static void set_ellipsoid(PJ *P) {
/* Check if there's any ellipsoid specification in the global params. */
/* If not, use WGS84 as default */
- if (pj_ell_set(P->ctx, P->params, &P->a, &P->es)) {
+ if (0 != pj_ellipsoid (P)) {
P->a = 6378137.0;
P->es = .00669438002290341575;
}
diff --git a/src/pj_ell_set.c b/src/pj_ell_set.c
index 8a9b39a1..8ba928b7 100644
--- a/src/pj_ell_set.c
+++ b/src/pj_ell_set.c
@@ -1,13 +1,424 @@
/* set ellipsoid parameters a and es */
#include <string.h>
+#include <errno.h>
+#include <proj.h>
#include "proj_internal.h"
#include "projects.h"
+
+/* series coefficients for calculating ellipsoid-equivalent spheres */
#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 */
+/* Prototypes of the pj_ellipsoid helper functions */
+static int pj_ellps_handler (PJ *P);
+static int pj_size (PJ *P);
+static int pj_shape (PJ *P);
+static int pj_spherification (PJ *P);
+
+static paralist *pj_get_param (paralist *list, char *key);
+static char *pj_param_value (paralist *list);
+static PJ_ELLPS *pj_find_ellps (char *name);
+
+
+/***************************************************************************************/
+int pj_ellipsoid (PJ *P) {
+/****************************************************************************************
+ This is a replacement for the clasic PROJ pj_ell_set function. The main difference
+ is that pj_ellipsoid augments the PJ object with a copy of the exact tags used to
+ define its related ellipsoid.
+
+ This makes it possible to let a new PJ object inherit the geometrical properties
+ of an existing one.
+
+ A complete ellipsoid definition comprises a size (primary) and a shape (secondary)
+ parameter.
+
+ Size parameters supported are:
+ R, defining the radius of a spherical planet
+ a, defining the semimajor axis of an ellipsoidal planet
+
+ Shape parameters supported are:
+ rf, the reverse flattening of the ellipsoid
+ f, the flattening of the ellipsoid
+ es, the eccentricity squared
+ e, the eccentricity
+ b, the semiminor axis
+
+ The ellps=xxx parameter provides both size and shape for a number of built in
+ ellipsoid definitions.
+
+ The ellipsoid definition may be augmented with a spherification flag, turning
+ the ellipsoid into a sphere with features defined by the ellipsoid.
+
+ Spherification parameters supported are:
+ R_A, which gives a sphere with the same surface area as the ellipsoid
+ R_A, which gives a sphere with the same volume as the ellipsoid
+
+ R_a, which gives a sphere with R = (a + b)/2 (arithmetic mean)
+ R_g, which gives a sphere with R = sqrt(a*b) (geometric mean)
+ R_h, which gives a sphere with R = 2*a*b/(a+b) (harmonic mean)
+
+ R_lat_a=phi, which gives a sphere with R being the arithmetic mean of
+ of the corresponding ellipsoid at latitude phi.
+ R_lat_g=phi, which gives a sphere with R being the geometric mean of
+ of the corresponding ellipsoid at latitude phi.
+
+ If R is given as size parameter, any shape and spherification parameters
+ given are ignored.
+
+ If size and shape is given as ellps=xxx, later shape and size parameters
+ are are taken into account as modifiers for the built in ellipsoid definition.
+
+ While this may seem strange, it is in accordance with historical PROJ
+ behaviour. It can e.g. be used to define coordinates on the ellipsoid
+ scaled to unit semimajor axis by specifying "+ellps=xxx +a=1"
+
+****************************************************************************************/
+ int last_errno;
+
+ last_errno = proj_errno_reset (P);
+
+ P->def_size = P->def_shape = P->def_spherification = P->def_ellps = 0;
+
+ /* If an ellps argument is specified, start by using that */
+ if (0 != pj_ellps_handler (P))
+ return proj_errno (P);
+
+ /* We may overwrite the size */
+ if (0 != pj_size (P))
+ return proj_errno (P);
+
+ /* We may also overwrite the shape */
+ if (0 != pj_shape (P))
+ return proj_errno (P);
+
+ /* When we're done with it, we compute all related ellipsoid parameters */
+ pj_calc_ellps_params (P, P->a, P->es);
+
+ /* And finally, we may turn it into a sphere */
+ if (0 != pj_spherification (P))
+ return proj_errno (P);
+
+ if (proj_errno (P))
+ return proj_errno (P);
+
+ /* success - restore previous error status */
+ return proj_errno_restore (P, last_errno);
+}
+
+
+/***************************************************************************************/
+static int pj_ellps_handler (PJ *P) {
+/***************************************************************************************/
+ PJ B;
+ PJ_ELLPS *ellps;
+ paralist *par = 0;
+ char *def, *name;
+
+ /* ellps specified? */
+ par = pj_get_param (P->params, "ellps");
+ if (0==par)
+ return 0;
+
+ /* This amounts to zeroing all ellipsoidal parameters in P */
+ memset (&B, 0, sizeof (B));
+ pj_inherit_ellipsoid_defs(&B, P);
+
+ /* move B into P's context */
+ B.ctx = P->ctx;
+
+ /* Store definition */
+ P->def_ellps = def = par->param;
+ par->used = 1;
+
+ if (strlen (def) < 7)
+ return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+ name = def + 6;
+ ellps = pj_find_ellps (name);
+ if (0==ellps)
+ return proj_errno_set (P, PJD_ERR_UNKNOWN_ELLP_PARAM);
+
+ /* Now, we have the right size and shape parameters and can produce */
+ /* a fake PJ to make pj_shape do the hard work for us */
+ B = *P;
+ B.params = pj_mkparam (ellps->major);
+ B.params->next = pj_mkparam (ellps->ell);
+
+ pj_size (&B);
+ pj_shape (&B);
+
+ P->a = B.a;
+ P->b = B.b;
+ P->f = B.f;
+ P->e = B.e;
+ P->es = B.es;
+
+ pj_dealloc (B.params->next);
+ pj_dealloc (B.params);
+ if (proj_errno (&B))
+ return proj_errno (&B);
+ return 0;
+}
+
+
+/***************************************************************************************/
+static int pj_size (PJ *P) {
+/***************************************************************************************/
+ char *keys[] = {"R", "a"};
+ paralist *par = 0;
+ size_t i, len;
+ len = sizeof (keys) / sizeof (char *);
+ /* Check which size key is specified */
+ for (i = 0; i < len; i++) {
+ par = pj_get_param (P->params, keys[i]);
+ if (par)
+ break;
+ }
+
+ /* A size parameter *must* be given, but may have been given as ellps prior */
+ if (0==par)
+ return 0!=P->a? 0: proj_errno_set (P, PJD_ERR_MAJOR_AXIS_NOT_GIVEN);
+
+ P->def_size = par->param;
+ par->used = 1;
+ P->a = pj_atof (pj_param_value (par));
+ if (P->a <= 0)
+ return proj_errno_set (P, PJD_ERR_MAJOR_AXIS_NOT_GIVEN);
+ if (HUGE_VAL==P->a)
+ return proj_errno_set (P, PJD_ERR_MAJOR_AXIS_NOT_GIVEN);
+
+ if ('R'==par->param[0]) {
+ P->es = P->f = P->e = P->rf = 0;
+ P->b = P->a;
+ }
+ return 0;
+}
+
+
+/***************************************************************************************/
+static int pj_shape (PJ *P) {
+/***************************************************************************************/
+ char *keys[] = {"rf", "f", "es", "e", "b"};
+ paralist *par = 0;
+ char *def = 0;
+ size_t i, len;
+
+ par = 0;
+ len = sizeof (keys) / sizeof (char *);
+
+ /* Check which shape key is specified */
+ for (i = 0; i < len; i++) {
+ par = pj_get_param (P->params, keys[i]);
+ if (par)
+ break;
+ }
+
+ /* Not giving a shape parameter means selecting a sphere, unless shape */
+ /* has been selected previously via ellps=xxx */
+ if (0==par && P->es != 0)
+ return 0;
+ if (0==par && P->es==0) {
+ P->es = P->f = 0;
+ P->b = P->a;
+ return 0;
+ }
+
+ P->def_shape = def = par->param;
+ par->used = 1;
+ P->es = P->f = P->b = P->e = P->rf = 0;
+
+ switch (i) {
+
+ /* reverse flattening, rf */
+ case 0:
+ P->rf = pj_atof (pj_param_value (par));
+ if (HUGE_VAL==P->rf)
+ return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+ P->f = 1 / P->rf;
+ P->es = 2*P->f - P->f*P->f;
+ break;
+
+ /* flattening, f */
+ case 1:
+ P->f = pj_atof (pj_param_value (par));
+ if (HUGE_VAL==P->f)
+ return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+ if (0==P->f)
+ return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+ P->rf = 1 / P->f;
+ P->es = 2*P->f - P->f*P->f;
+ break;
+
+ /* eccentricity squared, es */
+ case 2:
+ P->es = pj_atof (pj_param_value (par));
+ if (HUGE_VAL==P->es)
+ return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+ if (1==P->es)
+ return proj_errno_set (P, PJD_ERR_ECCENTRICITY_IS_ONE);
+ break;
+
+ /* eccentricity, e */
+ case 3:
+ P->e = pj_atof (pj_param_value (par));
+ if (HUGE_VAL==P->e)
+ return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+ if (0==P->e)
+ return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+ if (1==P->e)
+ return proj_errno_set (P, PJD_ERR_ECCENTRICITY_IS_ONE);
+ P->es = P->e * P->e;
+ break;
+
+ /* semiminor axis, b */
+ case 4:
+ P->b = pj_atof (pj_param_value (par));
+ if (HUGE_VAL==P->b)
+ return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+ if (0==P->b)
+ return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+ P->f = (P->a - P->b) / P->a;
+ P->es = 2*P->f - P->f*P->f;
+ break;
+ default:
+ return PJD_ERR_INVALID_ARG;
+
+ }
+ if (P->es < 0)
+ return proj_errno_set (P, PJD_ERR_ES_LESS_THAN_ZERO);
+ return 0;
+}
+
+
+/***************************************************************************************/
+static int pj_spherification (PJ *P) {
+/***************************************************************************************/
+ char *keys[] = {"R_A", "R_V", "R_a", "R_g", "R_h", "R_lat_a", "R_lat_g"};
+ size_t len, i;
+ paralist *par = 0;
+ char *def = 0;
+
+ double t;
+ char *v, *endp;
+
+ len = sizeof (keys) / sizeof (char *);
+ P->def_spherification = 0;
+
+ /* Check which spherification key is specified */
+ for (i = 0; i < len; i++) {
+ par = pj_get_param (P->params, keys[i]);
+ if (par)
+ break;
+ }
+
+ /* No spherification specified? Then we're done */
+ if (i==len)
+ return 0;
+
+ /* Store definition */
+ P->def_spherification = def = par->param;
+ par->used = 1;
+
+ switch (i) {
+
+ /* R_A - a sphere with same area as ellipsoid */
+ case 0:
+ P->a *= 1. - P->es * (SIXTH + P->es * (RA4 + P->es * RA6));
+ break;
+
+ /* R_V - a sphere with same volume as ellipsoid */
+ case 1:
+ P->a *= 1. - P->es * (SIXTH + P->es * (RV4 + P->es * RV6));
+ break;
+
+ /* R_a - a sphere with R = the arithmetic mean of the ellipsoid */
+ case 2:
+ P->a = (P->a + P->b) / 2;
+ break;
+
+ /* R_g - a sphere with R = the geometric mean of the ellipsoid */
+ case 3:
+ P->a = sqrt (P->a * P->b);
+ break;
+
+ /* R_h - a sphere with R = the harmonic mean of the ellipsoid */
+ case 4:
+ if (P->a + P->b == 0)
+ return proj_errno_set (P, PJD_ERR_TOLERANCE_CONDITION);
+ P->a = (2*P->a * P->b) / (P->a + P->b);
+ break;
+
+ /* R_lat_a - a sphere with R = the arithmetic mean of the ellipsoid at given latitude */
+ case 5:
+ /* R_lat_g - a sphere with R = the geometric mean of the ellipsoid at given latitude */
+ case 6:
+ v = pj_param_value (par);
+ t = proj_dmstor (v, &endp);
+ if (fabs (t) > M_HALFPI)
+ return proj_errno_set (P, PJD_ERR_REF_RAD_LARGER_THAN_90);
+ t = sin (t);
+ t = 1 - P->es * t * t;
+ if (i==5) /* arithmetic */
+ P->a *= (1. - P->es + t) / (2 * t * sqrt(t));
+ else /* geometric */
+ P->a *= sqrt (1 - P->es) / t;
+ break;
+ }
+
+ /* Clean up the ellipsoidal parameters to reflect the sphere */
+ P->es = P->e = P->rf = P->f = 0;
+ P->b = P->a;
+ pj_calc_ellps_params (P, P->a, 0);
+
+ return 0;
+}
+
+
+
+
+
+/* locate parameter in list */
+static paralist *pj_get_param (paralist *list, char *key) {
+ size_t l = strlen(key);
+ while (list && !(0==strncmp(list->param, key, l) && (0==list->param[l] || list->param[l] == '=') ) )
+ list = list->next;
+ return list;
+}
+
+
+static char *pj_param_value (paralist *list) {
+ char *key, *value;
+ if (0==list)
+ return 0;
+
+ key = list->param;
+ value = strchr (key, '=');
+
+ /* a flag (i.e. a key without value) has its own name (key) as value */
+ return value? value + 1: key;
+}
+
+
+static PJ_ELLPS *pj_find_ellps (char *name) {
+ int i;
+ char *s;
+ if (0==name)
+ return 0;
+
+ /* Search through internal ellipsoid list for name */
+ for (i = 0; (s = pj_ellps[i].id) && strcmp(name, s) ; ++i);
+ if (0==s)
+ return 0;
+ return pj_ellps + i;
+}
+
+
+
+
+
/* copy ellipsoidal parameters from src to dst */
void pj_inherit_ellipsoid_defs(const PJ *src, PJ *dst) {
@@ -44,13 +455,41 @@ void pj_inherit_ellipsoid_defs(const PJ *src, PJ *dst) {
dst->a_orig = src->a_orig;
}
+
+/***************************************************************************************/
int pj_calc_ellps_params(PJ *P, double a, double es) {
+/****************************************************************************************
+ Calculate a large number of ancillary ellipsoidal parameters, in addition to
+ the two traditional PROJ defining parameters: Semimajor axis, a, and the
+ eccentricity squared, es.
+
+ Most of these parameters are fairly cheap to compute in comparison to the overall
+ effort involved in initializing a PJ object. They may, however, take a substantial
+ part of the time taken in computing an individual point transformation.
+
+ So by providing them up front, we can amortize the (already modest) cost over all
+ transformations carried out over the entire lifetime of a PJ object, rather than
+ incur that cost for every single transformation.
+
+ Most of the parameter calculations here are based on the "angular eccentricity",
+ i.e. the angle, measured from the semiminor axis, of a line going from the north
+ pole to one of the foci of the ellipsoid - or in other words: The arc sine of the
+ eccentricity.
+
+ The formulae used are mostly taken from:
+
+ Richard H. Rapp: Geometric Geodesy, Part I, (178 pp, 1991).
+ Columbus, Ohio: Dept. of Geodetic Science
+ and Surveying, Ohio State University.
+
+****************************************************************************************/
P->a = a;
P->es = es;
/* Compute some ancillary ellipsoidal parameters */
- P->e = sqrt(P->es); /* eccentricity */
+ if (P->e==0)
+ P->e = sqrt(P->es); /* eccentricity */
P->alpha = asin (P->e); /* angular eccentricity */
/* second eccentricity */
@@ -62,7 +501,8 @@ int pj_calc_ellps_params(PJ *P, double a, double es) {
P->e3s = P->e3 * P->e3;
/* flattening */
- P->f = 1 - cos (P->alpha); /* = 1 - sqrt (1 - PIN->es); */
+ if (0==P->f)
+ P->f = 1 - cos (P->alpha); /* = 1 - sqrt (1 - PIN->es); */
P->rf = P->f != 0.0 ? 1.0/P->f: HUGE_VAL;
/* second flattening */
@@ -74,7 +514,8 @@ int pj_calc_ellps_params(PJ *P, double a, double es) {
P->rn = P->n != 0.0 ? 1/P->n: HUGE_VAL;
/* ...and a few more */
- P->b = (1 - P->f)*P->a;
+ if (0==P->b)
+ P->b = (1 - P->f)*P->a;
P->rb = 1. / P->b;
P->ra = 1. / P->a;
@@ -89,8 +530,45 @@ int pj_calc_ellps_params(PJ *P, double a, double es) {
return 0;
}
-/* initialize geographic shape parameters */
+
+
+#ifndef KEEP_ORIGINAL_PJ_ELL_SET
int pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
+ PJ B;
+ int ret;
+
+ memset (&B, 0, sizeof (B));
+ B.ctx = ctx;
+ B.params = pl;
+
+ ret = pj_ellipsoid (&B);
+ if (ret)
+ return ret;
+
+ *a = B.a;
+ *es = B.es;
+ return 0;
+}
+#else
+
+
+
+/**************************************************************************************/
+int pj_ell_set (projCtx ctx, paralist *pl, double *a, double *es) {
+/***************************************************************************************
+ Initialize ellipsoidal parameters: This is the original ellipsoid setup
+ function by Gerald Evenden - significantly more compact than pj_ellipsoid and
+ its many helper functions, and still quite readable.
+
+ It is, however, also so tight that it is hard to modify and add functionality,
+ and equally hard to find the right place to add further commentary for improved
+ future maintainability.
+
+ Hence, when the need to store in the PJ object, the parameters actually used to
+ define the ellipsoid came up, rather than modifying this little gem of
+ "economy of expression", a much more verbose reimplementation, pj_ellipsoid,
+ was written.
+***************************************************************************************/
int i;
double b=0.0, e;
char *name;
@@ -120,6 +598,7 @@ int pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
}
*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 */
@@ -142,6 +621,8 @@ int pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
} /* else *es == 0. and sphere of radius *a */
if (b == 0.0)
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));
@@ -192,3 +673,4 @@ bomb:
{ pj_ctx_set_errno( ctx, -13); return 1; }
return 0;
}
+#endif
diff --git a/src/pj_init.c b/src/pj_init.c
index 53fd1039..0cdfcc4d 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -536,7 +536,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
return pj_default_destructor (PIN, PJD_ERR_MISSING_ARGS);
if (PIN->need_ellps) {
- if (pj_ell_set(ctx, start, &PIN->a, &PIN->es)) {
+ if (0 != pj_ellipsoid (PIN)) {
pj_log (ctx, PJ_LOG_DEBUG_MINOR, "pj_init_ctx: Must specify ellipsoid or sphere");
return pj_default_destructor (PIN, PJD_ERR_MISSING_ARGS);
}
diff --git a/src/pj_malloc.c b/src/pj_malloc.c
index c003c717..c9275074 100644
--- a/src/pj_malloc.c
+++ b/src/pj_malloc.c
@@ -103,7 +103,7 @@ void pj_dalloc(void *ptr) {
/**********************************************************************/
free(ptr);
}
-
+
/**********************************************************************/
void *pj_dealloc (void *ptr) {
@@ -150,7 +150,7 @@ void *pj_dealloc_params (projCtx ctx, paralist *start, int errlev) {
Also called from pj_init_ctx when encountering errors before the PJ
proper is allocated.
-******************************************************************************/
+******************************************************************************/
paralist *t, *n;
for (t = start; t; t = n) {
n = t->next;
@@ -177,12 +177,12 @@ void *pj_default_destructor (PJ *P, int errlev) { /* Destructor */
if (0==P)
return 0;
-
+
/* free grid lists */
pj_dealloc( P->gridlist );
pj_dealloc( P->vgridlist_geoid );
pj_dealloc( P->catalog_name );
-
+
/* We used to call pj_dalloc( P->catalog ), but this will leak */
/* memory. The safe way to clear catalog and grid is to call */
/* pj_gc_unloadall(pj_get_default_ctx()); and pj_deallocate_grids(); */
@@ -191,7 +191,7 @@ void *pj_default_destructor (PJ *P, int errlev) { /* Destructor */
/* free the interface to Charles Karney's geodesic library */
pj_dealloc( P->geod );
-
+
/* free parameter list elements */
pj_dealloc_params (pj_get_ctx(P), P->params, errlev);
diff --git a/src/pj_param.c b/src/pj_param.c
index 93f5cf53..ee952eca 100644
--- a/src/pj_param.c
+++ b/src/pj_param.c
@@ -2,8 +2,9 @@
#include <projects.h>
#include <stdio.h>
#include <string.h>
- paralist * /* create parameter list entry */
-pj_mkparam(char *str) {
+
+/* create parameter list entry */
+paralist *pj_mkparam(char *str) {
paralist *newitem;
if((newitem = (paralist *)pj_malloc(sizeof(paralist) + strlen(str))) != NULL) {
@@ -16,6 +17,7 @@ pj_mkparam(char *str) {
return newitem;
}
+
/************************************************************************/
/* pj_param() */
/* */
diff --git a/src/proj.def b/src/proj.def
index c282525c..ca10489a 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -147,3 +147,6 @@ EXPORTS
proj_angular_input @133
proj_angular_output @134
+
+ pj_ellipsoid @135
+ pj_calc_ellps_params @136
diff --git a/src/proj.h b/src/proj.h
index fbd59b45..f4640e15 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -382,9 +382,9 @@ double proj_xyz_dist (XYZ a, XYZ b);
/* Set or read error level */
int proj_errno (PJ *P);
-void proj_errno_set (PJ *P, int err);
+int proj_errno_set (PJ *P, int err);
int proj_errno_reset (PJ *P);
-void proj_errno_restore (PJ *P, int err);
+int proj_errno_restore (PJ *P, int err);
PJ_DERIVS proj_derivatives(PJ *P, const LP lp);
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 7196b90c..4f044101 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -520,37 +520,33 @@ int proj_errno (PJ *P) {
}
/*****************************************************************************/
-void proj_errno_set (PJ *P, int err) {
+int proj_errno_set (PJ *P, int err) {
/******************************************************************************
- Sets errno at the context and bubble it up to the thread local errno
+ Set context-errno, bubble it up to the thread local errno, return err
******************************************************************************/
/* Use proj_errno_reset to explicitly clear the error status */
if (0==err)
- return;
+ return 0;
/* For P==0 err goes to the default context */
proj_context_errno_set (pj_get_ctx (P), err);
errno = err;
- return;
+ return err;
}
/*****************************************************************************/
-void proj_errno_restore (PJ *P, int err) {
+int 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.
+ Use proj_errno_restore when the current function succeeds, but the
+ error flag was set on entry, and stored/reset using proj_errno_reset
+ in order to monitor for new errors.
See usage example under proj_errno_reset ()
******************************************************************************/
if (0==err)
- return;
+ return 0;
proj_errno_set (P, err);
+ return 0;
}
/*****************************************************************************/
@@ -562,17 +558,24 @@ int proj_errno_reset (PJ *P) {
Returns the previous value of the errno, for convenient reset/restore
operations:
- void foo (PJ *P) {
+ int foo (PJ *P) {
+ // errno may be set on entry, but we need to reset it to be able to
+ // check for errors from "do_something_with_P(P)"
int last_errno = proj_errno_reset (P);
+ // local failure
+ if (0==P)
+ return proj_errno_set (P, 42);
+
+ // call to function that may fail
do_something_with_P (P);
- // failure - keep latest error status
+ // failure in do_something_with_P? - keep latest error status
if (proj_errno(P))
- return;
- // success - restore previous error status
- proj_errno_restore (P, last_errno);
- return;
+ return proj_errno (P);
+
+ // success - restore previous error status, return 0
+ return proj_errno_restore (P, last_errno);
}
******************************************************************************/
int last_errno;
diff --git a/src/projects.h b/src/projects.h
index 32b74b08..f8f1fa38 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -229,6 +229,11 @@ struct PJconsts {
projCtx_t *ctx;
const char *descr; /* From pj_list.h or individual PJ_*.c file */
paralist *params; /* Parameter list */
+ char *def_size; /* Shape and size parameters extracted from params */
+ char *def_shape;
+ char *def_spherification;
+ char *def_ellps;
+
struct geod_geodesic *geod; /* For geodesic computations */
struct pj_opaque *opaque; /* Projection specific parameters, Defined in PJ_*.c */
int inverted; /* Tell high level API functions to swap inv/fwd */
@@ -413,11 +418,6 @@ struct ARG_list {
typedef union { double f; int i; char *s; } PROJVALUE;
-struct PJ_SELFTEST_LIST {
- char *id; /* projection keyword */
- int (* testfunc)(void); /* projection entry point */
-};
-
struct PJ_ELLPS {
char *id; /* ellipse keyword name */
char *major; /* a= value */
@@ -534,6 +534,7 @@ enum deprecated_constants_for_now_dropped_analytical_factors {
#define PJD_ERR_LAT_0_IS_ZERO -55
#define PJD_ERR_ELLIPSOIDAL_UNSUPPORTED -56
#define PJD_ERR_TOO_MANY_INITS -57
+#define PJD_ERR_INVALID_ARG -58
struct projFileAPI_t;
@@ -676,6 +677,7 @@ double aacos(projCtx,double), aasin(projCtx,double), asqrt(double), aatan2(doubl
PROJVALUE pj_param(projCtx ctx, paralist *, const char *);
paralist *pj_mkparam(char *);
+int pj_ellipsoid (PJ *);
int pj_ell_set(projCtx ctx, paralist *, double *, double *);
int pj_datum_set(projCtx,paralist *, PJ *);
int pj_prime_meridian_set(paralist *, PJ *);
diff --git a/test/gie/ellipsoid.gie b/test/gie/ellipsoid.gie
new file mode 100644
index 00000000..120e407d
--- /dev/null
+++ b/test/gie/ellipsoid.gie
@@ -0,0 +1,121 @@
+===============================================================================
+
+Test pj_ellipsoid, the reimplementation of pj_ell_set
+
+===============================================================================
+
+
+BEGIN
+
+-------------------------------------------------------------------------------
+First a spherical example
+-------------------------------------------------------------------------------
+operation proj=merc R=6400000
+-------------------------------------------------------------------------------
+tolerance 10 nm
+accept 1 2
+expect 111701.0721276371 223447.5262032605
+
+accept 12 55
+expect 1340412.8655316452 7387101.1430967357
+-------------------------------------------------------------------------------
+
+
+-------------------------------------------------------------------------------
+Then an explicitly defined ellipsoidal example
+-------------------------------------------------------------------------------
+operation proj=merc a=6400000 rf=297
+-------------------------------------------------------------------------------
+tolerance 10 nm
+accept 1 2
+expect 111701.0721276371 221945.9681832088
+
+accept 12 55
+expect 1340412.8655316452 7351803.9151705895
+-------------------------------------------------------------------------------
+
+
+-------------------------------------------------------------------------------
+Then try using a built in ellipsoid
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 10 nm
+accept 1 2
+expect 111319.4907932736 221194.0771604237
+
+accept 12 55
+expect 1335833.8895192828 7326837.7148738774
+-------------------------------------------------------------------------------
+
+
+-------------------------------------------------------------------------------
+Then try to fail deliberately
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80000000000
+expect failure
+operation proj=merc +a=-1
+expect failure
+operation proj=merc +no_defs
+expect failure
+operation proj=merc +es=-1 +no_defs
+expect failure
+operation
+expect failure
+operation cobra
+expect failure
+-------------------------------------------------------------------------------
+
+
+-------------------------------------------------------------------------------
+Finally test the spherification functionality
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_A
+tolerance 10 nm
+accept 12 55
+expect 1334340.6237297705 7353636.6296552019
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_V
+tolerance 10 nm
+accept 12 55
+expect 1334339.2852675652 7353629.2533042720
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_a
+tolerance 10 nm
+accept 12 55
+expect 1333594.4904527504 7349524.6413825499
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_g
+tolerance 10 nm
+accept 12 55
+expect 1333592.6102291327 7349514.2793497816
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_h
+tolerance 10 nm
+accept 12 55
+expect 1333590.7300081658 7349503.9173316229
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_lat_a=60
+tolerance 10 nm
+accept 12 55
+expect 1338073.7436268919 7374210.0924803326
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_lat_g=60
+tolerance 10 nm
+accept 12 55
+expect 1338073.2696101593 7374207.4801437631
+-------------------------------------------------------------------------------
+
+
+-------------------------------------------------------------------------------
+This one from testvarious failed at first version of the pull request
+-------------------------------------------------------------------------------
+operation proj=healpix a=1 lon_0=0 ellps=WGS84
+-------------------------------------------------------------------------------
+accept 0 41.937853904844985
+expect 0 0.78452
+accept -90 0
+expect -1.56904 0
+-------------------------------------------------------------------------------
+
+END
diff --git a/travis/install.sh b/travis/install.sh
index 049669d4..f71108c6 100755
--- a/travis/install.sh
+++ b/travis/install.sh
@@ -84,6 +84,7 @@ PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/builtins.gie
PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/more_builtins.gie
PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/deformation.gie
PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/axisswap.gie
+PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/ellipsoid.gie
# install & run the working GIGS test
# create locations that pyproj understands
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 08:15:22 +0000