diff options
| -rw-r--r-- | src/PJ_aea.c | 1 | ||||
| -rw-r--r-- | src/PJ_natearth2.c | 146 |
2 files changed, 109 insertions, 38 deletions
diff --git a/src/PJ_aea.c b/src/PJ_aea.c index c5480d1e..aea92854 100644 --- a/src/PJ_aea.c +++ b/src/PJ_aea.c @@ -364,7 +364,6 @@ int pj_latlong_selftest (void) {return 10000;} int pj_lonlat_selftest (void) {return 10000;} int pj_longlat_selftest (void) {return 10000;} -int pj_natearth2_selftest (void) {return 10000;} int pj_nell_selftest (void) {return 10000;} int pj_nell_h_selftest (void) {return 10000;} int pj_nicol_selftest (void) {return 10000;} diff --git a/src/PJ_natearth2.c b/src/PJ_natearth2.c index c3776a8b..5a57fb14 100644 --- a/src/PJ_natearth2.c +++ b/src/PJ_natearth2.c @@ -1,13 +1,13 @@ /* -The Natural Earth II projection was designed by Tom Patterson, US National +The Natural Earth II projection was designed by Tom Patterson, US National Park Service, in 2012, using Flex Projector. The polynomial equation was -developed by Bojan Savric and Bernhard Jenny, College of Earth, Ocean, +developed by Bojan Savric and Bernhard Jenny, College of Earth, Ocean, and Atmospheric Sciences, Oregon State University. Port to PROJ.4 by Bojan Savric, 4 April 2016 */ - #define PJ_LIB__ -#include <projects.h> +#include <projects.h> + PROJ_HEAD(natearth2, "Natural Earth 2") "\n\tPCyl., Sph."; #define A0 0.84719 @@ -27,49 +27,121 @@ PROJ_HEAD(natearth2, "Natural Earth 2") "\n\tPCyl., Sph."; #define EPS 1e-11 #define MAX_Y (0.84719 * 0.535117535153096 * PI) -FORWARD(s_forward); /* spheroid */ - double phi2, phi4, phi6; - phi2 = lp.phi * lp.phi; - phi4 = phi2 * phi2; - phi6 = phi2 * phi4; +static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ + XY xy = {0.0,0.0}; + double phi2, phi4, phi6; + (void) P; - xy.x = lp.lam * (A0 + A1 * phi2 + phi6 * phi6 * (A2 + A3 * phi2 + A4 * phi4 + A5 * phi6)); - xy.y = lp.phi * (B0 + phi4 * phi4 * (B1 + B2 * phi2 + B3 * phi4)); - return (xy); + phi2 = lp.phi * lp.phi; + phi4 = phi2 * phi2; + phi6 = phi2 * phi4; + + xy.x = lp.lam * (A0 + A1 * phi2 + phi6 * phi6 * (A2 + A3 * phi2 + A4 * phi4 + A5 * phi6)); + xy.y = lp.phi * (B0 + phi4 * phi4 * (B1 + B2 * phi2 + B3 * phi4)); + return xy; } -INVERSE(s_inverse); /* spheroid */ - double yc, tol, y2, y4, y6, f, fder; + + +static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ + LP lp = {0.0,0.0}; + double yc, tol, y2, y4, y6, f, fder; + (void) P; /* make sure y is inside valid range */ - if (xy.y > MAX_Y) { - xy.y = MAX_Y; - } else if (xy.y < -MAX_Y) { - xy.y = -MAX_Y; - } + if (xy.y > MAX_Y) { + xy.y = MAX_Y; + } else if (xy.y < -MAX_Y) { + xy.y = -MAX_Y; + } /* latitude */ - yc = xy.y; + yc = xy.y; for (;;) { /* Newton-Raphson */ - y2 = yc * yc; - y4 = y2 * y2; - f = (yc * (B0 + y4 * y4 * (B1 + B2 * y2 + B3 * y4))) - xy.y; - fder = C0 + y4 * y4 * (C1 + C2 * y2 + C3 * y4); - yc -= tol = f / fder; - if (fabs(tol) < EPS) { - break; - } - } - lp.phi = yc; + y2 = yc * yc; + y4 = y2 * y2; + f = (yc * (B0 + y4 * y4 * (B1 + B2 * y2 + B3 * y4))) - xy.y; + fder = C0 + y4 * y4 * (C1 + C2 * y2 + C3 * y4); + yc -= tol = f / fder; + if (fabs(tol) < EPS) { + break; + } + } + lp.phi = yc; /* longitude */ - y2 = yc * yc; - y4 = y2 * y2; - y6 = y2 * y4; + y2 = yc * yc; + y4 = y2 * y2; + y6 = y2 * y4; + + lp.lam = xy.x / (A0 + A1 * y2 + y6 * y6 * (A2 + A3 * y2 + A4 * y4 + A5 * y6)); + + return lp; +} + - lp.lam = xy.x / (A0 + A1 * y2 + y6 * y6 * (A2 + A3 * y2 + A4 * y4 + A5 * y6)); +static void *freeup_new (PJ *P) { /* Destructor */ + if (0==P) + return 0; - return (lp); + return pj_dealloc(P); } -FREEUP; if (P) pj_dalloc(P); } -ENTRY0(natearth2) P->es = 0; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P) + + +static void freeup (PJ *P) { + freeup_new (P); + return; +} + + +PJ *PROJECTION(natearth2) { + P->es = 0; + P->inv = s_inverse; + P->fwd = s_forward; + + return P; +} + +#ifdef PJ_OMIT_SELFTEST +int pj_natearth2_selftest (void) {return 0;} +#else + +int pj_natearth2_selftest (void) { + double tolerance_lp = 1e-10; + double tolerance_xy = 1e-7; + + char s_args[] = {"+proj=natearth2 +a=6400000 +lat_1=0.5 +lat_2=2"}; + + LP fwd_in[] = { + { 2, 1}, + { 2,-1}, + {-2, 1}, + {-2,-1} + }; + + XY s_fwd_expect[] = { + { 189255.172934730799, 113022.495810907014}, + { 189255.172934730799, -113022.495810907014}, + {-189255.172934730799, 113022.495810907014}, + {-189255.172934730799, -113022.495810907014}, + }; + + XY inv_in[] = { + { 200, 100}, + { 200,-100}, + {-200, 100}, + {-200,-100} + }; + + LP s_inv_expect[] = { + { 0.00211344929691056112, 0.000884779612080993237}, + { 0.00211344929691056112, -0.000884779612080993237}, + {-0.00211344929691056112, 0.000884779612080993237}, + {-0.00211344929691056112, -0.000884779612080993237}, + }; + + return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect); +} + + +#endif |