Take into account +towgs84=0,0,0 in pipeline to still imply geodetic-…

…>cartesian->geodetic (fixes #881)
OSGeo · Mar 21, 2018 · c49fcf7 · c49fcf7
1 parent 11396b4
commit c49fcf7
Show file tree

Hide file tree

Showing 5 changed files with 36 additions and 11 deletions.
diff --git 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
@@ -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
@@ -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);

diff --git a/test/gie/4D-API_cs2cs-style.gie b/test/gie/4D-API_cs2cs-style.gie
@@ -227,4 +227,18 @@ accept      -100 40 0
 expect      -100.0004058367   40.0000058947        0.0000
 -------------------------------------------------------------------------------
 
+-------------------------------------------------------------------------------
+Test that +towgs84=0,0,0 parameter is handled as still implying cart
+transformation
+-------------------------------------------------------------------------------
+operation   +proj=pipeline
+            +step +proj=utm +zone=11 +ellps=clrk66 +towgs84=0,0,0 +inv
+            +step +proj=utm +zone=11 +datum=WGS84
+
+-------------------------------------------------------------------------------
+tolerance   20 cm
+accept      440720 3751320 0
+expect      440719.958709357 3751294.2109841 -4.44340920541435
+-------------------------------------------------------------------------------
+
 </gie>
diff --git a/test/gigs/5111.1.gie b/test/gigs/5111.1.gie
@@ -21,7 +21,7 @@ tolerance  0.05 m
 accept    100.0876483 77.6534822
 expect    2800000.0 15000000.0
 
-tolerance  0.05 m
+tolerance  0.055 m
 accept    100.0876483 73.1442856
 expect    2800000.0 13000000.0