Merge pull request #1093 from rouault/geoc_flag_fix

Assorted fixes related to +geoc flag handing

6 files changed, 68 insertions, 114 deletions

diff --git a/nad/testvarious b/nad/testvarious
index 30d1402f..c0a8ff29 100755
--- a/nad/testvarious
+++ b/nad/testvarious
@@ -139,20 +139,29 @@ $EXE +proj=geocent +datum=WGS84  \
 0.00    -0.00 6356752.31
 EOF
 #
-#echo "#############################################################" >> ${OUT}
-#echo Test conversion between geocentric latlong and geodetic latlong >> ${OUT}
-#
-# The +geoc flag does not currently work with pj_transform() so this is
-# disabled.
-#
-#$EXE +proj=latlong +datum=WGS84 +geoc \
-# +to +proj=latlong +datum=WGS84  \
-# -E >>${OUT} <<EOF
-#0d00'00.000"W 0d00'00.000"N 0.0
-#79d00'00.000"W 45d00'00.000"N 0.0
-#12d00'00.000"W 45d00'00.000"N 0.0
-#0d00'00.000"W 90d00'00.000"N 0.0
-#EOF
+echo "#############################################################" >> ${OUT}
+echo Test conversion from geocentric latlong to geodetic latlong >> ${OUT}
+#
+$EXE +proj=latlong +datum=WGS84 +geoc \
+ +to +proj=latlong +datum=WGS84  \
+ -E >>${OUT} <<EOF
+0d00'00.000"W 0d00'00.000"N 0.0
+79d00'00.000"W 45d00'00.000"N 0.0
+12d00'00.000"W 45d00'00.000"N 0.0
+0d00'00.000"W 90d00'00.000"N 0.0
+EOF
+#
+echo "#############################################################" >> ${OUT}
+echo Test conversion from geodetic latlong to geocentric latlong >> ${OUT}
+#
+$EXE +proj=latlong +datum=WGS84 \
+ +to +proj=latlong +datum=WGS84 +geoc  \
+ -E >>${OUT} <<EOF
+0d00'00.000"W 0d00'00.000"N 0.0
+79d00'00.000"W 44d48'27.276"N 0.000
+12d00'00.000"W 44d48'27.276"N 0.0
+0d00'00.000"W 90d00'00.000"N 0.0
+EOF
 #
 echo "##############################################################" >> ${OUT}
 echo "Test stere projection (re: win32 ticket 12)" >> ${OUT}
diff --git a/nad/tv_out.dist b/nad/tv_out.dist
index 743eeabe..f04e2c35 100644
--- a/nad/tv_out.dist
+++ b/nad/tv_out.dist
@@ -42,6 +42,18 @@ Test geocentric x/y/z consumption.
 6378147.00      -0.00 0.00	0dE	0dN 10.000
 861996.98       -4434590.01 4487348.41	79dW	45dN 0.001
 0.00    -0.00 6356752.31	0dE	90dN -0.004
+#############################################################
+Test conversion from geocentric latlong to geodetic latlong
+0d00'00.000"W 0d00'00.000"N 0.0	0dE	0dN 0.000
+79d00'00.000"W 45d00'00.000"N 0.0	79dW	44d48'27.276"N 0.000
+12d00'00.000"W 45d00'00.000"N 0.0	12dW	44d48'27.276"N 0.000
+0d00'00.000"W 90d00'00.000"N 0.0	0dE	90dN 0.000
+#############################################################
+Test conversion from geodetic latlong to geocentric latlong
+0d00'00.000"W 0d00'00.000"N 0.0	0dE	0dN 0.000
+79d00'00.000"W 44d48'27.276"N 0.000	79dW	45dN 0.000
+12d00'00.000"W 44d48'27.276"N 0.0	12dW	45dN 0.000
+0d00'00.000"W 90d00'00.000"N 0.0	0dE	90dN 0.000
 ##############################################################
 Test stere projection (re: win32 ticket 12)
 105 40	5577808.93	1494569.40 0.00
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index 2ed4c469..fa9cbbc8 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -137,39 +137,7 @@ static PJ_COORD fwd_finalize (PJ *P, PJ_COORD coo) {
         break;
 
     case PJ_IO_UNITS_ANGULAR:
-        if (INPUT_UNITS==PJ_IO_UNITS_ANGULAR)
-            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_FWD, coo); /* Go orthometric from geometric */
-        if (coo.lp.lam==HUGE_VAL)
-            return coo;
-        if (P->hgridshift)
-            coo = proj_trans (P->hgridshift, PJ_INV, coo);
-        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 */
-            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)
-            return coo;
-
-        /* If input latitude was geocentrical, convert back to geocentrical */
-        if (P->geoc)
-            coo = proj_geocentric_latitude (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);
+        break;
     }
 
     if (P->axisswap)
diff --git a/src/pj_inv.c b/src/pj_inv.c
index d1a02bca..ab9cd78f 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -50,49 +50,6 @@ static PJ_COORD inv_prepare (PJ *P, PJ_COORD coo) {
     if (P->axisswap)
         coo = proj_trans (P->axisswap, PJ_INV, coo);
 
-    /* Check validity of angular input coordinates */
-    if (INPUT_UNITS==PJ_IO_UNITS_ANGULAR) {
-        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_geocentric_latitude (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 || (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 (P->vgridshift)
-            coo = proj_trans (P->vgridshift, PJ_INV, coo); /* Go geometric from orthometric */
-        return coo;
-    }
-
     /* Handle remaining possible input types */
     switch (INPUT_UNITS) {
     case PJ_IO_UNITS_WHATEVER:
@@ -143,29 +100,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/src/pj_transform.c b/src/pj_transform.c
index b459e4a2..80294f6d 100644
--- a/src/pj_transform.c
+++ b/src/pj_transform.c
@@ -190,7 +190,7 @@ static int geographic_to_projected (PJ *P, long n, int dist, double *x, double *
     long i;
 
     /* Nothing to do? */
-    if (P->is_latlong)
+    if (P->is_latlong && !P->geoc)
         return 0;
     if (P->is_geocent)
         return 0;
@@ -290,7 +290,7 @@ static int projected_to_geographic (PJ *P, long n, int dist, double *x, double *
     long i;
 
     /* Nothing to do? */
-    if (P->is_latlong)
+    if (P->is_latlong && !P->geoc)
         return 0;
 
     /* Check first if projection is invertible. */
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index afd30be8..fd1ff10f 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -376,6 +376,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