Comparing `healpy.harmonic_ud_grade` with `skytools.change_resolution`

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.

#| echo: false
import numpy as np
import healpy as hp
import skytools as st
import matplotlib.pyplot as plt

plt.rcParams["figure.dpi"] = 120
plt.rcParams["font.size"] = 10
np.random.seed(42)


def planck_like_cls(lmax):
    """Very rough Planck-like TT, EE, TE, BB power spectra."""
    ell = np.arange(2, lmax + 1, dtype=float)
    tt = 3000.0 * np.exp(-ell / 1500) * (1 + (ell / 500) ** (-1.5))
    ee = 0.1 * tt
    te = 0.3 * np.sqrt(tt * ee)
    bb = ee * 0.01 + np.where(ell < 100, 0.01, 0.0)
    cls = np.zeros((6, lmax + 1))
    factor = 2 * np.pi / ell / (ell + 1)
    for arr, idx in [(tt, 0), (ee, 1), (bb, 2), (te, 3)]:
        cls[idx, 2:] = arr * factor
    return cls

print(f"healpy {hp.__version__}")

import importlib.metadata
import os
import pathlib
import subprocess
import skytools

skytools_dir = pathlib.Path(skytools.__file__).resolve().parent
try:
    skytools_version = importlib.metadata.version("skytools")
except importlib.metadata.PackageNotFoundError:
    skytools_version = "unknown"

try:
    repo_root = pathlib.Path(
        subprocess.check_output(
            ["git", "-C", str(skytools_dir), "rev-parse", "--show-toplevel"],
            stderr=subprocess.DEVNULL,
            text=True,
        ).strip()
    )
    pyproject = repo_root / "pyproject.toml"
    if pyproject.exists() and "skytools" in pyproject.read_text().lower():
        git_hash = subprocess.check_output(
            ["git", "-C", str(repo_root), "rev-parse", "--short", "HEAD"],
            text=True,
        ).strip()
        git_date = subprocess.check_output(
            ["git", "-C", str(repo_root), "log", "-1", "--format=%ci"],
            text=True,
        ).strip()
        print(f"skytools {skytools_version} git: {git_hash} ({git_date})")
    else:
        print(f"skytools {skytools_version} ({skytools_dir})")
except Exception:
    print(f"skytools {skytools_version} ({skytools_dir})")

healpy 1.19.1.dev89+g1b8b872ed.d20260529
skytools 0.0.1b7 (/home/zonca/p/software/healpy/.venv/lib/python3.12/site-packages/skytools)

1. Temperature-only downgrade (NSIDE 256 → 64)

Called with matching settings. The expected residual is numerical noise from different weighting strategies (healpy uses pixel-weights + optional iteration, skytools uses ring-weights).

NSIDE_IN = 256
NSIDE_OUT = 64
LMAX_IN = 3 * NSIDE_IN - 1
LMAX_OUT = 3 * NSIDE_OUT - 1

cls_in = planck_like_cls(LMAX_IN)
alm = hp.synalm(cls_in[0], lmax=LMAX_IN)
map_in = hp.alm2map(alm, nside=NSIDE_IN)

map_hp = hp.harmonic_ud_grade(
    map_in, nside_out=NSIDE_OUT, pol=False, pixwin=False, fwhm_out=0
)
map_st = st.change_resolution(
    map_in, nside_out=NSIDE_OUT, mode="i", lmax_sht=LMAX_OUT
)

res = map_hp - map_st

fig = plt.figure(figsize=(16, 3.5))
for i, (m, title) in enumerate([
    (map_in, f"Input T  (NSIDE={NSIDE_IN})"),
    (map_hp, "healpy harmonic_ud_grade"),
    (map_st, "skytools change_resolution"),
    (res, f"Residual  (σ={res.std():.2e})"),
]):
    hp.projview(m, title=title, min=-300 if i < 3 else -5,
                max=300 if i < 3 else 5, cmap="RdBu_r",
                sub=(1, 4, i + 1))
plt.show()

print(f"T downgrade residual: mean={res.mean():.2e} std={res.std():.2e} max|diff|={np.abs(res).max():.2e}")

T downgrade residual: mean=-1.47e-05 std=3.49e-03 max|diff|=2.56e-01

2. Polarization downgrade (IQU, NSIDE 128 → 32)

Three-component polarized maps test spin-2 transform agreement between the two libraries.

NSIDE_IN = 128
NSIDE_OUT = 32
LMAX_IN = 3 * NSIDE_IN - 1
LMAX_OUT = 3 * NSIDE_OUT - 1

cls_in = planck_like_cls(LMAX_IN)
alm_t, alm_e, alm_b = hp.synalm(cls_in[:4], lmax=LMAX_IN, new=True)
maps_in = hp.alm2map([alm_t, alm_e, alm_b], nside=NSIDE_IN)

maps_hp = hp.harmonic_ud_grade(
    maps_in, nside_out=NSIDE_OUT, pol=True, pixwin=False, fwhm_out=0
)
maps_st = st.change_resolution(
    maps_in, nside_out=NSIDE_OUT, mode="iqu", lmax_sht=LMAX_OUT
)

res_t = maps_hp[0] - maps_st[0]
res_q = maps_hp[1] - maps_st[1]
res_u = maps_hp[2] - maps_st[2]

fig = plt.figure(figsize=(16, 10))
for i, (comp, m_hp, m_st, r) in enumerate([
    ("T", maps_hp[0], maps_st[0], res_t),
    ("Q", maps_hp[1], maps_st[1], res_q),
    ("U", maps_hp[2], maps_st[2], res_u),
]):
    for j, (m, title, vmin, vmax) in enumerate([
        (m_hp, f"healpy {comp}", -200, 200),
        (m_st, f"skytools {comp}", -200, 200),
        (r, f"residual {comp}  (σ={r.std():.2e})", -2, 2),
    ]):
        idx = i * 3 + j + 1
        hp.projview(m, title=title, min=vmin, max=vmax,
                    cmap="RdBu_r", sub=(3, 3, idx))
plt.show()

print(f"\nPolarization residuals:")
print(f"  T std: {res_t.std():.2e}")
print(f"  Q std: {res_q.std():.2e}")
print(f"  U std: {res_u.std():.2e}")


Polarization residuals:
  T std: 1.18e-02
  Q std: 2.71e-03
  U std: 2.65e-03

3. Power-spectrum agreement after downgrade

After downgrading, the observed power spectrum changes. We compare the spectra of the two outputs to quantify statistical agreement.

cl_hp = hp.anafast(maps_hp, lmax=LMAX_OUT)
cl_st = hp.anafast(maps_st, lmax=LMAX_OUT)

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for i, (label, idx) in enumerate([("TT", 0), ("EE", 1), ("TE", 2)]):
    ax = axes[i]
    ells = np.arange(2, LMAX_OUT + 1)
    dls_hp = cl_hp[idx, 2:] * ells * (ells + 1) / (2 * np.pi)
    dls_st = cl_st[idx, 2:] * ells * (ells + 1) / (2 * np.pi)

    ax.plot(ells, dls_hp, label="healpy", lw=1.5)
    ax.plot(ells, dls_st, label="skytools", lw=1.0, ls="--")
    ax.set_xlabel(r"$\ell$")
    ax.set_ylabel(r"$D_\ell^{" + label + "}$ [μK²]")
    ax.set_title(label)
    ax.legend()
    ax.set_yscale("log")
plt.show()

fig, ax = plt.subplots(figsize=(8, 4))
for i, (label, idx) in enumerate([("TT", 0), ("EE", 1), ("TE", 2)]):
    with np.errstate(divide="ignore", invalid="ignore"):
        frac = np.abs(cl_hp[idx, 2:] - cl_st[idx, 2:]) / np.abs(cl_hp[idx, 2:])
    frac = frac[np.isfinite(frac)]
    print(f"  {label} median |ΔCl|/Cl: {np.median(frac):.2e}")
    ells = np.arange(2, LMAX_OUT + 1)
    frac_full = np.abs(cl_hp[idx, 2:] - cl_st[idx, 2:]) / np.clip(np.abs(cl_hp[idx, 2:]), 1e-30, np.inf)
    ax.semilogy(ells, frac_full, label=label, alpha=0.7)
ax.set_xlabel(r"$\ell$")
ax.set_ylabel(r"$|\Delta C_\ell| / C_\ell$")
ax.set_title("Fractional power-spectrum difference")
ax.legend()
ax.grid(alpha=0.3)
plt.show()

  TT median |ΔCl|/Cl: 3.84e-06
  EE median |ΔCl|/Cl: 2.15e-06
  TE median |ΔCl|/Cl: 2.32e-05

4. Beam re-convolution (0.5° → 1°)

Key API difference: healpy expects fwhm in radians, skytools in arcminutes.

We create an input map smoothed at 0.5°, then downgrade and apply a 1° output beam.

NSIDE_IN = 256
NSIDE_OUT = 128
LMAX_IN = 3 * NSIDE_IN - 1
LMAX_OUT = 3 * NSIDE_OUT - 1

fwhm_half_deg = 0.5   # degrees
fwhm_one_deg = 1.0   # degrees
fwhm_half_rad = np.radians(fwhm_half_deg)
fwhm_one_rad = np.radians(fwhm_one_deg)

cls_in = planck_like_cls(LMAX_IN)
alm = hp.synalm(cls_in[0], lmax=LMAX_IN)
map_smooth = hp.alm2map(alm, nside=NSIDE_IN, fwhm=fwhm_half_rad)

# healpy: fwhm in radians
map_hp_beam = hp.harmonic_ud_grade(
    map_smooth,
    nside_out=NSIDE_OUT,
    pol=False,
    pixwin=False,
    fwhm_in=fwhm_half_rad,
    fwhm_out=fwhm_one_rad,
)

# skytools: fwhm in arcminutes
map_st_beam = st.change_resolution(
    map_smooth,
    nside_out=NSIDE_OUT,
    mode="i",
    lmax_sht=LMAX_OUT,
    fwhm_in=fwhm_half_deg * 60,  # convert deg -> arcmin
    fwhm_out=fwhm_one_deg * 60,   # convert deg -> arcmin
)

res_beam = map_hp_beam - map_st_beam

fig = plt.figure(figsize=(16, 3.5))
for i, (m, title, vmin, vmax) in enumerate([
    (map_smooth, "Input (smoothed 0.5°)", -300, 300),
    (map_hp_beam, "healpy (1° out)", -300, 300),
    (map_st_beam, "skytools (1° out)", -300, 300),
    (res_beam, f"Residual (σ={res_beam.std():.2e})", -5, 5),
]):
    hp.projview(m, title=title, min=vmin, max=vmax,
                cmap="RdBu_r", sub=(1, 4, i + 1))
plt.show()

print(f"\nBeam comparison residual: std={res_beam.std():.2e} max|diff|={np.abs(res_beam).max():.2e}")


Beam comparison residual: std=3.26e-04 max|diff|=6.18e-02

5. Pixel-window function handling

Both libraries can deconvolve the input pixel window and apply the output one.

NSIDE_IN = 256
NSIDE_OUT = 64
LMAX_IN = 3 * NSIDE_IN - 1
LMAX_OUT = 3 * NSIDE_OUT - 1

cls_in = planck_like_cls(LMAX_IN)
alm = hp.synalm(cls_in[0], lmax=LMAX_IN)
map_in = hp.alm2map(alm, nside=NSIDE_IN)

# healpy
map_hp_pw = hp.harmonic_ud_grade(
    map_in, nside_out=NSIDE_OUT, pol=False, pixwin=True, fwhm_out=0
)

# skytools: pass nside integers
map_st_pw = st.change_resolution(
    map_in,
    nside_out=NSIDE_OUT,
    mode="i",
    lmax_sht=LMAX_OUT,
    pixwin_in=NSIDE_IN,
    pixwin_out=NSIDE_OUT,
)

res_pw = map_hp_pw - map_st_pw

fig = plt.figure(figsize=(16, 3.5))
for i, (m, title, vmin, vmax) in enumerate([
    (map_hp_pw, "healpy (with pixwin)", -300, 300),
    (map_st_pw, "skytools (with pixwin)", -300, 300),
    (res_pw, f"Residual (σ={res_pw.std():.2e})", -5, 5),
]):
    hp.projview(m, title=title, min=vmin, max=vmax,
                cmap="RdBu_r", sub=(1, 3, i + 1))
plt.show()

print(f"\nPixel-window residual: std={res_pw.std():.2e} max|diff|={np.abs(res_pw).max():.2e}")


Pixel-window residual: std=3.13e-03 max|diff|=2.67e-01

6. Custom beam transfer functions

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.

NSIDE_IN = 256
NSIDE_OUT = 64
LMAX_IN = 3 * NSIDE_IN - 1
LMAX_OUT = 3 * NSIDE_OUT - 1

# Cosine-tapered beam: smooth rolloff, non-Gaussian
ell = np.arange(LMAX_OUT + 1, dtype=float)
ell_knee = LMAX_OUT // 3
cosine_beam = np.where(ell <= ell_knee, 1.0,
              0.5 * (1 + np.cos(np.pi * (ell - ell_knee) / (LMAX_OUT - ell_knee))))
cosine_beam[LMAX_OUT:] = 0.0

cls_in = planck_like_cls(LMAX_IN)
alm = hp.synalm(cls_in[0], lmax=LMAX_IN)
map_in = hp.alm2map(alm, nside=NSIDE_IN)

# healpy: beam_window_out (1D for T-only)
map_hp_cos = hp.harmonic_ud_grade(
    map_in, nside_out=NSIDE_OUT, pol=False, pixwin=False,
    beam_window_out=cosine_beam,
)

# skytools: beam_out (1D for mode="i")
map_st_cos = st.change_resolution(
    map_in, nside_out=NSIDE_OUT, mode="i", lmax_sht=LMAX_OUT,
    beam_out=cosine_beam,
)

res_cos = map_hp_cos - map_st_cos

# Compare with a Gaussian beam of similar width
fwhm_gauss = np.radians(30.0 / 60.0)  # 30 arcmin in radians
gauss_beam_arr = hp.gauss_beam(fwhm_gauss, lmax=LMAX_OUT)

map_hp_gauss = hp.harmonic_ud_grade(
    map_in, nside_out=NSIDE_OUT, pol=False, pixwin=False,
    beam_window_out=gauss_beam_arr,
)

# Compare power spectra: cosine vs Gaussian beam output
cl_cos = hp.anafast(map_hp_cos, lmax=LMAX_OUT)
cl_gauss = hp.anafast(map_hp_gauss, lmax=LMAX_OUT)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
ell_plot = np.arange(LMAX_OUT + 1)
axes[0].plot(ell_plot, cosine_beam, label="Cosine-tapered", lw=2)
axes[0].plot(ell_plot, gauss_beam_arr, label="Gaussian (30′)", ls="--", lw=2)
axes[0].set_xlabel(r"$\ell$")
axes[0].set_ylabel(r"$b_\ell$")
axes[0].set_title("Beam transfer functions")
axes[0].legend()
axes[0].grid(alpha=0.3)

ells = np.arange(2, LMAX_OUT + 1)
axes[1].semilogy(ells, cl_cos[2:] * ells * (ells + 1) / (2*np.pi), label="Cosine beam", lw=2)
axes[1].semilogy(ells, cl_gauss[2:] * ells * (ells + 1) / (2*np.pi), label="Gaussian beam", ls="--", lw=1.5)
axes[1].set_xlabel(r"$\ell$")
axes[1].set_ylabel(r"$D_\ell^{TT}$")
axes[1].set_title("Output power spectra")
axes[1].legend()
axes[1].grid(alpha=0.3)
plt.show()

print(f"\nCustom beam residual (healpy vs skytools):")
print(f"  std={res_cos.std():.2e}  max|diff|={np.abs(res_cos).max():.2e}")


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.

NSIDE = 128
LMAX = 3 * NSIDE - 1

fwhm_10_arcmin = 10.0
fwhm_30_arcmin = 30.0
fwhm_10_rad = np.radians(fwhm_10_arcmin / 60.0)
fwhm_30_rad = np.radians(fwhm_30_arcmin / 60.0)

cls_in = planck_like_cls(LMAX)
alm = hp.synalm(cls_in[0], lmax=LMAX)
map_smooth10 = hp.alm2map(alm, nside=NSIDE, fwhm=fwhm_10_rad)

# healpy: use_pixel_weights=False avoids deconvolution artifacts
map_hp_reconv = hp.harmonic_ud_grade(
    map_smooth10,
    nside_out=NSIDE,
    pol=False,
    pixwin=False,
    fwhm_in=fwhm_10_rad,
    fwhm_out=fwhm_30_rad,
    use_pixel_weights=False,
)

# skytools: same (nside_out defaults to input NSIDE)
map_st_reconv = st.change_resolution(
    map_smooth10,
    nside_out=NSIDE,
    mode="i",
    lmax_sht=LMAX,
    fwhm_in=fwhm_10_arcmin,
    fwhm_out=fwhm_30_arcmin,
)

# Reference: directly smooth at 30' from the original alm
map_ref_30 = hp.alm2map(alm, nside=NSIDE, fwhm=fwhm_30_rad)

res_reconv = map_hp_reconv - map_st_reconv
res_vs_ref = map_hp_reconv - map_ref_30

fig = plt.figure(figsize=(16, 3.5))
for i, (m, title, vmin, vmax) in enumerate([
    (map_smooth10, "Input (10' beam)", -300, 300),
    (map_hp_reconv, "healpy reconv (30')", -300, 300),
    (map_st_reconv, "skytools reconv (30')", -300, 300),
    (res_reconv, f"healpy - skytools (std={res_reconv.std():.2e})", -0.1, 0.1),
]):
    hp.projview(m, title=title, min=vmin, max=vmax,
                cmap="RdBu_r", sub=(1, 4, i + 1))
plt.show()

# Check spectra
cl_input = hp.anafast(map_smooth10, lmax=LMAX)
cl_reconv = hp.anafast(map_hp_reconv, lmax=LMAX)
cl_ref = hp.anafast(map_ref_30, lmax=LMAX)

fig, ax = plt.subplots(figsize=(8, 4))
ells = np.arange(2, LMAX + 1)
ax.semilogy(ells, cl_input[2:] * ells * (ells + 1) / (2*np.pi), label="Input (10'')", alpha=0.5)
ax.semilogy(ells, cl_reconv[2:] * ells * (ells + 1) / (2*np.pi), label="Reconvolved (30'')", lw=2)
ax.semilogy(ells, cl_ref[2:] * ells * (ells + 1) / (2*np.pi), label="Direct 30''", ls="--", lw=1.5)
ax.set_xlabel(r"$\ell$")
ax.set_ylabel(r"$D_\ell^{TT}$")
ax.set_title("Reconvolution: power spectrum comparison")
ax.legend()
ax.grid(alpha=0.3)
plt.show()

print(f"\nReconvolution comparison:")
print(f"  healpy vs skytools: std={res_reconv.std():.2e}")
print(f"  healpy vs direct 30': std={res_vs_ref.std():.2e}")


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.

NSIDE_IN = 256
NSIDE_OUT = 64
LMAX_OUT = 3 * NSIDE_OUT - 1

np.random.seed(42)
cls_in = planck_like_cls(3 * NSIDE_IN - 1)
alm = hp.synalm(cls_in[0], lmax=3 * NSIDE_IN - 1)
map_in = hp.alm2map(alm, nside=NSIDE_IN)

# Map input (standard path)
map_hp_map = hp.harmonic_ud_grade(
    map_in, nside_out=NSIDE_OUT, pol=False, pixwin=False, fwhm_out=0
)

# alm input (skips map2alm)
map_hp_alm = hp.harmonic_ud_grade(
    alm, nside_out=NSIDE_OUT, nside_in=NSIDE_IN,
    input_type='alm', pol=False, pixwin=False, fwhm_out=0
)

diff = map_hp_map - map_hp_alm
print(f'Map  → harmonic_ud_grade:  std={np.std(map_hp_map):.6e}')
print(f'alm → harmonic_ud_grade:  std={np.std(map_hp_alm):.6e}')
print(f'Difference (map vs alm input): std={np.std(diff):.6e}')
print(f'The small difference comes from the map2alm round-trip in the map path.')

fig = plt.figure(figsize=(16, 3.5))
for i, (m, title, vmin, vmax) in enumerate([
    (map_hp_map, 'Map input path', None, None),
    (map_hp_alm, 'alm input path', None, None),
    (diff, 'Difference', -1e-4, 1e-4),
]):
    hp.projview(m, sub=(1, 3, i + 1), title=title,
                min=vmin, max=vmax)
plt.suptitle('input_type="map" vs "alm"', y=1.02)
plt.show()

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 arcminute fwhm convention, healpy.harmonic_ud_grade and skytools.change_resolution agree closely for the downgrade cases. Residuals are dominated by the different SHT weighting strategies (pixel weights with optional iteration in healpy vs ring weights in skytools) and by the exact transform path used in each package:

Test	Latest residual
Temperature downgrade	σ ≈ 3.5×10⁻³
Polarization (T/Q/U)	σ ≈ 1.2×10⁻² / 2.7×10⁻³ / 2.7×10⁻³
Power spectrum (median frac. TT/EE/TE)	≈ 3.8×10⁻⁶ / 2.2×10⁻⁶ / 2.3×10⁻⁵
Beam re-convolution	σ ≈ 3.3×10⁻⁴
Pixel-window correction	σ ≈ 3.1×10⁻³
Custom beam (cosine-tapered)	σ ≈ 2.9×10⁻³
Reconvolution at same NSIDE	σ ≈ 5.7×10⁻²

Key API differences

Feature	healpy `harmonic_ud_grade`	skytools `change_resolution`
FWHM units	radians	arcminutes
Custom beam params	`beam_window_in`, `beam_window_out`	`beam_in`, `beam_out`
Custom beam format (T)	1D `(lmax+1,)`	1D `(lmax+1,)`
Custom beam format (pol)	2D `(lmax+1, 4)` [T, E/B, T→E, T→B]	2D `(lmax+1, nmaps)` [I, Q, U]
Pixel window	`pixwin=True/False` (auto NSIDE)	`pixwin_in/out=NSIDE` (integer)
Reconvolution	`nside_out = nside_in`	`nside_out` defaults to input NSIDE
Default output beam	`fwhm_out=0` (no extra beam); `fwhm_out=None` for Planck-scaled FWHM	None
Direct alm input	`input_type='alm'` + `nside_in`	Not supported

healpy harmonic_ud_grade vs skytools change_resolution: A Comparison

Comparing healpy.harmonic_ud_grade with skytools.change_resolution