This notebook compares two SHT-based HEALPix resolution-change functions: - healpy.harmonic_ud_grade (new in healpy, now with custom beam and reconvolution support) - skytools.change_resolution (github.com/1cosmologist/skytools)
Both functions analyse a map into \(a_{\ell m}\), optionally correct for beam and pixel-window functions, and synthesise at a new NSIDE. Both support custom (non-Gaussian) beam transfer functions and reconvolution at the same NSIDE.
print(f"healpy {hp.__version__}")# skytools does not expose __version__, print git hash and dateimport skytoolsimport osfrom datetime import datetimeskytools_dir = os.path.dirname(skytools.__file__)try: git_hash =!git -C {skytools_dir} rev-parse --short HEAD git_date =!git -C {skytools_dir} log -1--format=%ciprint(f"skytools git: {git_hash[0]} ({git_date[0]})")exceptException:print("skytools (version unknown)")
Called with matching settings. The expected residual is numerical noise from different weighting strategies (healpy uses pixel-weights + optional iteration, skytools uses ring-weights).
Both libraries now support passing arbitrary beam transfer functions as arrays, not just Gaussian FWHM. This is essential for real instrumental beams that are often non-Gaussian (asymmetric, with sidelobes, etc.).
API difference in beam array format:
For temperature-only maps both use the same 1D shape (lmax+1,). For polarized maps the layouts differ:
healpy (beam_window_in / beam_window_out): matches the output of hp.gauss_beam(fwhm, lmax, pol=True) — a 2D array of shape (lmax+1, 4) with columns [T, E/B, T→E, T→B]. Column 0 is the temperature beam, column 1 is the polarization beam, and columns 2–3 are temperature–polarization cross-coupling terms. For simple unpolarized beams you can pass a 1D array or use only column 0.
skytools (beam_in / beam_out): a 2D array of shape (lmax+1, nmaps) where each column is the beam for the corresponding map component — column 0 = I, column 1 = Q, column 2 = U (no cross-coupling terms). For simple beams where T and polarization share the same beam, all columns are identical.
We use a cosine-tapered beam as a simple non-Gaussian example.
Custom beam residual (healpy vs skytools):
std=2.93e-03 max|diff|=2.73e-01
7. Reconvolution at same NSIDE
When nside_out == nside_in, both functions perform beam reconvolution: deconvolving the input beam and applying a new output beam at the same pixelisation. This is useful for changing the effective resolution of a map without changing its NSIDE — similar to the HEALPix process_alm facility.
skytools defaults nside_out to the input NSIDE when not specified, making reconvolution the default behaviour. healpy requires nside_out to be explicitly set.
Note on pixel weights: The beam deconvolution/reconvolution is applied in a single step as \(b_{\rm out}(\ell) / b_{\rm in}(\ell)\), so there is no amplification of high-\(\ell\) modes. However, healpy’s default use_pixel_weights=True produces \(a_{\ell m}\) that differ slightly from skytools’ ring-weighted SHT. At the same NSIDE these differences are not partially cancelled by a resolution change, so the residual shows a structured pixel-weight pattern. Using use_pixel_weights=False eliminates this pattern and matches skytools’ ring-weighted approach.
Reconvolution comparison:
healpy vs skytools: std=5.65e-02
healpy vs direct 30': std=4.15e+00
8. Direct \(a_{\ell m}\) input (input_type='alm')
healpy’s harmonic_ud_grade now accepts \(a_{\ell m}\) coefficients directly via input_type='alm', skipping the map2alm step. This is useful in simulation pipelines where you already have harmonics from a previous step. When using this mode, nside_in must be provided explicitly since it cannot be inferred from the \(a_{\ell m}\) array.
skytools does not have an equivalent — it always starts from a pixel-space map.
Map → harmonic_ud_grade: std=3.385528e+03
alm → harmonic_ud_grade: std=3.385528e+03
Difference (map vs alm input): std=5.988260e-03
The small difference comes from the map2alm round-trip in the map path.
Summary
After matching settings and correcting for the radian vs arcminutefwhm convention, healpy.harmonic_ud_grade and skytools.change_resolution agree to within a few ×10⁻³, with the residual dominated by the different weighting strategies (pixel weights vs ring weights):
Test
Residual σ
Temperature downgrade
~3.5×10⁻³
Polarization (T/Q/U)
~1–3×10⁻³
Power spectrum (median frac.)
~10⁻⁶
Beam re-convolution
~1.5×10⁻³
Pixel-window correction
~3.1×10⁻³
Custom beam (cosine-tapered)
~same order as Gaussian beam
Reconvolution at same NSIDE
~same order as downgrade
The dominant sources of difference are the weighting strategy (pixel weights with optional iteration in healpy vs ring weights with no iteration in skytools), which is well below the level that matters for real-world science.
