Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Enhance precision of inverse Spherical Mercator projection near the equator #931

Closed
wants to merge 2 commits into from
Closed

Enhance precision of inverse Spherical Mercator projection near the equator #931

Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
67c20c7
Enhance precision of inverse Spherical Mercator projection near the e…
jgoizueta Apr 12, 2018
38277fc
Fix syntax
jgoizueta Apr 16, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
18 changes: 17 additions & 1 deletion src/PJ_merc.c
View file
Open in desktop
Expand Up @@ -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
Expand All @@ -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)) */
Expand All @@ -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));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The extra ) at the end here causes the build failure.

}
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) {
Expand Down Expand Up @@ -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;
}
Expand Down