Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Assorted fixes related to +geoc flag handing #1093

Merged
merged 3 commits into from
Aug 24, 2018
Merged

Assorted fixes related to +geoc flag handing #1093

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8cfc813
pj_inv(): fix inverse of +proj=longlat +geoc
rouault Aug 19, 2018
bb8937e
Remove dead code in fwd_finalize() and inv_prepare()
rouault Aug 19, 2018
f05de97
pj_transform / cs2cs: honour +geoc flag
rouault Aug 19, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 23 additions & 14 deletions nad/testvarious
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down
12 changes: 12 additions & 0 deletions nad/tv_out.dist
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
34 changes: 1 addition & 33 deletions src/pj_fwd.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
85 changes: 20 additions & 65 deletions src/pj_inv.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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)
Expand Down
4 changes: 2 additions & 2 deletions src/pj_transform.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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. */
Expand Down
10 changes: 10 additions & 0 deletions test/gie/more_builtins.gie
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down