Pixell utils

The kszx.pixell_utils module contains wrappers around the pixell library.

References

kszx.pixell_utils.read_map(filename)

Reads pixell map in FITS format, and returns a pixell.enmap.

If the FITS file contains temperature + polarization, only temperature will be returned!

kszx.pixell_utils.write_map(filename, m)

Writes pixell map in FITS format.

kszx.pixell_utils.plot_map(m, downgrade, nolabels=True, filename=None, **kwds)

Plots pixell map ‘m’.

Thin wrapper around pixell.enplot(), just adding a few tweaks:

  • Make downgrade argument mandatory (usually needed in practice to avoid runaway heavyweight behavior). Suggest downgrade=20 for ACT.

  • Make nolabels=True the default (I haven’t figured out how to make pixell labels look reasonable).

  • Add ‘filename’ argument to select between show/write.

For more kwds, see:

https://pixell.readthedocs.io/en/latest/reference.html#pixell.enplot.plot

kszx.pixell_utils.map_from_catalog(shape, wcs, gcat, weights=None, rcat=None, rweights=None, normalized=True, allow_outliers=True)

Project a 3-d galaxy catalog to a 2-d pixell map.

Function args:

  • shape, wcs: these args define the pixell coordinate system.

  • gcat (Catalog): galaxy catalog. Must define columns ra_deg, dec_deg.

  • weights (1-d numpy array, optional): per-galaxy weights.

  • rcat, rweights (optional): random catalog to be subtracted. The rweights are renormalized so that sum(rweights) = -sum(weights).

  • normalized (boolean): If normalized=True, then the output map includes a factor 1 / (pixel area). This normalization best represents a sum of delta functions \(f(x) = \sum_j w_j \delta^2(\theta - \theta_j)\).

  • allow_outliers (boolean): If ‘allow_outliers’ is True, then “out-of-bounds” galaxies (i.e. galaxies outside the boundaries of the pixell map) will be silently discarded. If ‘allow_outliers’ is False, then out-of-bounds galaxies will raise an exception.

Returns a pixell.enmap.ndmap.

kszx.pixell_utils.eval_map_on_catalog(m, catalog, pad=None, return_mask=False)

Evaluates pixell map ‘m’ at (ra,dec) values from catalog, and returns a 1-d array.

Function arguments:

  • m (pixell.enmap.ndmap): 2-d map to be evaluated.

  • catalog (Catalog): galaxy catalog defining (ra, dec) values where map is evaluated. Must define columns ra_deg, dec_deg.

  • pad (boolean): Since pixell maps can be partial-sky, some of the (ra,dec) pairs may be “outliers”, i.e. outside the map footprint. If pad=None (the default), this raises an exception. If ‘pad’ is specified (e.g. pad=0.0) then outlier values are replaced by ‘pad’.

  • return_mask (boolean): Defines what quantity this function returns.

    If return_mask=False (the default), returns a 1-d array ‘map_vals’.

    If return_mask=True, returns a pair (map_vals, mask), where ‘mask’ is a 1-d boolean array which is True for non-outliers, False for outliers.

kszx.pixell_utils.alm2map(alm, shape, wcs)

Applies alm2map spherical transform, returns pixell.enmap.ndmap.

kszx.pixell_utils.ang2pix(shape, wcs, ra_deg, dec_deg, allow_outliers=False, outlier_errmsg=None)

Given (ra,dec) values on the sky, returns the corresponding pixel indices.

Note: if you are interested in this function in order to pixelize a Catalog, you may want to call map_from_catalog() instead.

Function arguments:

  • shape, wcs: these args define the pixell coordinate system.

  • ra_deg, dec_deg: numpy arrays with same shape, containing (ra, dec) values.

  • allow_outliers (boolean): since pixell maps can be partial-sky, pixel indices can be out of bounds. In this case:

    If allow_outliers=False, we raise an exception.

    If allow_outliers=True, then we “clip” the returned idec/ira arrays so that they are in bounds, and return a boolean mask to indicate which values are valid (see below).

  • outlier_errmsg (string, optional): error message if exception is thrown.

The return value depends on the value of allow_outliers:

  • If allow_outliers=False, returns a pair (idec, ira).

  • If allow_outliers=True, returns a triple (idec, ira, mask).

The (idec, ira) arrays are integer-valued, and contain pixel indices along each of the axes of the 2-d pixell map. The mask array is boolean-valued.

Roughly equivalent to healpy.ang2pix(..., latlon=True). Note that the argument ordering is (ra, dec), consistent with healpy but not pixell!

kszx.pixell_utils.uK_arcmin_from_ivar(ivar, max_uK_arcmin=1000000.0)

Given an inverse variance map, returns associated noise map (in uK-arcmin).

Based on Mat’s orphics library:

https://github.com/msyriac/orphics/blob/master/orphics/maps.py

Function args:

  • ivar (pixell map): inverse variance map, units (uK)^{-2}. (For example, the return value from read_ivar().)

  • max_uK_arcmin (float): max allowed value in the return map (prevents dividing by zero).

Return value is a pixell map with units uK-arcmin.

Equivalent to:

ps = ivar.pixsizemap()
uK_arcmin = sqrt(ps/ivar) * (60*180/np.pi)
return np.minimum(uK_arcmin, max_uK_arcmin)

but regulated to avoid dividing by zero.

kszx.pixell_utils.fkp_from_ivar(ivar, cl0, normalize=True, return_wvar=False)

Given an inverse variance map, returns associated FKP weighting.

Function args:

  • ivar (pixell map): inverse variance map, units (uK)^{-2}. (For example, the return value from read_ivar().)

  • cl0 (float): FKP weight function parameter \(C_l^{(0)}\).

    • Intuitively, cl0 = “fiducial signal \(C_l\) at wavenumber \(l\) of interest”.

    • cl0=0 corresponds to inverse noise weighting.

    • Large cl0 corresponds to uniform weighing (but cl0=np.inf won’t work).

    • I usually use cl0 = 0.01 for plotting all-sky CMB temperature maps.

    • For kSZ filtering, cl0=3e-5 is a reasonable choice.

  • normalize (boolean): if True, then we normalize the weight function so that \(\max(W(\theta))=1\).

  • return_wvar (boolean): if True, then we return \(W(\theta) / \mbox{ivar}(\theta)\), instead of returning \(W(\theta)\).

Returns a pixell map.

The FKP weighting is defined by:

\[\begin{split}\begin{align} W(\theta) &= \frac{1}{C_l^{(0)} + N(\theta)} \\ N(\theta) &\equiv \frac{\mbox{Pixel area}}{\mbox{ivar}(\theta)} \hspace{1cm} \mbox{"Local" noise power spectrum} \end{align}\end{split}\]

In implementation, in order to avoid divide-by-zero for ivar=0, we compute \(W(\theta)\) equivalently as:

\[W(\theta) = \frac{\mbox{ivar}(\theta)}{(\mbox{pixel area}) + C_l^{(0)} \mbox{ivar}(\theta)}\]

kszx.pixell_utils.sensitivity_curve(ivar, step, n)

Given an ivar map, computes the “sensitivity function” \(S(x)\) at x = [ step, …, n*step].

We define the “sensitivity function” \(S(x)\) to be the sky area (in deg^2) with sensitivity >= x, where x has units (uK-armcin)^{-2}. Note that \(S(x)\) is a decreasing function of x. This is useful for making a plot of sky sensitivity.

Function args:

  • ivar (pixell map): inverse variance map, units (uK)^{-2}. (For example, the return value from read_ivar().)

  • step: units (uK-arcmin)^{-2}, see above. Recommend step = 1.0e-5 for ACT.

  • n: number of returned points, see above.

Returns 1-d array [ S(step), S(2*step), ..., S(n*step) ].

Example usage:

step = 1.0e-5
nbins = 3000
ivar = kszx.act.read_ivar(150, dr=6, download=True)

x = step * np.arange(nbins)
Sx = kszx.pixell_utils.sensitivity_curve(ivar,step,nbins)

plt.plot(x, Sx)
plt.yscale('log')
plt.ylim(1, 10**5)
plt.xlabel(r'Threshold $x$ ($\mu$K-arcmin)$^{-2}$')
plt.ylabel(r'Sky area with sensitivity $> x$ (deg$^2$)')