Merge branch 'master' into kinematic-grids

16 files changed, 171 insertions, 150 deletions

diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index 99e1af7d..67154ab3 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -120,8 +120,13 @@ static double normal_radius_of_curvature (double a, double es, double phi) {
 
 /*********************************************************************/
 static double geocentric_radius (double a, double b, double phi) {
-/*********************************************************************/
-    /* This is from WP2, but uses hypot() for potentially better numerical robustness */
+/*********************************************************************
+    Return the geocentric radius at latitude phi, of an ellipsoid
+    with semimajor axis a and semiminor axis b.
+
+    This is from WP2, but uses hypot() for potentially better
+    numerical robustness
+***********************************************************************/
     return hypot (a*a*cos (phi), b*b*sin(phi)) / hypot (a*cos(phi), b*sin(phi));
 }
 
@@ -129,18 +134,15 @@ static double geocentric_radius (double a, double b, double phi) {
 /*********************************************************************/
 static XYZ cartesian (LPZ geodetic,  PJ *P) {
 /*********************************************************************/
-    double N, h, lam, phi, cosphi = cos(geodetic.phi);
+    double N, cosphi = cos(geodetic.phi);
     XYZ xyz;
 
     N   =  normal_radius_of_curvature(P->a, P->es, geodetic.phi);
-    phi =  geodetic.phi;
-    lam =  geodetic.lam;
-    h   =  geodetic.z;
 
     /* HM formula 5-27 (z formula follows WP) */
-    xyz.x = (N + h) * cosphi * cos(lam);
-    xyz.y = (N + h) * cosphi * sin(lam);
-    xyz.z = (N * (1 - P->es) + h) * sin(phi);
+    xyz.x = (N + geodetic.z) * cosphi      * cos(geodetic.lam);
+    xyz.y = (N + geodetic.z) * cosphi      * sin(geodetic.lam);
+    xyz.z = (N * (1 - P->es) + geodetic.z) * sin(geodetic.phi);
 
     return xyz;
 }
@@ -149,23 +151,19 @@ static XYZ cartesian (LPZ geodetic,  PJ *P) {
 /*********************************************************************/
 static LPZ geodetic (XYZ cartesian,  PJ *P) {
 /*********************************************************************/
-    double N, b, p, theta, c, s, e2s;
+    double N, p, theta, c, s;
     LPZ lpz;
 
     /* Perpendicular distance from point to Z-axis (HM eq. 5-28) */
     p = hypot (cartesian.x, cartesian.y);
 
-    /* Ancillary ellipsoidal parameters */
-    b   =  P->b;
-    e2s =  P->e2s;
-
     /* HM eq. (5-37) */
-    theta  =  atan2 (cartesian.z * P->a,  p * b);
+    theta  =  atan2 (cartesian.z * P->a,  p * P->b);
 
     /* HM eq. (5-36) (from BB, 1976) */
     c  =  cos(theta);
     s  =  sin(theta);
-    lpz.phi  =  atan2 (cartesian.z + e2s*b*s*s*s,  p - P->es*P->a*c*c*c);
+    lpz.phi  =  atan2 (cartesian.z + P->e2s*P->b*s*s*s,  p - P->es*P->a*c*c*c);
     lpz.lam  =  atan2 (cartesian.y, cartesian.x);
     N        =  normal_radius_of_curvature (P->a, P->es, lpz.phi);
 
@@ -176,7 +174,7 @@ static LPZ geodetic (XYZ cartesian,  PJ *P) {
         /* by computing the height as the cartesian z value      */
         /* minus the geocentric radius of the Earth at the given */
         /* latitude                                              */
-        double r = geocentric_radius (P->a, b, lpz.phi);
+        double r = geocentric_radius (P->a, P->b, lpz.phi);
         lpz.z = fabs (cartesian.z) - r;
     }
     else
diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c
index 503dc392..c9c30401 100644
--- a/src/PJ_helmert.c
+++ b/src/PJ_helmert.c
@@ -92,8 +92,9 @@ struct pj_opaque_helmert {
 #define R21 (Q->R[2][1])
 #define R22 (Q->R[2][2])
 
+/**************************************************************************/
 static void update_parameters(PJ *P) {
-    /**************************************************************************
+/***************************************************************************
 
     Update transformation parameters.
     ---------------------------------
@@ -114,7 +115,8 @@ static void update_parameters(PJ *P) {
     of that parameter.
 
     [0] http://itrf.ign.fr/doc_ITRF/Transfo-ITRF2008_ITRFs.txt
-    **************************************************************************/
+
+*******************************************************************************/
 
     struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque;
     double dt = Q->t_obs - Q->epoch;
@@ -146,8 +148,9 @@ static void update_parameters(PJ *P) {
     return;
 }
 
+/**************************************************************************/
 static void build_rot_matrix(PJ *P) {
-    /**************************************************************************
+/***************************************************************************
 
     Build rotation matrix.
     ----------------------
@@ -208,7 +211,7 @@ static void build_rot_matrix(PJ *P) {
     as published, the "transpose" option provides the ability to switch
     between the conventions.
 
-    ***************************************************************************/
+***************************************************************************/
     struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque;
 
     double  f,  t,  p;    /* phi/fi , theta, psi  */
diff --git a/src/PJ_horner.c b/src/PJ_horner.c
index d6d2c51c..a9a2c922 100644
--- a/src/PJ_horner.c
+++ b/src/PJ_horner.c
@@ -1,19 +1,3 @@
-#define PJ_LIB__
-#include "proj_internal.h"
-#include "projects.h"
-#include <assert.h>
-#include <stddef.h>
-#include <math.h>
-#include <errno.h>
-
-PROJ_HEAD(horner,    "Horner polynomial evaluation");
-
-/* make horner.h interface with proj's memory management */
-#define horner_dealloc(x) pj_dealloc(x)
-#define horner_calloc(n,x) pj_calloc(n,x)
-
-/* The next few hundred lines is mostly cut-and-paste from the horner.h header library */
-
 /***********************************************************************
 
     Interfacing to a classic piece of geodetic software
@@ -91,6 +75,20 @@ PROJ_HEAD(horner,    "Horner polynomial evaluation");
  *
  *****************************************************************************/
 
+#define PJ_LIB__
+#include "proj_internal.h"
+#include "projects.h"
+#include <assert.h>
+#include <stddef.h>
+#include <math.h>
+#include <errno.h>
+
+PROJ_HEAD(horner,    "Horner polynomial evaluation");
+
+/* make horner.h interface with proj's memory management */
+#define horner_dealloc(x) pj_dealloc(x)
+#define horner_calloc(n,x) pj_calloc(n,x)
+
 
 struct horner;
 typedef struct horner HORNER;
@@ -399,7 +397,7 @@ static int parse_coefs (PJ *P, double *coefs, char *param, int ncoefs) {
 
     buf = pj_calloc (strlen (param) + 2, sizeof(char));
     if (0==buf) {
-        proj_log_error (P, "Horner: Out of core");
+        proj_log_error (P, "Horner: No memory left");
         return 0;
     }
 
diff --git a/src/PJ_molodensky.c b/src/PJ_molodensky.c
index f09f07cd..50423fd8 100644
--- a/src/PJ_molodensky.c
+++ b/src/PJ_molodensky.c
@@ -90,6 +90,13 @@ static double RM (double a, double es, double phi) {
 
         E.J Krakiwsky & D.B. Thomson, 1974,
         GEODETIC POSITION COMPUTATIONS,
+
+        Fredericton NB, Canada:
+        University of New Brunswick,
+        Department of Geodesy and Geomatics Engineering,
+        Lecture Notes No. 39,
+        99 pp.
+
         http://www2.unb.ca/gge/Pubs/LN39.pdf
 
 **********************************************************/
diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c
index 137fdec8..f1411cb4 100644
--- a/src/PJ_pipeline.c
+++ b/src/PJ_pipeline.c
@@ -45,9 +45,10 @@
 
     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.
+    addition to the PROJ.4 code base: The implementation is by and large completed
+    by adding an extra projection called "pipeline" (i.e. this file), which
+    handles all business, and a small amount of added functionality in the
+    pj_init code, implementing support for multilevel, embedded pipelines.
 
     Syntactically, the pipeline system introduces the "+step" keyword (which
     indicates the start of each transformation step), the "+omit_fwd" and
@@ -63,10 +64,11 @@
 
     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,
+    Following geodetic conventions, the rotations are given in arcseconds,
     and the scale factor is given as parts-per-million.
 
-    [1] B. W. Kernighan & P. J. Plauger: Software tools. Addison-Wesley, 1976, 338 pp.
+    [1] B. W. Kernighan & P. J. Plauger: Software tools.
+        Reading, Massachusetts, Addison-Wesley, 1976, 338 pp.
 
 ********************************************************************************
 
@@ -156,10 +158,11 @@ static LP     pipeline_reverse (XY xyz, PJ *P);
 
     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.
+    For classic proj style projections, P->left has the value
+    PJ_IO_UNITS_RADIANS, while P->right has the value
+    PJ_IO_UNITS_CLASSIC, indicating that the forward driver expects
+    angular input coordinates, and provides linear output coordinates,
+    scaled by the P->a semimajor axis length.
 
     Newer projections may set P->left and P->right to either
     PJ_IO_UNITS_METERS, PJ_IO_UNITS_RADIANS or PJ_IO_UNITS_ANY,
@@ -377,7 +380,7 @@ PJ *PROJECTION(pipeline) {
 
         if (0==strcmp ("proj=pipeline", argv[i])) {
             if (-1 != i_pipeline) {
-                proj_log_error (P, "Pipeline: Nesting invalid");
+                proj_log_error (P, "Pipeline: Nesting only allowed when child pipelines are wrapped in +init's");
                 return destructor (P, PJD_ERR_MALFORMED_PIPELINE); /* ERROR: nested pipelines */
             }
             i_pipeline = i;
diff --git a/src/cct.c b/src/cct.c
index 72d8c6cc..90b84159 100644
--- a/src/cct.c
+++ b/src/cct.c
@@ -44,7 +44,7 @@ Hence, in honour of cct (the geodesist) this is cct (the program).
 
 ************************************************************************
 
-Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-09-19
+Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-10-26
 
 ************************************************************************
 
@@ -73,7 +73,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-09-19
 
 #include "optargpm.h"
 #include <proj.h>
-#include <projects.h>
+#include "projects.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <ctype.h>
@@ -85,8 +85,6 @@ double proj_atof(const char *str);
 
 char *column (char *buf, int n);
 PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double fixed_time);
-int print_output_line (FILE *fout, char *buf, PJ_COORD point);
-int main(int argc, char **argv);
 
 
 
@@ -237,16 +235,31 @@ int main(int argc, char **argv) {
     }
 
 
-    /* Loop over all lines of all input files */
+    /* Loop over all records of all input files */
     while (opt_input_loop (o, optargs_file_format_text)) {
         void *ret = fgets (buf, 10000, o->input);
-        int res;
         opt_eof_handler (o);
         if (0==ret) {
             fprintf (stderr, "Read error in record %d\n", (int) o->record_index);
             continue;
         }
         point = parse_input_line (buf, columns_xyzt, fixed_z, fixed_time);
+        if (HUGE_VAL==point.xyzt.x) {
+            char *c = column (buf, 1);
+
+            /* if it's a comment or blank line, we reflect it */
+            if (c && ((*c=='\0') || (*c=='#'))) {
+                fprintf (fout, "%s\n", buf);
+                continue;
+            }
+
+            /* otherwise, it must be a syntax error */
+            fprintf (fout, "# Record %d UNREADABLE: %s", (int) o->record_index, buf);
+            if (verbose)
+                fprintf (stderr, "%s: Could not parse file '%s' line %d\n", o->progname, opt_filename (o), opt_record (o));
+            continue;
+        }
+
         if (PJ_IO_UNITS_RADIANS==input_unit) {
             point.lpzt.lam = proj_torad (point.lpzt.lam);
             point.lpzt.phi = proj_torad (point.lpzt.phi);
@@ -256,13 +269,17 @@ int main(int argc, char **argv) {
             point.lpzt.lam = proj_todeg (point.lpzt.lam);
             point.lpzt.phi = proj_todeg (point.lpzt.phi);
         }
-        res = print_output_line (fout, buf, point);
-        if (0==res) {
-            fprintf (fout, "# UNREADABLE: %s", buf);
-            if (verbose)
-                fprintf (stderr, "%s: Could not parse file '%s' line %d\n", o->progname, opt_filename (o), opt_record (o));
+
+        if (HUGE_VAL==point.xyzt.x) {
+            /* transformation error (TODO provide existing internal errmsg here) */
+            fprintf (fout, "# Record %d TRANSFORMATION ERROR: %s", (int) o->record_index, buf);
+            continue;
         }
+
+        /* Time to print the result */
+        fprintf (fout, "%20.15f  %20.15f  %20.15f  %20.15f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t);
     }
+
     if (stdout != fout)
         fclose (fout);
     free (o);
@@ -322,15 +339,3 @@ PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double
     errno = prev_errno;
     return result;
 }
-
-
-int print_output_line (FILE *fout, char *buf, PJ_COORD point) {
-    char *c;
-    if (HUGE_VAL!=point.xyzt.x)
-        return fprintf (fout, "%20.15f  %20.15f  %20.15f  %20.15f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t);
-    c = column (buf, 1);
-    /* reflect comments and blanks */
-    if (c && ((*c=='\0') || (*c=='#')))
-        return fprintf (fout, "%s\n", buf);
-    return 0;
-}
diff --git a/src/gie.c b/src/gie.c
index f0850b61..269a7239 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -52,11 +52,12 @@ but all of them *just special functions*, and not particularly more
 special than the sin(), cos(), tan(), and hypot() already available in
 the C standard library.
 
-And hence, *they should not be particularly much harder to use*, for a
-programmer, than the sin()s, tan()s and hypot()s so readily available.
+And hence, according to Gerald, *they should not be particularly much
+harder to use*, for a programmer, than the sin()s, tan()s and hypot()s
+so readily available.
 
 Gerald's ingenuity also showed in the implementation of the vision,
-where he devised a highly comprehensible, yet simple, system of key-value
+where he devised a comprehensive, yet simple, system of key-value
 pairs for parameterising a map projection, and the highly flexible
 PJ struct, storing run-time compiled versions of those key-value pairs,
 hence making a map projection function call, pj_fwd(PJ, point), as easy
@@ -332,7 +333,7 @@ char *column (char *buf, int n) {
 }
 
 
-
+/* interpret <args> as a numeric followed by a linear decadal prefix - return the properly scaled numeric */
 static double strtod_scaled (char *args, double default_scale) {
     double s;
     char *endp = args;
diff --git a/src/makefile.vc b/src/makefile.vc
index 15fd6428..edb683b0 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -64,6 +64,7 @@ pipeline = \
 	PJ_vgridshift.obj PJ_hgridshift.obj PJ_unitconvert.obj PJ_molodensky.obj \
 	PJ_deformation.obj
 
+
 geodesic = geodesic.obj
 
 LIBOBJ	=	$(support) $(pseudo) $(azimuthal) $(conic) $(cylinder) $(misc) \
diff --git a/src/optargpm.h b/src/optargpm.h
index c47dfec1..23e375f8 100644
--- a/src/optargpm.h
+++ b/src/optargpm.h
@@ -10,8 +10,8 @@
 For PROJ.4 command line programs, we have a somewhat complex option
 decoding situation, since we have to navigate in a cocktail of classic
 single letter style options, prefixed by "-", GNU style long options
-prefixwd by "--", transformation specification elements prefixed by "+",
-and input file names prefixed by "" nothing.
+prefixed by "--", transformation specification elements prefixed by "+",
+and input file names prefixed by "" (i.e. nothing).
 
 Hence, classic getopt.h style decoding does not cut the mustard, so
 this is an attempt to catch up and chop the ketchup.
@@ -48,8 +48,8 @@ Operator specs are      +proj=utm +zone=32 +ellps=GRS80
 Operands are             bar baz
 
 
-While claiming neither to save the world, nor to hint at the "shape of
-jazz to come",  at least optargpm has shown useful in constructing cs2cs
+While neither claiming to save the world, nor to hint at the "shape of
+jazz to come", at least optargpm has shown useful in constructing cs2cs
 style transformation filters.
 
 Supporting a wide range of option syntax, the getoptpm API is somewhat
@@ -424,14 +424,6 @@ OPTARGS *opt_parse (int argc, char **argv, const char *flags, const char *keys,
     o->argc = argc;
     o->argv = argv;
     o->progname = opt_strip_path (argv[0]);
-/*    o->progname = argv[0];
-    last_path_delim = strrchr (argv[0], '\\');
-    if (last_path_delim > o->progname)
-        o->progname = last_path_delim;
-    last_path_delim = strrchr (argv[0], '/');
-    if (last_path_delim > o->progname)
-        o->progname = last_path_delim;
-*/
 
     /* Reset all flags */
     for (i = 0;  i < (int) strlen (flags);  i++)
diff --git a/src/pj_ellps.c b/src/pj_ellps.c
index 32fbbd42..2a2c5d3b 100644
--- a/src/pj_ellps.c
+++ b/src/pj_ellps.c
@@ -15,6 +15,7 @@ pj_ellps[] = {
 {"andrae",	"a=6377104.43",		"rf=300.0",		"Andrae 1876 (Den., Iclnd.)"},
 {"aust_SA",	"a=6378160.0",		"rf=298.25",		"Australian Natl & S. Amer. 1969"},
 {"GRS67",	"a=6378160.0",		"rf=298.2471674270",	"GRS 67(IUGG 1967)"},
+{"GSK2011",	"a=6378136.5",		"rf=298.2564151",	"GSK-2011"},
 {"bessel",	"a=6377397.155",	"rf=299.1528128",	"Bessel 1841"},
 {"bess_nam",	"a=6377483.865",	"rf=299.1528128",	"Bessel 1841 (Namibia)"},
 {"clrk66",	"a=6378206.4",		"b=6356583.8",		"Clarke 1866"},
@@ -40,6 +41,7 @@ pj_ellps[] = {
 {"mprts",	"a=6397300.",		"rf=191.",		"Maupertius 1738"},
 {"new_intl",	"a=6378157.5",		"b=6356772.2",		"New International 1967"},
 {"plessis",	"a=6376523.",		"b=6355863.", 		"Plessis 1817 (France)"},
+{"PZ90",	"a=6378136.0",		"rf=298.25784",		"PZ-90"},
 {"SEasia",	"a=6378155.0",		"b=6356773.3205",	"Southeast Asia"},
 {"walbeck",	"a=6376896.0",		"b=6355834.8467",	"Walbeck"},
 {"WGS60",	"a=6378165.0",		"rf=298.3",		"WGS 60"},
diff --git a/src/pj_init.c b/src/pj_init.c
index 704a8b55..2a945397 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -381,8 +381,9 @@ pj_init_plus_ctx( projCtx ctx, const char *definition )
 
                 if( argc+1 == MAX_ARG )
                 {
-                    pj_ctx_set_errno( ctx, -44 );
-                    goto bum_call;
+                    pj_dalloc( defn_copy );
+                    pj_ctx_set_errno( ctx, PJD_ERR_UNPARSEABLE_CS_DEF );
+                    return 0;
                 }
 
                 argv[argc++] = defn_copy + i + 1;
@@ -410,9 +411,7 @@ pj_init_plus_ctx( projCtx ctx, const char *definition )
     /* perform actual initialization */
     result = pj_init_ctx( ctx, argc, argv );
 
-bum_call:
     pj_dalloc( defn_copy );
-
     return result;
 }
 
@@ -441,6 +440,8 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
     int i;
     int found_def = 0;
     PJ *PIN = 0;
+    int n_pipelines = 0;
+    int n_inits = 0;
 
     if (0==ctx)
         ctx = pj_get_default_ctx ();
@@ -452,36 +453,45 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
         pj_ctx_set_errno (ctx, PJD_ERR_NO_ARGS);
         return 0;
     }
-    
+
+    /* count occurrences of pipelines and inits */
+    for (i = 0; i < argc; ++i) {
+        if (!strcmp (argv[i], "+proj=pipeline") || !strcmp(argv[i], "proj=pipeline") )
+                n_pipelines++;
+        if (!strncmp (argv[i], "+init=", 6) || !strncmp(argv[i], "init=", 5))
+            n_inits++;
+    }
+
+    /* can't have nested pipeline directly */
+    if (n_pipelines > 1) {
+        pj_ctx_set_errno (ctx, PJD_ERR_MALFORMED_PIPELINE);
+        return 0;
+    }
+
+    /* don't allow more than one +init in non-pipeline operations */
+    if (n_pipelines == 0 && n_inits > 1) {
+        pj_ctx_set_errno (ctx, PJD_ERR_TOO_MANY_INITS);
+        return 0;
+    }
+
     /* put arguments into internal linked list */
     start = curr = pj_mkparam(argv[0]);
     if (!curr)
         return pj_dealloc_params (ctx, start, ENOMEM);
 
-    /* build parameter list and expand +init's. Does not take care of a single +init. */
     for (i = 1; i < argc; ++i) {
         curr->next = pj_mkparam(argv[i]);
         if (!curr->next)
             return pj_dealloc_params (ctx, start, ENOMEM);
-
-        /* check if +init present */
-        if (pj_param(ctx, curr, "tinit").i) {
-            found_def = 0;
-            curr = get_init(ctx, &curr, curr->next, pj_param(ctx, curr, "sinit").s, &found_def);
-            if (!curr)
-                return pj_dealloc_params (ctx, start, PJD_ERR_NO_ARGS);
-
-            if (!found_def)
-                return pj_dealloc_params (ctx, start, PJD_ERR_NO_OPTION_IN_INIT_FILE);
-                
-        } else {
-            curr = curr->next;
-        }
+        curr = curr->next;
     }
 
-    /* in case the parameter list only consist of a +init parameter
-       it is expanded here (will not be handled in the above loop). */
-    if (pj_param(ctx, start, "tinit").i && argc == 1) {
+
+    /* 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     */
+    /* problem when +inits are expanded as late as possible.                     */
+    if (pj_param(ctx, start, "tinit").i && n_pipelines == 0) {
         found_def = 0;
         curr = get_init(ctx, &start, curr, pj_param(ctx, start, "sinit").s, &found_def);
         if (!curr)
@@ -489,10 +499,10 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
         if (!found_def)
             return pj_dealloc_params (ctx, start, PJD_ERR_NO_OPTION_IN_INIT_FILE);
     }
-    
+
     if (ctx->last_errno)
         return pj_dealloc_params (ctx, start, ctx->last_errno);
-        
+
     /* find projection selection */
     if (!(name = pj_param(ctx, start, "sproj").s))
         return pj_default_destructor (PIN, PJD_ERR_PROJ_NOT_NAMED);
@@ -500,10 +510,11 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
 
     if (!s)
         return pj_dealloc_params (ctx, start, PJD_ERR_UNKNOWN_PROJECTION_ID);
-        
-    /* set defaults, unless inhibited */
-    if (!(pj_param(ctx, start, "bno_defs").i))
-        curr = get_defaults(ctx,&start, curr, name);
+
+    /* set defaults, unless inhibited or we are initializing a pipeline */
+    if (!(pj_param(ctx, start, "bno_defs").i) && n_pipelines == 0)
+            curr = get_defaults(ctx,&start, curr, name);
+
     proj = (PJ *(*)(PJ *)) pj_list[i].proj;
 
     /* allocate projection structure */
@@ -511,6 +522,7 @@ 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;
@@ -526,13 +538,21 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
     PIN->vgridlist_geoid_count = 0;
 
     /* set datum parameters */
-    if (pj_datum_set(ctx, start, PIN)) 
+    if (pj_datum_set(ctx, start, PIN))
         return pj_default_destructor (PIN, PJD_ERR_MISSING_ARGS);
 
-    /* set ellipsoid/sphere parameters */
-    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);
+    /* 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 (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;
@@ -606,7 +626,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
         if( !(fabs(PIN->long_wrap_center) < 10 * M_TWOPI) )
             return pj_default_destructor (PIN, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
     }
-    
+
     /* axis orientation */
     if( (pj_param(ctx, start,"saxis").s) != NULL )
     {
@@ -619,7 +639,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
             || strchr( axis_legal, axis_arg[1] ) == NULL
             || strchr( axis_legal, axis_arg[2] ) == NULL)
             return pj_default_destructor (PIN, PJD_ERR_AXIS);
-            
+
         /* it would be nice to validate we don't have on axis repeated */
         strcpy( PIN->axis, axis_arg );
     }
@@ -643,7 +663,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
         PIN->k0 = pj_param(ctx, start, "dk").f;
     else
         PIN->k0 = 1.;
-    if (PIN->k0 <= 0.) 
+    if (PIN->k0 <= 0.)
         return pj_default_destructor (PIN, PJD_ERR_K_LESS_THAN_ZERO);
 
     /* set units */
@@ -704,7 +724,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
             && *next_str == '\0' )
             value = name;
 
-        if (!value) 
+        if (!value)
             return pj_default_destructor (PIN, PJD_ERR_UNKNOWN_PRIME_MERIDIAN);
         PIN->from_greenwich = dmstor_ctx(ctx,value,NULL);
     }
diff --git a/src/pj_strerrno.c b/src/pj_strerrno.c
index 89a5a025..fe137b87 100644
--- a/src/pj_strerrno.c
+++ b/src/pj_strerrno.c
@@ -62,7 +62,8 @@ pj_err_list[] = {
     "missing required arguments",                                      /* -54 */
     "lat_0 = 0",                                                       /* -55 */
     "ellipsoidal usage unsupported",                                   /* -56 */
-    
+    "only one +init allowed for non-pipeline operations",              /* -57 */
+
     /* When adding error messages, remember to update ID defines in
        projects.h, and transient_error array in pj_transform                  */
 };
diff --git a/src/proj.h b/src/proj.h
index ed885091..dd0672f9 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -53,26 +53,21 @@
  *           it remains as-is for now).
  *
  *           THIRD, I try to eliminate implicit type punning. Hence this API
- *           introduces the PJ_OBS ("observation") data type, for generic
- *           coordinate and handling of ancillary data.
+ *           introduces the PJ_COORD union data type, for generic 4D coordinate
+ *           handling.
  *
- *           It includes the PJ_COORD and PJ_TRIPLET unions making it possible
- *           to make explicit the previously used "implicit type punning", where
- *           a XY is turned into a LP by re#defining both as UV, behind the back
- *           of the user.
+ *           PJ_COORD makes it possible to make explicit the previously used
+ *           "implicit type punning", where a XY is turned into a LP by
+ *           re#defining both as UV, behind the back of the user.
  *
  *           The PJ_COORD union is used for storing 1D, 2D, 3D and 4D coordinates.
- *           The PJ_TRIPLET union is used for storing any set of up to 3 related
- *           observations. At the application code level, the names of these
- *           unions will usually not be used - they will only be accessed via
- *           their tag names in the PJ_OBS data type.
  *
  *           The bare essentials API presented here follows the PROJ.4
  *           convention of sailing the coordinate to be reprojected, up on
  *           the stack ("call by value"), and symmetrically returning the
- *           result on the stack. Although the PJ_OBS object is 4 times
- *           as large as the traditional XY and LP objects, timing results
- *           have shown the overhead to be very reasonable.
+ *           result on the stack. Although the PJ_COORD object is twice as large
+ *           as the traditional XY and LP objects, timing results have shown the
+ *           overhead to be very reasonable.
  *
  *           Contexts and thread safety
  *           --------------------------
@@ -82,8 +77,8 @@
  *           context subsystem is unavoidable in a multi-threaded world.
  *           Hence, instead of hiding it away, we move it into the limelight,
  *           highly recommending (but not formally requiring) the bracketing
- *           with calls to proj_context_create(...)/proj_context_destroy() of
- *           any code block calling PROJ.4 functions.
+ *           of any code block calling PROJ.4 functions with calls to
+ *           proj_context_create(...)/proj_context_destroy()
  *
  *           Legacy single threaded code need not do anything, but *may*
  *           implement a bit of future compatibility by using the backward
@@ -426,7 +421,7 @@ const PJ_ELLPS            *proj_list_ellps(void);
 const PJ_UNITS            *proj_list_units(void);
 const PJ_PRIME_MERIDIANS  *proj_list_prime_meridians(void);
 
-/* These are trivial, and while occasionaly useful in real code, primarily here to       */
+/* These are trivial, and while occasionally useful in real code, primarily here to      */
 /* simplify demo code, and in acknowledgement of the proj-internal discrepancy between   */
 /* angular units expected by classical proj, and by Charles Karney's geodesics subsystem */
 double proj_torad (double angle_in_degrees);
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index b1aa3883..5484b650 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -3,12 +3,6 @@
  * Purpose:  Implement a (currently minimalistic) proj API based primarily
  *           on the PJ_COORD 4D geodetic spatiotemporal data type.
  *
- *           proj thread contexts have not seen widespread use, so one of the
- *           intentions with this new API is to make them less visible on the
- *           API surface: Contexts do not have a life by themselves, they are
- *           visible only through their associated PJs, and the number of
- *           functions supporting them is limited.
- *
  * Author:   Thomas Knudsen,  thokn@sdfe.dk,  2016-06-09/2016-11-06
  *
  ******************************************************************************
diff --git a/src/proj_strtod.c b/src/proj_strtod.c
index 757dfaf6..fa683465 100644
--- a/src/proj_strtod.c
+++ b/src/proj_strtod.c
@@ -17,7 +17,7 @@ In the early versions of proj, iirc, a gnu version of strtod was used,
 mostly to work around cases where the same system library was used for
 C and Fortran linking, hence making strtod accept "D" and "d" as
 exponentiation indicators, following Fortran Double Precision constant
-syntax. This broke the proj angular syntax accepting a "d" to mean
+syntax. This broke the proj angular syntax, accepting a "d" to mean
 "degree": 12d34'56", meaning 12 degrees 34 minutes and 56 seconds.
 
 With an explicit MIT licence, PROJ.4 could not include GPL code any
diff --git a/src/projects.h b/src/projects.h
index b9d88cf3..6ea0c914 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -186,7 +186,7 @@ struct PJ_REGION_S;
 typedef struct PJ_REGION_S  PJ_Region;
 typedef struct ARG_list paralist;   /* parameter list */
 enum pj_io_units {
-    PJ_IO_UNITS_CLASSIC = 0,   /* LEFT: Radians     RIGHT: Scaled meters */
+    PJ_IO_UNITS_CLASSIC = 0,   /* Scaled meters (right) */
     PJ_IO_UNITS_METERS  = 1,   /* Meters  */
     PJ_IO_UNITS_RADIANS = 2    /* Radians */
 };
@@ -532,6 +532,7 @@ struct FACTORS {
 #define PJD_ERR_MISSING_ARGS            -54
 #define PJD_ERR_LAT_0_IS_ZERO           -55
 #define PJD_ERR_ELLIPSOIDAL_UNSUPPORTED -56
+#define PJD_ERR_TOO_MANY_INITS          -57
 
 struct projFileAPI_t;