From f1ea5041d5049498e3cd88f9a215ebb0f4b86a34 Mon Sep 17 00:00:00 2001 From: Javier Goizueta Date: Wed, 11 Apr 2018 17:56:08 +0200 Subject: Enhance the precision of Spherical Mercator projection near the equator This uses a linear approximation of tan(x+pi/4) for better precision at small latitudes. As a result points of latitude 0 maintain a 0 Y coordinate and 0,0 is transformed to 0,0 --- src/PJ_merc.c | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/PJ_merc.c b/src/PJ_merc.c index a3e5e846..3801b828 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -8,6 +8,13 @@ PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; #define EPS10 1.e-10 +static double _tan_near_fort_pi(double x) { + if (fabs(x) <= __DBL_EPSILON__) { + return 2*x + 1.0; + } + return tan(M_FORTPI + x); +} + static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ XY xy = {0.0,0.0}; if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) { @@ -27,7 +34,7 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ return xy; } xy.x = P->k0 * lp.lam; - xy.y = P->k0 * log(tan(M_FORTPI + .5 * lp.phi)); + xy.y = P->k0 * log(_tan_near_fort_pi(.5 * lp.phi)); return xy; } -- cgit v1.2.3 From bd2e1363b021e1a1c6a4ebbc0aa55e8dc98e2587 Mon Sep 17 00:00:00 2001 From: Javier Goizueta Date: Thu, 12 Apr 2018 08:50:59 +0200 Subject: Fix: use proper double epsilon declaration --- src/PJ_merc.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/PJ_merc.c b/src/PJ_merc.c index 3801b828..2c89fbc6 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -2,6 +2,7 @@ #include "proj_internal.h" #include "proj.h" #include "projects.h" +#include PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts="; PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; @@ -9,7 +10,7 @@ PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; #define EPS10 1.e-10 static double _tan_near_fort_pi(double x) { - if (fabs(x) <= __DBL_EPSILON__) { + if (fabs(x) <= DBL_EPSILON) { return 2*x + 1.0; } return tan(M_FORTPI + x); -- cgit v1.2.3 From f2b3604a94349697828baa7de63a30a97ac15b2b Mon Sep 17 00:00:00 2001 From: Javier Goizueta Date: Thu, 12 Apr 2018 09:47:43 +0200 Subject: Use log1p in forward spherical mercator --- src/PJ_merc.c | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/PJ_merc.c b/src/PJ_merc.c index 2c89fbc6..0bf98625 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -9,11 +9,31 @@ PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; #define EPS10 1.e-10 -static double _tan_near_fort_pi(double x) { +#if !defined(HAVE_C99_MATH) +#define HAVE_C99_MATH 0 +#endif + +#if HAVE_C99_MATH +#define log1px log1p +#else +static double log1px(double x) { + volatile double + y = 1 + x, + z = y - 1; + /* Here's the explanation for this magic: y = 1 + z, exactly, and z + * approx x, thus log(y)/z (which is nearly constant near z = 0) returns + * a good approximation to the true log(1 + x)/x. The multiplication x * + * (log(y)/z) introduces little additional error. */ + return z == 0 ? x : x * log(y) / z; +} +#endif + +static double logtanpfpim1(double x) { /* log(tan(x/2 + M_FORTPI)) */ if (fabs(x) <= DBL_EPSILON) { - return 2*x + 1.0; + /* tan(M_FORTPI + .5 * x) can be approximated by 1.0 + x */ + return log1px(x); } - return tan(M_FORTPI + x); + return log(tan(M_FORTPI + .5 * x)); } static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ @@ -35,7 +55,7 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ return xy; } xy.x = P->k0 * lp.lam; - xy.y = P->k0 * log(_tan_near_fort_pi(.5 * lp.phi)); + xy.y = P->k0 * logtanpfpim1(lp.phi); return xy; } -- cgit v1.2.3