Merge pull request #882 from rouault/fix_881

Take into account +towgs84=0,0,0 in pipeline to still imply geodetic->cartesian->geodetic (fixes #881)

3 files changed, 21 insertions, 10 deletions

diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index ef05962d..2a064e58 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -71,9 +71,10 @@ static PJ_COORD fwd_prepare (PJ *P, PJ_COORD coo) {
 
         if (P->hgridshift)
             coo = proj_trans (P->hgridshift, PJ_INV, coo);
-        else if (P->helmert) {
+        else if (P->helmert || (P->cart_wgs84 != 0 && P->cart != 0)) {
             coo = proj_trans (P->cart_wgs84, PJ_FWD, coo); /* Go cartesian in WGS84 frame */
-            coo = proj_trans (P->helmert,    PJ_INV, coo); /* Step into local frame */
+            if( P->helmert )
+                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)
@@ -146,9 +147,10 @@ static PJ_COORD fwd_finalize (PJ *P, PJ_COORD coo) {
             return coo;
         if (P->hgridshift)
             coo = proj_trans (P->hgridshift, PJ_INV, coo);
-        else if (P->helmert) {
+        else if (P->helmert || (P->cart_wgs84 != 0 && P->cart != 0)) {
             coo = proj_trans (P->cart_wgs84, PJ_FWD, coo); /* Go cartesian in WGS84 frame */
-            coo = proj_trans (P->helmert,    PJ_INV, coo); /* Step into local frame */
+            if( P->helmert )
+                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)
diff --git a/src/pj_inv.c b/src/pj_inv.c
index 4c2266fb..f66fd7d1 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -78,9 +78,10 @@ static PJ_COORD inv_prepare (PJ *P, PJ_COORD coo) {
 
         if (P->hgridshift)
             coo = proj_trans (P->hgridshift, PJ_FWD, coo);
-        else if (P->helmert) {
+        else if (P->helmert || (P->cart_wgs84 != 0 && P->cart != 0)) {
             coo = proj_trans (P->cart,       PJ_FWD, coo); /* Go cartesian in local frame */
-            coo = proj_trans (P->helmert,    PJ_FWD, coo); /* Step into WGS84 */
+            if( P->helmert )
+                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)
@@ -154,9 +155,10 @@ static PJ_COORD inv_finalize (PJ *P, PJ_COORD coo) {
                 return coo;
             if (P->hgridshift)
                 coo = proj_trans (P->hgridshift, PJ_FWD, coo);
-            else if (P->helmert) {
+            else if (P->helmert || (P->cart_wgs84 != 0 && P->cart != 0)) {
                 coo = proj_trans (P->cart,       PJ_FWD, coo); /* Go cartesian in local frame */
-                coo = proj_trans (P->helmert,    PJ_FWD, coo); /* Step into WGS84 */
+                if( P->helmert )
+                    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)
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 05f030e5..11a56ac0 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -419,6 +419,7 @@ Returns 1 on success, 0 on failure
 **************************************************************************************/
     PJ *Q;
     paralist *p;
+    int do_cart = 0;
     if (0==P)
         return 0;
 
@@ -481,8 +482,14 @@ Returns 1 on success, 0 on failure
         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])
+        if (0==d[0] && 0==d[1] && 0==d[2] && 0==d[3] && 0==d[4] && 0==d[5] && 0==d[6]) {
+            /* If the current ellipsoid is not WGS84, then make sure the */
+            /* change in ellipsoid is still done. */
+            if (!(fabs(P->a - 6378137.0) < 1e-8 && fabs(P->f - 1./ 298.257223563) < 1e-15)) {
+                do_cart = 1;
+            }
             break;
+        }
 
         if (n <= 8) /* 8==strlen ("towgs84=") */
             return 0;
@@ -503,7 +510,7 @@ Returns 1 on success, 0 on failure
 
     /* 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) {
+    if (P->is_geocent || P->helmert || do_cart) {
         char def[150];
         sprintf (def, "break_cs2cs_recursion     proj=cart   a=%40.20g  f=%40.20g", P->a, P->f);
         Q = proj_create (P->ctx, def);