Auxiliary latitudes feature?

[allixender]

Hi,
I am working with DGGS (Discrete Global Grid Systems), like HEALPIX, and many if not most are defined on the sphere. I would like to do the authalic adjustment/conversion between a latitude on the authalic sphere (sphere with the same surface like the Earth) and the ellipsoid (say WGS84). I have found the package Geodesy.jl which makes heavy reference to Charles Karney's GeographicLib C++ library, but doesn't implement AuxLatitude... and also says that the package is in maintenance mode.

Now I found this package that deals with similar geodesic calculations.

There is a Python package that is very explicit on the implementation of Karney's geographiclib features in Python, because Karney doesn't expose much from his C++ to his Python bindings (https://github.com/mrJean1/PyGeodesy/tree/master/pygeodesy/auxilats)

I wonder if we can try this in Julia?
Best,
Alex

[anowacki]

Hi @allixender. This geodesics package is pretty basic at present. Although at one time I planned to incorporate the implementation of Karney's algorithm in GeographicLib.jl into this, I'm afraid it hasn't been my priority, and probably the same will be true of exposing any authalic adjustment code. (I admit in fact that I don't really know what that is exactly.)

I suggest engaging with the JuliaGeo folk such as @asinghvi17 who are more active in this area. You might also want to open a topic about this on the Julia Discourse or on the Slack.

All that said, if you are keen to implement anything, I'd be happy for that to be included and would welcome a PR.

[allixender]

Hi @anowacki thanks for the response. In fact I managed to do a basic port from geographiclib's implementation in the Rust geodesy.rs crate (which seemed to be easier to understand). There are obviously many more functions and I focused only on my use case.

So a perfect sphere with the same surface like our Earth would be called an authalic sphere. Now, many Discrete Global Grid Systems (DGGS) do all their magic on the sphere. The longitudes stay the same, but using the latitudes directly unchanged from the authalic sphere on the ellipsoid (WGS84) is actually introducing a "small" bias. This is where the conversion from the authalic latitude to the geodetic/geographic latitude is need. DGGS are a super interesting topic 🤓 but this authalic sphere vs ellipsoid is just recently becoming more important. Historically, these grids have not been used at sub-kilometre resolutions, because global climate models and and indexing didn't need that. With the advant of H3 etc this has changed and now we want global but high resolution equal-areas cells.

This is a fairly short module. Do you want me to put it somewhere and propose to incorporate into the package? Happy to engage with @asinghvi17

"""
Authalic ↔ Geodetic latitude conversion for WGS84 ellipsoid.

This implementation follows https://github.com/busstoptaktik/geodesy Rust crate which follows GeographicLib's logic

References:
- GeographicLib: https://geographiclib.sourceforge.io/
"""

module AuthalicLatitude

using LinearAlgebra
using StaticArrays # For fixed-size arrays, which are more performant here

struct Ellipsoid
    a::Float64  # Semimajor axis
    f::Float64  # Flattening
end

# The order of the series is a global constant
const POLYNOMIAL_ORDER = 6

# The FourierCoefficients struct. This is mutable because we construct it piece by piece.
# The PolynomialCoefficients struct from your previous request.
# For clarity, I'm redefining it here.
# We'll use SMatrix from StaticArrays for maximum performance.
struct PolynomialCoefficients
    fwd::SMatrix{POLYNOMIAL_ORDER, POLYNOMIAL_ORDER, Float64}
    inv::SMatrix{POLYNOMIAL_ORDER, POLYNOMIAL_ORDER, Float64}
end

mutable struct FourierCoefficients
    fwd::MVector{POLYNOMIAL_ORDER, Float64} # Mutable static vector
    inv::MVector{POLYNOMIAL_ORDER, Float64} # Mutable static vector
    etc::MVector{2, Float64}

    # A convenient constructor to create an empty (zeroed) instance
    FourierCoefficients() = new(
        MVector{POLYNOMIAL_ORDER, Float64}(zeros(POLYNOMIAL_ORDER)),
        MVector{POLYNOMIAL_ORDER, Float64}(zeros(POLYNOMIAL_ORDER)),
        MVector{2, Float64}(zeros(2))
    )
end




# Coefficients for expansion of the normalized meridian arc unit in terms
# of *n²*, the square of the third flattening.
# See [Karney 2010](crate::Bibliography::Kar10) eq. (29)
const MERIDIAN_ARC_COEFFICIENTS = [
    1.0,
    1.0/4.0,
    1.0/64.0,
    1.0/256.0,
    25.0/16384.0,
    49.0/65536.0,
    441.0/1048576.0
]

# WGS84 ellipsoid parameters
const WGS84_A = 6378137.0              # Semi-major axis (meters)
const WGS84_F = 1.0 / 298.257223563    # Flattening
# const WGS84_B = WGS84_A * (1.0 - WGS84_F)  # Semi-minor axis
# const WGS84_E2 = WGS84_F * (2.0 - WGS84_F)  # First eccentricity squared
# const WGS84_E = sqrt(WGS84_E2)         # First eccentricity

const WGS84 = Ellipsoid(WGS84_A, WGS84_F)

# These coefficients are for high-precision latitude conversions in Float64.
# https://github.com/busstoptaktik/geodesy/blob/main/src/ellipsoid/constants.rs
const AUTHALIC = PolynomialCoefficients(
    # Geodetic to authalic: Coefficients for converting 𝜙 to 𝜉.
    # Eq. A19 in Karney (2022)
    [
        -4/3  -4/45   88/315   538/4725   20824/467775   -44732/2837835;
        0     34/45    8/105  -2482/14175  -37192/467775  -12467764/212837625;
        0     0   -1532/2835   -898/14175   54968/467775   100320856/1915538625;
        0     0       0      6007/14175   24496/467775   -5884124/70945875;
        0     0       0          0       -23356/66825   -839792/19348875;
        0     0       0          0           0         570284222/1915538625
    ],

    # Authalic to geodetic: Coefficients for converting 𝜉 to 𝜙.
    # Eq. A20 in Karney (2022)
    [
        4/3   4/45   -16/35   -2582/14175   60136/467775   28112932/212837625;
        0   46/45   152/945   -11966/14175  -21016/51975   251310128/638512875;
        0     0   3044/2835    3802/14175   -94388/66825   -8797648/10945935;
        0     0       0       6059/4725    41072/93555  -1472637812/638512875;
        0     0       0           0       768272/467775   455935736/638512875;
        0     0       0           0            0        4210684958/1915538625
    ]
)

# --- Base "Getter" Functions ---
# These are the direct equivalent of the abstract Rust methods. They just
# access the fields of the struct.

"The semimajor axis, `a`."
semimajor_axis(e::Ellipsoid) = e.a

"The flattening, `f = (a - b) / a`."
flattening(e::Ellipsoid) = e.f


# --- Synonyms ---
# These are convenience functions that call the more descriptive ones.

"Synonym for `semimajor_axis`."
a(e::Ellipsoid) = semimajor_axis(e)

"Synonym for `flattening`."
f(e::Ellipsoid) = flattening(e)


# --- Derived Geometric Properties ---
# These functions are composed from the base functions above.
# Notice how the code is clean, readable, and purely mathematical.

"The semimedian axis, `ay`. For an ellipsoid of revolution, this is equal to `a`."
semimedian_axis(e::Ellipsoid) = semimajor_axis(e)

"The semiminor axis, `b = a * (1 - f)`."
semiminor_axis(e::Ellipsoid) = a(e) * (1.0 - f(e))


# --- Derived Flattening Types ---

"The second flattening, `g = (a - b) / b`."
function second_flattening(e::Ellipsoid)
    # This could be a one-liner, but using an intermediate variable
    # like the Rust code can sometimes be clearer.
    b = semiminor_axis(e)
    return (a(e) - b) / b
end

"The third flattening, `n = (a - b) / (a + b) = f / (2 - f)`."
third_flattening(e::Ellipsoid) = f(e) / (2.0 - f(e))

"""
    horner(arg, coefficients)

Evaluate a polynomial Σ cᵢ ⋅ xⁱ using the efficient Horner's scheme.
`coefficients` should be ordered from the lowest power (c₀, c₁, ...) to the highest.

https://github.com/busstoptaktik/geodesy/blob/main/src/math/series.rs
"""
function horner(arg::Float64, coefficients)
    # The Rust code iterates from the back. Julia's muladd makes this clean.
    # Start with the highest-order coefficient
    val = coefficients[end]
    # Loop downwards from the second to last coefficient
    for i in (length(coefficients) - 1):-1:1
        val = muladd(arg, val, coefficients[i])
    end
    return val
end


"""
    fourier_coefficients(arg, poly_coeffs::PolynomialCoefficients)

Compute Fourier coefficients by evaluating their corresponding Taylor polynomials.
The `arg` is typically the third flattening `n` of an ellipsoid.
"""
function fourier_coefficients(arg::Float64, poly_coeffs::PolynomialCoefficients)
    result = FourierCoefficients()
    for i in 1:POLYNOMIAL_ORDER # Julia loops are 1-based
        # Note: We pass the i-th row of the matrix to horner
        row_fwd = poly_coeffs.fwd[i, :]
        row_inv = poly_coeffs.inv[i, :]
        result.fwd[i] = arg * horner(arg, row_fwd)
        result.inv[i] = arg * horner(arg, row_inv)
    end
    return result
end


"""
    fourier_sin(angle, coeffs)

Evaluate a Fourier sine series: Σ cᵢ * sin(i * angle).
This is the Julia implementation of the implied `fourier::sin` from Rust.
"""
function fourier_sin(angle::Float64, coeffs)
    val = 0.0
    for i in 1:length(coeffs)
        val = muladd(coeffs[i], sin(i * angle), val)
    end
    return val
end


# We need the Ellipsoid struct and associated functions
f(e::Ellipsoid) = e.f
third_flattening(e::Ellipsoid) = f(e) / (2.0 - f(e))

# latitude function from geodesy rust crate
# https://github.com/busstoptaktik/geodesy/blob/main/src/ellipsoid/latitudes.rs

# --- Placeholder for the missing function ---
"""
The Normalized Meridian Arc Unit, *Qn*, is the mean length of one radian
of the meridian. "Normalized", because we measure it in units of the
semimajor axis, *a*.

[König und Weise](crate::Bibliography::KuW51) p.50 (96), p.19 (38b), p.5 (2), here using the extended
version from [Karney 2010](crate::Bibliography::Kar10) eq. (29)
"""
function normalized_meridian_arc_unit(e::Ellipsoid)
    n = third_flattening(e);
    return horner(n * n, MERIDIAN_ARC_COEFFICIENTS) / (1. + n)
end


"""
    coefficients_for_authalic_latitude(e::Ellipsoid, poly_coeffs::PolynomialCoefficients)

Computes the full set of Fourier coefficients needed for authalic latitude conversions
for a given ellipsoid `e`.
"""
function coefficients_for_authalic_latitude(e::Ellipsoid, poly_coeffs::PolynomialCoefficients)
    n = third_flattening(e)
    result = fourier_coefficients(n, poly_coeffs)
    # The Rust `etc[0]` is the first element, which is `etc[1]` in Julia.
    result.etc[1] = normalized_meridian_arc_unit(e)
    return result
end


"""
    latitude_geographic_to_authalic(geographic_latitude, coeffs::FourierCoefficients)

Convert geographic latitude (ϕ) to authalic latitude (ξ).
Note: The ellipsoid `e` is not needed here as its properties are already
baked into the `coeffs`.
"""
function latitude_geographic_to_authalic(geographic_latitude::Float64, coeffs::FourierCoefficients)
    return geographic_latitude + fourier_sin(2.0 * geographic_latitude, coeffs.fwd)
end


"""
    latitude_authalic_to_geographic(authalic_latitude, coeffs::FourierCoefficients)

Convert authalic latitude (ξ) to geographic latitude (ϕ).
"""
function latitude_authalic_to_geographic(authalic_latitude::Float64, coeffs::FourierCoefficients)
    return authalic_latitude + fourier_sin(2.0 * authalic_latitude, coeffs.inv)
end



# --- the conversion functions with ellipsoid and coefficients ---
function _authalic_to_geodetic(beta, ellipsoid, authalic_coeffs; degrees=true)
    # Convert to radians if needed
    beta_rad = degrees ? deg2rad(beta) : beta

    # Compute Fourier coefficients for the given ellipsoid
    fourier_coeffs = coefficients_for_authalic_latitude(ellipsoid, authalic_coeffs)
    # Convert authalic to geodetic
    phi_rad = latitude_authalic_to_geographic(beta_rad, fourier_coeffs)

    # Convert back to degrees if needed
    return degrees ? rad2deg(phi_rad) : phi_rad
end

function _geodetic_to_authalic(phi, ellipsoid, authalic_coeffs; degrees=true)
    # Convert to radians if needed
    phi_rad = degrees ? deg2rad(phi) : phi

    # Compute Fourier coefficients for the given ellipsoid
    fourier_coeffs = coefficients_for_authalic_latitude(ellipsoid, authalic_coeffs)
    # Convert geodetic to authalic
    beta_rad = latitude_geographic_to_authalic(phi_rad, fourier_coeffs)

    # Convert back to degrees if needed
    return degrees ? rad2deg(beta_rad) : beta_rad
    
end

# --- convenience functions ---
function authalic_to_geodetic(beta; degrees=true)
    # Convert to radians if needed
    beta_rad = degrees ? deg2rad(beta) : beta

    # Compute Fourier coefficients for the given ellipsoid
    fourier_coeffs = coefficients_for_authalic_latitude(WGS84, AUTHALIC)
    # Convert authalic to geodetic
    phi_rad = latitude_authalic_to_geographic(beta_rad, fourier_coeffs)

    # Convert back to degrees if needed
    return degrees ? rad2deg(phi_rad) : phi_rad
end

function geodetic_to_authalic(phi; degrees=true)
    # Convert to radians if needed
    phi_rad = degrees ? deg2rad(phi) : phi

    # Compute Fourier coefficients for the given ellipsoid
    fourier_coeffs = coefficients_for_authalic_latitude(WGS84, AUTHALIC)
    # Convert geodetic to authalic
    beta_rad = latitude_geographic_to_authalic(phi_rad, fourier_coeffs)

    # Convert back to degrees if needed
    return degrees ? rad2deg(beta_rad) : beta_rad
    
end

# --- Exports ---
# default conversion functions for authalic <-> WGS84
export authalic_to_geodetic, geodetic_to_authalic
# For manually using ellippsoid and coefficients
export _authalic_to_geodetic, _geodetic_to_authalic, WGS84, AUTHALIC

end  # module AuthalicLatitude

[danlooo]

@allixender thanks for your code on authalic latitudes! Would it also make sense to add this directly to e.g.
ISEA in PROJ and Proj.jl? We are currently working on implementing native ISEA in Julia. This would help improving the accuracy of DGGS implementations in Julia.

[asinghvi17]

As another collaborator on https://github.com/meggart/SphericalSpatialTrees.jl, I think we'd also be happy to have more such projections implemented in Julia where you get x, y, face_index out instead of just projecting all the faces onto a plane. That allows a lot more flexibility when getting datasets from other CRS into a DGGS.

[anowacki]

Thank you for the example code, @allixender, and for the explanation. I'm sure the PROJ people would welcome this, as would SphericalSpatialTrees.jl. I wonder where the best home for such a basic/fundamental transformation would be.

[asinghvi17]

I guess the right place is probably Geodesy.jl? But we could keep it in SphericalSpatialTrees for now while we iterate on it. It strikes me that you don't actually need mutable vectors in the fourier coefficients, for instance (though the code is probably clearer with the MVectors).

This should also probably be packaged up into a CoordinateTransformation so your coefficients get precomputed once for the ellipsoid, and then you're off to the races. Might also want to make the number type variable (e.g. MultiFloats.Float64x2) such that it becomes possible to compute in more precise numbers (I have some applications on the meter scale, for instance).

[allixender]

Hi @asinghvi17, thanks for the feedback and directions. I am maybe a little overwhelmed by the options.

I have seen @danlooo opening an issue on SphericalSpatialTrees. Shall I try to make a pull request?

[danlooo]

@allixender PRs are always welcome :)

[asinghvi17]

Yeah, I think let's get a PR up and then workshop the performance / framing a bit once that's up. I can also commit to adding some tests against Proj etc.