foscat.alm_latlon#

Spherical-harmonic transform for maps defined on an arbitrary colatitude / longitude grid organised into rings.

Differences with respect to foscat.alm.alm#

  • No dependency on HEALPix for pixel positioning.

  • Rings may have arbitrary colatitudes (not HEALPix colatitudes) and arbitrary longitudes (not necessarily uniform).

  • Longitude step: direct DFT when ring φ values are irregular; FFT + phase shift (like alm.comp_tf) when φ are uniformly spaced.

  • Colatitude step: same Legendre recurrence as alm.compute_legendre_m, evaluated at the cosines of the provided colatitudes.

  • Quadrature weights: trapezoidal in θ (sin θ · Δθ) × uniform in φ (2π/N_r) by default, or user-supplied weight array.

Main API#

build_rings_from_latlon(lat, lon, atol)

Groups a flat list of (lat, lon) into rings of equal colatitude.

compute_weights(ring_theta, ring_phi_list, ring_counts)

Computes the quadrature weights per pixel (steradians).

comp_tf_latlon(im, ring_phi_list, ring_counts, pixel_weights, mmax)

Ring-weighted DFT → ft[…, R, mmax+1].

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

Map → alm.

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

Map → Cl.

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

alm → map (synthesis).

Minimal example#

import numpy as np from foscat.alm_latlon import alm_latlon

# Regular grid (ntheta=64 rings × nphi=128 pixels per ring) ntheta, nphi = 64, 128 theta_1d = np.linspace(np.pi / (2*ntheta), np.pi - np.pi/(2*ntheta), ntheta) phi_1d = np.linspace(0, 2*np.pi*(1 - 1/nphi), nphi)

lat = np.repeat(theta_1d, nphi) # colatitude of each pixel lon = np.tile(phi_1d, ntheta) # longitude of each pixel

ring_theta, ring_phi_list, ring_counts, sort_idx =

alm_latlon.build_rings_from_latlon(lat, lon)

obj = alm_latlon(lmax=32) im = np.random.randn(ntheta * nphi)

alm_coeffs = obj.map2alm_latlon(

im[sort_idx], ring_theta, ring_phi_list, ring_counts

) cl = obj.anafast_latlon(

im[sort_idx], ring_theta, ring_phi_list, ring_counts

)

Classes#

alm_latlon

Spherical-harmonic transform on an arbitrary lat/lon grid organised into rings.

Module Contents#

class foscat.alm_latlon.alm_latlon(backend=None, lmax=24, limit_range=10000000000.0)[source]#

Bases: foscat.alm.alm

Spherical-harmonic transform on an arbitrary lat/lon grid organised into rings.

static build_rings_from_latlon(lat, lon, atol=1e-10, convention='colatitude_rad')[source]#

Group a flat list of pixels into rings of equal colatitude.

Parameters:

  • lat (array [N]  angular coordinate of each pixel (see convention).)

  • lon (array [N]  longitudinal coordinate of each pixel (see convention).)

  • atol (tolerance in radians for grouping two pixels into the same ring.)

  • convention (str  format of the input coordinates.) –

    ‘colatitude_rad’ (default)

    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°

  • All conventions are converted internally to colatitude + longitude

  • in radians before processing.

Returns:

  • ring_theta (ndarray [R]            colatitude θ (radians) per ring)

  • ring_phi_list (list[ndarray [N_r]]    longitudes  φ (radians) per ring)

  • ring_counts (ndarray int64 [R]      number of pixels per ring)

  • sort_idx (ndarray int64 [N]      permutation  im_sorted = im[sort_idx])

static compute_weights(ring_theta, ring_phi_list, ring_counts, quadrature='trapeze')[source]#

Computes the quadrature weights per pixel (steradians).

Parameters:

  • ring_theta ([R] colatitudes in radians)

  • ring_phi_list (list[ndarray]  longitudes per ring)

  • ring_counts ([R] number of pixels per ring)

  • quadrature (str  quadrature method in θ.) –

    ‘trapeze’ (default)

    Trapezoidal rule: w_θ = sin(θ_r) × Δθ_r Suitable for regular θ grids.

    ‘gauss_legendre’

    Exact Gauss-Legendre weights. Required for Gaussian grids (ERA5, ECMWF, IFS, ARPEGE…) where the colatitudes are the zeros of P_R(cos θ). The integral ∫f dΩ = ∫f dφ dx (x = cos θ) is then exact up to ℓ ≈ 2R-1.

    ‘equal_area’

    Equal weights: 4π / N_total. For equal-area grids (HEALPix).

Returns:

weights (ndarray float64 [N_total])

comp_tf_latlon(im, ring_phi_list, ring_counts, pixel_weights, mmax)[source]#

Pixel-weighted DFT for each ring.

For a ring with uniformly spaced φ values, uses the FFT + phase shift (same logic as alm.comp_tf). For an irregular ring, performs a direct DFT.

Parameters:

  • im ([..., N_total]  map (backend tensor or ndarray))

  • ring_phi_list (list[ndarray]   longitudes per ring)

  • ring_counts (ndarray int64   number of pixels per ring)

  • pixel_weights (ndarray float64 [N_total]  quadrature weights)

  • mmax (int  maximum frequency)

Returns:

ft (tensor [..., R, mmax+1]  complex)

map2alm_latlon(im, ring_theta, ring_phi_list, ring_counts, lmax=None, weights=None, quadrature='trapeze')[source]#

Compute the alm coefficients of a map defined on an arbitrary grid organised into rings.

Parameters:

  • im ([..., N_total]  map values, ordered ring by ring) – (use sort_idx from build_rings_from_latlon if the map is initially in arbitrary order)

  • ring_theta ([R] colatitude (radians) of each ring)

  • ring_phi_list (list[ndarray] or ndarray [N_total]) – longitudes (radians) per ring (or flat array)

  • ring_counts ([R] number of pixels per ring)

  • lmax (maximum multipole (default: self.lmax))

  • weights ([N_total] per-pixel weights (steradians).) – If None, trapezoidal quadrature is computed automatically. Pass weights=’uniform’ to use 1/N_total (alm.map2alm convention).

Returns:

alm_out ([..., n_alm]  complex, n_alm = Σ_``{m=0}``:py:class:^```{lmax}` (lmax-m+1)) – Layout: [m=0: l=0..lmax | m=1: l=1..lmax | …]

alm2map_latlon(alm, ring_theta, ring_phi_list, ring_counts, lmax=None)[source]#

Synthesis: reconstruct the map from alm coefficients.

Parameters:

  • alm ([..., n_alm]  harmonic coefficients)

  • ring_theta ([R] colatitudes (radians))

  • ring_phi_list (list[ndarray] or ndarray [N_total] longitudes)

  • ring_counts ([R] number of pixels per ring)

  • lmax (maximum multipole used when computing the alm)

Returns:

im_out ([..., N_total]  reconstructed map)

anafast_latlon(im, ring_theta, ring_phi_list, ring_counts, lmax=None, weights=None, quadrature='trapeze')[source]#

Estimate the power spectrum Cl of a map on an arbitrary grid.

Returns:

cl (tensor [..., lmax+1])

grid_summary(ring_theta, ring_phi_list, ring_counts, lmax=None)[source]#

Print a summary of the grid and the estimated computation cost.