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@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/proj
_______________________________________________
PROJ mailing list
PROJ@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/proj
