From 67c20c7a6eb9fa45c8ee4a8b759821eef82d39ac Mon Sep 17 00:00:00 2001 From: Javier Goizueta Date: Thu, 12 Apr 2018 18:12:26 +0200 Subject: [PATCH 1/2] Enhance precision of inverse Spherical Mercator projection near the equator --- src/PJ_merc.c | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/src/PJ_merc.c b/src/PJ_merc.c index 0bf98625d4..9697de6720 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -15,6 +15,7 @@ PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; #if HAVE_C99_MATH #define log1px log1p +#define expm1x expm1 #else static double log1px(double x) { volatile double @@ -26,6 +27,13 @@ static double log1px(double x) { * (log(y)/z) introduces little additional error. */ return z == 0 ? x : x * log(y) / z; } +static double expm1x(double x) { + /* very simplistic exmp1 implementation */ + if (fabs(x) <= DBL_EPSILON) { + return x; + } + return exp(x) - 1.0; +} #endif static double logtanpfpim1(double x) { /* log(tan(x/2 + M_FORTPI)) */ @@ -36,6 +44,14 @@ static double logtanpfpim1(double x) { /* log(tan(x/2 + M_FORTPI)) */ return log(tan(M_FORTPI + .5 * x)); } +static double hpip2atanexp(double x) { /* M_HALFPI - 2. * atan(exp(x)) */ + if (fabs(x) <= DBL_EPSILON) { + /* atan(exp(x)) can be approximated by (exp(x)-1)/2 + pi/4 */ + return -expm1x(x)); + } + return M_HALFPI - 2. * atan(exp(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) { @@ -73,7 +89,7 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ LP lp = {0.0,0.0}; - lp.phi = M_HALFPI - 2. * atan(exp(-xy.y / P->k0)); + lp.phi = hpip2atanexp(xy.y / P->k0); lp.lam = xy.x / P->k0; return lp; } From 38277fc2abb68479d2ec499d72d1b6dd0c3a9dbc Mon Sep 17 00:00:00 2001 From: Javier Goizueta Date: Mon, 16 Apr 2018 15:27:12 +0200 Subject: [PATCH 2/2] Fix syntax --- src/PJ_merc.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/PJ_merc.c b/src/PJ_merc.c index 9697de6720..c4923c92eb 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -47,7 +47,7 @@ static double logtanpfpim1(double x) { /* log(tan(x/2 + M_FORTPI)) */ static double hpip2atanexp(double x) { /* M_HALFPI - 2. * atan(exp(x)) */ if (fabs(x) <= DBL_EPSILON) { /* atan(exp(x)) can be approximated by (exp(x)-1)/2 + pi/4 */ - return -expm1x(x)); + return -expm1x(x); } return M_HALFPI - 2. * atan(exp(x)); }