Flat-sky power spectra: powerspectra, powerspectra_lonlat, ps#

Three helpers compute an isotropic 1D power spectrum by binning a 2D FFT radially. They are not spherical-harmonic transforms and should not be confused with sht.anafast (HEALPixSHT — spherical harmonics on the full sky) or LocalFFT.ps (Local two-dimensional FFT on HEALPix patches), which are the production-ready spectra for, respectively, a full-sky HEALPix map and a local gnomonic patch.

from healpix_analyse.powerspectra import powerspectra
from healpix_analyse.powerspectra_lonlat import powerspectra_lonlat
from healpix_analyse.ps import ps

ps(data, data_cross=None, plot_2D_fft=False)#

A plain 2D-FFT radial power spectrum on an already-regular 2D array (ny, nx). It has no notion of HEALPix, the sphere, or physical units — frequency axes are in cycles per array index. Use it on data you have already gridded yourself, for example the output of LocalFFT.project. This function is self-contained and has no dependency on AlmTransform.

powerspectra(cell_ids, level, data, ...) and powerspectra_lonlat(lon, lat, data, ...)#

These compute a similar radial spectrum starting from HEALPix cells (powerspectra) or from arbitrary (lon, lat) samples on iso-latitude rings (powerspectra_lonlat). Both build an AlmTransform (healpix_analyse/alm.py) internally and call its .fft() method.

Current status: under development#

AlmTransform is explicitly a work in progress, and powerspectra / powerspectra_lonlat inherit its limitations:

  • AlmTransform.ifft() is not implemented. It raises NotImplementedError rather than silently returning a wrong result — the Legendre-projection stage it needs, and the ring bookkeeping it depends on, are not built yet. powerspectra/powerspectra_lonlat only ever call .fft(), so this does not affect them directly, but any code trying to reconstruct a map from their intermediate AlmTransform will hit it.

  • The second transform stage is a plain FFT along latitude, not a Legendre sum. AlmTransform.fft() does one FFT along each ring’s longitude, then a second FFT along the ring (latitude) axis. That second step is only mathematically equivalent to a spherical-harmonic projection when rings are uniformly spaced in colatitude, which is not guaranteed in general — see the module’s own TODO: NUFFT or direct Legendre analysis for the second stage, currently only FFT stage implemented. Treat the result as a flat-sky-style spectrum, not a calibrated angular power spectrum.

  • Only periodic boundary conditions are supported. fft(..., pbc=False) raises NotImplementedError; pbc=True (the default) assumes the domain wraps, which is an approximation for a regional patch.

  • No handling of the 0/360° longitude wrap in the pixel-shift correction (noted as a TODO in the source) — patches crossing the date line are not validated.

  • method="alm" is not implemented — only method="fft" (the default) actually runs; indexing_scheme only accepts "2D_array" despite the docstring mentioning "ring" as a default.

  • Debug plotting is now gated, not removed. Earlier versions of AlmTransform.fft() (and the constructor, through compute_phase_shift) called matplotlib.pyplot.show() unconditionally on every call, which opened blocking plot windows even in batch/headless use. These diagnostic plots are now only shown when the transform is constructed with debug=True; with the default debug=False (used internally by powerspectra/powerspectra_lonlat), no plot window opens.

Recommendation#

For a spectrum you can rely on today:

Use powerspectra / powerspectra_lonlat only where their specific AlmTransform-based code path is what you are testing or extending, and keep the limitations above in mind when interpreting the result.