Spherical Harmonics

The nuclei of a resolution lie on isolatitude rings with equally spaced longitudes (see Isolatitude Rings), which is the layout that ducc0’s ring-based spherical harmonic routines accept. rhealpixdggs.harmonics feeds them the ring table and handles the bookkeeping: it needs the optional harmonics extra, which brings ducc0 (wheels are published for Linux, macOS and Windows), and the rest of the package works without it.

Conventions

Coefficients follow healpy: a real field band-limited at degree \(\ell_{max}\) has one complex coefficient \(a_{\ell m}\) for each \(0 \le m \le \ell \le \ell_{max}\), stored m-major with \((\ell, m)\) at alm_index(lmax, l, m) in an array of alm_size(lmax) entries, and its value is

\[f(\theta, \lambda) = \sum_{\ell} a_{\ell 0} Y_{\ell 0}(\theta, \lambda) + 2 \operatorname{Re} \sum_{\ell} \sum_{m > 0} a_{\ell m} Y_{\ell m}(\theta, \lambda),\]

so the \(m = 0\) coefficients of a real field are real. Values are attached to nuclei on the authalic sphere, where the cells are equal-area and the rings sit at their closed-form latitudes, and they travel as a flat array in ring_order: ring by ring from the north pole, each ring in cells_on_ring order. to_ring_order and from_ring_order convert between that array and a dictionary keyed by cell index.

Synthesis and analysis

synthesis evaluates a field given by its coefficients at every nucleus; analysis finds the coefficients that best fit values at the nuclei. The nuclei are not a quadrature (their rings are not at Gauss-Legendre latitudes), so analysis is ducc0’s iterative least-squares solve, pseudo_analysis, and it returns how the fit stopped alongside the coefficients. A round trip on resolution 2:

>>> import numpy as np
>>> from rhealpixdggs.dggs import WGS84_003
>>> from rhealpixdggs import harmonics
>>> rdggs, resolution, lmax = WGS84_003, 2, 12
>>> geometry = harmonics.ring_geometry(rdggs, resolution)
>>> sorted(geometry), int(geometry['nphi'].sum())
(['nphi', 'phi0', 'ringstart', 'theta'], 486)
>>> rng = np.random.default_rng(0)
>>> size = harmonics.alm_size(lmax)
>>> alm = rng.normal(size=size) + 1j * rng.normal(size=size)
>>> alm[:lmax + 1] = alm[:lmax + 1].real
>>> values = harmonics.synthesis(rdggs, resolution, alm, lmax)
>>> by_cell = harmonics.from_ring_order(rdggs, resolution, values)
>>> len(by_cell), harmonics.ring_order(rdggs, resolution)[0]
(486, np.str_('N44'))
>>> recovered, info = harmonics.analysis(rdggs, resolution, by_cell, lmax)
>>> info.stop_reason, info.iterations < 20
(1, True)
>>> bool(np.abs(recovered - alm).max() < 1e-8)
True

ring_geometry is the description ducc0 takes as keyword arguments, straight from ring_table: colatitude per ring on the authalic sphere, cells per ring, longitude of the first nucleus, and each ring’s offset in the flat array. It is exposed so that other ducc0 routines, such as adjoint_synthesis or the spin-weighted transforms, can be called on the same layout.

Analysing a field sampled at the nuclei works the same way. The field \(\sin\theta \cos\lambda\) of the Isolatitude Rings page is a single spherical harmonic, and the fit finds exactly that:

>>> ids = harmonics.ring_order(rdggs, resolution)
>>> theta = np.repeat(geometry['theta'], geometry['nphi'].astype(int))
>>> lon = np.radians(rdggs.nuclei(ids)[:, 0])
>>> alm, info = harmonics.analysis(rdggs, resolution, np.sin(theta) * np.cos(lon), lmax=4)
>>> nonzero = np.abs(alm) > 1e-8
>>> [(l, m) for m in range(5) for l in range(m, 5) if nonzero[harmonics.alm_index(4, l, m)]]
[(1, 1)]

How high a degree a resolution supports

An equatorial ring of \(4n\) cells (\(n = N_{side}^k\)) resolves longitudinal wavenumbers up to \(2n\), and the polar rings far fewer, so a band-limited field is only recoverable below degree \(2n\), and the least-squares problem grows ill-conditioned as the degree approaches that limit. Measured on WGS84 with \(N_{side} = 3\), synthesising random coefficients and analysing them back with the default 100 iterations:

Resolution	Cells	\(n\)	Degree	Degree / \(n\)	Iterations and relative error
2	486	9	9 / 13 / 17 / 18	1.0 / 1.4 / 1.9 / 2.0	6, 7, 20 iterations at \(10^{-10}\) or better; at degree 18 the solver reports convergence but the coefficients are wrong (error 0.8)
3	4374	27	27 / 40 / 48 / 53 / 54	1.0 / 1.5 / 1.8 / 2.0 / 2.0	5, 6, 12 iterations at \(10^{-10}\); 77 iterations at \(10^{-9}\); at degree 54 the iteration limit is hit and the coefficients are wrong (error 0.3)
4	39366	81	81 / 121 / 145 / 153 / 157	1.0 / 1.5 / 1.8 / 1.9 / 1.9	5, 5, 12 iterations at \(10^{-9}\) or better; 93 iterations at \(10^{-8}\); degree 157 needs 666 iterations for \(10^{-7}\)

At degree \(2n\) several different sets of coefficients produce the same values at the nuclei, so the solver returns one of them, and it is not the one the field was made from: the values are reproduced, the residual is tiny, and the stop reason can still say converged. The output does not warn you, so the degree has to be chosen in advance. A safe rule is degrees up to about \(1.8 n\), where a dozen iterations suffice; between there and \(2n\), raise maxiter and watch info.stop_reason (1 means converged, 7 means the iteration limit was hit) and info.residual, knowing that a clean stop is necessary but not sufficient. These numbers are for exact band-limited input; a real field sampled at the nuclei has power above any chosen degree, and that power aliases into the fit, so choose the degree from the field’s smoothness, not from this table alone.

healpy’s map2alm and alm2map cannot be used instead of this module: they assume HEALPix’s ring populations, which differ from rHEALPix’s.

Module reference

Spherical harmonic transforms on the nuclei of an rHEALPix grid.

The nuclei of a resolution lie on isolatitude rings with equally spaced longitudes (see Isolatitude Rings), which is the layout ducc0’s ring-based spherical harmonic routines accept. This module builds ducc0’s description of the rings from RHEALPixDGGS.ring_table, lays cell values out in cells_on_ring order, and wraps ducc0.sht.synthesis (coefficients to nucleus values) and ducc0.sht.pseudo_analysis (nucleus values to coefficients, a least-squares fit, since the nuclei are not a quadrature).

Coefficients follow the healpy convention: a real field has complex coefficients \(a_{\ell m}\) for \(0 \le m \le \ell \le \ell_{max}\), stored m-major with \((\ell, m)\) at alm_index(lmax, l, m), and its value is \(\sum_\ell a_{\ell 0} Y_{\ell 0} + 2 \operatorname{Re} \sum_{m > 0} a_{\ell m} Y_{\ell m}\).

The transforms work on the authalic sphere, where the cells are equal-area and the rings sit at the closed-form latitudes. With n cells per base-cell side, a field band-limited up to about degree 1.8 n is recovered by analysis in a dozen iterations; towards degree 2 n the fit is ill-conditioned and slow, and at 2 n the short polar rings alias: the solver returns coefficients that reproduce the values but are not the ones the field was made from, and may still report convergence.

Requires ducc0: pip install rhealpixdggs[harmonics].

class rhealpixdggs.harmonics.AnalysisInfo(stop_reason: int, iterations: int, residual: float, quality: float)

Bases: NamedTuple

How analysis stopped: stop_reason as ducc0 reports it (1: the solution converged to the requested tolerance, 2: a least-squares solution was found, 3: the problem is too ill-conditioned, 7: maxiter was reached), the iterations taken, the residual norm relative to the norm of the values, and ducc0’s quality measure of the least-squares solution.

iterations: int

Alias for field number 1

quality: float

Alias for field number 3

residual: float

Alias for field number 2

stop_reason: int

Alias for field number 0

rhealpixdggs.harmonics.alm_index(lmax: int, l: int, m: int) → int

Return the position of the coefficient of degree l and order m in the m-major layout healpy and ducc0 use for a field band-limited at degree lmax.

EXAMPLES:

>>> [alm_index(2, l, m) for m in range(3) for l in range(m, 3)]
[0, 1, 2, 3, 4, 5]

rhealpixdggs.harmonics.alm_size(lmax: int) → int

Return the number of coefficients of a real field band-limited at degree lmax: one for each 0 <= m <= l <= lmax.

EXAMPLES:

>>> alm_size(2)
6

rhealpixdggs.harmonics.analysis(rdggs: RHEALPixDGGS, resolution: int, values: Mapping[str, float] | NDArray[float64], lmax: int, maxiter: int = 100, epsilon: float = 1e-10, nthreads: int = 1) → tuple[ndarray, AnalysisInfo]

Return the coefficients (healpy layout, alm_size(lmax) complex numbers) of the real field band-limited at degree lmax that best fits values at the resolution resolution nuclei, and how the fit stopped. values is either a flat array in ring_order or a mapping from index string to value, which goes through to_ring_order.

The nuclei are not a quadrature, so this is ducc0.sht.pseudo_analysis, an iterative least-squares solve with at most maxiter iterations to relative tolerance epsilon. With n = N_side ** resolution, a field band-limited up to about degree 1.8 n is recovered to rounding in a dozen iterations; between that and 2 n the fit is ill-conditioned and may need maxiter raised; at 2 n and above the short polar rings alias and the result is not the coefficients you started from. Check info.stop_reason and info.residual rather than trusting the result.

EXAMPLES:

>>> from rhealpixdggs.dggs import WGS84_003
>>> import numpy as np
>>> alm = np.zeros(alm_size(2), dtype=complex)
>>> alm[alm_index(2, 2, 1)] = 0.5 - 0.25j
>>> values = synthesis(WGS84_003, 1, alm, lmax=2)
>>> recovered, info = analysis(WGS84_003, 1, values, lmax=2)
>>> bool(np.abs(recovered - alm).max() < 1e-9), info.stop_reason
(True, 1)

rhealpixdggs.harmonics.from_ring_order(rdggs: RHEALPixDGGS, resolution: int, values: NDArray[float64]) → dict[str, float]

Return the flat array values, laid out as ring_order describes, as a dictionary keyed by index string: the inverse of to_ring_order.

EXAMPLES:

>>> from rhealpixdggs.dggs import WGS84_003
>>> import numpy as np
>>> back = from_ring_order(WGS84_003, 1, np.arange(54.0))
>>> back['N4'], back['O0']
(0.0, 9.0)

rhealpixdggs.harmonics.ring_geometry(rdggs: RHEALPixDGGS, resolution: int) → dict[str, ndarray]

Return the description of the resolution resolution nuclei that ducc0’s ring-based routines take as keyword arguments: theta, the colatitude of each ring on the authalic sphere in radians; nphi, the cells on it; phi0, the longitude in radians of its first nucleus in cells_on_ring order (0 on a single-cell pole ring); and ringstart, its offset in the flat array to_ring_order builds.

EXAMPLES:

>>> from rhealpixdggs.dggs import WGS84_003
>>> geometry = ring_geometry(WGS84_003, 1)
>>> geometry['nphi'].tolist()
[1, 8, 12, 12, 12, 8, 1]
>>> geometry['ringstart'].tolist()
[0, 1, 9, 21, 33, 45, 53]

rhealpixdggs.harmonics.ring_order(rdggs: RHEALPixDGGS, resolution: int) → ndarray

Return the index strings of every resolution resolution cell in the order of the flat array the transforms use: ring by ring from the north pole, each ring in cells_on_ring order.

EXAMPLES:

>>> from rhealpixdggs.dggs import WGS84_003
>>> ring_order(WGS84_003, 1)[:10].tolist()
['N4', 'N6', 'N7', 'N8', 'N5', 'N2', 'N1', 'N0', 'N3', 'O0']

rhealpixdggs.harmonics.synthesis(rdggs: RHEALPixDGGS, resolution: int, alm: Any, lmax: int, nthreads: int = 1) → NDArray[float64]

Return the value at every resolution resolution nucleus, in ring_order, of the real field with coefficients alm (healpy layout, alm_size(lmax) complex numbers), evaluated on the authalic sphere by ducc0.sht.synthesis.

EXAMPLES:

>>> from rhealpixdggs.dggs import WGS84_003
>>> import numpy as np
>>> alm = np.zeros(alm_size(1), dtype=complex)
>>> alm[alm_index(1, 1, 0)] = 1.0  # Y_10, proportional to cos(colatitude)
>>> values = synthesis(WGS84_003, 1, alm, lmax=1)
>>> values.round(4)[[0, 9, 21, 33]].tolist()
[0.4886, 0.2172, 0.0, -0.2172]

rhealpixdggs.harmonics.to_ring_order(rdggs: RHEALPixDGGS, resolution: int, values: Mapping[str, float]) → NDArray[float64]

Return the values of every resolution resolution cell, keyed by index string in values, as the flat float64 array the transforms use (see ring_order). Every cell of the resolution must be present: raise a ValueError naming the first missing one otherwise.

EXAMPLES:

>>> from rhealpixdggs.dggs import WGS84_003
>>> values = {str(c): c.nucleus(plane=False)[1] for c in WGS84_003.grid(1)}
>>> to_ring_order(WGS84_003, 1, values).round(3)[:9].tolist()
[90.0, 58.528, 58.528, 58.528, 58.528, 58.528, 58.528, 58.528, 58.528]