healpix_analyse.alm_latlon#

Spherical harmonic transform (map → alm → Cl) for maps defined on arbitrary iso-latitude ring grids (ERA5/ECMWF Gaussian grids, regular lat/lon grids, HEALPix pixel coordinates, or any custom ring-structured grid).

This module is entirely self-contained: it depends only on numpy and torch. It does not depend on the foscat package.

Public API#

build_rings_from_latlon(lat, lon, atol, convention)

Sort a flat list of pixels into iso-latitude rings.

compute_weights(ring_theta, ring_phi_list, ring_counts, quadrature)

Compute per-pixel quadrature weights (steradians).

map2alm_latlon(im, ring_theta, ring_phi_list, ring_counts, lmax, weights, quadrature)

Map → spherical harmonic coefficients a_lm.

alm2map_latlon(alm, ring_theta, ring_phi_list, ring_counts, lmax)

Spherical harmonic coefficients → reconstructed map (synthesis).

anafast_latlon(im, ring_theta, ring_phi_list, ring_counts, lmax, weights, quadrature)

Map → angular power spectrum Cl.

grid_summary(ring_theta, ring_phi_list, ring_counts, lmax)

Print a diagnostic summary of the grid geometry.

Mathematical convention#

The transform follows the standard orthonormal convention (same as healpy):

a_lm = ∫ f(θ,φ) Y_lm*(θ,φ) dΩ

Y_lm(θ,φ) = sqrt((2l+1)/(4π)) * P̃_lm(cosθ) * exp(imφ)

C_l = 1/(2l+1) * [|a_l0|² + 2 Σ_{m=1}^{l} |a_lm|²]

where P̃_lm are the normalised associated Legendre polynomials.

Attributes#

ArrayLike

Functions#

build_rings_from_latlon(→ Tuple[numpy.ndarray, ...)

Group a flat list of N pixels into iso-latitude rings.

compute_weights(→ numpy.ndarray)

Compute per-pixel quadrature weights in steradians.

comp_tf_latlon(→ torch.Tensor)

Compute the weighted longitude DFT of a map, ring by ring.

map2alm_latlon(→ torch.Tensor)

Compute spherical harmonic coefficients a_lm from a map.

alm2map_latlon(→ torch.Tensor)

Synthesise a map from spherical harmonic coefficients (inverse transform).

anafast_latlon(→ torch.Tensor)

Estimate the angular power spectrum C_l from a map.

grid_summary(→ dict)

Print a human-readable summary of the grid and estimated compute cost.

Module Contents#

healpix_analyse.alm_latlon.ArrayLike#
healpix_analyse.alm_latlon.build_rings_from_latlon(lat: ArrayLike, lon: ArrayLike, atol: float = 1e-10, convention: str = 'colatitude_rad') → Tuple[numpy.ndarray, List[numpy.ndarray], numpy.ndarray, numpy.ndarray][source]#

Group a flat list of N pixels into iso-latitude rings.

This is the mandatory first step before any transform call. The function sorts pixels by colatitude (then longitude) and identifies contiguous groups sharing the same colatitude within tolerance atol.

Parameters:

  • lat (array-like, shape [N]) – Latitude or colatitude of each pixel. Unit and sign depend on convention (see below).

  • lon (array-like, shape [N]) – Longitude of each pixel. Unit depends on convention.

  • atol (float, default 1e-10) – Angular tolerance in radians to merge two pixels into the same ring. Increase to ~1e-6 when coordinates come from single-precision data.

  • convention (str, default "colatitude_rad") – Coordinate convention. Accepted values:

    "colatitude_rad"

    lat = colatitude θ in radians, 0 → π lon = longitude φ in radians, 0 → 2π

    "colatitude_deg"

    lat = colatitude θ in degrees, 0° → 180° lon = longitude φ in degrees, 0° → 360°

    "geographic_rad"

    lat = geographic latitude in radians, -π/2 → +π/2 lon = longitude in radians, -π → +π or 0 → 2π

    "geographic_deg"

    lat = geographic latitude in degrees, -90° → +90° lon = longitude in degrees, -180° → +180° or 0° → 360°

Returns:

  • ring_theta (np.ndarray, shape [R]) – Colatitude θ in radians of each ring, sorted in ascending order.

  • ring_phi_list (list of np.ndarray) – Longitudes φ in radians for each ring.

  • ring_counts (np.ndarray, shape [R], dtype int64) – Number of pixels N_r in each ring.

  • sort_idx (np.ndarray, shape [N], dtype int64) – Permutation index. Apply as im_sorted = im[sort_idx] before passing a map to any transform function.

healpix_analyse.alm_latlon.compute_weights(ring_theta: numpy.ndarray, ring_phi_list: List[numpy.ndarray], ring_counts: numpy.ndarray, quadrature: str = 'trapeze') → numpy.ndarray[source]#

Compute per-pixel quadrature weights in steradians.

The weights satisfy ∫ f dΩ ≈ Σ_i f_i * w_i when summed over all pixels in ring order.

Parameters:

  • ring_theta (np.ndarray, shape [R]) – Colatitudes in radians (output of build_rings_from_latlon).

  • ring_phi_list (list of np.ndarray) – Longitudes in radians per ring.

  • ring_counts (np.ndarray, shape [R]) – Number of pixels per ring.

  • quadrature (str, default "trapeze") – Quadrature rule applied in the colatitude direction.

    "trapeze"

    Trapezoidal rule: w_θ = sin(θ_r) * Δθ_r. Accurate for grids with regular θ spacing.

    "gauss_legendre"

    Exact Gauss-Legendre weights for nodes x_r = cos(θ_r). Mandatory for Gaussian grids (ERA5, ECMWF, IFS, ARPEGE). The quadrature is exact up to ℓ ≈ 2R - 1. A UserWarning is raised if the provided colatitudes do not match Gauss-Legendre nodes within a tolerance of 1e-6.

    "equal_area"

    Uniform weight 4π / N_total. Appropriate for equal-area pixelisations such as HEALPix.

Returns:

weights (np.ndarray, shape [N_total], dtype float64) – Per-pixel quadrature weights in steradians, in ring order.

healpix_analyse.alm_latlon.comp_tf_latlon(im: torch.Tensor, ring_phi_list: List[numpy.ndarray], ring_counts: numpy.ndarray, pixel_weights: numpy.ndarray, mmax: int, cdtype: torch.dtype = torch.complex128) → torch.Tensor[source]#

Compute the weighted longitude DFT of a map, ring by ring.

For each ring r with N_r pixels at longitudes φ_j:

F_r(m) = Σ_j f(θ_r, φ_j) * w_j * exp(-i m φ_j)

For uniformly-spaced rings the computation uses torch.fft.rfft followed by a phase-shift exp(-i m φ_0), matching the FOSCAT comp_tf convention. For irregular rings a direct DFT is used.

Parameters:

  • im (torch.Tensor, shape [..., N_total]) – Map values in ring order (apply sort_idx before calling).

  • ring_phi_list (list of np.ndarray) – Longitudes in radians per ring.

  • ring_counts (np.ndarray, shape [R]) – Number of pixels per ring.

  • pixel_weights (np.ndarray, shape [N_total]) – Per-pixel quadrature weights (output of compute_weights).

  • mmax (int) – Maximum azimuthal order to compute.

  • cdtype (torch.dtype, default torch.complex128) – Complex dtype for the output tensor.

Returns:

ft (torch.Tensor, shape [..., R, mmax+1], complex) – Ring-by-ring Fourier coefficients.

healpix_analyse.alm_latlon.map2alm_latlon(im: numpy.ndarray | torch.Tensor, ring_theta: numpy.ndarray, ring_phi_list: List[numpy.ndarray] | numpy.ndarray, ring_counts: numpy.ndarray, lmax: int = 64, weights: numpy.ndarray | str | None = None, quadrature: str = 'trapeze', limit_range: float = 10000000000.0) → torch.Tensor[source]#

Compute spherical harmonic coefficients a_lm from a map.

The transform follows the standard orthonormal convention (same as healpy):

a_lm = ∫ f(θ,φ) Y_lm*(θ,φ) dΩ

Parameters:

  • im (array-like or torch.Tensor, shape [..., N_total]) – Map values in ring order (apply sort_idx from build_rings_from_latlon first). An optional leading batch dimension is supported.

  • ring_theta (np.ndarray, shape [R]) – Colatitudes in radians (output of build_rings_from_latlon).

  • ring_phi_list (list of np.ndarray or flat np.ndarray [N_total]) – Longitudes in radians per ring.

  • ring_counts (np.ndarray, shape [R]) – Number of pixels per ring.

  • lmax (int, default 24) – Maximum multipole ℓ_max.

  • weights (np.ndarray [N_total], None, or "uniform") – Per-pixel quadrature weights.

    • None : computed automatically via compute_weights(quadrature=quadrature).

    • "uniform" : w_i = 1/N_total (comparable to healpy.map2alm).

    • array : user-supplied weights (quadrature is then ignored).

  • quadrature (str, default "trapeze") – Quadrature rule passed to compute_weights. Accepted values: "trapeze", "gauss_legendre", "equal_area". Ignored when weights is not None.

    Note

    Use "gauss_legendre" for ERA5/ECMWF Gaussian grids to avoid oscillations at high ℓ.

  • limit_range (float, default 1e10) – Overflow guard in the Legendre recurrence.

Returns:

alm_out (torch.Tensor, shape [..., n_alm], complex128) – Spherical harmonic coefficients where n_alm = Σ_{m=0}^{lmax} (lmax - m + 1).

Memory layout (same as healpy):

[m=0: ℓ=0..lmax | m=1: ℓ=1..lmax | … | m=lmax: ℓ=lmax]

healpix_analyse.alm_latlon.alm2map_latlon(alm: torch.Tensor, ring_theta: numpy.ndarray, ring_phi_list: List[numpy.ndarray] | numpy.ndarray, ring_counts: numpy.ndarray, lmax: int = 24, limit_range: float = 10000000000.0) → torch.Tensor[source]#

Synthesise a map from spherical harmonic coefficients (inverse transform).

This is the adjoint of map2alm_latlon. The output map is evaluated at the same grid positions as the input to the analysis, but it can also be used to synthesise on a different target grid.

Parameters:

  • alm (torch.Tensor, shape [..., n_alm], complex) – Spherical harmonic coefficients in the layout produced by map2alm_latlon.

  • ring_theta (np.ndarray, shape [R]) – Colatitudes in radians of the target grid.

  • ring_phi_list (list of np.ndarray or flat np.ndarray [N_total]) – Longitudes in radians of the target grid.

  • ring_counts (np.ndarray, shape [R]) – Number of pixels per ring of the target grid.

  • lmax (int, default 24) – Maximum multipole used when the alm were computed.

  • limit_range (float, default 1e10) – Overflow guard in the Legendre recurrence.

Returns:

im_out (torch.Tensor, shape [..., N_total], float64) – Reconstructed map values in ring order.

healpix_analyse.alm_latlon.anafast_latlon(im: numpy.ndarray | torch.Tensor, ring_theta: numpy.ndarray, ring_phi_list: List[numpy.ndarray] | numpy.ndarray, ring_counts: numpy.ndarray, lmax: int = 64, weights: numpy.ndarray | str | None = None, quadrature: str = 'trapeze', limit_range: float = 10000000000.0) → torch.Tensor[source]#

Estimate the angular power spectrum C_l from a map.

Internally calls map2alm_latlon then accumulates:

C_l = 1/(2l+1) * [|a_l0|² + 2 Σ_{m=1}^{l} |a_lm|²]

All parameters are identical to map2alm_latlon.

Returns:

cl (torch.Tensor, shape [..., lmax+1], float64) – Angular power spectrum.

healpix_analyse.alm_latlon.grid_summary(ring_theta: numpy.ndarray, ring_phi_list: List[numpy.ndarray] | numpy.ndarray, ring_counts: numpy.ndarray, lmax: int = 24) → dict[source]#

Print a human-readable summary of the grid and estimated compute cost.

Parameters:

  • ring_theta (np.ndarray, shape [R])

  • ring_phi_list (list of np.ndarray)

  • ring_counts (np.ndarray, shape [R])

  • lmax (int)

Returns:

info (dict) – Dictionary with the same fields printed to stdout.