Compare WebSky Radio Galaxies maps from PySM to Planck

healpy
pysm
Published

December 5, 2022

I need to check that the amplitude of the maps of Radio Galaxies produced by WebSky are reasonable if compared to experimental results. The main worry is that being point sources, when we apply a beam using Spherical Harmonics transform we get spurious results due to the input maps not being band-limited.

pip install pysm3 --pre
Requirement already satisfied: pysm3 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (3.4.0b4.dev138+g864a61e.d20221114)
Requirement already satisfied: toml in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from pysm3) (0.10.2)
Requirement already satisfied: astropy in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from pysm3) (5.1.1)
Requirement already satisfied: numba in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from pysm3) (0.56.3)
Requirement already satisfied: healpy in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from pysm3) (1.16.1)
Requirement already satisfied: numpy>=1.18 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from astropy->pysm3) (1.23.4)
Requirement already satisfied: packaging>=19.0 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from astropy->pysm3) (21.3)
Requirement already satisfied: pyerfa>=2.0 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from astropy->pysm3) (2.0.0.1)
Requirement already satisfied: PyYAML>=3.13 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from astropy->pysm3) (6.0)
Requirement already satisfied: scipy in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from healpy->pysm3) (1.9.3)
Requirement already satisfied: matplotlib in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from healpy->pysm3) (3.6.2)
Requirement already satisfied: llvmlite<0.40,>=0.39.0dev0 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from numba->pysm3) (0.39.1)
Requirement already satisfied: setuptools in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from numba->pysm3) (65.5.1)
Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from packaging>=19.0->astropy->pysm3) (3.0.9)
Requirement already satisfied: contourpy>=1.0.1 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from matplotlib->healpy->pysm3) (1.0.6)
Requirement already satisfied: kiwisolver>=1.0.1 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from matplotlib->healpy->pysm3) (1.4.4)
Requirement already satisfied: python-dateutil>=2.7 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from matplotlib->healpy->pysm3) (2.8.2)
Requirement already satisfied: pillow>=6.2.0 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from matplotlib->healpy->pysm3) (9.2.0)
Requirement already satisfied: cycler>=0.10 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from matplotlib->healpy->pysm3) (0.11.0)
Requirement already satisfied: fonttools>=4.22.0 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from matplotlib->healpy->pysm3) (4.38.0)
Requirement already satisfied: six>=1.5 in /global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages (from python-dateutil>=2.7->matplotlib->healpy->pysm3) (1.16.0)
Note: you may need to restart the kernel to use updated packages.
import healpy as hp
import pysm3
from pysm3 import units as u
%matplotlib inline

Get the Planck beam widths

!wget -nc https://irsa.ipac.caltech.edu/data/Planck/release_3/ancillary-data/LFI_RIMO_R3.31.fits
File ‘LFI_RIMO_R3.31.fits’ already there; not retrieving.

!wget -nc https://irsa.ipac.caltech.edu/data/Planck/release_3/ancillary-data/HFI_RIMO_R3.00.fits
File ‘HFI_RIMO_R3.00.fits’ already there; not retrieving.

from astropy.io import fits
fwhm = {}

with fits.open("LFI_RIMO_R3.31.fits") as f:
  for ch, ch_fwhm in zip(f[2].data["FREQUENCY"], f[2].data["FWHM"]):
    fwhm[ch[0]] = ch_fwhm * u.arcmin
with fits.open("HFI_RIMO_R3.00.fits") as f:
  for ch, ch_fwhm in zip(f[1].data["FREQUENCY"], f[1].data["FWHM"]):
    fwhm[ch[0]] = ch_fwhm * u.arcmin

Select a channel

“030” or “143”

ch = "030"
freq = float(ch) * u.GHz
instrument = "LFI" if freq < 100 * u.GHz else "HFI"
nside = 1024 if instrument == "LFI" else 2048
rel = "R3" if instrument == "LFI" else "R4"

Load Planck maps

!wget -nc https://irsa.ipac.caltech.edu/data/Planck/release_3/all-sky-maps/maps/LFI_SkyMap_030_1024_R3.00_full.fits
File ‘LFI_SkyMap_030_1024_R3.00_full.fits’ already there; not retrieving.

!wget -nc https://irsa.ipac.caltech.edu/data/Planck/release_3/all-sky-maps/maps/HFI_SkyMap_143_2048_R4.00_full.fits
File ‘HFI_SkyMap_143_2048_R4.00_full.fits’ already there; not retrieving.

m_planck = hp.read_map(f"{instrument}_SkyMap_{ch}_{nside}_{rel}.00_full.fits")
m_planck = (m_planck * u.K_CMB).to(u.uK_CMB)
hp.mollview(hp.remove_dipole(m_planck), max=1e3)

Create maps with PySM

sky = pysm3.Sky(nside=nside, preset_strings=["rg1"])
m_pysm = sky.get_emission(freq)
m_pysm_smoothed = pysm3.apply_smoothing_and_coord_transform(m_pysm, fwhm=fwhm[ch])
hp.map2alm_lsq did not converge in 10 iterations, residual relative error is 0.69
hp.mollview(m_pysm_smoothed[0], max=1e3)

Visualize a point source

Check for ringing around point sources

max_pix = m_pysm_smoothed[0].argmax()
lon,lat= hp.pix2ang(nside, max_pix, lonlat=True)
lon
275.2809917355372
lat
-62.0850844513906
hp.gnomview(m_pysm_smoothed[0], rot=(lon,lat))

m_pysm_smoothed[0].max()

\(156551.45 \; \mathrm{\mu K_{{RJ}}}\)

m_planck.max()

\(215075.31 \; \mathrm{\mu K_{{CMB}}}\)