Merge pull request #731 from busstoptaktik/4D-API_cs2cs-style

[WIP] 4D API cs2cs style

31 files changed, 838 insertions, 358 deletions

diff --git a/src/PJ_axisswap.c b/src/PJ_axisswap.c
index d57834fc..44446d9c 100644
--- a/src/PJ_axisswap.c
+++ b/src/PJ_axisswap.c
@@ -274,8 +274,8 @@ PJ *CONVERSION(axisswap,0) {
     }
 
     if (pj_param(P->ctx, P->params, "tangularunits").i) {
-        P->left  = PJ_IO_UNITS_RADIANS;
-        P->right = PJ_IO_UNITS_RADIANS;
+        P->left  = PJ_IO_UNITS_ANGULAR;
+        P->right = PJ_IO_UNITS_ANGULAR;
     } else {
         P->left  = PJ_IO_UNITS_PROJECTED;
         P->right = PJ_IO_UNITS_PROJECTED;
diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index 8e11c030..0746ec08 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -213,7 +213,7 @@ PJ *CONVERSION(cart,1) {
     P->inv3d  =  geodetic;
     P->fwd    =  cart_forward;
     P->inv    =  cart_reverse;
-    P->left   =  PJ_IO_UNITS_RADIANS;
+    P->left   =  PJ_IO_UNITS_ANGULAR;
     P->right  =  PJ_IO_UNITS_CARTESIAN;
     return P;
 }
diff --git a/src/PJ_geoc.c b/src/PJ_geoc.c
index 865b1089..20064f65 100644
--- a/src/PJ_geoc.c
+++ b/src/PJ_geoc.c
@@ -48,8 +48,8 @@ static PJ *CONVERSION(geoc, 1) {
     P->inv4d = inverse;
     P->fwd4d = forward;
 
-    P->left   =  PJ_IO_UNITS_RADIANS;
-    P->right  =  PJ_IO_UNITS_RADIANS;
+    P->left   =  PJ_IO_UNITS_ANGULAR;
+    P->right  =  PJ_IO_UNITS_ANGULAR;
 
     P->is_latlong = 1;
     return P;
diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c
index f726d946..229e30c2 100644
--- a/src/PJ_helmert.c
+++ b/src/PJ_helmert.c
@@ -486,6 +486,25 @@ PJ *TRANSFORMATION(helmert, 0) {
         P->right = PJ_IO_UNITS_PROJECTED;
     }
 
+    /* Support the classic PROJ towgs84 parameter, but allow later overrides.*/
+    /* Note that if towgs84 is specified, the datum_params array is set up   */
+    /* for us automagically by the pj_datum_set call in pj_init_ctx */
+    if (pj_param_exists (P->params, "towgs84")) {
+        Q->xyz_0.x = P->datum_params[0];
+        Q->xyz_0.y = P->datum_params[1];
+        Q->xyz_0.z = P->datum_params[2];
+
+        Q->opk_0.o = P->datum_params[3];
+        Q->opk_0.p = P->datum_params[4];
+        Q->opk_0.k = P->datum_params[5];
+
+        /* We must undo conversion to absolute scale from pj_datum_set */
+        if (0==P->datum_params[6])
+            Q->scale_0 = 0;
+        else
+            Q->scale_0 = (P->datum_params[6] - 1) * 1e6;
+    }
+
     /* Translations */
     if (pj_param (P->ctx, P->params, "tx").i)
         Q->xyz_0.x = pj_param (P->ctx, P->params, "dx").f;
@@ -570,14 +589,13 @@ PJ *TRANSFORMATION(helmert, 0) {
     /* Let's help with debugging */
     if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) {
         proj_log_debug(P, "Helmert parameters:");
-        proj_log_debug(P, "x=  % 3.5f  y=  % 3.5f  z=  % 3.5f", Q->xyz.x, Q->xyz.y, Q->xyz.z);
-        proj_log_debug(P, "rx= % 3.5f  ry= % 3.5f  rz= % 3.5f",
+        proj_log_debug(P, "x=  %8.5f  y=  %8.5f  z=  %8.5f", Q->xyz.x, Q->xyz.y, Q->xyz.z);
+        proj_log_debug(P, "rx= %8.5f  ry= %8.5f  rz= %8.5f",
                 Q->opk.o / ARCSEC_TO_RAD, Q->opk.p / ARCSEC_TO_RAD, Q->opk.k / ARCSEC_TO_RAD);
-        proj_log_debug(P, "s=% 3.5f  exact=% d  transpose=% d",
-                Q->scale, Q->exact, Q->transpose);
-        proj_log_debug(P, "dx= % 3.5f  dy= % 3.5f  dz= % 3.5f", Q->dxyz.x, Q->dxyz.y, Q->dxyz.z);
-        proj_log_debug(P, "drx=% 3.5f  dry=% 3.5f  drz=% 3.5f", Q->dopk.o, Q->dopk.p, Q->dopk.k);
-        proj_log_debug(P, "ds=% 3.5f  t_epoch=% 5.5f  t_obs=% 5.5f", Q->dscale, Q->t_epoch, Q->t_obs);
+        proj_log_debug(P, "s=  %8.5f  exact=%d  transpose=%d", Q->scale, Q->exact, Q->transpose);
+        proj_log_debug(P, "dx= %8.5f  dy= %8.5f  dz= %8.5f",   Q->dxyz.x, Q->dxyz.y, Q->dxyz.z);
+        proj_log_debug(P, "drx=%8.5f  dry=%8.5f  drz=%8.5f",   Q->dopk.o, Q->dopk.p, Q->dopk.k);
+        proj_log_debug(P, "ds= %8.5f  t_epoch=%8.5f  t_obs=%8.5f", Q->dscale, Q->t_epoch, Q->t_obs);
     }
 
     if ((Q->opk.o==0)  && (Q->opk.p==0)  && (Q->opk.k==0) && (Q->scale==0) &&
diff --git a/src/PJ_hgridshift.c b/src/PJ_hgridshift.c
index 5c2b944d..54440822 100644
--- a/src/PJ_hgridshift.c
+++ b/src/PJ_hgridshift.c
@@ -54,8 +54,8 @@ PJ *TRANSFORMATION(hgridshift,0) {
     P->fwd    = 0;
     P->inv    = 0;
 
-    P->left  = PJ_IO_UNITS_RADIANS;
-    P->right = PJ_IO_UNITS_RADIANS;
+    P->left  = PJ_IO_UNITS_ANGULAR;
+    P->right = PJ_IO_UNITS_ANGULAR;
 
     if (0==pj_param(P->ctx, P->params, "tgrids").i) {
         proj_log_error(P, "hgridshift: +grids parameter missing.");
diff --git a/src/PJ_latlong.c b/src/PJ_latlong.c
index b1909954..1331d59a 100644
--- a/src/PJ_latlong.c
+++ b/src/PJ_latlong.c
@@ -98,8 +98,8 @@ static PJ *latlong_setup (PJ *P) {
     P->fwd3d = latlong_forward_3d;
     P->inv4d = latlong_inverse_4d;
     P->fwd4d = latlong_forward_4d;
-    P->left = PJ_IO_UNITS_RADIANS;
-    P->right = PJ_IO_UNITS_RADIANS;
+    P->left = PJ_IO_UNITS_ANGULAR;
+    P->right = PJ_IO_UNITS_ANGULAR;
     return P;
 }
 
diff --git a/src/PJ_molodensky.c b/src/PJ_molodensky.c
index 765e1d50..b537e802 100644
--- a/src/PJ_molodensky.c
+++ b/src/PJ_molodensky.c
@@ -281,8 +281,8 @@ PJ *TRANSFORMATION(molodensky,1) {
     P->fwd    = forward_2d;
     P->inv    = reverse_2d;
 
-    P->left   = PJ_IO_UNITS_RADIANS;
-    P->right  = PJ_IO_UNITS_RADIANS;
+    P->left   = PJ_IO_UNITS_ANGULAR;
+    P->right  = PJ_IO_UNITS_ANGULAR;
 
     /* read args */
     if (pj_param(P->ctx, P->params, "tdx").i)
diff --git a/src/PJ_ob_tran.c b/src/PJ_ob_tran.c
index debb211c..c447ac08 100644
--- a/src/PJ_ob_tran.c
+++ b/src/PJ_ob_tran.c
@@ -232,7 +232,7 @@ PJ *PROJECTION(ob_tran) {
 
     /* Support some rather speculative test cases, where the rotated projection */
     /* is actually latlong. We do not want scaling in that case... */
-    if (Q->link->right==PJ_IO_UNITS_RADIANS)
+    if (Q->link->right==PJ_IO_UNITS_ANGULAR)
         P->right = PJ_IO_UNITS_PROJECTED;
 
 
diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c
index 850734fc..bb8b4d4d 100644
--- a/src/PJ_vgridshift.c
+++ b/src/PJ_vgridshift.c
@@ -69,8 +69,8 @@ PJ *TRANSFORMATION(vgridshift,0) {
     P->fwd    = 0;
     P->inv    = 0;
 
-    P->left  = PJ_IO_UNITS_RADIANS;
-    P->right = PJ_IO_UNITS_RADIANS;
+    P->left  = PJ_IO_UNITS_ANGULAR;
+    P->right = PJ_IO_UNITS_ANGULAR;
 
     return P;
 }
diff --git a/src/gie.c b/src/gie.c
index e0045ebd..b2673d13 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -387,6 +387,7 @@ Interpret <args> as a numeric followed by a linear decadal prefix.
 Return the properly scaled numeric
 ******************************************************************************/
     double s;
+    const double GRS80_DEG = 111319.4908; /* deg-to-m at equator of GRS80 */
     const char *endp = args;
     s = proj_strtod (args, (char **) &endp);
     if (args==endp)
@@ -408,6 +409,10 @@ Return the properly scaled numeric
         s /= 1e6;
     else if (0==strcmp(endp, "nm"))
         s /= 1e9;
+    else if (0==strcmp(endp, "rad"))
+        s = GRS80_DEG * proj_todeg (s);
+    else if (0==strcmp(endp, "deg"))
+        s = GRS80_DEG * s;
     else
         s *= default_scale;
     return s;
@@ -553,19 +558,29 @@ using the "builtins" command verb.
 }
 
 
-static PJ_COORD torad_coord (PJ_COORD a) {
-    PJ_COORD c = a;
-    c.lpz.lam = proj_torad (a.lpz.lam);
-    c.lpz.phi = proj_torad (a.lpz.phi);
-    return c;
+static PJ_COORD torad_coord (PJ *P, PJ_DIRECTION dir, PJ_COORD a) {
+    size_t i, n;
+    char *axis = "enut";
+    paralist *l = pj_param_exists (P->params, "axis");
+    if (l && dir==PJ_INV)
+        axis = l->param + strlen ("axis=");
+    for (i = 0,  n = strlen (axis);  i < n;  i++)
+        if (strchr ("news", axis[i]))
+            a.v[i] = proj_torad (a.v[i]);
+    return a;
 }
 
 
-static PJ_COORD todeg_coord (PJ_COORD a) {
-    PJ_COORD c = a;
-    c.lpz.lam = proj_todeg (a.lpz.lam);
-    c.lpz.phi = proj_todeg (a.lpz.phi);
-    return c;
+static PJ_COORD todeg_coord (PJ *P, PJ_DIRECTION dir, PJ_COORD a) {
+    size_t i, n;
+    char *axis = "enut";
+    paralist *l = pj_param_exists (P->params, "axis");
+    if (l && dir==PJ_FWD)
+        axis = l->param + strlen ("axis=");
+    for (i = 0,  n = strlen (axis);  i < n;  i++)
+        if (strchr ("news", axis[i]))
+            a.v[i] = proj_todeg (a.v[i]);
+    return a;
 }
 
 
@@ -627,7 +642,7 @@ back/forward transformation pairs.
     d = d==HUGE_VAL?  T.tolerance:  d;
 
     /* input ("accepted") values - probably in degrees */
-    coo = proj_angular_input  (T.P, T.dir)? torad_coord (T.a):  T.a;
+    coo = proj_angular_input  (T.P, T.dir)? torad_coord (T.P, T.dir, T.a):  T.a;
 
     r = proj_roundtrip (T.P, T.dir, ntrips, &coo);
     if (r <= d)
@@ -747,7 +762,7 @@ Tell GIE what to expect, when transforming the ACCEPTed input
         proj_errno_reset (T.P);
 
         /* Try to carry out the operation - and expect failure */
-        ci = proj_angular_input (T.P, T.dir)? torad_coord (T.a): T.a;
+        ci = proj_angular_input (T.P, T.dir)? torad_coord (T.P, T.dir, T.a): T.a;
         co = expect_trans_n_dim (ci);
 
         /* Failed to fail? - that's a failure */
@@ -780,20 +795,20 @@ Tell GIE what to expect, when transforming the ACCEPTed input
         return expect_message_cannot_parse (args);
 
     /* expected angular values, probably in degrees */
-    ce = proj_angular_output (T.P, T.dir)? torad_coord (T.e): T.e;
+    ce = proj_angular_output (T.P, T.dir)? torad_coord (T.P, T.dir, T.e): T.e;
     if (T.verbosity > 3)
         printf ("EXPECTS  %.4f  %.4f  %.4f  %.4f\n", ce.v[0],ce.v[1],ce.v[2],ce.v[3]);
 
     /* input ("accepted") values, also probably in degrees */
-    ci = proj_angular_input (T.P, T.dir)? torad_coord (T.a): T.a;
+    ci = proj_angular_input (T.P, T.dir)? torad_coord (T.P, T.dir, T.a): T.a;
     if (T.verbosity > 3)
         printf ("ACCEPTS  %.4f  %.4f  %.4f  %.4f\n", ci.v[0],ci.v[1],ci.v[2],ci.v[3]);
 
     /* angular output from proj_trans comes in radians */
     co = expect_trans_n_dim (ci);
-    T.b = proj_angular_output (T.P, T.dir)? todeg_coord (co): co;
+    T.b = proj_angular_output (T.P, T.dir)? todeg_coord (T.P, T.dir, co): co;
     if (T.verbosity > 3)
-        printf ("GOT      %.4f  %.4f  %.4f  %.4f\n", ci.v[0],ci.v[1],ci.v[2],ci.v[3]);
+        printf ("GOT      %.4f  %.4f  %.4f  %.4f\n", co.v[0],co.v[1],co.v[2],co.v[3]);
 
     /* but there are a few more possible input conventions... */
     if (proj_angular_output (T.P, T.dir)) {
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index a1025aac..0509858a 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -32,58 +32,137 @@
 #include "proj_internal.h"
 #include "projects.h"
 
-PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD coo) {
+#define INPUT_UNITS  P->left
+#define OUTPUT_UNITS P->right
+
+
+static PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD coo) {
     if (HUGE_VAL==coo.v[0])
         return proj_coord_error ();
 
+    /* The helmert datum shift will choke unless it gets a sensible 4D coordinate */
+    if (HUGE_VAL==coo.v[2] && P->helmert) coo.v[2] = 0.0;
+    if (HUGE_VAL==coo.v[3] && P->helmert) coo.v[3] = 0.0;
+
     /* Check validity of angular input coordinates */
-    if (P->left==PJ_IO_UNITS_RADIANS) {
+    if (INPUT_UNITS==PJ_IO_UNITS_ANGULAR) {
         double t;
-        LP lp = coo.lp;
 
         /* check for latitude or longitude over-range */
-        t = (lp.phi < 0  ?  -lp.phi  :  lp.phi) - M_HALFPI;
-        if (t > PJ_EPS_LAT  ||  lp.lam > 10  ||  lp.lam < -10) {
+        t = (coo.lp.phi < 0  ?  -coo.lp.phi  :  coo.lp.phi) - M_HALFPI;
+        if (t > PJ_EPS_LAT  ||  coo.lp.lam > 10  ||  coo.lp.lam < -10) {
             proj_errno_set (P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
             return proj_coord_error ();
         }
 
         /* Clamp latitude to -90..90 degree range */
-        if (lp.phi > M_HALFPI)
-            lp.phi = M_HALFPI;
-        if (lp.phi < -M_HALFPI)
-            lp.phi = -M_HALFPI;
+        if (coo.lp.phi > M_HALFPI)
+            coo.lp.phi = M_HALFPI;
+        if (coo.lp.phi < -M_HALFPI)
+            coo.lp.phi = -M_HALFPI;
 
         /* If input latitude is geocentrical, convert to geographical */
         if (P->geoc)
             coo = proj_geoc_lat (P, PJ_INV, coo);
 
+        /* Ensure longitude is in the -pi:pi range */
+        if (0==P->over)
+            coo.lp.lam = adjlon(coo.lp.lam);
+
+        if (P->hgridshift)
+            coo = proj_trans (P->hgridshift, PJ_INV, coo);
+        else if (P->helmert) {
+            coo = proj_trans (P->cart_wgs84, PJ_FWD, coo); /* Go cartesian in WGS84 frame */
+            coo = proj_trans (P->helmert,    PJ_FWD, coo); /* Step into local frame */
+            coo = proj_trans (P->cart,       PJ_INV, coo); /* Go back to angular using local ellps */
+        }
+        if (P->vgridshift)
+            coo = proj_trans (P->vgridshift, PJ_FWD, coo);
+
         /* Distance from central meridian, taking system zero meridian into account */
-        lp.lam = (lp.lam - P->from_greenwich) - P->lam0;
+        coo.lp.lam = (coo.lp.lam - P->from_greenwich) - P->lam0;
 
         /* Ensure longitude is in the -pi:pi range */
         if (0==P->over)
-            lp.lam = adjlon(lp.lam);
+            coo.lp.lam = adjlon(coo.lp.lam);
 
-        coo.lp = lp;
+        return coo;
     }
 
+
+    /* We do not support gridshifts on cartesian input */
+    if (INPUT_UNITS==PJ_IO_UNITS_CARTESIAN && P->helmert)
+            return proj_trans (P->helmert, PJ_FWD, coo);
     return coo;
 }
 
 
-PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo) {
+
+static PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo) {
+
+    switch (OUTPUT_UNITS) {
+
+    /* Handle false eastings/northings and non-metric linear units */
+    case PJ_IO_UNITS_CARTESIAN:
+
+        if (P->is_geocent) {
+            coo = proj_trans (P->cart, PJ_FWD, coo);
+        }
+
+        coo.xyz.x = P->fr_meter * (coo.xyz.x + P->x0);
+        coo.xyz.y = P->fr_meter * (coo.xyz.y + P->y0);
+        coo.xyz.z = P->fr_meter * (coo.xyz.z + P->z0);
+        break;
 
     /* Classic proj.4 functions return plane coordinates in units of the semimajor axis */
-    if (P->right==PJ_IO_UNITS_CLASSIC) {
+    case PJ_IO_UNITS_CLASSIC:
         coo.xy.x *= P->a;
         coo.xy.y *= P->a;
+
+    /* Falls through */ /* (<-- GCC warning silencer) */
+    /* to continue processing in common with PJ_IO_UNITS_PROJECTED */
+    case PJ_IO_UNITS_PROJECTED:
+        coo.xyz.x = P->fr_meter  * (coo.xyz.x + P->x0);
+        coo.xyz.y = P->fr_meter  * (coo.xyz.y + P->y0);
+        coo.xyz.z = P->vfr_meter * (coo.xyz.z + P->z0);
+        break;
+
+    case PJ_IO_UNITS_WHATEVER:
+        break;
+
+    case PJ_IO_UNITS_ANGULAR:
+        if (INPUT_UNITS==PJ_IO_UNITS_ANGULAR)
+            break;
+
+        /* adjust longitude to central meridian */
+        if (0==P->over)
+            coo.lpz.lam = adjlon(coo.lpz.lam);
+
+        if (P->vgridshift)
+            coo = proj_trans (P->vgridshift, PJ_INV, coo);
+        if (P->hgridshift)
+            coo = proj_trans (P->hgridshift, PJ_FWD, coo);
+        else if (P->helmert) {
+            coo = proj_trans (P->cart_wgs84, PJ_FWD, coo); /* Go cartesian in WGS84 frame */
+            coo = proj_trans (P->helmert,    PJ_FWD, coo); /* Step into local frame */
+            coo = proj_trans (P->cart,       PJ_INV, coo); /* Go back to angular using local ellps */
+        }
+
+        /* If input latitude was geocentrical, convert back to geocentrical */
+        if (P->geoc)
+            coo = proj_geoc_lat (P, PJ_FWD, coo);
+
+
+        /* Distance from central meridian, taking system zero meridian into account */
+        coo.lp.lam = coo.lp.lam + P->from_greenwich + P->lam0;
+
+        /* adjust longitude to central meridian */
+        if (0==P->over)
+            coo.lpz.lam = adjlon(coo.lpz.lam);
     }
 
-    /* Handle false eastings/northings and non-metric linear units */
-    coo.xyz.x = P->fr_meter  * (coo.xyz.x + P->x0);
-    coo.xyz.y = P->fr_meter  * (coo.xyz.y + P->y0);
-    coo.xyz.z = P->vfr_meter * (coo.xyz.z + P->z0);
+    if (P->axisswap)
+        coo = proj_trans (P->axisswap, PJ_FWD, coo);
 
     return coo;
 }
diff --git a/src/pj_geocent.c b/src/pj_geocent.c
index eca62080..6530b103 100644
--- a/src/pj_geocent.c
+++ b/src/pj_geocent.c
@@ -54,6 +54,8 @@ PJ *PROJECTION(geocent) {
     P->y0 = 0.0;
     P->inv = inverse;
     P->fwd = forward;
+    P->left = PJ_IO_UNITS_ANGULAR;
+    P->right = PJ_IO_UNITS_CARTESIAN;
 
     return P;
 }
diff --git a/src/pj_init.c b/src/pj_init.c
index f25c0344..51a742c7 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -769,26 +769,3 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
     proj_errno_restore (PIN, err);
     return PIN;
 }
-
-
-
-/************************************************************************/
-/*                              pj_free()                               */
-/*                                                                      */
-/*      This is the application callable entry point for destroying     */
-/*      a projection definition.  It does work generic to all           */
-/*      projection types, and then calls the projection specific        */
-/*      free function, P->destructor(), to do local work.               */
-/*      In most cases P->destructor()==pj_default_destructor.           */
-/************************************************************************/
-
-void pj_free(PJ *P) {
-    if (0==P)
-        return;
-    /* free projection parameters - all the hard work is done by */
-    /* pj_default_destructor (in pj_malloc.c), which is supposed */
-    /* to be called as the last step of the local destructor     */
-    /* pointed to by P->destructor. In most cases,               */
-    /* pj_default_destructor actually *is* what is pointed to    */
-    P->destructor (P, 0);
-}
diff --git a/src/pj_inv.c b/src/pj_inv.c
index e1fb05f6..45e43820 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -32,46 +32,131 @@
 #include "proj_internal.h"
 #include "projects.h"
 
+#define INPUT_UNITS  P->right
+#define OUTPUT_UNITS P->left
 
-
-PJ_COORD pj_inv_prepare (PJ *P, PJ_COORD coo) {
+static PJ_COORD pj_inv_prepare (PJ *P, PJ_COORD coo) {
     if (coo.xyz.x == HUGE_VAL) {
         proj_errno_set (P, PJD_ERR_INVALID_X_OR_Y);
         return proj_coord_error ();
     }
 
+    /* The helmert datum shift will choke unless it gets a sensible 4D coordinate */
+    if (HUGE_VAL==coo.v[2] && P->helmert) coo.v[2] = 0.0;
+    if (HUGE_VAL==coo.v[3] && P->helmert) coo.v[3] = 0.0;
+
+    if (P->axisswap)
+        coo = proj_trans (P->axisswap, PJ_INV, coo);
+
+    /* Check validity of angular input coordinates */
+    if (INPUT_UNITS==PJ_IO_UNITS_ANGULAR) {
+        double t;
+
+        /* check for latitude or longitude over-range */
+        t = (coo.lp.phi < 0  ?  -coo.lp.phi  :  coo.lp.phi) - M_HALFPI;
+        if (t > PJ_EPS_LAT  ||  coo.lp.lam > 10  ||  coo.lp.lam < -10) {
+            proj_errno_set (P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
+            return proj_coord_error ();
+        }
+
+        /* Clamp latitude to -90..90 degree range */
+        if (coo.lp.phi > M_HALFPI)
+            coo.lp.phi = M_HALFPI;
+        if (coo.lp.phi < -M_HALFPI)
+            coo.lp.phi = -M_HALFPI;
+
+        /* If input latitude is geocentrical, convert to geographical */
+        if (P->geoc)
+            coo = proj_geoc_lat (P, PJ_INV, coo);
+
+        /* Distance from central meridian, taking system zero meridian into account */
+        coo.lp.lam = (coo.lp.lam + P->from_greenwich) - P->lam0;
+
+        /* Ensure longitude is in the -pi:pi range */
+        if (0==P->over)
+            coo.lp.lam = adjlon(coo.lp.lam);
+
+        if (P->hgridshift)
+            coo = proj_trans (P->hgridshift, PJ_FWD, coo);
+        else if (P->helmert) {
+            coo = proj_trans (P->cart,       PJ_FWD, coo); /* Go cartesian in local frame */
+            coo = proj_trans (P->helmert,    PJ_INV, coo); /* Step into WGS84 */
+            coo = proj_trans (P->cart_wgs84, PJ_INV, coo); /* Go back to angular using WGS84 ellps */
+        }
+        if (P->vgridshift)
+            coo = proj_trans (P->vgridshift, PJ_INV, coo);
+        return coo;
+    }
+
+    /* Handle remaining possible input types */
+    switch (INPUT_UNITS) {
+    case PJ_IO_UNITS_WHATEVER:
+        return coo;
+
     /* de-scale and de-offset */
-    coo.xyz.x = (coo.xyz.x * P->to_meter  - P->x0);
-    coo.xyz.y = (coo.xyz.y * P->to_meter  - P->y0);
-    coo.xyz.z = (coo.xyz.z * P->vto_meter - P->z0);
-
-    /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis  */
-    /* Multiplying by ra, rather than dividing by a because the CALCOFI projection       */
-    /* stomps on a and hence (apparently) depends on this to roundtrip correctly         */
-    /* (CALCOFI avoids further scaling by stomping - a better solution must be possible) */
-    if (P->right==PJ_IO_UNITS_CLASSIC) {
+    case PJ_IO_UNITS_CARTESIAN:
+        coo.xyz.x = P->to_meter * coo.xyz.x - P->x0;
+        coo.xyz.y = P->to_meter * coo.xyz.y - P->y0;
+        coo.xyz.z = P->to_meter * coo.xyz.z - P->z0;
+
+        if (P->is_geocent) {
+            coo = proj_trans (P->cart, PJ_INV, coo);
+        }
+
+        return coo;
+
+    case PJ_IO_UNITS_PROJECTED:
+    case PJ_IO_UNITS_CLASSIC:
+        coo.xyz.x = P->to_meter  * coo.xyz.x - P->x0;
+        coo.xyz.y = P->to_meter  * coo.xyz.y - P->y0;
+        coo.xyz.z = P->vto_meter * coo.xyz.z - P->z0;
+        if (INPUT_UNITS==PJ_IO_UNITS_PROJECTED)
+            return coo;
+
+        /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis  */
+        /* Multiplying by ra, rather than dividing by a because the CALCOFI projection       */
+        /* stomps on a and hence (apparently) depends on this to roundtrip correctly         */
+        /* (CALCOFI avoids further scaling by stomping - but a better solution is possible)  */
         coo.xyz.x *= P->ra;
         coo.xyz.y *= P->ra;
+        return coo;
+    /* Silence some compiler warnings about PJ_IO_UNITS_ANGULAR not handled */
+    default:
+        break;
     }
 
+    /* Should not happen, so we could return pj_coord_err here */
     return coo;
 }
 
 
 
-PJ_COORD pj_inv_finalize (PJ *P, PJ_COORD coo) {
+static PJ_COORD pj_inv_finalize (PJ *P, PJ_COORD coo) {
     if (coo.xyz.x == HUGE_VAL) {
         proj_errno_set (P, PJD_ERR_INVALID_X_OR_Y);
         return proj_coord_error ();
     }
 
-    if (P->left==PJ_IO_UNITS_RADIANS) {
-        /* Distance from central meridian, taking system zero meridian into account */
-        coo.lp.lam = coo.lp.lam + P->from_greenwich + P->lam0;
-
-         /* adjust longitude to central meridian */
-        if (0==P->over)
-            coo.lpz.lam = adjlon(coo.lpz.lam);
+    if (OUTPUT_UNITS==PJ_IO_UNITS_ANGULAR) {
+
+        if (INPUT_UNITS!=PJ_IO_UNITS_ANGULAR) {
+            /* Distance from central meridian, taking system zero meridian into account */
+            coo.lp.lam = coo.lp.lam + P->from_greenwich + P->lam0;
+
+            /* adjust longitude to central meridian */
+            if (0==P->over)
+                coo.lpz.lam = adjlon(coo.lpz.lam);
+
+            if (P->vgridshift)
+                coo = proj_trans (P->vgridshift, PJ_INV, coo);
+            if (P->hgridshift)
+                coo = proj_trans (P->hgridshift, PJ_FWD, coo);
+            else if (P->helmert) {
+                coo = proj_trans (P->cart,       PJ_FWD, coo); /* Go cartesian in local frame */
+                coo = proj_trans (P->helmert,    PJ_INV, coo); /* Step into WGS84 */
+                coo = proj_trans (P->cart_wgs84, PJ_INV, coo); /* Go back to angular using WGS84 ellps */
+            }
+        }
 
         /* If input latitude was geocentrical, convert back to geocentrical */
         if (P->geoc)
diff --git a/src/pj_malloc.c b/src/pj_malloc.c
index 63278a56..cfb7f658 100644
--- a/src/pj_malloc.c
+++ b/src/pj_malloc.c
@@ -162,6 +162,32 @@ void *pj_dealloc_params (PJ_CONTEXT *ctx, paralist *start, int errlev) {
 }
 
 
+
+
+/************************************************************************/
+/*                              pj_free()                               */
+/*                                                                      */
+/*      This is the application callable entry point for destroying     */
+/*      a projection definition.  It does work generic to all           */
+/*      projection types, and then calls the projection specific        */
+/*      free function, P->destructor(), to do local work.               */
+/*      In most cases P->destructor()==pj_default_destructor.           */
+/************************************************************************/
+
+void pj_free(PJ *P) {
+    if (0==P)
+        return;
+    /* free projection parameters - all the hard work is done by */
+    /* pj_default_destructor, which is supposed */
+    /* to be called as the last step of the local destructor     */
+    /* pointed to by P->destructor. In most cases,               */
+    /* pj_default_destructor actually *is* what is pointed to    */
+    P->destructor (P, 0);
+}
+
+
+
+
 /*****************************************************************************/
 void *pj_default_destructor (PJ *P, int errlev) {   /* Destructor */
 /*****************************************************************************
@@ -196,6 +222,14 @@ void *pj_default_destructor (PJ *P, int errlev) {   /* Destructor */
     /* free parameter list elements */
     pj_dealloc_params (pj_get_ctx(P), P->params, errlev);
 
+    /* free the cs2cs emulation elements */
+    pj_free (P->axisswap);
+    pj_free (P->helmert);
+    pj_free (P->cart);
+    pj_free (P->cart_wgs84);
+    pj_free (P->hgridshift);
+    pj_free (P->vgridshift);
+
     pj_dealloc (P->opaque);
     return pj_dealloc(P);
 }
diff --git a/src/proj.c b/src/proj.c
index 1d702ea4..1c86854e 100644
--- a/src/proj.c
+++ b/src/proj.c
@@ -152,12 +152,12 @@ static void process(FILE *fid) {
             }
         } else {    /* x-y or decimal degree ascii output, scale if warranted by output units */
             if (inverse) {
-                if (Proj->left == PJ_IO_UNITS_RADIANS) {
+                if (Proj->left == PJ_IO_UNITS_ANGULAR) {
                     data.uv.v *= RAD_TO_DEG;
                     data.uv.u *= RAD_TO_DEG;
                 }
             } else {
-                if (Proj->right == PJ_IO_UNITS_RADIANS) {
+                if (Proj->right == PJ_IO_UNITS_ANGULAR) {
                     data.uv.v *= RAD_TO_DEG;
                     data.uv.u *= RAD_TO_DEG;
                 }
@@ -262,7 +262,7 @@ static void vprocess(FILE *fid) {
         }
 
         /* apply rad->deg scaling in case the output from a pipeline has degrees as units */
-        if (!inverse && Proj->right == PJ_IO_UNITS_RADIANS) {
+        if (!inverse && Proj->right == PJ_IO_UNITS_ANGULAR) {
             dat_xy.x *= RAD_TO_DEG;
             dat_xy.y *= RAD_TO_DEG;
         }
diff --git a/src/proj.def b/src/proj.def
index 6ab0c007..099a78a5 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -155,3 +155,5 @@ EXPORTS
     pj_approx_2D_trans             @138
     pj_approx_3D_trans             @139
     pj_has_inverse                 @140
+    pj_param_exists                @141
+
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 2c6bac2d..2b250398 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -52,8 +52,8 @@ int proj_angular_input (PJ *P, enum PJ_DIRECTION dir) {
     dir: {PJ_FWD, PJ_INV}
 ******************************************************************************/
     if (PJ_FWD==dir)
-        return pj_left (P)==PJ_IO_UNITS_RADIANS;
-    return pj_right (P)==PJ_IO_UNITS_RADIANS;
+        return pj_left (P)==PJ_IO_UNITS_ANGULAR;
+    return pj_right (P)==PJ_IO_UNITS_ANGULAR;
 }
 
 /*****************************************************************************/
@@ -100,7 +100,7 @@ double proj_xyz_dist (XYZ a, XYZ b) {
 /* Measure numerical deviation after n roundtrips fwd-inv (or inv-fwd) */
 double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo) {
     int i;
-    PJ_COORD o, u, org;
+    PJ_COORD t, org;
 
     if (0==P)
         return HUGE_VAL;
@@ -111,23 +111,23 @@ double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo) {
     }
 
     /* in the first half-step, we generate the output value */
-    u = org = *coo;
-    o = *coo = proj_trans (P, direction, u);
+    org  = *coo;
+    *coo = proj_trans (P, direction, org);
+    t = *coo;
 
-    /* now we take n-1 full steps */
-    for (i = 0;  i < n - 1;  i++) {
-        u = proj_trans (P, -direction, o);
-        o = proj_trans (P,  direction, u);
-    }
+    /* now we take n-1 full steps in inverse direction: We are */
+    /* out of phase due to the half step already taken */
+    for (i = 0;  i < n - 1;  i++)
+        t = proj_trans (P,  direction,  proj_trans (P, -direction, t) );
 
     /* finally, we take the last half-step */
-    u = proj_trans (P, -direction, o);
+    t = proj_trans (P, -direction, t);
 
     /* checking for angular *input* since we do a roundtrip, and end where we begin */
     if (proj_angular_input (P, direction))
-        return proj_lpz_dist (P, org.lpz, u.lpz);
+        return proj_lpz_dist (P, org.lpz, t.lpz);
 
-    return proj_xyz_dist (org.xyz, u.xyz);
+    return proj_xyz_dist (org.xyz, t.xyz);
 }
 
 
@@ -380,7 +380,114 @@ char*  proj_rtodms(char *s, double r, int pos, int neg) {
     return rtodms(s, r, pos, neg);
 }
 
+/*************************************************************************************/
+static PJ* skip_prep_fin(PJ *P) {
+/**************************************************************************************
+Skip prepare and finalize function for the various "helper operations" added to P when
+in cs2cs compatibility mode.
+**************************************************************************************/
+    P->skip_fwd_prepare  = 1;
+    P->skip_fwd_finalize = 1;
+    P->skip_inv_prepare  = 1;
+    P->skip_inv_finalize = 1;
+    return P;
+}
 
+/*************************************************************************************/
+static int pj_cs2cs_emulation_setup (PJ *P) {
+/**************************************************************************************
+If any cs2cs style modifiers are given (axis=..., towgs84=..., ) create the 4D API
+equivalent operations, so the preparation and finalization steps in the pj_inv/pj_fwd
+invocators can emulate the behaviour of pj_transform and the cs2cs app.
+**************************************************************************************/
+    PJ *Q;
+    paralist *p;
+    char def[1000];
+    if (0==P)
+        return 0;
+
+    /* Don't recurse when calling proj_create (which calls us back) */
+    if (pj_param_exists (P->params, "break_cs2cs_recursion"))
+        return 1;
+
+    /* Swap axes? */
+    p = pj_param_exists (P->params, "axis");
+
+    /* Don't axisswap if data are already in "enu" order */
+    if (p && (0!=strcmp ("enu", p->param))) {
+        sprintf (def, "break_cs2cs_recursion     proj=axisswap  axis=%s", P->axis);
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        P->axisswap = skip_prep_fin(Q);
+    }
+
+    /* Geoid grid(s) given? */
+    p = pj_param_exists (P->params, "geoidgrids");
+    if (p  &&  strlen (p->param) > strlen ("geoidgrids=")) {
+        char *gridnames = p->param + strlen ("geoidgrids=");
+        sprintf (def, "break_cs2cs_recursion     proj=vgridshift  grids=%s", gridnames);
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        P->vgridshift = skip_prep_fin(Q);
+    }
+
+    /* Datum shift grid(s) given? */
+    p = pj_param_exists (P->params, "nadgrids");
+    if (p  &&  strlen (p->param) > strlen ("nadgrids=")) {
+        char *gridnames = p->param + strlen ("nadgrids=");
+        sprintf (def, "break_cs2cs_recursion     proj=hgridshift  grids=%s", gridnames);
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        P->hgridshift = skip_prep_fin(Q);
+    }
+
+    /* We ignore helmert if we have grid shift */
+    p = P->hgridshift ? 0 : pj_param_exists (P->params, "towgs84");
+    while (p) {
+        char *s = p->param;
+        double *d = P->datum_params;
+        size_t n = strlen (s);
+
+        /* We ignore null helmert shifts (common in auto-translated resource files, e.g. epsg) */
+        if (0==d[0] && 0==d[1] && 0==d[2] && 0==d[3] && 0==d[4] && 0==d[5] && 0==d[6])
+            break;
+
+        if (n > 900)
+            return 0;
+        if (n <= 8) /* 8==strlen ("towgs84=") */
+            return 0;
+        sprintf (def, "break_cs2cs_recursion     proj=helmert %s", s);
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        P->helmert = skip_prep_fin(Q);
+
+        break;
+    }
+
+    /* We also need cartesian/geographical transformations if we are working in */
+    /* geocentric/cartesian space or we need to do a Helmert transform.         */
+    if (P->is_geocent || P->helmert) {
+        char *wgs84 = "ellps=WGS84";
+        sprintf (def, "break_cs2cs_recursion     proj=cart");
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        pj_inherit_ellipsoid_def(P, Q);
+        P->cart = skip_prep_fin(Q);
+
+        sprintf (def, "break_cs2cs_recursion     proj=cart  %s", wgs84);
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        P->cart_wgs84 = skip_prep_fin(Q);
+    }
+
+    return 1;
+}
 
 
 /*************************************************************************************/
@@ -394,9 +501,10 @@ PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) {
     It may even use free formatting "proj  =  utm;  zone  =32  ellps= GRS80".
     Note that the semicolon separator is allowed, but not required.
 **************************************************************************************/
-    PJ   *P;
-    char *args, **argv;
+    PJ    *P;
+    char  *args, **argv;
     size_t argc, n;
+    int    ret;
 
     if (0==ctx)
         ctx = pj_get_default_ctx ();
@@ -421,6 +529,12 @@ PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) {
 
     pj_dealloc (argv);
     pj_dealloc (args);
+
+    /* Support cs2cs-style modifiers */
+    ret = pj_cs2cs_emulation_setup  (P);
+    if (0==ret)
+        return proj_destroy (P);
+
     return P;
 }
 
@@ -449,6 +563,7 @@ indicator, as in {"+proj=utm", "+zone=32"}, or leave it out, as in {"proj=utm",
         return 0;
 
     P = proj_create (ctx, c);
+
     pj_dealloc ((char *) c);
     return P;
 }
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 41f95c07..1b98c346 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -68,7 +68,7 @@ enum pj_io_units {
     PJ_IO_UNITS_CLASSIC   = 1,  /* Scaled meters (right), projected system */
     PJ_IO_UNITS_PROJECTED = 2,  /* Meters, projected system */
     PJ_IO_UNITS_CARTESIAN = 3,  /* Meters, 3D cartesian system */
-    PJ_IO_UNITS_RADIANS   = 4   /* Radians */
+    PJ_IO_UNITS_ANGULAR   = 4   /* Radians */
 };
 enum pj_io_units pj_left (PJ *P);
 enum pj_io_units pj_right (PJ *P);
@@ -84,10 +84,6 @@ void proj_context_inherit (PJ *parent, PJ *child);
 PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P);
 PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P);
 
-PJ_COORD pj_fwd_prepare  (PJ *P, PJ_COORD coo);
-PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo);
-PJ_COORD pj_inv_prepare  (PJ *P, PJ_COORD coo);
-PJ_COORD pj_inv_finalize (PJ *P, PJ_COORD coo);
 PJ_COORD pj_approx_2D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
 PJ_COORD pj_approx_3D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
 
diff --git a/src/projects.h b/src/projects.h
index 86ca2524..e4bdca4f 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -196,7 +196,7 @@ enum pj_io_units {
     PJ_IO_UNITS_CLASSIC   = 1,  /* Scaled meters (right), projected system */
     PJ_IO_UNITS_PROJECTED = 2,  /* Meters, projected system */
     PJ_IO_UNITS_CARTESIAN = 3,  /* Meters, 3D cartesian system */
-    PJ_IO_UNITS_RADIANS   = 4   /* Radians */
+    PJ_IO_UNITS_ANGULAR   = 4   /* Radians */
 };
 #endif
 #ifndef PROJ_H
@@ -349,6 +349,14 @@ struct PJconsts {
     enum pj_io_units left;          /* Flags for input/output coordinate types */
     enum pj_io_units right;
 
+    /* These PJs are used for implementing cs2cs style coordinate handling in the 4D API */
+    PJ *axisswap;
+    PJ *cart;
+    PJ *cart_wgs84;
+    PJ *helmert;
+    PJ *hgridshift;
+    PJ *vgridshift;
+
 
     /*************************************************************************************
 
@@ -394,7 +402,7 @@ struct PJconsts {
     double  from_greenwich;            /* prime meridian offset (in radians) */
     double  long_wrap_center;          /* 0.0 for -180 to 180, actually in radians*/
     int     is_long_wrap_set;
-    char    axis[4];                   /* TODO: Description needed */
+    char    axis[4];                   /* Axis order, pj_transform/pj_adjust_axis */
 
     /* New Datum Shift Grid Catalogs */
     char   *catalog_name;
@@ -615,7 +623,7 @@ C_NAMESPACE PJ *pj_##name (PJ *P) {                          \
     P->destructor = pj_default_destructor;                   \
     P->descr = des_##name;                                   \
     P->need_ellps = NEED_ELLPS;                              \
-    P->left  = PJ_IO_UNITS_RADIANS;                          \
+    P->left  = PJ_IO_UNITS_ANGULAR;                          \
     P->right = PJ_IO_UNITS_CLASSIC;                          \
     return P;                                                \
 }                                                            \
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index a91492da..7475f715 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -7,18 +7,21 @@ proj_add_gie_test("Axisswap" "gie/axisswap.gie")
 proj_add_gie_test("Deformation" "gie/deformation.gie")
 proj_add_gie_test("Ellipsoid" "gie/ellipsoid.gie")
 proj_add_gie_test("GDA" "gie/GDA.gie")
+proj_add_gie_test("4D-API-cs2cs-style" "gie/4D-API_cs2cs-style.gie")
 
+# GIGS tests. Uncommented tests are expected to fail due to issues with
+# various projections. Should be investigated further and fixed.
 proj_add_gie_test("GIGS-5101.1-jhs" "gigs/5101.1-jhs.gie")
 proj_add_gie_test("GIGS-5101.2-jhs" "gigs/5101.2-jhs.gie")
 proj_add_gie_test("GIGS-5101.3-jhs" "gigs/5101.3-jhs.gie")
 proj_add_gie_test("GIGS-5101.4-jhs-etmerc" "gigs/5101.4-jhs-etmerc.gie")
-#proj_add_gie_test("GIGS-5101.4-jhs" "gigs/5101.4-jhs.gie")
+#proj_add_gie_test("GIGS-5101.4-jhs" "gigs/5101.4-jhs.gie") # Same as above, but using etmerc instead of tmerc
 proj_add_gie_test("GIGS-5102.1" "gigs/5102.1.gie")
 #proj_add_gie_test("GIGS-5102.2" "gigs/5102.2.gie")
 proj_add_gie_test("GIGS-5103.1" "gigs/5103.1.gie")
 proj_add_gie_test("GIGS-5103.2" "gigs/5103.2.gie")
 proj_add_gie_test("GIGS-5103.3" "gigs/5103.3.gie")
-proj_add_gie_test("GIGS-5104" "gigs/5104.gie")
+#proj_add_gie_test("GIGS-5104" "gigs/5104.gie")
 #proj_add_gie_test("GIGS-5105.1" "gigs/5105.1.gie")
 proj_add_gie_test("GIGS-5105.2" "gigs/5105.2.gie")
 proj_add_gie_test("GIGS-5106" "gigs/5106.gie")
@@ -29,12 +32,15 @@ proj_add_gie_test("GIGS-5109" "gigs/5109.gie")
 proj_add_gie_test("GIGS-5111.1" "gigs/5111.1.gie")
 #proj_add_gie_test("GIGS-5111.2" "gigs/5111.2.gie")
 proj_add_gie_test("GIGS-5112" "gigs/5112.gie")
-#proj_add_gie_test("GIGS-5113" "gigs/5113.gie")
-#proj_add_gie_test("GIGS-5201" "gigs/5201.gie")
+proj_add_gie_test("GIGS-5113" "gigs/5113.gie")
+proj_add_gie_test("GIGS-5201" "gigs/5201.gie")
 #proj_add_gie_test("GIGS-5203" "gigs/5203.1.gie")
 #proj_add_gie_test("GIGS-5204.1" "gigs/5204.1.gie")
 #proj_add_gie_test("GIGS-5205.1" "gigs/5205.1.gie")
 #proj_add_gie_test("GIGS-5206" "gigs/5206.gie")
 #proj_add_gie_test("GIGS-5207.1" "gigs/5207.1.gie")
 #proj_add_gie_test("GIGS-5207.2" "gigs/5207.2.gie")
-#proj_add_gie_test("GIGS-5208" "gigs/5208.gie")
+proj_add_gie_test("GIGS-5208" "gigs/5208.gie")
+
+
+
diff --git a/test/gie/4D-API_cs2cs-style.gie b/test/gie/4D-API_cs2cs-style.gie
new file mode 100644
index 00000000..853df993
--- /dev/null
+++ b/test/gie/4D-API_cs2cs-style.gie
@@ -0,0 +1,161 @@
+===============================================================================
+
+Test the 4D API handling of cs2cs style transformation options.
+
+These tests are mostly based on the same material as those in
+more_builtins.gie, since we are testing the same kinds of things,
+but provided through a different interface.
+
+===============================================================================
+
+
+<gie>
+
+-------------------------------------------------------------------------------
+Test the handling of the +towgs84 parameter
+-------------------------------------------------------------------------------
+This example is from Lotti Jivall: "Simplified transformations from
+ITRF2008/IGS08 to ETRS89 for maritime applications". The original
+XYZ data (cf. more_builtins.gie) have been transformed to LPZ using
+this command: echo x_val y_val z_val | cct -It0 proj=cart ellps=GRS80 --
+
+NOTE: Here, the ellipsoid has been swapped to WGS84, to align with
+      the WGS84 ellipsoid used in the cs2cs emulation introduced by
+      pj_cs2cs_emulation_setup()
+-------------------------------------------------------------------------------
+operation proj=latlong ellps=WGS84
+          towgs84 =  0.676780,  0.654950, -0.528270,
+                    -0.022742,  0.012667,  0.022704, -0.010700
+-------------------------------------------------------------------------------
+tolerance 0.05 mm
+direction forward
+
+accept    13.4999969828397   54.9999995493206   -0.6034
+expect    13.4999906103972   54.9999957469562   -0.6374
+
+direction inverse
+
+accept    13.4999906103972   54.9999957469562   -0.6374
+expect    13.4999969828397   54.9999995493206   -0.6034
+-------------------------------------------------------------------------------
+
+
+
+-------------------------------------------------------------------------------
+This example is a random point, transformed from ED50 to ETRS89 using KMStrans2.
+
+NOTE: Signs swapped wrt. KMStrans2, which apparently uses frame rotation.
+-------------------------------------------------------------------------------
+operation proj=latlong ellps=intl
+          towgs84 =  081.07030, 089.36030, 115.75260,
+                     000.48488, 000.02436, 000.41321, 0.540645
+-------------------------------------------------------------------------------
+tolerance  0.1 mm
+# echo SWAPPED and trimmed - expect succcess
+# accept     16.819999997    55.170000002        61.0
+accept     16.82    55.17        61.0
+expect     16.8210462130   55.1705688946       29.0317
+-------------------------------------------------------------------------------
+This commented-out version of the example above was used to detect the sign
+discrepancy between KMSTrans2 and PROJ.
+-------------------------------------------------------------------------------
+operation proj=latlong ellps=intl
+          towgs84 =  -081.07030, -089.36030, -115.75260,
+                     -000.48488, -000.02436, -000.41321, -0.540645
+-------------------------------------------------------------------------------
+# echo NOT SWAPPED - Expect failure
+# accept     16.8210462130   55.1705688946       29.0317
+# expect     16.819999997    55.170000002        61.0
+-------------------------------------------------------------------------------
+
+
+
+-------------------------------------------------------------------------------
+operation  proj=latlong nadgrids=nzgd2kgrid0005.gsb  ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 1 nm
+accept    173                 -45                  0 0
+expect    172.999892181021551 -45.001620431954613  0 0
+direction inverse
+accept    172.999892181021551 -45.001620431954613  0 0
+expect    173                 -45                  0 0
+-------------------------------------------------------------------------------
+
+
+
+-------------------------------------------------------------------------------
+operation  proj=latlong geoidgrids=egm96_15.gtx  ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 15 cm      # lax tolerance due to widespread bad egm96 file
+accept    12.5 55.5   0
+expect    12.5 55.5 -36.0213
+direction inverse
+accept    12.5 55.5 -36.0213
+expect    12.5 55.5  0
+-------------------------------------------------------------------------------
+operation  proj=merc geoidgrids=egm96_15.gtx  ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept    12.5 55.5 0
+expect    1391493.63492   7424275.19462  -36.0213
+direction inverse
+accept    1391493.63492   7424275.19462  -36.0213
+expect    12.5 55.5 0
+-------------------------------------------------------------------------------
+
+
+
+-------------------------------------------------------------------------------
+Same as the two above, but also do axis swapping.
+-------------------------------------------------------------------------------
+operation  proj=latlong geoidgrids=egm96_15.gtx  axis=neu ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 15 cm      # lax tolerance due to widely distributed, bad egm96 file
+accept    12.5 55.5  0
+expect    55.5 12.5 -36.0213
+direction inverse
+accept    55.5 12.5 -36.0213
+expect    12.5 55.5  0
+-------------------------------------------------------------------------------
+operation  proj=latlong geoidgrids=egm96_15.gtx  axis=dne ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 15 cm      # lax tolerance due to widely distributed, bad egm96 file
+accept    12.5 55.5  0
+expect    36.0213 55.5 12.5
+direction inverse
+accept    36.0213 55.5 12.5
+expect    12.5 55.5  0
+-------------------------------------------------------------------------------
+operation  proj=merc geoidgrids=egm96_15.gtx  ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept    12.5 55.5 0
+expect    1391493.63492   7424275.19462  -36.0213
+direction inverse
+accept    1391493.63492   7424275.19462  -36.0213
+expect    12.5 55.5 0
+-------------------------------------------------------------------------------
+
+
+-------------------------------------------------------------------------------
+Some more complex axis swapping.
+-------------------------------------------------------------------------------
+operation  proj=latlong geoidgrids=egm96_15.gtx  axis=nue ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 15 cm      # lax tolerance due to widely distributed, bad egm96 file
+accept    12.5 55.5 0
+expect    55.5 -36.0213 12.5
+direction inverse
+accept    55.5 -36.0213 12.5
+expect    12.5 55.5 0
+-------------------------------------------------------------------------------
+operation  proj=merc geoidgrids=egm96_15.gtx axis=sue ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 15 cm
+accept    12.5 55.5 0
+expect    -7424275.1946 -36.0213  1391493.6349 0.0000
+direction inverse
+accept    -7424275.1946 -36.0213  1391493.6349 0.0000
+expect    12.5 55.5 0
+-------------------------------------------------------------------------------
+</gie>
diff --git a/test/gie/axisswap.gie b/test/gie/axisswap.gie
index ac148a92..5fe86200 100644
--- a/test/gie/axisswap.gie
+++ b/test/gie/axisswap.gie
@@ -53,9 +53,17 @@ expect      3 -2 1 4
 roundtrip   100
 
 operation   proj=axisswap axis=neu
-tolerance   0.000001 m
-accept      1 2 3 4
-expect      2 1 3 4
+tolerance   0 m
+accept      1 2 3
+expect      2 1 3
+
+# when using the +axis parameter we specify the order of the INPUT coordinate,
+# as opposed to +order which relates to the OUTPUT coordinate. Here we test
+# that n(1), u(2) and e(3) are swapped correctly to enu ordering.
+operation   proj=axisswap axis=nue
+tolerance   0 m
+accept      1 2 3
+expect      2 3 1
 
 operation   proj=axisswap axis=swd
 tolerance   0.000001 m
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index 3877b58a..e8f1bfe1 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -1233,14 +1233,14 @@ Geocentric
 operation +proj=geocent   +ellps=GRS80  +lat_1=0.5 +lat_2=2
 -------------------------------------------------------------------------------
 tolerance 0.00010 mm
-accept  2 1
-expect  222638.981586547 111319.490793274
-accept  2 -1
-expect  222638.981586547 -111319.490793274
-accept  -2 1
-expect  -222638.981586547 111319.490793274
-accept  -2 -1
-expect  -222638.981586547 -111319.490793274
+accept  2 1 0
+expect  6373287.27950247        222560.09599219 110568.77482092
+accept  2 -1 0
+expect  6373287.27950247        222560.09599219 -110568.77482092
+accept  -2 1 0
+expect  6373287.27950247        -222560.09599219 110568.77482092
+accept  -2 -1 0
+expect  6373287.27950247        -222560.09599219 -110568.77482092
 
 direction inverse
 accept  200 100
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index 20bffbd9..e21f0321 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -131,11 +131,11 @@ expect    12 55 0 0
 -------------------------------------------------------------------------------
 Finally test a pipeline with more than one init step
 -------------------------------------------------------------------------------
-operation proj=pipeline        step
-          init=epsg:25832 inv  step
-          init=epsg:25833      step
-          init=epsg:25833 inv  step
-          init=epsg:25832
+operation proj=pipeline
+          step init=epsg:25832 inv
+          step init=epsg:25833
+          step init=epsg:25833 inv
+          step init=epsg:25832
 -------------------------------------------------------------------------------
 accept    691875.63214  6098907.82501  0  0
 expect    691875.63214  6098907.82501  0  0
diff --git a/test/gigs/5103.1.gie b/test/gigs/5103.1.gie
index f00ed019..c972b738 100644
--- a/test/gigs/5103.1.gie
+++ b/test/gigs/5103.1.gie
@@ -10,254 +10,201 @@ Test 5103 (part 1), Lambert Conic Conformal (2SP), v2-0_2011-06-28.
 operation  +proj=pipeline
            +step +init=epsg:4313 +inv
            +step +init=epsg:31370
+
+tolerance  30 mm
 --------------------------------------------------------------------------------
-tolerance  0.03 m
 accept    5 58
 expect    187742.7 969521.653
 
-tolerance  0.03 m
 accept    5 57
 expect    188698.877 857277.135
 
-tolerance  0.03 m
 accept    5 56
 expect    189652.853 745291.184
 
-tolerance  0.03 m
 accept    5 55
 expect    190604.967 633523.672
 
-tolerance  0.03 m
 accept    5 54
 expect    191555.55 521935.9
 
-tolerance  0.03 m
 accept    5 53
 expect    192504.921 410490.433
 
-tolerance  0.03 m
 accept    5.3876389 52.1561606
 expect    219843.841 316827.604
 
-tolerance  0.03 m
 accept    4 51
 expect    124202.936 187756.876
 
-tolerance  0.03 m
 accept    4 50
 expect    123652.406 76521.628
 
-tolerance  0.03 m
 accept    4 49
 expect    123101.889 -34711.068
 
-tolerance  0.03 m
 accept    3.3137281 47.9752611
 expect    71254.553 -148236.592
 
-tolerance  0.03 m
 accept    3 53
 expect    58108.966 411155.591
 
-tolerance  0.03 m
 accept    4 53
 expect    125304.704 410370.504
 
-tolerance  0.03 m
 accept    5 53
 expect    192504.921 410490.433
 
-tolerance  0.03 m
 accept    6 53
 expect    259697.429 411515.356
 
-tolerance  0.03 m
 accept    7 53
 expect    326870.04 413445.087
 
-tolerance  0.03 m
 accept    8 53
 expect    394010.571 416279.276
 
-tolerance  0.03 m
 accept    9 53
 expect    461106.844 420017.408
 
-tolerance  0.03 m
 accept    10 53
 expect    528146.69 424658.807
 
-tolerance  0.03 m
 accept    11 53
 expect    595117.95 430202.63
 
+
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
-           +step +init=epsg:31370 +inv
-           +step +init=epsg:4313
+operation  proj=pipeline
+           step init=epsg:31370 inv
+           step init=epsg:4313
+
+tolerance  30 mm
 --------------------------------------------------------------------------------
-tolerance  0.03 m
 accept    187742.7 969521.653
 expect    5 58
 
-tolerance  0.03 m
 accept    188698.877 857277.135
 expect    5 57
 
-tolerance  0.03 m
 accept    189652.853 745291.184
 expect    5 56
 
-tolerance  0.03 m
 accept    190604.967 633523.672
 expect    5 55
 
-tolerance  0.03 m
 accept    191555.55 521935.9
 expect    5 54
 
-tolerance  0.03 m
 accept    192504.921 410490.433
 expect    5 53
 
-tolerance  0.03 m
 accept    219843.841 316827.604
 expect    5.3876389 52.1561606
 
-tolerance  0.03 m
 accept    124202.936 187756.876
 expect    4 51
 
-tolerance  0.03 m
 accept    123652.406 76521.628
 expect    4 50
 
-tolerance  0.03 m
 accept    123101.889 -34711.068
 expect    4 49
 
-tolerance  0.03 m
 accept    71254.553 -148236.592
 expect    3.3137281 47.9752611
 
-tolerance  0.03 m
 accept    58108.966 411155.591
 expect    3 53
 
-tolerance  0.03 m
 accept    125304.704 410370.504
 expect    4 53
 
-tolerance  0.03 m
 accept    192504.921 410490.433
 expect    5 53
 
-tolerance  0.03 m
 accept    259697.429 411515.356
 expect    6 53
 
-tolerance  0.03 m
 accept    326870.04 413445.087
 expect    7 53
 
-tolerance  0.03 m
 accept    394010.571 416279.276
 expect    8 53
 
-tolerance  0.03 m
 accept    461106.844 420017.408
 expect    9 53
 
-tolerance  0.03 m
 accept    528146.69 424658.807
 expect    10 53
 
-tolerance  0.03 m
 accept    595117.95 430202.63
 expect    11 53
 
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
+operation  +proj=pipeline   towgs84=0,0,0   # turn off dual datum shift
            +step +init=epsg:4313 +inv
            +step +init=epsg:31370
+
+tolerance  6 mm
 --------------------------------------------------------------------------------
-tolerance  0.006 m
 accept    5 58
 roundtrip 1000
 
-tolerance  0.006 m
 accept    5 57
 roundtrip 1000
 
-tolerance  0.006 m
 accept    5 56
 roundtrip 1000
 
-tolerance  0.006 m
 accept    5 55
 roundtrip 1000
 
-tolerance  0.006 m
 accept    5 54
 roundtrip 1000
 
-tolerance  0.006 m
 accept    5 53
 roundtrip 1000
 
-tolerance  0.006 m
 accept    5.3876389 52.1561606
 roundtrip 1000
 
-tolerance  0.006 m
 accept    4 51
 roundtrip 1000
 
-tolerance  0.006 m
 accept    4 50
 roundtrip 1000
 
-tolerance  0.006 m
 accept    4 49
 roundtrip 1000
 
-tolerance  0.006 m
 accept    3.3137281 47.9752611
 roundtrip 1000
 
-tolerance  0.006 m
 accept    3 53
 roundtrip 1000
 
-tolerance  0.006 m
 accept    4 53
 roundtrip 1000
 
-tolerance  0.006 m
 accept    5 53
 roundtrip 1000
 
-tolerance  0.006 m
 accept    6 53
 roundtrip 1000
 
-tolerance  0.006 m
 accept    7 53
 roundtrip 1000
 
-tolerance  0.006 m
 accept    8 53
 roundtrip 1000
 
-tolerance  0.006 m
 accept    9 53
 roundtrip 1000
 
-tolerance  0.006 m
 accept    10 53
 roundtrip 1000
 
-tolerance  0.006 m
 accept    11 53
 roundtrip 1000
 
diff --git a/test/gigs/5110.gie b/test/gigs/5110.gie
index 031dfac9..2d26d1ab 100644
--- a/test/gigs/5110.gie
+++ b/test/gigs/5110.gie
@@ -4,6 +4,14 @@ Test 5110, Lambert Azimuthal Equal Area, v2-0_2011-06-28.
 
 --------------------------------------------------------------------------------
 
+# ETRS89
+<4258> +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs  <>
+
+# ETRS89 / LAEA Europe
+<3035> +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
+       +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs  <>
+
+
 <gie>
 
 --------------------------------------------------------------------------------
diff --git a/test/gigs/5111.1.gie b/test/gigs/5111.1.gie
index 2517b185..335e1d68 100644
--- a/test/gigs/5111.1.gie
+++ b/test/gigs/5111.1.gie
@@ -4,10 +4,16 @@ Test 5111 (part 1), Mercator (variant A), v2-0_2011-06-28.
 
 --------------------------------------------------------------------------------
 
+# Batavia
+<4211> +proj=longlat +ellps=bessel +towgs84=-377,681,-50,0,0,0,0 +no_defs  <>
+# Batavia / NEIEZ
+<3001> +proj=merc +lon_0=110 +k=0.997 +x_0=3900000 +y_0=900000 +ellps=bessel +towgs84=-377,681,-50,0,0,0,0 +units=m +no_defs  <>
+
+
 <gie>
 
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
+operation  +proj=pipeline towgs84=0,0,0
            +step +init=epsg:4211 +inv
            +step +init=epsg:3001
 --------------------------------------------------------------------------------
@@ -152,7 +158,7 @@ accept    -69.0 -2.0
 expect    -15964105.84 679490.65
 
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
+operation  +proj=pipeline towgs84=0,0,0
            +step +init=epsg:3001 +inv
            +step +init=epsg:4211
 --------------------------------------------------------------------------------
@@ -297,7 +303,7 @@ accept    -15964105.84 679490.65
 expect    -69.0 -2.0
 
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
+operation  +proj=pipeline towgs84=0,0,0
            +step +init=epsg:4211 +inv
            +step +init=epsg:3001
 --------------------------------------------------------------------------------
diff --git a/test/gigs/5112.gie b/test/gigs/5112.gie
index 3c4bb367..207632d2 100644
--- a/test/gigs/5112.gie
+++ b/test/gigs/5112.gie
@@ -7,77 +7,68 @@ Test 5112, Mercator (variant B),  v2-0_2011-06-28.
 <gie>
 
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
-           +step +init=epsg:4284 +inv
-           +step +init=epsg:3388
+operation  proj=pipeline
+           step init=epsg:4284 inv
+           step init=epsg:3388
+
+tolerance 50 mm
 --------------------------------------------------------------------------------
-tolerance  0.05 m
 accept    51.0 42.0
 expect    0.0 3819897.85
 
-tolerance  0.05 m
 accept    51.0 0.0
 expect    0.0 0.0
 
-tolerance  0.05 m
 accept    57.0 0.0
 expect    497112.88 0.0
 
-tolerance  0.05 m
 accept    54.0 20.5
 expect    248556.44 1724781.5
 
-tolerance  0.05 m
 accept    67.0 -41.0
 expect    1325634.35 -3709687.25
 
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
-           +step +init=epsg:3388 +inv
-           +step +init=epsg:4284
+operation  proj=pipeline
+           step init=epsg:3388 inv
+           step init=epsg:4284
+
+tolerance 50 mm
 --------------------------------------------------------------------------------
-tolerance  0.05 m
 accept    0.0 3819897.85
 expect    51.0 42.0
 
-tolerance  0.05 m
 accept    0.0 0.0
 expect    51.0 0.0
 
-tolerance  0.05 m
 accept    497112.88 0.0
 expect    57.0 0.0
 
-tolerance  0.05 m
 accept    248556.44 1724781.5
 expect    54.0 20.5
 
-tolerance  0.05 m
 accept    1325634.35 -3709687.25
 expect    67.0 -41.0
 
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
-           +step +init=epsg:4284 +inv
-           +step +init=epsg:3388
+operation proj=pipeline towgs84=0,0,0
+          step init=epsg:4284 inv
+          step init=epsg:3388
+
+tolerance 6 mm
 --------------------------------------------------------------------------------
-tolerance  0.006 m
 accept    51.0 42.0
 roundtrip 1000
 
-tolerance  0.006 m
 accept    51.0 0.0
 roundtrip 1000
 
-tolerance  0.006 m
 accept    57.0 0.0
 roundtrip 1000
 
-tolerance  0.006 m
 accept    54.0 20.5
 roundtrip 1000
 
-tolerance  0.006 m
 accept    67.0 -41.0
 roundtrip 1000
 
diff --git a/test/gigs/5201.gie b/test/gigs/5201.gie
index e7bdcb1e..532443a9 100644
--- a/test/gigs/5201.gie
+++ b/test/gigs/5201.gie
@@ -4,6 +4,12 @@ Test 5201, Geographic Geocentric conversions, v2.0_2011-09-28. (EPSG 4979 - WGS8
 
 --------------------------------------------------------------------------------
 
+# WGS 84
+<4978> +proj=geocent +datum=WGS84 +units=m +no_defs  <>
+# WGS 84
+<4326> +proj=longlat +datum=WGS84 +no_defs  <>
+
+
 <gie>
 
 --------------------------------------------------------------------------------
@@ -11,111 +17,111 @@ operation  +proj=pipeline
            +step +init=epsg:4978 +inv
            +step +init=epsg:4326
 --------------------------------------------------------------------------------
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01 m
 accept    -962479.5924 555687.8517 6260738.6526
 expect    150 80 1214.137
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01 m
 accept    -962297.0059 555582.4354 6259542.961
 expect    150 80 0
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -1598248.169 2768777.623 5501278.468
 expect    119.99524538 60.00475191 619.6317
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -1598023.169 2768387.912 5500499.045
 expect    119.9952447 60.00475258 -280.3683
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    2764210.4054 4787752.865 3170468.5199
 expect    60 30 189.569
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    2764128.3196 4787610.6883 3170373.7354
 expect    60 30 0
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    6377934.396 -112 434
 expect    -0.00100615 0.00392509 -202.5882
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    6374934.396 -112 434
 expect    -0.00100662 0.00392695 -3202.5881
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    6367934.396 -112 434
 expect    -0.00100773 0.00393129 -10202.5881
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    2764128.3196 -4787610.6883 -3170373.7354
 expect    -60 -30 0
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    2763900.3489 -4787215.8313 -3170110.4974
 expect    -60 -30 -526.476
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    2763880.8633 -4787182.0813 -3170087.9974
 expect    -60 -30 -571.476
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -1598023.169 -2768611.912 -5499631.045
 expect    -119.99323757 -59.99934884 -935.0995
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -1597798.169 -2768222.201 -5498851.622
 expect    -119.99323663 -59.99934874 -1835.0995
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -962297.0059 -555582.4354 -6259542.961
 expect    -150 -80 0
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -962150.945 -555498.1071 -6258586.4616
 expect    -150 -80 -971.255
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -961798.2951 -555294.5046 -6256277.0874
 expect    -150 -80 -3316.255
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -2187336.719 -112 5971017.093
 expect    -179.99706624 70.00490733 -223.6178
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -2904698.5551 -2904698.5551 4862789.0377
 expect    -135 50 0
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    371 -5783593.614 2679326.11
 expect    -89.99632465 25.00366329 -274.7286
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    6378137 0 0
 expect    0 0 0
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -4087095.478 2977467.559 -3875457.429
 expect    143.92649252 -37.65282217 737.7182
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -4085919.959 2976611.233 -3874335.274
 expect    143.92649211 -37.65282206 -1099.2288
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -4084000.165 2975212.729 -3872502.631
 expect    143.92649143 -37.65282187 -4099.2288
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -4079520.647 2971949.553 -3868226.465
 expect    143.92648984 -37.65282143 -11099.2288
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -2904698.5551 2904698.5551 -4862789.0377
 expect    135 -50 0
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -2187336.719 -112 -5970149.093
 expect    -179.99706624 -70.00224647 -1039.2896
 
@@ -124,111 +130,111 @@ operation  +proj=pipeline
            +step +init=epsg:4326 +inv
            +step +init=epsg:4978
 --------------------------------------------------------------------------------
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    150 80 1214.137
 expect    -962479.5924 555687.8517 6260738.6526
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    150 80 0
 expect    -962297.0059 555582.4354 6259542.961
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    119.99524538 60.00475191 619.6317
 expect    -1598248.169 2768777.623 5501278.468
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    119.9952447 60.00475258 -280.3683
 expect    -1598023.169 2768387.912 5500499.045
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    60 30 189.569
 expect    2764210.4054 4787752.865 3170468.5199
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    60 30 0
 expect    2764128.3196 4787610.6883 3170373.7354
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -0.00100615 0.00392509 -202.5882
 expect    6377934.396 -112 434
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -0.00100662 0.00392695 -3202.5881
 expect    6374934.396 -112 434
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -0.00100773 0.00393129 -10202.5881
 expect    6367934.396 -112 434
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -60 -30 0
 expect    2764128.3196 -4787610.6883 -3170373.7354
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -60 -30 -526.476
 expect    2763900.3489 -4787215.8313 -3170110.4974
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -60 -30 -571.476
 expect    2763880.8633 -4787182.0813 -3170087.9974
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -119.99323757 -59.99934884 -935.0995
 expect    -1598023.169 -2768611.912 -5499631.045
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -119.99323663 -59.99934874 -1835.0995
 expect    -1597798.169 -2768222.201 -5498851.622
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -150 -80 0
 expect    -962297.0059 -555582.4354 -6259542.961
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -150 -80 -971.255
 expect    -962150.945 -555498.1071 -6258586.4616
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -150 -80 -3316.255
 expect    -961798.2951 -555294.5046 -6256277.0874
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -179.99706624 70.00490733 -223.6178
 expect    -2187336.719 -112 5971017.093
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -135 50 0
 expect    -2904698.5551 -2904698.5551 4862789.0377
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -89.99632465 25.00366329 -274.7286
 expect    371 -5783593.614 2679326.11
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    0 0 0
 expect    6378137 0 0
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    143.92649252 -37.65282217 737.7182
 expect    -4087095.478 2977467.559 -3875457.429
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    143.92649211 -37.65282206 -1099.2288
 expect    -4085919.959 2976611.233 -3874335.274
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    143.92649143 -37.65282187 -4099.2288
 expect    -4084000.165 2975212.729 -3872502.631
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    143.92648984 -37.65282143 -11099.2288
 expect    -4079520.647 2971949.553 -3868226.465
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    135 -50 0
 expect    -2904698.5551 2904698.5551 -4862789.0377
 
-tolerance  [8.333333333333333e-08, 8.333333333333333e-08, 0.01] m
+tolerance  0.01  m
 accept    -179.99706624 -70.00224647 -1039.2896
 expect    -2187336.719 -112 -5970149.093
 
@@ -237,111 +243,111 @@ operation  +proj=pipeline
            +step +init=epsg:4978 +inv
            +step +init=epsg:4326
 --------------------------------------------------------------------------------
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -962479.5924 555687.8517 6260738.6526
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -962297.0059 555582.4354 6259542.961
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -1598248.169 2768777.623 5501278.468
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -1598023.169 2768387.912 5500499.045
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    2764210.4054 4787752.865 3170468.5199
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    2764128.3196 4787610.6883 3170373.7354
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    6377934.396 -112 434
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    6374934.396 -112 434
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    6367934.396 -112 434
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    2764128.3196 -4787610.6883 -3170373.7354
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    2763900.3489 -4787215.8313 -3170110.4974
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    2763880.8633 -4787182.0813 -3170087.9974
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -1598023.169 -2768611.912 -5499631.045
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -1597798.169 -2768222.201 -5498851.622
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -962297.0059 -555582.4354 -6259542.961
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -962150.945 -555498.1071 -6258586.4616
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -961798.2951 -555294.5046 -6256277.0874
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -2187336.719 -112 5971017.093
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -2904698.5551 -2904698.5551 4862789.0377
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    371 -5783593.614 2679326.11
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    6378137 0 0
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -4087095.478 2977467.559 -3875457.429
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -4085919.959 2976611.233 -3874335.274
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -4084000.165 2975212.729 -3872502.631
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -4079520.647 2971949.553 -3868226.465
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -2904698.5551 2904698.5551 -4862789.0377
 roundtrip 1000
 
-tolerance  [5.555555555555556e-08, 5.555555555555556e-08, 0.006] m
+tolerance  0.01  m
 accept    -2187336.719 -112 -5970149.093
 roundtrip 1000
 
diff --git a/test/gigs/5208.gie b/test/gigs/5208.gie
index 5d94474a..25a10b96 100644
--- a/test/gigs/5208.gie
+++ b/test/gigs/5208.gie
@@ -2,8 +2,24 @@
 
 Test 5208, Longitude Rotation, v2.0_2011-06-28.
 
+The test tolerance is 0.01". Since gie can only use linear tolerances we
+convert that to an approximate liniar distance instead, by multiplying with
+111km:
+
+    0.01" * 111 km = 2.777777778-7 * 111000 m = 0.03 m
+
+To be on the safe side we, use 0.01 m as the tolerance.
+
 --------------------------------------------------------------------------------
 
+# NTF
+<4275> +proj=longlat +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +no_defs  <>
+
+# NTF (Paris)
+<4807> +proj=longlat +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +no_defs  <>
+
+
+
 <gie>
 
 --------------------------------------------------------------------------------
@@ -11,59 +27,59 @@ operation  +proj=pipeline
            +step +init=epsg:4275 +inv
            +step +init=epsg:4807
 --------------------------------------------------------------------------------
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    5 58
 expect    2.66277083 58
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    5 56
 expect    2.66277083 56
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    5 55
 expect    2.66277083 55
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    5 53
 expect    2.66277083 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    4 51
 expect    1.66277083 51
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    4 49
 expect    1.66277083 49
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    2.33722917 46.8
 expect    0 46.8
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    3 53
 expect    0.66277083 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    4 53
 expect    1.66277083 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    6 53
 expect    3.66277083 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    7 53
 expect    4.66277083 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    9 53
 expect    6.66277083 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    10 53
 expect    7.66277083 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    11 53
 expect    8.66277083 53
 
@@ -72,59 +88,59 @@ operation  +proj=pipeline
            +step +init=epsg:4807 +inv
            +step +init=epsg:4275
 --------------------------------------------------------------------------------
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    2.66277083 58
 expect    5 58
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    2.66277083 56
 expect    5 56
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    2.66277083 55
 expect    5 55
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    2.66277083 53
 expect    5 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    1.66277083 51
 expect    4 51
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    1.66277083 49
 expect    4 49
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    0 46.8
 expect    2.33722917 46.8
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    0.66277083 53
 expect    3 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    1.66277083 53
 expect    4 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    3.66277083 53
 expect    6 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    4.66277083 53
 expect    7 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    6.66277083 53
 expect    9 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    7.66277083 53
 expect    10 53
 
-tolerance  2.77777777778e-07 m
+tolerance 0.01 m
 accept    8.66277083 53
 expect    11 53
 
@@ -133,59 +149,59 @@ operation  +proj=pipeline
            +step +init=epsg:4275 +inv
            +step +init=epsg:4807
 --------------------------------------------------------------------------------
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    5 58
 roundtrip 1000
 
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    5 56
 roundtrip 1000
 
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    5 55
 roundtrip 1000
 
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    5 53
 roundtrip 1000
 
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    4 51
 roundtrip 1000
 
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    4 49
 roundtrip 1000
 
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    2.33722917 46.8
 roundtrip 1000
 
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    3 53
 roundtrip 1000
 
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    4 53
 roundtrip 1000
 
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    6 53
 roundtrip 1000
 
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    7 53
 roundtrip 1000
 
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    9 53
 roundtrip 1000
 
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    10 53
 roundtrip 1000
 
-tolerance  5.55555555556e-08 m
+tolerance 0.01 m
 accept    11 53
 roundtrip 1000