From 8cfc81380617ff4a17a06a97635f77c5e9ed7d5b Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 19 Aug 2018 14:13:06 +0200 Subject: 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 +++++++++++++++++++++----------------------- 1 file changed, 21 insertions(+), 23 deletions(-) (limited to 'src') diff --git a/src/pj_inv.c b/src/pj_inv.c index d1a02bca..285d8b85 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) -- cgit v1.2.3 From bb8937e047afef45bbb790613106165c71746c49 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 19 Aug 2018 14:18:40 +0200 Subject: 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(-) (limited to 'src') 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 285d8b85..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 && 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: -- cgit v1.2.3 From f05de97c210ece96da6f230465e483e0cfff2c7d Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 19 Aug 2018 14:31:01 +0200 Subject: pj_transform / cs2cs: honour +geoc flag --- src/pj_transform.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src') 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. */ -- cgit v1.2.3