Implement 'exact' authalic latitude->geographic latitude and use it in +proj=aea/cea/laea/eqearth/healpix/rhealpix

[rouault]

Fixes #4440

$ echo 1020000 1562000 | bin/cs2cs epsg:3035 epsg:4937 -d 10 | bin/cs2cs epsg:4937 epsg:3035 -d 4
1020000.0000 1562000.0000 0.0000

[rouault]

@kbevers I'm not sure if this is appropriate for backport. This is somewhere in between a new feature and a bugfix. This could cause downstream software to potentially have failures in their test suite if they test coordinates with sufficient precision, so perhaps not appropriate for a bugfix release.

[jjimenezshaw]

@cffk was suggesting something in 2023 about the authalic latitude in https://lists.osgeo.org/pipermail/proj/2023-March/010956.html

[kbevers]

I'm not sure if this is appropriate for backport.

Let's call it a new feature, this is a changed approach that delivers better accuracy in some cases.

@cffk since you've written the book (well, sort of) on this topic it would be foolish not request your input here :-)

[jjimenezshaw]

I was going to add the reference in docs/source/references.bib in a PR, but I think that it does not appear if it is not referenced anywhere.
So here it is

@Article{Karney2023,
  author                   = {Charles F. F. Karney},
  title                    = {On auxiliary latitudes},
  journal                  = {Survey Review},
  year                     = {2024},
  volume                   = {56},
  number                   = {395},
  pages                    = {165-180},
  eprint                   = {2212.05818},
  doi                      = {10.1080/00396265.2023.2217604}
}

[cffk]

I'm not sure if this is appropriate for backport.

Let's call it a new feature, this is a changed approach that delivers better accuracy in some cases.

@cffk since you've written the book (well, sort of) on this topic it would be foolish not request your input here :-)

I'll take a look and give feedback tomorrow (03-21).

[cffk]

A. Addressing this PR is the narrowest way...

1. Presumably the last line of pj_authalic_lat_inverse_exact should be

return phi;

2. When implementing Newton's method (the same function), it's helpful
   to include a comment as to how the problem is being formulated. How about

// Given beta, solve
//   f(phi) = qp*sin(beta)/(1-e^2) - q(phi)/(1-e^2)
// for phi, where
//   q(phi)/(1-e^2) = sin(phi)/(1 - e^2*sin(phi)^2) + atanh(e*sin(phi))/e
// and
//   df(phi)/dphi = - dq(phi)/dphi / (1-e^2)
//                = - 2 * (1-e^2) * cos(phi) / (1 - e^2*sinphi^2)^2

3. Given that pj_authalic_lat_inverse_exact replaces the logarithm
   expression with atanh, then so too should pj_authalic_lat_q_coeff.
   This way the function is exactly odd and maintains high relative
   accuracy near the equator.

B. However, widening the scope of this comment ...

The method for computing the auxiliary latitude pj_authalic_lat is
badly ill-conditioned near the poles because of the properties of
arcsine. (And also near the equator, unless atanh is used.) The
spatial resolution of asin(q/qp) at the equator is about 90 mm.

Hence, pj_authalic_lat_inverse_exact is also inaccurate near the
poles. A better fix would therefore be to upgrade
pj_authalic_lat_inverse_approx so that it returns an accurate result.
The current implementation uses the series of Adams (1921) which is 3rd
order in e^2. I recommend using the 6th order series in n with the
coefficients given by (A20) of my paper. (Here is the link to the
preprint: https://arxiv.org/abs/2212.05818.) This is accurate to
roundoff for |f| < 1/150. The series can be evaluated with the
existing machinery in mlfn.cpp which evaluates similar series (for the
conversions between geographic and rectifying latitude) with Horner and
Clenshaw. The result would then be more accurate than
pj_authalic_lat_inverse_exact

C. Finally ...

A longer range goal would be to replace the ill-conditioned formula for
pj_authatic_lat by a similar series, Eq. (A19) in my paper.

Indeed latitude.cpp might be fixed also to provide the series results
for pj_conformal_lat_inverse (currently using a simple-minded and slow
iterative method). (As an aside, in pj_conformal_lat, Snyder's (3-1)
should be replaced with my better-conditioned Eqs. (8) and (10).)

And it might be extended to support the conversions in mlfn.cpp and in
tmerc.cpp. (The latter would require the implementation of complex
Clenshaw from that file.)

In summary, I recommend A. Now also would be a good time to implement
B. C should be considered later with the replacement of
pj_authalic_lat being at the top of the queue.

[rouault]

@cffk Thanks for the review. I've added a commit using your comments A. B and C are a bit beyond my understanding / available time unfortunately. If we were going to do this, I believe the most straightforward way would probably to just borrow your implementation AuxLatitude.hpp, AuxLatitude.cpp from https://doi.org/10.5281/zenodo.7382666 .

[cffk]

Let's go with this. To address Even's remark about absorbing my AuxLatitude implementation: I fear this may be overkill. Perhaps @jjimenezshaw can be persuaded to start a cleanup of latitude.cpp. Or maybe I'll take a crack at doing the two authalic conversions.

[rouault]

hum it looks like using atanh() in the forward path degrades the accuracy on i386 at the pole. Maybe we should revert that part ?

[cffk]

hum it looks like using atanh() in the forward path degrades the accuracy on i386 at the pole. Maybe we should revert that part ?

How are you measuring this? The log expression entails significant loss of precision (you're evalutating log(1 + small) with the result being small -- but you've lost significant bits with 1 + small); so don't compare the atanh expression to this unless you evaluation the log expression with long double.

My recommendation: keep the atanh formulation. The accuracy at the pole is already shot with the asin.

[rouault]

How are you measuring this?

well, https://github.com/OSGeo/PROJ/actions/runs/13995771272/job/39190443832?pr=4441 now fails, whereas it didn't before my additional commit

I can comment out the problematic test points / relax their precision requirements for now