diff options
| author | Even Rouault <even.rouault@spatialys.com> | 2018-08-19 14:13:06 +0200 |
|---|---|---|
| committer | Even Rouault <even.rouault@spatialys.com> | 2018-08-19 14:15:28 +0200 |
| commit | 8cfc81380617ff4a17a06a97635f77c5e9ed7d5b (patch) | |
| tree | b4f8ffdd48eeff98df545269723293804da5b57c /src | |
| parent | 121deeb3ff5c44aa1ae67944be0ad8c3f38d3986 (diff) | |
| download | PROJ-8cfc81380617ff4a17a06a97635f77c5e9ed7d5b.tar.gz PROJ-8cfc81380617ff4a17a06a97635f77c5e9ed7d5b.zip | |
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
Diffstat (limited to 'src')
| -rw-r--r-- | src/pj_inv.c | 44 |
1 files changed, 21 insertions, 23 deletions
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) |