[PROJ] Geodetic to Authalic latitude conversions

[Jérôme St-Louis]

Sorry I should correct a few confusing mistakes in my previous e-mail:

     - I was talking mainly about authalic latitude ==> geodetic 
latitude conversion, not the other way around.
/pj_authlat()/ returns the geodetic latitude for a given authalic 
latitude (the name of the function is quite confusing -- I would propose 
to rename to something less ambiguous if it is only used internally)

     - The function currently uses 3 /sin()/ calls for order 3 (not 2, 
made a typo)

     - The ISEA projection does not currently make use of /pj_authlat()/ 
or /pj_qsfn() /for automatically converting authalic / geodetic 
latitudes, but I am proposing to introduce this when using non-spherical 
ellipsoid (just like LAEA of which it is a modified form).

Best,

-Jerome

On 9/10/24 9:19 PM, Jérôme St-Louis via PROJ wrote:
>
> 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);
>> }
>
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20240910/f5effbd3/attachment-0001.htm>