Fix numerous bugs in the cs2cs emulation

3 files changed, 27 insertions, 15 deletions

diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index 66b86aab..dbb7748e 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -73,9 +73,11 @@ static PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD coo) {
             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->helmert,    PJ_INV, coo); /* Step into local frame */
             coo = proj_trans (P->cart,       PJ_INV, coo); /* Go back to angular using local ellps */
         }
+        if (coo.lp.lam==HUGE_VAL)
+            return coo;
         if (P->vgridshift)
             coo = proj_trans (P->vgridshift, PJ_FWD, coo);
 
@@ -92,7 +94,7 @@ static PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD 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 proj_trans (P->helmert, PJ_INV, coo);
     return coo;
 }
 
@@ -140,13 +142,17 @@ static PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo) {
 
         if (P->vgridshift)
             coo = proj_trans (P->vgridshift, PJ_INV, coo);
+        if (coo.lp.lam==HUGE_VAL)
+            return coo;
         if (P->hgridshift)
-            coo = proj_trans (P->hgridshift, PJ_FWD, coo);
+            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->helmert,    PJ_INV, coo); /* Step into local frame */
             coo = proj_trans (P->cart,       PJ_INV, coo); /* Go back to angular using local ellps */
         }
+        if (coo.lp.lam==HUGE_VAL)
+            return coo;
 
         /* If input latitude was geocentrical, convert back to geocentrical */
         if (P->geoc)
diff --git a/src/pj_inv.c b/src/pj_inv.c
index 5e44eb34..7b47f5d7 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -80,11 +80,13 @@ static PJ_COORD pj_inv_prepare (PJ *P, PJ_COORD coo) {
             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->helmert,    PJ_FWD, coo); /* Step into WGS84 */
             coo = proj_trans (P->cart_wgs84, PJ_INV, coo); /* Go back to angular using WGS84 ellps */
         }
+        if (coo.lp.lam==HUGE_VAL)
+            return coo;
         if (P->vgridshift)
-            coo = proj_trans (P->vgridshift, PJ_INV, coo);
+            coo = proj_trans (P->vgridshift, PJ_FWD, coo);
         return coo;
     }
 
@@ -99,9 +101,8 @@ static PJ_COORD pj_inv_prepare (PJ *P, PJ_COORD coo) {
         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) {
+        if (P->is_geocent)
             coo = proj_trans (P->cart, PJ_INV, coo);
-        }
 
         return coo;
 
@@ -148,14 +149,18 @@ static PJ_COORD pj_inv_finalize (PJ *P, PJ_COORD coo) {
                 coo.lpz.lam = adjlon(coo.lpz.lam);
 
             if (P->vgridshift)
-                coo = proj_trans (P->vgridshift, PJ_INV, coo);
+                coo = proj_trans (P->vgridshift, PJ_FWD, coo);
+            if (coo.lp.lam==HUGE_VAL)
+                return 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->helmert,    PJ_FWD, coo); /* Step into WGS84 */
                 coo = proj_trans (P->cart_wgs84, PJ_INV, coo); /* Go back to angular using WGS84 ellps */
             }
+            if (coo.lp.lam==HUGE_VAL)
+                return coo;
         }
 
         /* If input latitude was geocentrical, convert back to geocentrical */
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 42e3cbf9..8f6305a4 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -473,12 +473,13 @@ invocators can emulate the behaviour of pj_transform and the cs2cs app.
         def = malloc (100+n);
         if (0==def)
             return 0;
-        sprintf (def, "break_cs2cs_recursion     proj=helmert %s", s);
+        sprintf (def, "break_cs2cs_recursion     proj=helmert %s transpose", s);
         Q = proj_create (P->ctx, def);
+        pj_inherit_ellipsoid_def (P, Q);
         free (def);
         if (0==Q)
             return 0;
-        P->helmert = skip_prep_fin(Q);
+        P->helmert = skip_prep_fin (Q);
 
         break;
     }
@@ -491,14 +492,14 @@ invocators can emulate the behaviour of pj_transform and the cs2cs app.
         Q = proj_create (P->ctx, def);
         if (0==Q)
             return 0;
-        pj_inherit_ellipsoid_def(P, Q);
-        P->cart = skip_prep_fin(Q);
+        pj_inherit_ellipsoid_def (P, Q);
+        P->cart = skip_prep_fin (Q);
 
         sprintf (def, "break_cs2cs_recursion     proj=cart  ellps=WGS84");
         Q = proj_create (P->ctx, def);
         if (0==Q)
             return 0;
-        P->cart_wgs84 = skip_prep_fin(Q);
+        P->cart_wgs84 = skip_prep_fin (Q);
     }
 
     return 1;