Move pipeline initialization logic to PJ_pipeline.c and decrease the number of special cases to handle in pj_init.c

4 files changed, 97 insertions, 59 deletions

diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c
index c89917a9..0729b1da 100644
--- a/src/PJ_pipeline.c
+++ b/src/PJ_pipeline.c
@@ -98,8 +98,9 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-20
 ********************************************************************************/
 
 #define PJ_LIB__
+#include <geodesic.h>
 #include "proj_internal.h"
-#include <projects.h>
+#include "projects.h"
 
 #include <assert.h>
 #include <stddef.h>
@@ -291,9 +292,45 @@ static char **argv_params (paralist *params, size_t argc) {
     return argv;
 }
 
+/* Being the special operator that the pipeline is, we have to handle the    */
+/* ellipsoid differently than usual. In general, the pipeline operation does */
+/* not need an ellipsoid, but in some cases it is beneficial nonetheless.    */
+/* Unfortunately we can't use the normal ellipsoid setter in pj_init, since  */
+/* it adds a +ellps parameter to the global args if nothing else is specified*/
+/* This is problematic since that ellipsoid spec is then passed on to the    */
+/* pipeline children. This is rarely what we want, so here we implement our  */
+/* own logic instead. If an ellipsoid is set in the global args, it is used  */
+/* as the pipeline ellipsoid. Otherwise we use WGS84 parameters as default.  */
+/* At last we calculate the rest of the ellipsoid parameters and             */
+/* re-initialize P->geod.                                                    */
+static void set_ellipsoid(PJ *P) {
+    paralist *cur, *attachment;
+
+    /* Break the linked list after the global args */
+    attachment = 0;
+    for (cur = P->params; cur != 0; cur = cur->next)
+        if (strcmp("step", cur->next->param) == 0) {
+            attachment = cur->next;
+            cur->next = 0;
+            break;
+        }
+
+    /* 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)) {
+        P->a  = 6378137.0;
+        P->es = .00669438002290341575;
+    }
 
+    pj_calc_ellps_params(P, P->a, P->es);
+    geod_init(P->geod, P->a,  (1 - sqrt (1 - P->es)));
 
-PJ *PROJECTION(pipeline) {
+    /* Re-attach the dangling list */
+    cur->next = attachment;
+}
+
+
+PJ *OPERATION(pipeline,0) {
     int i, nsteps = 0, argc;
     int i_pipeline = -1, i_first_step = -1, i_current_step;
     char **argv, **current_argv;
@@ -353,6 +390,8 @@ PJ *PROJECTION(pipeline) {
     if (0==pj_create_pipeline (P, nsteps))
         return destructor (P, ENOMEM);
 
+    set_ellipsoid(P);
+
     /* 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;
diff --git a/src/pj_ell_set.c b/src/pj_ell_set.c
index d43feecb..8a9b39a1 100644
--- a/src/pj_ell_set.c
+++ b/src/pj_ell_set.c
@@ -44,6 +44,51 @@ 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) {
+
+    P->a = a;
+    P->es = es;
+
+    /* Compute some ancillary ellipsoidal parameters */
+    P->e = sqrt(P->es);      /* eccentricity */
+    P->alpha = asin (P->e);  /* angular eccentricity */
+
+    /* second eccentricity */
+    P->e2  = tan (P->alpha);
+    P->e2s = P->e2 * P->e2;
+
+    /* third eccentricity */
+    P->e3    = (0!=P->alpha)? sin (P->alpha) / sqrt(2 - sin (P->alpha)*sin (P->alpha)): 0;
+    P->e3s = P->e3 * P->e3;
+
+    /* flattening */
+    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 */
+    P->f2  = (cos(P->alpha)!=0)? 1/cos (P->alpha) - 1: 0;
+    P->rf2 = P->f2 != 0.0 ? 1/P->f2: HUGE_VAL;
+
+    /* third flattening */
+    P->n  = pow (tan (P->alpha/2), 2);
+    P->rn = P->n != 0.0 ? 1/P->n: HUGE_VAL;
+
+    /* ...and a few more */
+    P->b  = (1 - P->f)*P->a;
+    P->rb = 1. / P->b;
+    P->ra = 1. / P->a;
+
+    P->one_es = 1. - P->es;
+    if (P->one_es == 0.) {
+        pj_ctx_set_errno( P->ctx, PJD_ERR_ECCENTRICITY_IS_ONE);
+        return PJD_ERR_ECCENTRICITY_IS_ONE;
+    }
+
+    P->rone_es = 1./P->one_es;
+
+    return 0;
+}
+
 /* initialize geographic shape parameters */
 int pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
     int i;
diff --git a/src/pj_init.c b/src/pj_init.c
index 2a945397..f69d9eae 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -28,12 +28,13 @@
  *****************************************************************************/
 
 #define PJ_LIB__
-#include <projects.h>
 #include <geodesic.h>
 #include <stdio.h>
 #include <string.h>
 #include <errno.h>
 #include <ctype.h>
+#include "proj_internal.h"
+#include "projects.h"
 
 /* Maximum size of files using the "escape carriage return" feature */
 #define MAX_CR_ESCAPE 65537
@@ -129,7 +130,6 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
         int *found_def) {
     pj_read_state *state = (pj_read_state*) calloc(1,sizeof(pj_read_state));
     char sword[MAX_CR_ESCAPE];
-    char *pipeline;
     int len;
     int in_target = 0;
     const char *next_char = NULL;
@@ -220,10 +220,7 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
             sword[word_len+1] = '\0';
 
             /* do not override existing parameter value of same name - unless in pipeline definition */
-            pipeline = pj_param(ctx, *start, "sproj").s;
-            if (pipeline && strcmp (pipeline, "pipeline"))
-                pipeline = 0;
-            if (pipeline || !pj_param(ctx, *start, sword).i) {
+            if (!pj_param(ctx, *start, sword).i) {
                 /* don't default ellipse if datum, ellps or any earth model
                    information is set. */
                 if( strncmp(sword+1,"ellps=",6) != 0
@@ -260,8 +257,7 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
 /************************************************************************/
 /*                            get_defaults()                            */
 /************************************************************************/
-static paralist *
-get_defaults(projCtx ctx, paralist **start, paralist *next, char *name) {
+static paralist *get_defaults(projCtx ctx, paralist **start, paralist *next, char *name) {
     PAFile fid;
 
     if ( (fid = pj_open_lib(ctx,"proj_def.dat", "rt")) != NULL) {
@@ -486,7 +482,6 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
         curr = curr->next;
     }
 
-
     /* Only expand +init's in non-pipeline operations. +init's in pipelines are  */
     /* expanded in the individual pipeline steps during pipeline initialization. */
     /* Potentially this leads to many nested pipelines, which shouldn't be a     */
@@ -522,7 +517,6 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
     if (0==PIN)
         return pj_dealloc_params (ctx, start, ENOMEM);
 
-
     PIN->ctx = ctx;
     PIN->params = start;
     PIN->is_latlong = 0;
@@ -541,58 +535,17 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
     if (pj_datum_set(ctx, start, PIN))
         return pj_default_destructor (PIN, PJD_ERR_MISSING_ARGS);
 
-    /* set ellipsoid/sphere parameters. If we are initializing a pipeline we */
-    /* override ellipsoid-setter, since it also adds a "+ellps=WGS84" to the */
-    /* parameter list which creates problems when the pipeline driver passes */
-    /* the global parameters on the it's children.                           */
-    if (n_pipelines > 0) {
-        PIN->a = 6378137.0;
-        PIN->es = .00669438002290341575;
-    } else {
+    if (PIN->need_ellps) {
         if (pj_ell_set(ctx, start, &PIN->a, &PIN->es)) {
             pj_log (ctx, PJ_LOG_DEBUG_MINOR, "pj_init_ctx: Must specify ellipsoid or sphere");
             return pj_default_destructor (PIN, PJD_ERR_MISSING_ARGS);
         }
-    }
-
-    PIN->a_orig = PIN->a;
-    PIN->es_orig = 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    = (0!=PIN->alpha)? sin (PIN->alpha) / sqrt(2 - sin (PIN->alpha)*sin (PIN->alpha)): 0;
-    PIN->e3s = PIN->e3 * PIN->e3;
-
-    /* flattening */
-    PIN->f  = 1 - cos (PIN->alpha);   /* = 1 - sqrt (1 - PIN->es); */
-    PIN->rf = PIN->f != 0.0 ? 1.0/PIN->f: HUGE_VAL;
-
-    /* second flattening */
-    PIN->f2  = (cos(PIN->alpha)!=0)? 1/cos (PIN->alpha) - 1: 0;
-    PIN->rf2 = PIN->f2 != 0.0 ? 1/PIN->f2: HUGE_VAL;
-
-    /* third flattening */
-    PIN->n  = pow (tan (PIN->alpha/2), 2);
-    PIN->rn = PIN->n != 0.0 ? 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.)
-        return pj_default_destructor (PIN, PJD_ERR_ECCENTRICITY_IS_ONE);
-    PIN->rone_es = 1./PIN->one_es;
-
 
+        PIN->a_orig = PIN->a;
+        PIN->es_orig = PIN->es;
+        if (pj_calc_ellps_params(PIN, PIN->a, PIN->es))
+            return pj_default_destructor (PIN, PJD_ERR_ECCENTRICITY_IS_ONE);
+    }
 
     /* Now that we have ellipse information check for WGS84 datum */
     if( PIN->datum_type == PJD_3PARAM
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 33cfbdfd..2d6a2406 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -99,6 +99,7 @@ void proj_log_trace (PJ *P, const char *fmt, ...);
 void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION log);
 
 void pj_inherit_ellipsoid_defs(const PJ *src, PJ *dst);
+int pj_calc_ellps_params(PJ *P, double a, double es);
 
 /* Lowest level: Minimum support for fileapi */
 void proj_fileapi_set (PJ *P, void *fileapi);