<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Thanks Even.</p>
<p>
<blockquote type="cite">
<p>I was wondering if we could not save trigonometric
computations noting that</p>
<p>sin(a + b) = sin(a) * cos(b) + cos(a) * sin(b)</p>
<p>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.</p>
</blockquote>
That sounds brilliant. I will certainly investigate that!<br>
</p>
<div class="moz-signature">
<blockquote type="cite">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... ?</blockquote>
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:</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">n^1: 0.0016792203863837205<br>
n^2: 0.0000028197811060<br>
n^3: 0.0000000047350339<br>
n^4: 0.0000000000079512</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">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).<br>
<br>
I am guessing that improved accuracy for adding the 4th order is
somewhere between millimeters and centimeters.<br>
We will surely find out when we run the unit tests for those
projections?</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">
<blockquote type="cite">
<p>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. </p>
</blockquote>
</div>
<div class="moz-signature">It's quite straightforward, but optimized
for the precalculations like PROJ was doing in <i>auth.cpp</i>.
There is no "interpretation" involved.<br>
</div>
<div class="moz-signature">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):<br>
</div>
<div class="moz-signature"><font size="4"><br>
</font></div>
<div class="moz-signature">
<blockquote type="cite"><font size="4" face="monospace">double
latAuthalicToGeodetic(double auth) <br>
{ <br>
// <a class="moz-txt-link-freetext"
href="https://arxiv.org/pdf/2212.05818">https://arxiv.org/pdf/2212.05818</a>
<br>
// ∆η(ζ) = S^(L)(ζ) · Cηζ · P^(M) (n) + O(n^L+1) -- (20)
<br>
static const double n = (wgs84Major - wgs84Minor) /
(wgs84Major + wgs84Minor); // Third flattening <br>
#define ORDER 6 <br>
static const double C[21] = // (A20) <br>
{ <br>
4 / 3.0, 4 / 45.0, -16/35.0, -2582 /14175.0, 60136
/467775.0, 28112932/ 212837625.0, <br>
46 / 45.0, 152/945.0, -11966 /14175.0, -21016
/ 51975.0, 251310128/ 638512875.0, <br>
3044/2835.0, 3802 /14175.0, -94388
/ 66825.0, -8797648/ 10945935.0, <br>
6059 / 4725.0, 41072
/ 93555.0, -1472637812/ 638512875.0, <br>
768272
/467775.0, -455935736/ 638512875.0, <br>
4210684958/1915538625.0 <br>
}; <br>
double r = 0; <br>
int l, m, i = 0, j; <br>
<br>
for(l = 1; l <= ORDER; l++) <br>
{ <br>
double s = sin(2 * l * auth); <br>
for(m = l; m <= ORDER; m++, i++) <br>
{ <br>
double c = C[i], p = n; <br>
for(j = 1; j < m; j++) p *= n; <br>
r += s * c * p; <br>
} <br>
} <br>
return auth + r; <br>
} </font></blockquote>
</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">The nice thing about all this is that
this really allows to convert between most of the auxiliary
latitudes (in both directions).<br>
So actually the proposed <i>pj_authlat() </i>I shared earlier
could equally apply to parameteric or geocentric or conformal
latitudes.<br>
And the only change for the setup in <i>pj_authset()</i> would be
to use a different C matrix based on the type of latitudes
involved.<br>
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.<br>
</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">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 <b>actually</b> equal
area :) (with the sole exception of the 12 pentagons each being
5/6th the area of the hexagons).<br>
<br>
</div>
<div class="moz-signature">e.g. you can see
<a class="moz-txt-link-freetext" href="https://maps.gnosis.earth/ogcapi/dggs/ISEA3H/zones?zone-level=5&compact-zones=false">https://maps.gnosis.earth/ogcapi/dggs/ISEA3H/zones?zone-level=5&compact-zones=false</a>
(or <a
href="http://geojson.io/#data=data:text/x-url,https://maps.gnosis.earth/ogcapi/dggs/ISEA3H/zones.geojson%3Fzone-level%3D5%26compact-zones%3DFalse">see
it on the geojson.io 3D globe</a>)<br>
</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">
<blockquote type="cite">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) ?<br>
</blockquote>
</div>
<div class="moz-signature">That is certainly possible...</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">The order 6 Cφξ (C[phi,xi]) A20 matrix
from page 13 is here in <i>AuxLatitude.cpp</i> (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):<br>
</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">
<a class="moz-txt-link-freetext" href="https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp#L631">https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp#L631</a></div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">The code itself is a bit overwhelming,
which is why I did not actually look at it much.<br>
My version follows exactly the existing pattern of <i>pj_auth()</i>
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.<br>
</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">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.<br>
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.<br>
I believe the summation is happening from here with <i>exact =
false</i>:</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">
<a class="moz-txt-link-freetext" href="https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp#L312">https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp#L312</a></div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">which invokes <i>Clenshaw()</i> here:</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">
<a class="moz-txt-link-freetext" href="https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp#L1319">https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp#L1319</a><br>
</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">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.<br>
</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">
<blockquote type="cite">
<p>No, we don't use the C++ GeographicLib, just the C port of it
from <a class="moz-txt-link-freetext"
href="https://github.com/geographiclib/geographiclib-c/tree/main/src">https://github.com/geographiclib/geographiclib-c/tree/main/src</a></p>
</blockquote>
</div>
<div class="moz-signature">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).</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">Thank you!</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">Kind regards,</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">-Jerome</div>
<div class="moz-signature"><br>
</div>
<div class="moz-cite-prefix">On 9/11/24 3:27 AM, Even Rouault wrote:<br>
</div>
<blockquote type="cite"
cite="mid:ac408a8f-9044-40ad-bdbc-debd10e162f6@spatialys.com">
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<p>Hi Jérôme,</p>
<p>I was wondering if we could not save trigonometric computations
noting that</p>
<p>sin(a + b) = sin(a) * cos(b) + cos(a) * sin(b)</p>
<p>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.</p>
<p>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... ?</p>
<p>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) ?<br>
</p>
<p>Besides those remarks, your proposal generally looks good to my
non-expert eyes on the topic of authalic latitudes :-)</p>
<p>No, we don't use the C++ GeographicLib, just the C port of it
from <a class="moz-txt-link-freetext"
href="https://github.com/geographiclib/geographiclib-c/tree/main/src"
moz-do-not-send="true">https://github.com/geographiclib/geographiclib-c/tree/main/src</a><br>
</p>
<p>Even<br>
</p>
<div class="moz-cite-prefix">Le 11/09/2024 à 03:36, Jérôme
St-Louis a écrit :<br>
</div>
<blockquote type="cite"
cite="mid:450be707-b926-4e03-9025-f1b66badcd64@ecere.com">
<meta http-equiv="Content-Type"
content="text/html; charset=UTF-8">
<p>Sorry I should correct a few confusing mistakes in my
previous e-mail:</p>
<p> - I was talking mainly about authalic latitude ==>
geodetic latitude conversion, not the other way around.<br>
<i>pj_authlat()</i> 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)<br>
<br>
- The function currently uses 3 <i>sin()</i> calls for
order 3 (not 2, made a typo)<br>
<br>
- The ISEA projection does not currently make use of <i>pj_authlat()</i>
or <i>pj_qsfn() </i>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).</p>
<p>Best,</p>
<p>-Jerome</p>
<div class="moz-cite-prefix">On 9/10/24 9:19 PM, Jérôme St-Louis
via PROJ wrote:<br>
</div>
<blockquote type="cite"
cite="mid:9ed4372a-f6a3-4420-b32b-232c3c4638b0@ecere.com">
<meta http-equiv="content-type"
content="text/html; charset=UTF-8">
<p>Hi Even, all,</p>
<p>As suggested in <a class="moz-txt-link-freetext"
href="https://github.com/OSGeo/PROJ/pull/4211"
moz-do-not-send="true">https://github.com/OSGeo/PROJ/pull/4211</a>
, 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 (<i>pj_authlat()</i> for
the geodetic ==> authalic), which is used for ISEA, but
also by:</p>
<p>- Equal Area Cylindrical (cea.cpp)<br>
- EqualEarth (eqearth.cpp)<br>
- HEALPix (healpix.cpp) -- which might also imply rHEALPix<br>
- Lambert Azimuthal Equal Area (laea.cpp)</p>
<p>The current version of <i>pj_authlat()</i> uses a series
expansion based on the eccentricity squared (I believe this
is based on a <a
href="https://en.wikipedia.org/wiki/Lagrange_inversion_theorem"
moz-do-not-send="true">Lagrange reversion</a> as <a
href="https://en.wikipedia.org/wiki/Latitude#Inverse_formulae_and_series"
moz-do-not-send="true">mentioned here</a>), and Charles
Karney in <a href="https://arxiv.org/pdf/2212.05818"
moz-do-not-send="true">this paper</a> establishes that:<br>
</p>
<p> </p>
<blockquote type="cite"><i>series expansions using the third
flattening as the small parameter uniformly result in
smaller truncation errors compared to using the
eccentricity.</i></blockquote>
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.<br>
The current order 3, which requires 2 <i>sin()</i> calls,
gives precision up to roughly 8 decimals.<br>
Order 6 would mean twice the number of <i>sin()</i> calls,
and the 8 decimals should be enough to <a
href="https://xkcd.com/2170/" moz-do-not-send="true">point
to Waldo on a page</a>, however C.K. argues that:
<p> </p>
<blockquote type="cite"><i>Modern geodetic libraries should
strive to achieve full double-precision accuracy </i>(page
4)<i><br>
</i></blockquote>
<p>and I would much agree with that.<br>
</p>
<p>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 <a
href="https://github.com/OSGeo/PROJ/pull/4211#pullrequestreview-2233903246"
moz-do-not-send="true">discussed here</a>), and
adjustments for tests etc.</p>
<p>PROJ currently does not use the full GeographicLib which
already includes code for conversion between different
auxiliary latitudes, does it?</p>
<p> ( <a class="moz-txt-link-freetext"
href="https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp"
moz-do-not-send="true">https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp</a>
)<br>
</p>
<p>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 <i>pj_authlat()</i>.</p>
<p>Kind regards,</p>
<p>-Jerome<br>
</p>
<blockquote type="cite"><font size="5" face="monospace">void
pj_authset(double a, double b, double cp[6]) <br>
{ <br>
// <a class="moz-txt-link-freetext"
href="https://arxiv.org/pdf/2212.05818"
moz-do-not-send="true">https://arxiv.org/pdf/2212.05818</a>
<br>
// ∆η(ζ) = S^(L)(ζ) · Cηζ · P^(M) (n) + O(n^L+1) --
(20) <br>
double n = (a - b) / (a + b); // Third flattening <br>
static const double C[21] = // (A20) <br>
{ <br>
4 / 3.0, 4 / 45.0, -16/35.0, -2582 /14175.0,
60136 /467775.0, 28112932/ 212837625.0, <br>
46 / 45.0, 152/945.0, -11966 /14175.0,
-21016 / 51975.0, 251310128/ 638512875.0, <br>
3044/2835.0, 3802 /14175.0,
-94388 / 66825.0, -8797648/ 10945935.0, <br>
6059 / 4725.0,
41072 / 93555.0, -1472637812/ 638512875.0, <br>
768272 /467775.0, -455935736/ 638512875.0, <br>
4210684958/1915538625.0 <br>
}; <br>
double p; <br>
<br>
cp[0] = C[ 0] * n; <br>
<br>
p = n * n; <br>
cp[0] += C[ 1] * p; <br>
cp[1] = C[ 6] * p; <br>
<br>
p *= n; <br>
cp[0] += C[ 2] * p; <br>
cp[1] += C[ 7] * p; <br>
cp[2] = C[11] * p; <br>
<br>
p *= n; <br>
cp[0] += C[ 3] * p; <br>
cp[1] += C[ 8] * p; <br>
cp[2] += C[12] * p; <br>
cp[3] = C[15] * p; <br>
<br>
p *= n; <br>
cp[0] += C[ 4] * p; <br>
cp[1] += C[ 9] * p; <br>
cp[2] += C[13] * p; <br>
cp[3] += C[16] * p; <br>
cp[4] = C[18] * p; <br>
<br>
p *= n; <br>
cp[0] += C[ 5] * p; <br>
cp[1] += C[10] * p; <br>
cp[2] += C[14] * p; <br>
cp[3] += C[17] * p; <br>
cp[4] += C[19] * p; <br>
cp[5] = C[20] * p; <br>
} <br>
<br>
double pj_authlat(double auth, const double * cp) <br>
{ <br>
// <a class="moz-txt-link-freetext"
href="https://arxiv.org/pdf/2212.05818"
moz-do-not-send="true">https://arxiv.org/pdf/2212.05818</a>
<br>
// ∆η(ζ) = S^(L)(ζ) · Cηζ · P^(M) (n) + O(n^L+1) --
(20) <br>
double a2x = auth + auth, a4x = a2x + a2x, a8x = a4x +
a4x; <br>
return auth + <br>
cp[0] * sin(a2x) + <br>
cp[1] * sin(a4x) + <br>
cp[2] * sin(a4x + a2x) + <br>
cp[3] * sin(a8x) + <br>
cp[4] * sin(a8x + a2x) + <br>
cp[5] * sin(a8x + a4x); <br>
} </font><br>
</blockquote>
<br>
<fieldset class="moz-mime-attachment-header"></fieldset>
<pre class="moz-quote-pre" wrap="">_______________________________________________
PROJ mailing list
<a class="moz-txt-link-abbreviated moz-txt-link-freetext"
href="mailto:PROJ@lists.osgeo.org" moz-do-not-send="true">PROJ@lists.osgeo.org</a>
<a class="moz-txt-link-freetext"
href="https://lists.osgeo.org/mailman/listinfo/proj"
moz-do-not-send="true">https://lists.osgeo.org/mailman/listinfo/proj</a>
</pre>
</blockquote>
</blockquote>
<pre class="moz-signature" cols="72">--
<a class="moz-txt-link-freetext" href="http://www.spatialys.com"
moz-do-not-send="true">http://www.spatialys.com</a>
My software is free, but my time generally not.</pre>
</blockquote>
</body>
</html>