foscat.alm_latlon ================= .. py:module:: foscat.alm_latlon .. autoapi-nested-parse:: alm_latlon.py ============= 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 ------- .. autoapisummary:: foscat.alm_latlon.alm_latlon Module Contents --------------- .. py:class:: alm_latlon(backend=None, lmax=24, limit_range=10000000000.0) Bases: :py:obj:`foscat.alm.alm` Spherical-harmonic transform on an arbitrary lat/lon grid organised into rings. .. py:method:: build_rings_from_latlon(lat, lon, atol=1e-10, convention='colatitude_rad') :staticmethod: Group a flat list of pixels into rings of equal colatitude. :Parameters: * **lat** (:py:class:`array [N] angular coordinate` of :py:class:`each pixel (see convention).`) * **lon** (:py:class:`array [N] longitudinal coordinate` of :py:class:`each pixel (see convention).`) * **atol** (:py:class:`tolerance in radians for grouping two pixels into the same ring.`) * **convention** (:py:class:`str format` of :py:class:`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** (:py:class:`ndarray [R] colatitude θ (radians) per ring`) * **ring_phi_list** (:py:class:`list[ndarray [N_r]] longitudes φ (radians) per ring`) * **ring_counts** (:py:class:`ndarray int64 [R] number` of :py:class:`pixels per ring`) * **sort_idx** (:py:class:`ndarray int64 [N] permutation im_sorted = im[sort_idx]`) .. py:method:: compute_weights(ring_theta, ring_phi_list, ring_counts, quadrature='trapeze') :staticmethod: Computes the quadrature weights per pixel (steradians). :Parameters: * **ring_theta** (:py:class:`[R] colatitudes in radians`) * **ring_phi_list** (:py:class:`list[ndarray] longitudes per ring`) * **ring_counts** (:py:class:`[R] number` of :py:class:`pixels per ring`) * **quadrature** (:py:class:`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** (:py:class:`ndarray float64 [N_total]`) .. py:method:: comp_tf_latlon(im, ring_phi_list, ring_counts, pixel_weights, mmax) 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** (:py:class:`[...`, :py:class:`N_total] map (backend tensor` or :py:class:`ndarray)`) * **ring_phi_list** (:py:class:`list[ndarray] longitudes per ring`) * **ring_counts** (:py:class:`ndarray int64 number` of :py:class:`pixels per ring`) * **pixel_weights** (:py:class:`ndarray float64 [N_total] quadrature weights`) * **mmax** (:py:class:`int maximum frequency`) :returns: **ft** (:py:class:`tensor [...`, :py:class:`R`, :py:class:`mmax+1] complex`) .. py:method:: map2alm_latlon(im, ring_theta, ring_phi_list, ring_counts, lmax=None, weights=None, quadrature='trapeze') Compute the alm coefficients of a map defined on an arbitrary grid organised into rings. :Parameters: * **im** (:py:class:`[...`, :py:class:`N_total] map values`, :py:class:`ordered ring by ring`) -- (use sort_idx from build_rings_from_latlon if the map is initially in arbitrary order) * **ring_theta** (:py:class:`[R] colatitude (radians)` of :py:class:`each ring`) * **ring_phi_list** (:py:class:`list[ndarray]` or :py:class:`ndarray [N_total]`) -- longitudes (radians) per ring (or flat array) * **ring_counts** (:py:class:`[R] number` of :py:class:`pixels per ring`) * **lmax** (:py:class:`maximum multipole (default`: :py:class:`self.lmax)`) * **weights** (:py:class:`[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** (:py:class:`[...`, :py:class:`n_alm] complex`, n_alm = Σ_``{m=0}``:py:class:`^```{lmax}`` (lmax-m+1)) -- Layout: [m=0: l=0..lmax | m=1: l=1..lmax | …] .. py:method:: alm2map_latlon(alm, ring_theta, ring_phi_list, ring_counts, lmax=None) Synthesis: reconstruct the map from alm coefficients. :Parameters: * **alm** (:py:class:`[...`, :py:class:`n_alm] harmonic coefficients`) * **ring_theta** (:py:class:`[R] colatitudes (radians)`) * **ring_phi_list** (:py:class:`list[ndarray]` or :py:class:`ndarray [N_total] longitudes`) * **ring_counts** (:py:class:`[R] number` of :py:class:`pixels per ring`) * **lmax** (:py:class:`maximum multipole used when computing the alm`) :returns: **im_out** (:py:class:`[...`, :py:class:`N_total] reconstructed map`) .. py:method:: anafast_latlon(im, ring_theta, ring_phi_list, ring_counts, lmax=None, weights=None, quadrature='trapeze') Estimate the power spectrum Cl of a map on an arbitrary grid. :returns: **cl** (:py:class:`tensor [...`, :py:class:`lmax+1]`) .. py:method:: grid_summary(ring_theta, ring_phi_list, ring_counts, lmax=None) Print a summary of the grid and the estimated computation cost.