Support additional cylindrical (aka equiangular) sampling schemes

[Magwos]

In CMB data analysis (aiming at probing small angular scales), an interest has been growing towards cylindrical griding, named CAR in https://phy-act1.princeton.edu/public/saiola/act_dr4_C20.pdf or ECP in https://www.aanda.org/articles/aa/full_html/2013/08/aa21575-13/aa21575-13.html. It corresponds to a rectangular projection, equi-spaced in latitude, and not equal-area projection, but with the same number of pixels per ring, and relies on Clenshaw-Curtis quadrature weightings for the associated transforms.

This sampling scheme is the baseline one used for the ACT (Atacama Cosmology Telescope) collaboration and central in many tools developed for the LAT analysis in the Simons Observatory collaboration (in which I am part of).

I wanted to ask if it would be possible to have s2fft support for such sampling scheme.

[jasonmcewen]

Thanks @Magwos, sounds like a very useful feature. We can certainly add this. Note that MW/DH equiangular sampling is already supported but I appreciate you probably want exactly the same sampling as used elsewhere.

Do you have a reference (paper and/or code) with the precise details of the sampling? I had a super quick look at the references you posted but didn't immediately see the details there.

[Magwos]

Thanks a lot for your positive answer!

So, currently the implementation that is most used is done through the package https://github.com/simonsobs/pixell (in particular https://github.com/simonsobs/pixell/blob/master/pixell/curvedsky.py), which itself calls the ducc package (https://mtr.pages.mpcdf.de/ducc/sht.html#module-ducc0.sht) for their transforms.

An example of CMB-science oriented package calling pixell is this one: https://github.com/simonsobs/pspy
Pspy is mostly used towards small-scale science, and it is usually applied to partial sky maps with very high resolution (typically nside=4096, the SO LAT will have a resolution of 1 arcmin), storing only the relevant pixels and not the 0 of the unobserved area (and of course the coordinates of the patch observed).

If you need more details about all of this, please let me know.

[jasonmcewen]

Thanks @Magwos !

Adding further details from your email that you sent separately so we have everything here in one place.

there is no paper detailing this CAR pixelisation which was adopted by ACT (named after the Plate-Carré projection, or Equirectangular one), as it corresponds in fact to the ECP (equidistant cylindrical projection) grid also detailed here: https://www.aanda.org/articles/aa/pdf/2013/06/aa21494-13.pdf (4.4.1)

[jasonmcewen]

@Magwos Could you let us know the precise sample positions that are used for ACT, i.e. (theta, phi) coordinates for a given L or $\ell_\text{max}$, so we can be sure we support what you need?

[Magwos]

Hi @jasonmcewen, thanks a lot for your message!
My most demanding needs regarding the sky area and $\ell_{\rm max}$ will be more regarding Simons Observatory than ACT (the SO LAT analysis will use the same sampling scheme), and are basically given by (from the SO forecast paper here https://arxiv.org/pdf/1808.07445):

• $\ell_{\rm max} = 8000$
• $f_{\rm sky}$ possibly reaching 40% of the sky, as represented in Figure 3 of the article (with the galaxy plane avoided)

Here is an example script using healpy to reproduce the large area represented in the Figure, spanning phi between 0 and 360 and theta between 55 and 165

nside = 512

deg_to_rad = np.pi/180

theta_max = 35 + 130
theta_min = -75 + 130
phi_max = 360 
phi_min = -0 

map_hp_template = np.ones(12*nside**2)

theta, phi = hp.pix2ang(nside, np.arange(12*nside**2))

cond_1 = np.logical_or(theta>theta_max*deg_to_rad, theta<theta_min*deg_to_rad)
cond_2 = np.logical_or(phi>phi_max*deg_to_rad, phi<phi_min*deg_to_rad)
cond = np.logical_or(cond_1, cond_2)


map_hp_template[cond] = 0

hp.mollview(map_hp_template, title='HEALPix template', cbar=True)
print(map_hp_template.sum()/map_hp_template.size)

plt.show()

And here another script mimicking a cut of the galaxy plane

nside = 512

deg_to_rad = np.pi/180

theta_max = 35 + 130
theta_min = -75 + 130
phi_max = 270 
phi_min = 190 
phi_max_2 = 160
phi_min_2 = 80

map_hp_template = np.ones(12*nside**2)

theta, phi = hp.pix2ang(nside, np.arange(12*nside**2))

cond_1 = np.logical_or(theta>theta_max*deg_to_rad, theta<theta_min*deg_to_rad)
cond_phi_1 = np.logical_and(phi<phi_max*deg_to_rad, phi>phi_min*deg_to_rad)
cond_phi_2 = np.logical_and(phi<phi_max_2*deg_to_rad, phi>phi_min_2*deg_to_rad)
cond_phi = np.logical_or(cond_phi_1, cond_phi_2)
cond = np.logical_or(cond_1, cond_phi)


map_hp_template[cond] = 0

hp.mollview(map_hp_template, title='HEALPix template', cbar=True)

print(map_hp_template.sum()/map_hp_template.size)
plt.show()

[jasonmcewen]

Thanks @Magwos. I just had a chance to look into this. I was actually after the precise (theta,phi) sampling positions for the grid you're interested in, rather than coverage.

In any case, after looking into the standard Clenshaw-Curtis sample positions (i.e. the Fejer type 2 sample positions with poles), as described in, e.g., Waldvogel (2006), it turns out the sample positions are identical to what we call mwss sampling (as first introduced this old paper and also discussed in our recent paper).

So Clenshaw-Curtis samples are already supported, you just need to set sampling=mwss.

At present we aren't using Clenshaw-Curtis quadrature but rather the exact quadrature outlined in our papers. While Clenshaw-Curtis quadrature is not exact for the numbers of samples used for the mwss grid, it's nevertheless highly accurate. So we will still try to support Clenshaw-Curtis quadrature in furture.