pj_inv(): fix inverse of +proj=longlat +geoc

There is a regression in PROJ 5 regarding the handling of the +geoc flag, specific to the case of +proj=longlat, and the inverse transformation, which makes it a no-op: echo "12 55 0 " | src/cct +proj=pipeline +step +proj=longlat +ellps=GRS80 +geoc +inv 12.0000000000 55.0000000000 0.0000 inf With this fix, we now get: echo "12 55 0 " | src/cct +proj=pipeline +step +proj=longlat +ellps=GRS80 +geoc +inv 12.0000000000 54.8189733083 0.0000 inf The fix consists in making inv_prepare() do symetrically as fwd_finalize(), ie skip a number of transforms when left and right units are angular, and in inv_finalize() apply the (OUTPUT_UNITS==PJ_IO_UNITS_ANGULAR) processing even if INPUT_UNITS == PJ_IO_UNITS_ANGULAR
OSGeo · Aug 19, 2018 · 8cfc813 · 8cfc813
1 parent 121deeb
commit 8cfc813
Show file tree

Hide file tree

Showing 2 changed files with 31 additions and 23 deletions.
diff --git a/src/pj_inv.c b/src/pj_inv.c
@@ -51,7 +51,7 @@ static PJ_COORD inv_prepare (PJ *P, PJ_COORD coo) {
         coo = proj_trans (P->axisswap, PJ_INV, coo);
 
     /* Check validity of angular input coordinates */
-    if (INPUT_UNITS==PJ_IO_UNITS_ANGULAR) {
+    if (INPUT_UNITS==PJ_IO_UNITS_ANGULAR && OUTPUT_UNITS!=PJ_IO_UNITS_ANGULAR) {
         double t;
 
         /* check for latitude or longitude over-range */
@@ -143,29 +143,27 @@ static PJ_COORD inv_finalize (PJ *P, PJ_COORD coo) {
 
     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); /* Go geometric from orthometric */
-            if (coo.lp.lam==HUGE_VAL)
-                return coo;
-            if (P->hgridshift)
-                coo = proj_trans (P->hgridshift, PJ_FWD, coo);
-            else if (P->helmert || (P->cart_wgs84 != 0 && P->cart != 0)) {
-                coo = proj_trans (P->cart,       PJ_FWD, coo); /* Go cartesian in local frame */
-                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)
-                return 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);
+
+        if (P->vgridshift)
+            coo = proj_trans (P->vgridshift, PJ_INV, coo); /* Go geometric from orthometric */
+        if (coo.lp.lam==HUGE_VAL)
+            return coo;
+        if (P->hgridshift)
+            coo = proj_trans (P->hgridshift, PJ_FWD, coo);
+        else if (P->helmert || (P->cart_wgs84 != 0 && P->cart != 0)) {
+            coo = proj_trans (P->cart,       PJ_FWD, coo); /* Go cartesian in local frame */
+            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)
+            return coo;
 
         /* If input latitude was geocentrical, convert back to geocentrical */
         if (P->geoc)

diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
@@ -356,6 +356,16 @@ accept     12   89.99999999999      0 0
 expect     12   89.999999999989996  0 0
 -------------------------------------------------------------------------------
 
+-------------------------------------------------------------------------------
+geocentric latitude using old +geoc flag
+-------------------------------------------------------------------------------
+operation  proj=pipeline step proj=longlat ellps=GRS80 geoc inv
+accept     12   55                  0 0
+expect     12   54.818973308324573  0 0
+roundtrip  1
+-------------------------------------------------------------------------------
+
+
 
 -------------------------------------------------------------------------------
 some less used options