From 5229741fe5cdf2e5fa4772b9141298417331fa98 Mon Sep 17 00:00:00 2001 From: Frank Warmerdam Date: Sun, 4 Mar 2012 00:15:57 +0000 Subject: added the Natural Earth projection git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@2182 4e78687f-474d-0410-85f9-8d5e500ac6b2 --- NEWS | 2 ++ nad/testvarious | 31 ++++++++++++++++++++++ nad/tv_out.dist | 27 +++++++++++++++++++ src/Makefile.am | 2 +- src/Makefile.in | 4 +-- src/PJ_natearth.c | 77 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/makefile.vc | 2 +- src/pj_list.h | 1 + 8 files changed, 142 insertions(+), 4 deletions(-) create mode 100644 src/PJ_natearth.c diff --git a/NEWS b/NEWS index 4ad2f9f6..99ecdc5f 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,8 @@ 4.8.0 Release Notes ------------------- + o Added the Natural Earth projection. + o Added HEALPIX, rHEALPIX and Icosahedral Snyder Equal Area projections. o nad2bin now produces "CTable2" format grid shift files by default which diff --git a/nad/testvarious b/nad/testvarious index aadb8796..3886de56 100755 --- a/nad/testvarious +++ b/nad/testvarious @@ -371,6 +371,37 @@ $EXE +proj=latlong +ellps=sphere \ 3820138.08 2888664.15 EOF # +echo "##############################################################" >> ${OUT} +echo "Test the natural earth projection" >> ${OUT} +$EXE +proj=latlong +a=1 +lon_0=0 \ + +to +proj=natearth +a=6371008.7714 +b=6371008.7714 -f '%.'7'f' \ + -E >>${OUT} < +PROJ_HEAD(natearth, "Natural Earth") "\n\tPCyl., Sph."; +#define A0 0.8707 +#define A1 -0.131979 +#define A2 -0.013791 +#define A3 0.003971 +#define A4 -0.001529 +#define B0 1.007226 +#define B1 0.015085 +#define B2 -0.044475 +#define B3 0.028874 +#define B4 -0.005916 +#define C0 B0 +#define C1 (3 * B1) +#define C2 (7 * B2) +#define C3 (9 * B3) +#define C4 (11 * B4) +#define EPS 1e-11 +#define MAX_Y (0.8707 * 0.52 * PI) + +FORWARD(s_forward); /* spheroid */ + double phi2, phi4; + + phi2 = lp.phi * lp.phi; + phi4 = phi2 * phi2; + xy.x = lp.lam * (A0 + phi2 * (A1 + phi2 * (A2 + phi4 * phi2 * (A3 + phi2 * A4)))); + xy.y = lp.phi * (B0 + phi2 * (B1 + phi4 * (B2 + B3 * phi2 + B4 * phi4))); + return (xy); +} +INVERSE(s_inverse); /* spheroid */ + double yc, tol, y2, y4, f, fder; + + /* 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; + } + + /* latitude */ + yc = xy.y; + for (;;) { /* Newton-Raphson */ + y2 = yc * yc; + y4 = y2 * y2; + f = (yc * (B0 + y2 * (B1 + y4 * (B2 + B3 * y2 + B4 * y4)))) - xy.y; + fder = C0 + y2 * (C1 + y4 * (C2 + C3 * y2 + C4 * y4)); + yc -= tol = f / fder; + if (fabs(tol) < EPS) { + break; + } + } + lp.phi = yc; + + /* longitude */ + y2 = yc * yc; + lp.lam = xy.x / (A0 + y2 * (A1 + y2 * (A2 + y2 * y2 * y2 * (A3 + y2 * A4)))); + + return (lp); +} +FREEUP; if (P) pj_dalloc(P); } +ENTRY0(natearth) P->es = 0; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P) diff --git a/src/makefile.vc b/src/makefile.vc index 220468ab..58c348f6 100644 --- a/src/makefile.vc +++ b/src/makefile.vc @@ -27,7 +27,7 @@ misc = \ PJ_lask.obj PJ_nocol.obj PJ_ob_tran.obj PJ_oea.obj \ PJ_tpeqd.obj PJ_vandg.obj PJ_vandg2.obj PJ_vandg4.obj \ PJ_wag7.obj pj_latlong.obj PJ_krovak.obj pj_geocent.obj \ - PJ_healpix.obj + PJ_healpix.obj PJ_natearth.obj pseudo = \ PJ_boggs.obj PJ_collg.obj PJ_crast.obj PJ_denoy.obj \ diff --git a/src/pj_list.h b/src/pj_list.h index ac88f9d8..8230f92f 100644 --- a/src/pj_list.h +++ b/src/pj_list.h @@ -82,6 +82,7 @@ PROJ_HEAD(moll, "Mollweide") PROJ_HEAD(murd1, "Murdoch I") PROJ_HEAD(murd2, "Murdoch II") PROJ_HEAD(murd3, "Murdoch III") +PROJ_HEAD(natearth, "Natural Earth") PROJ_HEAD(nell, "Nell") PROJ_HEAD(nell_h, "Nell-Hammer") PROJ_HEAD(nicol, "Nicolosi Globular") -- cgit v1.2.3