[PROJ] Geodetic to Authalic latitude conversions
Even Rouault
even.rouault at spatialys.com
Wed Sep 11 00:27:55 PDT 2024
Hi Jérôme,
I was wondering if we could not save trigonometric computations noting that
sin(a + b) = sin(a) * cos(b) + cos(a) * sin(b)
and applying that recursively. At the end we would just to have to
compute sin(auth) and cos(auth). But I've no idea about the impact on
numerical stability of using this trigonometric trick.
My other question would be the "backward" compatibility of this change,
like if we would need to provide a "+approx_auth_lat" flag to replicate
current results to use current formulas (like the +approx flag of
+proj=tmerc when more precise formulas were introduced). But I've no
specific opinion if it is needed here. It would be good to have some
examples using let's say WGS84 to see if we're talking about differences
in results of the order of meters, centimeters, micrometers... ?
It is also not immediately obvious to me to correlate your proposed code
with the formulas in the paper by just starring at both at the same
time. Looks like some "interpretation" of the paper has been done. It
might be straightforward for someone who is familiar with the topic
enough, but an error could also be easy to make when transcribing in
code. I'm wondering if it wouldn't be possible to compare numeric
results with the implementation in AuxLatitude.cpp to have some
confidence (I've also tried to compare your formulas with the ones of
that file and can't directly correlate them) ?
Besides those remarks, your proposal generally looks good to my
non-expert eyes on the topic of authalic latitudes :-)
No, we don't use the C++ GeographicLib, just the C port of it from
https://github.com/geographiclib/geographiclib-c/tree/main/src
Even
Le 11/09/2024 à 03:36, Jérôme St-Louis a écrit :
>
> 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
--
http://www.spatialys.com
My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20240911/5b51c0dc/attachment.htm>
More information about the PROJ
mailing list