From d33de1c5bb193740796092ffd5641d11b731a4e0 Mon Sep 17 00:00:00 2001 From: Frank Warmerdam Date: Mon, 27 Feb 2012 07:56:32 +0000 Subject: added +sweep for +proj=geos (#146) git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@2176 4e78687f-474d-0410-85f9-8d5e500ac6b2 --- ChangeLog | 3 +++ nad/testvarious | 66 +++++++++++++++++++++++++++++++++++++++++++++++++++ nad/tv_out.dist | 42 +++++++++++++++++++++++++++++++++ src/PJ_geos.c | 70 +++++++++++++++++++++++++++++++++++++++++++++++-------- src/pj_strerrno.c | 1 + 5 files changed, 172 insertions(+), 10 deletions(-) diff --git a/ChangeLog b/ChangeLog index 05ff7f32..a80a9111 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,8 @@ 2012-02-26 Frank Warmerdam + * src/PJ_geos.c, nad/testvarious: Added GEOS +sweep and add GEOS + to the test suite (#146) + * nad/CH: added swiss datum related definitions from strk (#145) * src/Makefile.am, src/mutltistresstest.c: provide for building diff --git a/nad/testvarious b/nad/testvarious index 3505dd16..7a753404 100755 --- a/nad/testvarious +++ b/nad/testvarious @@ -295,6 +295,72 @@ $EXE +proj=latlong +a=1 +lon_0=0 +ellps=WGS84 \ 0 0.7853981633974483 -1.5707963267948966 0 EOF +echo "##############################################################" >> ${OUT} +echo "Test geos projection" >> ${OUT} +echo "Test geos on a sphere" >> ${OUT} +$EXE +proj=latlong +ellps=sphere \ + +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=sphere -E >>${OUT} <> ${OUT} +$EXE +proj=latlong +ellps=sphere \ + +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=WGS84 -E >>${OUT} <> ${OUT} +$EXE +proj=latlong +ellps=sphere \ + +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=sphere -I -E >>${OUT} <> ${OUT} +$EXE +proj=latlong +ellps=sphere \ + +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=WGS84 -I -E >>${OUT} <> ${OUT} +$EXE +proj=latlong +ellps=sphere \ + +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=sphere +sweep=y -E >>${OUT} <> ${OUT} +$EXE +proj=latlong +ellps=sphere \ + +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=WGS84 +sweep=y -E >>${OUT} <> ${OUT} +$EXE +proj=latlong +ellps=sphere \ + +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=sphere +sweep=y -I -E >>${OUT} <> ${OUT} +$EXE +proj=latlong +ellps=sphere \ + +to +proj=geos +h=35785831.0 +lon_0=0 +ellps=WGS84 +sweep=y -I -E >>${OUT} < @@ -54,8 +57,16 @@ FORWARD(s_forward); /* spheroid */ if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz) < 0.) F_ERROR; /* Calculation based on view angles from satellite.*/ tmp = P->radius_g - Vx; - xy.x = P->radius_g_1 * atan(Vy / tmp); - xy.y = P->radius_g_1 * atan(Vz / hypot(Vy, tmp)); + if(P->flip_axis) + { + xy.x = P->radius_g_1 * atan(Vy / hypot(Vz, tmp)); + xy.y = P->radius_g_1 * atan(Vz / tmp); + } + else + { + xy.x = P->radius_g_1 * atan(Vy / tmp); + xy.y = P->radius_g_1 * atan(Vz / hypot(Vy, tmp)); + } return (xy); } FORWARD(e_forward); /* ellipsoid */ @@ -74,8 +85,16 @@ FORWARD(e_forward); /* ellipsoid */ F_ERROR; /* Calculation based on view angles from satellite. */ tmp = P->radius_g - Vx; - xy.x = P->radius_g_1 * atan (Vy / tmp); - xy.y = P->radius_g_1 * atan (Vz / hypot (Vy, tmp)); + if(P->flip_axis) + { + xy.x = P->radius_g_1 * atan (Vy / hypot (Vz, tmp)); + xy.y = P->radius_g_1 * atan (Vz / tmp); + } + else + { + xy.x = P->radius_g_1 * atan (Vy / tmp); + xy.y = P->radius_g_1 * atan (Vz / hypot (Vy, tmp)); + } return (xy); } INVERSE(s_inverse); /* spheroid */ @@ -83,8 +102,16 @@ INVERSE(s_inverse); /* spheroid */ /* Setting three components of vector from satellite to position.*/ Vx = -1.0; - Vy = tan (xy.x / (P->radius_g - 1.0)); - Vz = tan (xy.y / (P->radius_g - 1.0)) * sqrt (1.0 + Vy * Vy); + if(P->flip_axis) + { + Vz = tan (xy.y / (P->radius_g - 1.0)); + Vy = tan (xy.x / (P->radius_g - 1.0)) * sqrt (1.0 + Vz * Vz); + } + else + { + Vy = tan (xy.x / (P->radius_g - 1.0)); + Vz = tan (xy.y / (P->radius_g - 1.0)) * sqrt (1.0 + Vy * Vy); + } /* Calculation of terms in cubic equation and determinant.*/ a = Vy * Vy + Vz * Vz + Vx * Vx; b = 2 * P->radius_g * Vx; @@ -104,8 +131,16 @@ INVERSE(e_inverse); /* ellipsoid */ /* Setting three components of vector from satellite to position.*/ Vx = -1.0; - Vy = tan (xy.x / P->radius_g_1); - Vz = tan (xy.y / P->radius_g_1) * hypot(1.0, Vy); + if(P->flip_axis) + { + Vz = tan (xy.y / P->radius_g_1); + Vy = tan (xy.x / P->radius_g_1) * hypot(1.0, Vz); + } + else + { + Vy = tan (xy.x / P->radius_g_1); + Vz = tan (xy.y / P->radius_g_1) * hypot(1.0, Vy); + } /* Calculation of terms in cubic equation and determinant.*/ a = Vz / P->radius_p; a = Vy * Vy + a * a + Vx * Vx; @@ -126,7 +161,22 @@ FREEUP; if (P) free(P); } ENTRY0(geos) if ((P->h = pj_param(P->ctx, P->params, "dh").f) <= 0.) E_ERROR(-30); if (P->phi0) E_ERROR(-46); - P->radius_g = 1. + (P->radius_g_1 = P->h / P->a); + P->sweep_axis = pj_param(P->ctx, P->params, "ssweep").s; + if (P->sweep_axis == NULL) + P->flip_axis = 0; + else + { + if (P->sweep_axis[1] != '\0' || + (P->sweep_axis[0] != 'x' && + P->sweep_axis[0] != 'y')) + E_ERROR(-49); + if (P->sweep_axis[0] == 'y') + P->flip_axis = 1; + else + P->flip_axis = 0; + } + P->radius_g_1 = P->h / P->a; + P->radius_g = 1. + P->radius_g_1; P->C = P->radius_g * P->radius_g - 1.0; if (P->es) { P->radius_p = sqrt (P->one_es); diff --git a/src/pj_strerrno.c b/src/pj_strerrno.c index 186a50e8..9d237039 100644 --- a/src/pj_strerrno.c +++ b/src/pj_strerrno.c @@ -54,6 +54,7 @@ pj_err_list[] = { "unknown prime meridian conversion id", /* -46 */ "illegal axis orientation combination", /* -47 */ "point not within available datum shift grids", /* -48 */ + "invalid sweep axis, choose x or y", /* -49 */ }; char * pj_strerrno(int err) -- cgit v1.2.3