Introduce compatibility for cs2cs-style proj-strings into the 4D API.

Parameters such as towgs84, nadgrids and geoidgrids was previously only
handled by pj_transform(). This commit add a compatibility layer in
proj_create() by calling the pj_cs2cs_emulation_setup() function. This
function sets up a handful of predefined transformation objects on the
PJ object that is being created. Each of these transformation objects
are related to the cs2cs-style parameters we are trying to emulate in
the 4D API. That is, if the +towgs84 parameters is used we create
P->helmert with the parameters specified in +towgs84. Similarly for
+axis, +nadgrids and +geoidgrids. When these transformation objects
exists we use them in the prepare and finalize functions in pj_fwd/
pj_inv. If no cs2cs-style parametes are specified we skip those
parts of the prepare and finalizing steps.

Co-authored-by:Thomas Knudsen <thokn@sdfe.dk>
Co-authored-by:Kristian Evers <kristianevers@gmail.com>

11 files changed, 421 insertions, 90 deletions

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/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..2293cf2e 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_RADIANS) {
         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_RADIANS:
+        if (INPUT_UNITS==PJ_IO_UNITS_RADIANS)
+            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..e676b735 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_RADIANS;
+    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..02a2ca53 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_RADIANS) {
+        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_RADIANS 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_RADIANS) {
+
+        if (INPUT_UNITS!=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 (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.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..fa24c430 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -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..2c6cbee2 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -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..aeaf09ae 100644
--- a/src/projects.h
+++ b/src/projects.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;