From 8cfc81380617ff4a17a06a97635f77c5e9ed7d5b Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 19 Aug 2018 14:13:06 +0200 Subject: [PATCH 1/3] 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 --- src/pj_inv.c | 44 ++++++++++++++++++-------------------- test/gie/more_builtins.gie | 10 +++++++++ 2 files changed, 31 insertions(+), 23 deletions(-) diff --git a/src/pj_inv.c b/src/pj_inv.c index d1a02bca9a..285d8b85f3 100644 --- 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 index 2908fd633e..15964a37ad 100644 --- 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 From bb8937e047afef45bbb790613106165c71746c49 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 19 Aug 2018 14:18:40 +0200 Subject: [PATCH 2/3] Remove dead code in fwd_finalize() and inv_prepare() Those functions contain code specific of input != angular and output = angular which to the best of my knowledge never happens in PROJ. Looking at that code, I feel there are a number of errors in it, and anyway removing it shows absolutely no change in the test suite, which shows it is unused. I've also added exit(1) to verify that those code paths are never taken. This was found while investigating the fix for 8cfc81380617ff4a17a06a97635f77c5e9ed7d5b --- src/pj_fwd.c | 34 +--------------------------------- src/pj_inv.c | 43 ------------------------------------------- 2 files changed, 1 insertion(+), 76 deletions(-) diff --git a/src/pj_fwd.c b/src/pj_fwd.c index 2ed4c46994..fa9cbbc833 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 285d8b85f3..ab9cd78ffc 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 && OUTPUT_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: From f05de97c210ece96da6f230465e483e0cfff2c7d Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 19 Aug 2018 14:31:01 +0200 Subject: [PATCH 3/3] pj_transform / cs2cs: honour +geoc flag --- nad/testvarious | 37 +++++++++++++++++++++++-------------- nad/tv_out.dist | 12 ++++++++++++ src/pj_transform.c | 4 ++-- 3 files changed, 37 insertions(+), 16 deletions(-) diff --git a/nad/testvarious b/nad/testvarious index 30d1402fd3..c0a8ff2954 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} <> ${OUT} +echo Test conversion from geocentric latlong to geodetic latlong >> ${OUT} +# +$EXE +proj=latlong +datum=WGS84 +geoc \ + +to +proj=latlong +datum=WGS84 \ + -E >>${OUT} <> ${OUT} +echo Test conversion from geodetic latlong to geocentric latlong >> ${OUT} +# +$EXE +proj=latlong +datum=WGS84 \ + +to +proj=latlong +datum=WGS84 +geoc \ + -E >>${OUT} <> ${OUT} echo "Test stere projection (re: win32 ticket 12)" >> ${OUT} diff --git a/nad/tv_out.dist b/nad/tv_out.dist index 743eeabe16..f04e2c3507 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_transform.c b/src/pj_transform.c index b459e4a2b2..80294f6dcd 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. */