[PROJ] Geodetic to Authalic latitude conversions
Jérôme St-Louis
jerome at ecere.com
Wed Sep 11 01:47:18 PDT 2024
Thanks Even.
> 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.
>
That sounds brilliant. I will certainly investigate that!
> 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... ?
We should be able to get an idea of the order of magnitude from the
powers of the third flattening which are multiplying the series terms:
n^1: 0.0016792203863837205
n^2: 0.0000028197811060
n^3: 0.0000000047350339
n^4: 0.0000000000079512
With my earlier comparison with the order 3 using the eccentricity
squared used in PROJ vs this method, the PROJ method was accurate to 8
decimals in decimal degrees at ~58 degrees (the greatest delta is around
45 degrees).
I am guessing that improved accuracy for adding the 4th order is
somewhere between millimeters and centimeters.
We will surely find out when we run the unit tests for those projections?
> 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's quite straightforward, but optimized for the precalculations like
PROJ was doing in /auth.cpp/. There is no "interpretation" involved.
Here is my earlier more "direct" mapping of the paper (specifically
Equation 20 on page 3 and matrix Cφξ A20 on page 13), with the loops
still rolled up and no pre-computation, as a single function, which you
might find easier to correlate with both the optimized version as well
as to the summation function in equation 20 on page 3 (and the related
equations 15-19 above on the same page):
> double latAuthalicToGeodetic(double auth)
> {
> // https://arxiv.org/pdf/2212.05818
> // ∆η(ζ) = S^(L)(ζ) · Cηζ · P^(M) (n) + O(n^L+1) -- (20)
> static const double n = (wgs84Major - wgs84Minor) / (wgs84Major +
> wgs84Minor); // Third flattening
> #define ORDER 6
> 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 r = 0;
> int l, m, i = 0, j;
>
> for(l = 1; l <= ORDER; l++)
> {
> double s = sin(2 * l * auth);
> for(m = l; m <= ORDER; m++, i++)
> {
> double c = C[i], p = n;
> for(j = 1; j < m; j++) p *= n;
> r += s * c * p;
> }
> }
> return auth + r;
> }
The nice thing about all this is that this really allows to convert
between most of the auxiliary latitudes (in both directions).
So actually the proposed /pj_authlat() /I shared earlier could equally
apply to parameteric or geocentric or conformal latitudes.
And the only change for the setup in /pj_authset()/ would be to use a
different C matrix based on the type of latitudes involved.
So I would probably make the functions generic so that they could
eventually be used for that purpose as well if needed, with an enum to
pick the right conversion matrix.
I've also ran this through quite a few checks and verified that as a
result, our ISEA3H and ISEA9R discrete global grid hierarchies are now
*actually* equal area :) (with the sole exception of the 12 pentagons
each being 5/6th the area of the hexagons).
e.g. you can see
https://maps.gnosis.earth/ogcapi/dggs/ISEA3H/zones?zone-level=5&compact-zones=false
(or see it on the geojson.io 3D globe
<http://geojson.io/#data=data:text/x-url,https://maps.gnosis.earth/ogcapi/dggs/ISEA3H/zones.geojson%3Fzone-level%3D5%26compact-zones%3DFalse>)
> 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) ?
That is certainly possible...
The order 6 Cφξ (C[phi,xi]) A20 matrix from page 13 is here in
/AuxLatitude.cpp/ (but you need to read it from right-to-left compared
to the rows in my version, which is organized as shown in the paper on
page 13):
https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp#L631
The code itself is a bit overwhelming, which is why I did not actually
look at it much.
My version follows exactly the existing pattern of /pj_auth()/ in PROJ
as well as the functions on page 3 of the paper, which I could
understand much easier, so that should be easy to correlate.
The GeographicLib code may be doing some more fancy things to improve
numerical stability. There might be a lot of C++ boilerplate with the
options to convert to and from any type of auxiliary latitudes, and
compiler options for order 4, 6 and 8.
There's also code in there I believe to iterate to the "exact" value
without using the series, which may be useful either to figure out those
coefficients, and/or for a huge flattening value.
I believe the summation is happening from here with /exact = false/:
https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp#L312
which invokes /Clenshaw()/ here:
https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp#L1319
It is of course still possible that I made a mistake with all this, but
I imagine that updating this and running the tests should help to be
confident about the results.
> No, we don't use the C++ GeographicLib, just the C port of it from
> https://github.com/geographiclib/geographiclib-c/tree/main/src
>
Right, I believe that version predates that whole AuxLatitude module
(the paper is from 2023 and it was first added to GeographicLib 2.2 I
believe).
Thank you!
Kind regards,
-Jerome
On 9/11/24 3:27 AM, Even Rouault wrote:
>
> 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/d1258f14/attachment-0001.htm>
More information about the PROJ
mailing list