[PROJ] Geodetic to Authalic latitude conversions

[Jérôme St-Louis]

Hi Even, all,

As suggested in https://github.com/OSGeo/PROJ/pull/4211 , I've just 
written an implementation of the authalic / geodetic latitude conversion 
with the goal to automatically perform this conversion when using ISEA 
with an non-spherical ellipsoid, and was wondering whether we would want 
to improve the PROJ version (/pj_authlat()/ for the geodetic ==> 
authalic), which is used for ISEA, but also by:

- Equal Area Cylindrical (cea.cpp)
- EqualEarth (eqearth.cpp)
- HEALPix (healpix.cpp) -- which might also imply rHEALPix
- Lambert Azimuthal Equal Area (laea.cpp)

The current version of /pj_authlat()/ uses a series expansion based on 
the eccentricity squared (I believe this is based on a Lagrange 
reversion <https://en.wikipedia.org/wiki/Lagrange_inversion_theorem> as 
mentioned here 
<https://en.wikipedia.org/wiki/Latitude#Inverse_formulae_and_series>), 
and Charles Karney in this paper <https://arxiv.org/pdf/2212.05818> 
establishes that:

> /series expansions using the third flattening as the small parameter 
> uniformly result in smaller truncation errors compared to using the 
> eccentricity./
The current version also only uses order 3, whereas GeographicLib would 
normally use order 6 for double precision, with order 4 being the 
minimum configuration.
The current order 3, which requires 2 /sin()/ calls, gives precision up 
to roughly 8 decimals.
Order 6 would mean twice the number of /sin()/ calls, and the 8 decimals 
should be enough to point to Waldo on a page <https://xkcd.com/2170/>, 
however C.K. argues that:

> /Modern geodetic libraries should strive to achieve full 
> double-precision accuracy /(page 4)/
> /

and I would much agree with that.

Based on the direction we want to go, I would submit a PR including 
adjustments to avoid the separate heap allocation for the doubles (as 
discussed here 
<https://github.com/OSGeo/PROJ/pull/4211#pullrequestreview-2233903246>), 
and adjustments for tests etc.

PROJ currently does not use the full GeographicLib which already 
includes code for conversion between different auxiliary latitudes, does it?

     ( 
https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp 
)

See my proposed new version below, which I wrote based on equation 20 on 
page 3, and matrix A20 on page 13 from the paper, applying optimizations 
similar to the current /pj_authlat()/.

Kind regards,

-Jerome

> void pj_authset(double a, double b, double cp[6])
> {
>    // https://arxiv.org/pdf/2212.05818
>    // ∆η(ζ) = S^(L)(ζ) · Cηζ · P^(M) (n) + O(n^L+1)    -- (20)
>    double n = (a - b) / (a + b);  // Third flattening
>    static const double C[21] = // (A20)
>    {
>       4 / 3.0,  4 / 45.0,   -16/35.0,  -2582 /14175.0,  60136 
> /467775.0,    28112932/ 212837625.0,
>                46 / 45.0,  152/945.0, -11966 /14175.0, -21016 / 
> 51975.0,   251310128/ 638512875.0,
>                          3044/2835.0,   3802 /14175.0, -94388 / 
> 66825.0,    -8797648/  10945935.0,
>                                         6059 / 4725.0,  41072 / 
> 93555.0, -1472637812/ 638512875.0,
>                                                        768272 
> /467775.0,  -455935736/ 638512875.0,
> 4210684958/1915538625.0
>    };
>    double p;
>
>    cp[0]  = C[ 0] * n;
>
>    p = n * n;
>    cp[0] += C[ 1] * p;
>    cp[1]  = C[ 6] * p;
>
>    p *= n;
>    cp[0] += C[ 2] * p;
>    cp[1] += C[ 7] * p;
>    cp[2]  = C[11] * p;
>
>    p *= n;
>    cp[0] += C[ 3] * p;
>    cp[1] += C[ 8] * p;
>    cp[2] += C[12] * p;
>    cp[3]  = C[15] * p;
>
>    p *= n;
>    cp[0] += C[ 4] * p;
>    cp[1] += C[ 9] * p;
>    cp[2] += C[13] * p;
>    cp[3] += C[16] * p;
>    cp[4]  = C[18] * p;
>
>    p *= n;
>    cp[0] += C[ 5] * p;
>    cp[1] += C[10] * p;
>    cp[2] += C[14] * p;
>    cp[3] += C[17] * p;
>    cp[4] += C[19] * p;
>    cp[5]  = C[20] * p;
> }
>
> double pj_authlat(double auth, const double * cp)
> {
>    // https://arxiv.org/pdf/2212.05818
>    // ∆η(ζ) = S^(L)(ζ) · Cηζ · P^(M) (n) + O(n^L+1)    -- (20)
>    double a2x = auth + auth, a4x = a2x + a2x, a8x = a4x + a4x;
>    return auth +
>       cp[0] * sin(a2x) +
>       cp[1] * sin(a4x) +
>       cp[2] * sin(a4x + a2x) +
>       cp[3] * sin(a8x) +
>       cp[4] * sin(a8x + a2x) +
>       cp[5] * sin(a8x + a4x);
> }
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20240910/8866f84c/attachment.htm>