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.

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

Enhance the precision of Spherical Mercator projection near the equator #928

Merged
merged 3 commits into from Apr 12, 2018

Conversation

Projects
None yet
3 participants
@jgoizueta
Contributor

jgoizueta commented Apr 11, 2018

This uses a linear approximation of tan(x+pi/4) for better precision at small latitudes.

Before this patch:

cs2cs +init=epsg:4326 +to +init=epsg:3857  -v -f "%.9e" <<EOF
0 0
0 1E-14
0 1E-13
0 1E-12
0 1E-11
0 1E-10
0 1E-9
0 1
EOF
0.000000000e+00	-7.081154552e-10 0.000000000e+00
0.000000000e+00	1.416230910e-09 0.000000000e+00
0.000000000e+00	9.913616372e-09 0.000000000e+00
0.000000000e+00	1.104660110e-07 0.000000000e+00
0.000000000e+00	1.113157496e-06 0.000000000e+00
0.000000000e+00	1.113015872e-05 0.000000000e+00
0.000000000e+00	1.113199982e-04 0.000000000e+00
0.000000000e+00	1.113251429e+05 0.000000000e+00

After:

./src/cs2cs +init=epsg:4326 +to +init=epsg:3857  -v -f "%.9e" <<EOF
0 0
0 1E-14
0 1E-13
0 1E-12
0 1E-11
0 1E-10
0 1E-9
0 1
EOF
0.000000000e+00	0.000000000e+00 0.000000000e+00
0.000000000e+00	1.416230910e-09 0.000000000e+00
0.000000000e+00	1.132984728e-08 0.000000000e+00
0.000000000e+00	1.118822419e-07 0.000000000e+00
0.000000000e+00	1.113157496e-06 0.000000000e+00
0.000000000e+00	1.113157496e-05 0.000000000e+00
0.000000000e+00	1.113199982e-04 0.000000000e+00
0.000000000e+00	1.113251429e+05 0.000000000e+00

All values are now more accurate (I checked using arbitrary precision arithmetic), and we have the nice property of preserving the (0,0) coordinates of the origin point.

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
@cffk

This comment has been minimized.

Contributor

cffk commented Apr 11, 2018

Perhaps now is a good time to revisit the Mercator formulas. Approximating tan(pi/4 + phi/2) by 1 + phi, still leaves the problem that you have to take the log of a quantity close to unity.

Much better, from a numerical point of view, than log(tan(pi/4 + phi/2)) is the mathematically equivalent expression asinh(tan(x)). This requires that you have a good implementation of asinh(x) which you do have if HAVE_C99_MATH. Otherwise, you can use

y = fabs(x);
y = log1p(y * (1 + y/(hypot(1.0, y) + 1)));
return x < 0 ? -y : (x > 0 ? y : x);  // asinh(-0.0) = -0.0

and log1p can be implemented as in geodesic.c.

@kbevers

This comment has been minimized.

Member

kbevers commented Apr 11, 2018

Perhaps now is a good time to revisit the Mercator formulas.

Now is as good a time as any. What about the inverse formulas? There's an atan(exp(x)) thing going on. Presumably that can also be done better numerically.

@jgoizueta

This comment has been minimized.

Contributor

jgoizueta commented Apr 12, 2018

I won't get into larger changes like using atanh at the moment (I think it would required further analysis) but I think using logp1 is worth and straightforward to include here.

Looking at the inverse formula would also be interesting, but similarly I won't get
into analyzing tan(exp(x)) now, but will look instead at the cancellation near Y=0 for consistency with the other changes.

Thanks @cffk & @kbevers!

@kbevers

This comment has been minimized.

Member

kbevers commented Apr 12, 2018

Thanks for giving this a go @jgoizueta. We're at the limit of my numerical math capabilities here so I'll defer to @cffk's judgement on this. Are you happy with this approach, Charles?

@kbevers

This comment has been minimized.

Member

kbevers commented Apr 12, 2018

Looking a bit more in to this it turns out that we have a log1px defined already in the extended transverse mercator:
https://github.com/OSGeo/proj.4/blob/master/src/proj_etmerc.c#L69-L78

There's even a asinh function as well:
https://github.com/OSGeo/proj.4/blob/master/src/proj_etmerc.c#L84-L88

We are getting close to the point were a dedicated non-C99-math function collection is needed.

@cffk

This comment has been minimized.

Contributor

cffk commented Apr 12, 2018

I agree... Perhaps we should let this modification go through but queue of a bigger survey of these hygene + accuracy issues.

Note well that several conformal projections use similar formulas (converting between geographic latitude to conformal latitude so that a spherical conformal projection can be used). Off the top of my head, I suspect that

  • Mercator
  • transverse Mercator
  • polar stereographic
  • Lambert conformal conic

are all candidates for a uniform treatment.

Resources are my paper on transverse Mercator, eqs (3), (4). Here the treatment works for arbitrary eccentricity. Engsager and Poder (and etmerc) use a trigonometric series for the transformation.

In addition to a C99 math library collection, perhaps it would be good to consider a centralized location for Clenshaw and Horner sums, and the basic ellipsoidal Mercator formulas (for use by etmerc etc.)

The various equal-area projections should also use the same machinery (if they don't already) to do the transformation between geographic latitude and authalic latitude. See also my table of trigonometric series for transformations between all the various latitudes (E+K's series appear here φ − χ and χ − φ).

@kbevers kbevers merged commit 8e1cdad into OSGeo:master Apr 12, 2018

3 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
coverage/coveralls Coverage increased (+0.2%) to 77.504%
Details
@kbevers

This comment has been minimized.

Member

kbevers commented Apr 12, 2018

I agree... Perhaps we should let this modification go through but queue of a bigger survey of these hygene + accuracy issues.

I have merged this (thanks @jgoizueta) and opened for #930 for further discussion on the generic solution to the math functions problem.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment