alm_latlon — Part 1: Quickstart

Module healpix_analyse.alm_latlon
Dependencies numpy, torch only — no HEALPix, no FOSCAT
Parts 1 · Quickstart · 2 · Mathematics · 3 · API reference

---

Installation

alm_latlon ships as a single file inside the healpix_analyse package. No extra dependency beyond numpy and torch is required.

from healpix_analyse.alm_latlon import (
    build_rings_from_latlon,
    compute_weights,
    map2alm_latlon,
    alm2map_latlon,
    anafast_latlon,
    grid_summary,
)

---

The four-line workflow

Every use of alm_latlon follows the same pattern regardless of the grid type:

1.  build_rings_from_latlon   →  ring_theta, ring_phi_list, ring_counts, sort_idx
2.  im_sorted = im[sort_idx]  →  reorder your map into ring order
3.  map2alm_latlon  /  anafast_latlon  →  alm  /  Cl
4.  alm2map_latlon  (optional)  →  reconstructed map

---

Example 1 — HEALPix pixels (validation against healpy)

Note

healpy appears here, and only here, because this example exists to cross-check our spectrum against a reference implementation. It is not a dependency of healpix-analyse: every other example takes its geometry from healpix-geo.

The simplest way to check correctness is to use HEALPix pixel coordinates.

import numpy as np
import healpy as hp
from healpix_analyse.alm_latlon import build_rings_from_latlon, anafast_latlon

nside = 64
npix  = 12 * nside**2
lmax  = 3 * nside   # same default as healpy

# Simulate a random map
np.random.seed(0)
im = np.random.randn(npix)

# Step 1 — build rings from HEALPix pixel coordinates (ring order)
theta, phi = hp.pix2ang(nside, np.arange(npix))   # colatitude, longitude in radians
ring_theta, ring_phi_list, ring_counts, sort_idx = build_rings_from_latlon(
    theta, phi,
    convention="colatitude_rad"   # default, no conversion needed
)

# Step 2 — compute Cl
# HEALPix is an equal-area grid → use quadrature="equal_area"
cl = anafast_latlon(
    im[sort_idx],
    ring_theta, ring_phi_list, ring_counts,
    lmax=lmax,
    quadrature="equal_area",
)

# Reference
cl_hp = hp.anafast(im, lmax=lmax)

# Both Cl arrays should agree to within numerical precision
import matplotlib.pyplot as plt
ell = np.arange(lmax + 1)
plt.loglog(ell[2:], cl.numpy()[2:],  label="anafast_latlon")
plt.loglog(ell[2:], cl_hp[2:],       label="healpy.anafast", ls="--")
plt.xlabel("ℓ"); plt.ylabel("Cℓ"); plt.legend(); plt.show()

---

Example 2 — ERA5 / ECMWF Gaussian grid

ERA5 data comes on a Gaussian grid where latitudes are the zeros of a Legendre polynomial. The "gauss_legendre" quadrature rule must be used — any other rule introduces oscillations at high ℓ (see Part 2 for the reason).

import numpy as np
import xarray as xr
from healpix_analyse.alm_latlon import build_rings_from_latlon, anafast_latlon

# Load an ERA5 field (e.g. surface temperature, single level)
ds  = xr.open_dataset("era5_t2m.nc")
t2m = ds["t2m"].values.ravel()          # shape [N]
lat = ds["latitude"].values             # 1-D, geographic degrees [-90, +90]
lon = ds["longitude"].values            # 1-D, degrees [0, 360]

# Broadcast to pixel arrays if the field is on a 2-D (lat × lon) grid
lat2d, lon2d = np.meshgrid(lat, lon, indexing="ij")
lat_flat = lat2d.ravel()
lon_flat = lon2d.ravel()

# Step 1 — build rings
ring_theta, ring_phi_list, ring_counts, sort_idx = build_rings_from_latlon(
    lat_flat, lon_flat,
    convention="geographic_deg"   # ERA5 uses geographic degrees
)

# Step 2 — compute Cl with Gauss-Legendre quadrature
lmax = 2 * len(ring_theta)   # safe upper bound for a Gaussian grid
cl = anafast_latlon(
    t2m[sort_idx],
    ring_theta, ring_phi_list, ring_counts,
    lmax=lmax,
    quadrature="gauss_legendre",   # mandatory for Gaussian grids
)

ell = np.arange(lmax + 1)
import matplotlib.pyplot as plt
plt.loglog(ell[2:], cl.numpy()[2:])
plt.xlabel("ℓ"); plt.ylabel("Cℓ  [K²]"); plt.title("ERA5 T2m power spectrum")
plt.show()

---

Example 3 — Regular lat/lon grid

A grid with equally-spaced latitudes and longitudes (e.g. from a global atmospheric model output on a regular 0.5° grid).

import numpy as np
from healpix_analyse.alm_latlon import build_rings_from_latlon, anafast_latlon

# Build a regular 1° global grid
dlat, dlon = 1.0, 1.0
lat_1d = np.arange(-90 + dlat/2,  90, dlat)   # geographic degrees, N→S
lon_1d = np.arange(0,            360, dlon)   # degrees, 0→360

lat2d, lon2d = np.meshgrid(lat_1d, lon_1d, indexing="ij")
lat_flat = lat2d.ravel()
lon_flat = lon2d.ravel()

np.random.seed(1)
im = np.random.randn(len(lat_flat))

ring_theta, ring_phi_list, ring_counts, sort_idx = build_rings_from_latlon(
    lat_flat, lon_flat,
    convention="geographic_deg"
)

# Trapezoidal rule is fine for a regular grid
lmax = 180
cl = anafast_latlon(
    im[sort_idx],
    ring_theta, ring_phi_list, ring_counts,
    lmax=lmax,
    quadrature="trapeze",
)
print(cl.shape)   # torch.Size([181])

---

Example 4 — Getting the alm coefficients

When you need the individual alm coefficients (e.g. for filtering, cross-spectra, or further processing):

from healpix_analyse.alm_latlon import map2alm_latlon
import numpy as np
from healpix_geo import nested

depth = 5                                  # nside = 32
npix  = 12 * 4 ** depth
im    = np.random.randn(npix)
lmax  = 3 * 2 ** depth

lon, lat = nested.healpix_to_lonlat(np.arange(npix, dtype=np.uint64), depth)
theta, phi = np.radians(90.0 - lat), np.radians(lon)
ring_theta, ring_phi_list, ring_counts, sort_idx = build_rings_from_latlon(
    theta, phi, convention="colatitude_rad"
)

alm = map2alm_latlon(
    im[sort_idx],
    ring_theta, ring_phi_list, ring_counts,
    lmax=lmax,
    quadrature="equal_area",
)

# alm is a complex128 torch.Tensor of shape [n_alm]
# n_alm = sum(lmax - m + 1 for m in range(lmax+1))
print(alm.shape)        # torch.Size([n_alm])
print(alm.dtype)        # torch.complex128
print(alm[:3])          # a_00, a_10, a_20  (m=0 block)

---

Example 5 — Round-trip analysis / synthesis

alm2map_latlon reconstructs the map from its alm coefficients. The residual measures the truncation error introduced by the finite lmax.

from healpix_analyse.alm_latlon import map2alm_latlon, alm2map_latlon
import numpy as np
from healpix_geo import nested

depth = 5                                  # nside = 32
npix  = 12 * 4 ** depth
lmax  = 3 * 2 ** depth
im    = np.random.randn(npix)

lon, lat = nested.healpix_to_lonlat(np.arange(npix, dtype=np.uint64), depth)
theta, phi = np.radians(90.0 - lat), np.radians(lon)
ring_theta, ring_phi_list, ring_counts, sort_idx = build_rings_from_latlon(
    theta, phi, convention="colatitude_rad"
)

im_sorted = im[sort_idx]

# Analysis
alm = map2alm_latlon(
    im_sorted, ring_theta, ring_phi_list, ring_counts,
    lmax=lmax, quadrature="equal_area"
)

# Synthesis (same grid)
im_rec = alm2map_latlon(
    alm, ring_theta, ring_phi_list, ring_counts, lmax=lmax
)

residual = (im_rec.numpy() - im_sorted)
print(f"Max |residual| : {np.abs(residual).max():.2e}")
print(f"RMS  residual  : {np.std(residual):.2e}")

---

Example 6 — Batch maps

All transform functions support a leading batch dimension, so you can process many maps in a single call:

import numpy as np
from healpix_geo import nested
from healpix_analyse.alm_latlon import build_rings_from_latlon, anafast_latlon

depth = 6                                  # nside = 64
B     = 50                                 # number of maps
lmax  = 3 * 2 ** depth
npix  = 12 * 4 ** depth

# Batch of B independent maps, shape [B, npix]
im_batch = np.random.randn(B, npix)

lon, lat = nested.healpix_to_lonlat(np.arange(npix, dtype=np.uint64), depth)
theta, phi = np.radians(90.0 - lat), np.radians(lon)
ring_theta, ring_phi_list, ring_counts, sort_idx = build_rings_from_latlon(
    theta, phi, convention="colatitude_rad"
)

# Apply sort_idx along the pixel axis
cl_batch = anafast_latlon(
    im_batch[:, sort_idx],
    ring_theta, ring_phi_list, ring_counts,
    lmax=lmax,
    quadrature="equal_area",
)
print(cl_batch.shape)   # torch.Size([50, lmax+1])

# Mean spectrum across realisations
cl_mean = cl_batch.mean(dim=0)

---

Example 7 — Inspecting the grid

grid_summary prints a diagnostic before committing to a long computation:

from healpix_analyse.alm_latlon import build_rings_from_latlon, grid_summary
import numpy as np
from healpix_geo import nested

depth = 7                                  # nside = 128
npix  = 12 * 4 ** depth
lon, lat = nested.healpix_to_lonlat(np.arange(npix, dtype=np.uint64), depth)
theta, phi = np.radians(90.0 - lat), np.radians(lon)
ring_theta, ring_phi_list, ring_counts, _ = build_rings_from_latlon(
    theta, phi, convention="colatitude_rad"
)

info = grid_summary(ring_theta, ring_phi_list, ring_counts, lmax=3 * 2 ** depth)

Example output:

=== Grid summary  (lmax=384) ===
  Total pixels      : 196608
  Latitude rings    : 511
  Uniform rings     : 511/511  (FFT acceleration)
  θ range           : [0.110°, 179.890°]
  Pixels/ring       : min=4, max=512, mean=384.8
  n_alm             : 74305
  Cost estimate     : 1.50e+06 ops  (FFT: 1.50e+06, DFT: 0.00e+00)

---

Choosing the right quadrature

Grid type	quadrature	Notes
HEALPix (any nside)	"equal_area"	All pixels have the same area.
Regular lat/lon (equidistant θ)	"trapeze"	Trapezoidal rule works well.
ERA5, ECMWF, IFS, ARPEGE, GFS	"gauss_legendre"	Mandatory — eliminates Cℓ oscillations at high ℓ.
Custom explicit weights	pass weights= array	quadrature is ignored.

Symptom of wrong quadrature: regular oscillations in Cℓ that grow towards high ℓ. If you see this pattern, switch to "gauss_legendre".

---

Coordinate convention quick reference

build_rings_from_latlon accepts four input conventions via the convention argument:

convention	lat	lon
"colatitude_rad" (default)	colatitude θ, radians, 0 → π	longitude φ, radians, 0 → 2π
"colatitude_deg"	colatitude θ, degrees, 0° → 180°	longitude φ, degrees, 0° → 360°
"geographic_rad"	geographic latitude, radians, −π/2 → +π/2	longitude, radians, −π → +π or 0 → 2π
"geographic_deg"	geographic latitude, degrees, −90° → +90°	longitude, degrees, −180° → +180° or 0° → 360°

All conventions convert internally to colatitude + longitude in radians.