diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Makefile.am | 2 | ||||
| -rw-r--r-- | src/Makefile.in | 4 | ||||
| -rw-r--r-- | src/PJ_natearth.c | 77 | ||||
| -rw-r--r-- | src/makefile.vc | 2 | ||||
| -rw-r--r-- | src/pj_list.h | 1 |
5 files changed, 82 insertions, 4 deletions
diff --git a/src/Makefile.am b/src/Makefile.am index 97e9e833..621b0b4a 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -57,7 +57,7 @@ libproj_la_SOURCES = \ pj_open_lib.c pj_param.c pj_phi2.c pj_pr_list.c \ pj_qsfn.c pj_strerrno.c pj_tsfn.c pj_units.c pj_ctx.c pj_log.c \ pj_zpoly1.c rtodms.c vector1.c pj_release.c pj_gauss.c \ - PJ_healpix.c \ + PJ_healpix.c PJ_natearth.c \ \ nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \ pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \ diff --git a/src/Makefile.in b/src/Makefile.in index 412bd0bd..61ef536d 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -102,7 +102,7 @@ am_libproj_la_OBJECTS = PJ_aeqd.lo PJ_gnom.lo PJ_laea.lo \ pj_open_lib.lo pj_param.lo pj_phi2.lo pj_pr_list.lo pj_qsfn.lo \ pj_strerrno.lo pj_tsfn.lo pj_units.lo pj_ctx.lo pj_log.lo \ pj_zpoly1.lo rtodms.lo vector1.lo pj_release.lo pj_gauss.lo \ - PJ_healpix.lo nad_cvt.lo nad_init.lo nad_intr.lo emess.lo \ + PJ_healpix.lo PJ_natearth.lo nad_cvt.lo nad_init.lo nad_intr.lo emess.lo \ pj_apply_gridshift.lo pj_datums.lo pj_datum_set.lo \ pj_transform.lo geocent.lo pj_utils.lo pj_gridinfo.lo \ pj_gridlist.lo jniproj.lo pj_mutex.lo pj_initcache.lo \ @@ -316,7 +316,7 @@ libproj_la_SOURCES = \ pj_open_lib.c pj_param.c pj_phi2.c pj_pr_list.c \ pj_qsfn.c pj_strerrno.c pj_tsfn.c pj_units.c pj_ctx.c pj_log.c \ pj_zpoly1.c rtodms.c vector1.c pj_release.c pj_gauss.c \ - PJ_healpix.c \ + PJ_healpix.c PJ_natearth.c \ \ nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \ pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \ diff --git a/src/PJ_natearth.c b/src/PJ_natearth.c new file mode 100644 index 00000000..f83a3ffd --- /dev/null +++ b/src/PJ_natearth.c @@ -0,0 +1,77 @@ +/* +The Natural Earth projection was designed by Tom Patterson, US National Park +Service, in 2007, using Flex Projector. The shape of the original projection +was defined at every 5 degrees and piece-wise cubic spline interpolation was +used to compute the complete graticule. +The code here uses polynomial functions instead of cubic splines and +is therefore much simpler to program. The polynomial approximation was +developed by Bojan Savric, in collaboration with Tom Patterson and Bernhard +Jenny, Institute of Cartography, ETH Zurich. It slightly deviates from +Patterson's original projection by adding additional curvature to meridians +where they meet the horizontal pole line. This improvement is by intention +and designed in collaboration with Tom Patterson. +Port to PROJ.4 by Bernhard Jenny, 6 June 2011 +*/ + +#define PJ_LIB__ +#include <projects.h> +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") |