healpix_analyse.alm_latlon ========================== .. py:module:: healpix_analyse.alm_latlon .. autoapi-nested-parse:: alm_latlon.py ============= 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 ---------- .. autoapisummary:: healpix_analyse.alm_latlon.ArrayLike Functions --------- .. autoapisummary:: healpix_analyse.alm_latlon.build_rings_from_latlon healpix_analyse.alm_latlon.compute_weights healpix_analyse.alm_latlon.comp_tf_latlon healpix_analyse.alm_latlon.map2alm_latlon healpix_analyse.alm_latlon.alm2map_latlon healpix_analyse.alm_latlon.anafast_latlon healpix_analyse.alm_latlon.grid_summary Module Contents --------------- .. py:data:: ArrayLike .. py:function:: 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] 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** (:py:class:`array-like`, :py:class:`shape [N]`) -- Latitude or colatitude of each pixel. Unit and sign depend on ``convention`` (see below). * **lon** (:py:class:`array-like`, :py:class:`shape [N]`) -- Longitude of each pixel. Unit depends on ``convention``. * **atol** (:py:class:`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** (:py:class:`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** (:py:class:`np.ndarray`, :py:class:`shape [R]`) -- Colatitude θ in **radians** of each ring, sorted in ascending order. * **ring_phi_list** (:py:class:`list` of :py:class:`np.ndarray`) -- Longitudes φ in **radians** for each ring. * **ring_counts** (:py:class:`np.ndarray`, :py:class:`shape [R]`, :py:class:`dtype int64`) -- Number of pixels N_r in each ring. * **sort_idx** (:py:class:`np.ndarray`, :py:class:`shape [N]`, :py:class:`dtype int64`) -- Permutation index. Apply as ``im_sorted = im[sort_idx]`` before passing a map to any transform function. .. py:function:: compute_weights(ring_theta: numpy.ndarray, ring_phi_list: List[numpy.ndarray], ring_counts: numpy.ndarray, quadrature: str = 'trapeze') -> numpy.ndarray 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** (:py:class:`np.ndarray`, :py:class:`shape [R]`) -- Colatitudes in radians (output of ``build_rings_from_latlon``). * **ring_phi_list** (:py:class:`list` of :py:class:`np.ndarray`) -- Longitudes in radians per ring. * **ring_counts** (:py:class:`np.ndarray`, :py:class:`shape [R]`) -- Number of pixels per ring. * **quadrature** (:py:class:`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** (:py:class:`np.ndarray`, :py:class:`shape [N_total]`, :py:class:`dtype float64`) -- Per-pixel quadrature weights in steradians, in ring order. .. py:function:: 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 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** (:py:class:`torch.Tensor`, :py:class:`shape [...`, :py:class:`N_total]`) -- Map values in ring order (apply sort_idx before calling). * **ring_phi_list** (:py:class:`list` of :py:class:`np.ndarray`) -- Longitudes in radians per ring. * **ring_counts** (:py:class:`np.ndarray`, :py:class:`shape [R]`) -- Number of pixels per ring. * **pixel_weights** (:py:class:`np.ndarray`, :py:class:`shape [N_total]`) -- Per-pixel quadrature weights (output of ``compute_weights``). * **mmax** (:py:class:`int`) -- Maximum azimuthal order to compute. * **cdtype** (:py:class:`torch.dtype`, *default* :py:class:`torch.complex128`) -- Complex dtype for the output tensor. :returns: **ft** (:py:class:`torch.Tensor`, :py:class:`shape [...`, :py:class:`R`, :py:class:`mmax+1]`, :py:class:`complex`) -- Ring-by-ring Fourier coefficients. .. py:function:: map2alm_latlon(im: Union[numpy.ndarray, torch.Tensor], ring_theta: numpy.ndarray, ring_phi_list: Union[List[numpy.ndarray], numpy.ndarray], ring_counts: numpy.ndarray, lmax: int = 64, weights: Optional[Union[numpy.ndarray, str]] = None, quadrature: str = 'trapeze', limit_range: float = 10000000000.0) -> torch.Tensor 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** (:py:class:`array-like` or :py:class:`torch.Tensor`, :py:class:`shape [...`, :py:class:`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** (:py:class:`np.ndarray`, :py:class:`shape [R]`) -- Colatitudes in radians (output of ``build_rings_from_latlon``). * **ring_phi_list** (:py:class:`list` of :py:class:`np.ndarray` or :py:class:`flat np.ndarray [N_total]`) -- Longitudes in radians per ring. * **ring_counts** (:py:class:`np.ndarray`, :py:class:`shape [R]`) -- Number of pixels per ring. * **lmax** (:py:class:`int`, *default* ``24``) -- Maximum multipole ℓ_max. * **weights** (:py:class:`np.ndarray [N_total]`, :py:obj:`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** (:py:class:`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** (:py:class:`float`, *default* ``1e10``) -- Overflow guard in the Legendre recurrence. :returns: **alm_out** (:py:class:`torch.Tensor`, :py:class:`shape [...`, :py:class:`n_alm]`, :py:class:`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]`` .. py:function:: alm2map_latlon(alm: torch.Tensor, ring_theta: numpy.ndarray, ring_phi_list: Union[List[numpy.ndarray], numpy.ndarray], ring_counts: numpy.ndarray, lmax: int = 24, limit_range: float = 10000000000.0) -> torch.Tensor 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** (:py:class:`torch.Tensor`, :py:class:`shape [...`, :py:class:`n_alm]`, :py:class:`complex`) -- Spherical harmonic coefficients in the layout produced by ``map2alm_latlon``. * **ring_theta** (:py:class:`np.ndarray`, :py:class:`shape [R]`) -- Colatitudes in radians of the target grid. * **ring_phi_list** (:py:class:`list` of :py:class:`np.ndarray` or :py:class:`flat np.ndarray [N_total]`) -- Longitudes in radians of the target grid. * **ring_counts** (:py:class:`np.ndarray`, :py:class:`shape [R]`) -- Number of pixels per ring of the target grid. * **lmax** (:py:class:`int`, *default* ``24``) -- Maximum multipole used when the alm were computed. * **limit_range** (:py:class:`float`, *default* ``1e10``) -- Overflow guard in the Legendre recurrence. :returns: **im_out** (:py:class:`torch.Tensor`, :py:class:`shape [...`, :py:class:`N_total]`, :py:class:`float64`) -- Reconstructed map values in ring order. .. py:function:: anafast_latlon(im: Union[numpy.ndarray, torch.Tensor], ring_theta: numpy.ndarray, ring_phi_list: Union[List[numpy.ndarray], numpy.ndarray], ring_counts: numpy.ndarray, lmax: int = 64, weights: Optional[Union[numpy.ndarray, str]] = None, quadrature: str = 'trapeze', limit_range: float = 10000000000.0) -> torch.Tensor 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** (:py:class:`torch.Tensor`, :py:class:`shape [...`, :py:class:`lmax+1]`, :py:class:`float64`) -- Angular power spectrum. .. py:function:: grid_summary(ring_theta: numpy.ndarray, ring_phi_list: Union[List[numpy.ndarray], numpy.ndarray], ring_counts: numpy.ndarray, lmax: int = 24) -> dict Print a human-readable summary of the grid and estimated compute cost. :Parameters: * **ring_theta** (:py:class:`np.ndarray`, :py:class:`shape [R]`) * **ring_phi_list** (:py:class:`list` of :py:class:`np.ndarray`) * **ring_counts** (:py:class:`np.ndarray`, :py:class:`shape [R]`) * **lmax** (:py:class:`int`) :returns: **info** (:py:class:`dict`) -- Dictionary with the same fields printed to stdout.