<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<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">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">http://www.spatialys.com</a>
My software is free, but my time generally not.</pre>
</body>
</html>