import mapsims
import numpy as np
import healpy as hp
import pysm3.units as u
from pathlib import Path
%matplotlib inline
The most important repository for Simons Observatory map based simulations is:
It contains all the documentation about available released map and software, in particular:
mapsims
¶channels = mapsims.parse_channels("tube:ST0")
print(channels)
[(Channel ST0_UHF1, Channel ST0_UHF2)]
channel = channels[0][0]
channel.beam
Check the documentation of astropy.units
for details about handling units.
channel.center_frequency
channel.beam.to(u.deg)
for ch in mapsims.parse_channels("telescope:SA"):
print(ch)
Channel ST0_UHF1 Channel ST0_UHF2 Channel ST1_MFF1 Channel ST1_MFF2 Channel ST2_MFS1 Channel ST2_MFS2 Channel ST3_LF1 Channel ST3_LF2
mapsims.parse_channels("030", instrument_parameters="planck_deltabandpass")
[Channel 030]
for ch in mapsims.parse_channels("all") + mapsims.parse_channels("all", instrument_parameters="planck_deltabandpass"):
print(ch.telescope, ch, ch.beam, ch.center_frequency)
LA Channel LT0_UHF1 1.0 arcmin 225.7 GHz LA Channel LT0_UHF2 0.9 arcmin 285.4 GHz LA Channel LT1_UHF1 1.0 arcmin 225.7 GHz LA Channel LT1_UHF2 0.9 arcmin 285.4 GHz LA Channel LT2_MFF1 2.2 arcmin 92.0 GHz LA Channel LT2_MFF2 1.4 arcmin 147.5 GHz LA Channel LT3_MFF1 2.2 arcmin 92.0 GHz LA Channel LT3_MFF2 1.4 arcmin 147.5 GHz LA Channel LT4_MFS1 2.2 arcmin 88.6 GHz LA Channel LT4_MFS2 1.4 arcmin 146.5 GHz LA Channel LT5_MFS1 2.2 arcmin 88.6 GHz LA Channel LT5_MFS2 1.4 arcmin 146.5 GHz LA Channel LT6_LF1 7.4 arcmin 25.7 GHz LA Channel LT6_LF2 5.1 arcmin 38.9 GHz SA Channel ST0_UHF1 19.0 arcmin 225.7 GHz SA Channel ST0_UHF2 17.0 arcmin 285.4 GHz SA Channel ST1_MFF1 42.0 arcmin 92.0 GHz SA Channel ST1_MFF2 27.0 arcmin 147.5 GHz SA Channel ST2_MFS1 42.0 arcmin 88.6 GHz SA Channel ST2_MFS2 27.0 arcmin 146.5 GHz SA Channel ST3_LF1 144.0 arcmin 25.7 GHz SA Channel ST3_LF2 99.0 arcmin 38.9 GHz planck Channel 030 33.102652125 arcmin 28.4 GHz planck Channel 044 27.94348615 arcmin 44.1 GHz planck Channel 070 13.07645961 arcmin 70.4 GHz planck Channel 100 9.682 arcmin 100.89 GHz planck Channel 143 7.303 arcmin 142.876 GHz planck Channel 217 5.021 arcmin 221.156 GHz planck Channel 353 4.944 arcmin 357.5 GHz planck Channel 545 4.831 arcmin 555.2 GHz planck Channel 857 4.638 arcmin 866.8 GHz
PySM 3 also adds the capability of converting between $K_{CMB}$, $K_{RJ}$ and $MJ/sr$
(100 * u.uK_CMB).to(u.uK_RJ, equivalencies=u.cmb_equivalencies(channel.center_frequency))
(100 * u.uK_CMB).to(u.MJy/u.sr, equivalencies=u.cmb_equivalencies(channel.center_frequency))
folder = "/project/projectdirs/sobs/v4_sims/mbs/201906_highres_foregrounds_extragalactic_tophat"
!ls $folder/512
ame cmb_lensed_solardipole combined_solardip ksz cib cmb_unlensed dust synchrotron cmb combined freefree tsz
filename_template = "{nside}/{content}/{num:04d}/simonsobs_{content}_uKCMB_{telescope}{band}_nside{nside}_{num:04d}.fits"
filename = Path(folder) / filename_template.format(nside=512, content="combined", num=0,
telescope=channel.telescope.lower(), band=channel.band)
print(filename)
/project/projectdirs/sobs/v4_sims/mbs/201906_highres_foregrounds_extragalactic_tophat/512/combined/0000/simonsobs_combined_uKCMB_saUHF1_nside512_0000.fits
m = hp.read_map(filename, (0,1,2))
/global/common/software/sobs/mbs/lib/python3.7/site-packages/healpy/fitsfunc.py:352: UserWarning: If you are not specifying the input dtype and using the default np.float64 dtype of read_map(), please consider that it will change in a future version to None as to keep the same dtype of the input file: please explicitly set the dtype if it is important to you. "If you are not specifying the input dtype and using the default "
NSIDE = 512 ORDERING = NESTED in fits file INDXSCHM = IMPLICIT Ordering converted to RING Ordering converted to RING Ordering converted to RING
m_P = np.sqrt(m[1]**2 + m[2]**2)
hp.mollview(m_P, min=0, max=300, title="Polarization amplitude: " + channel.tag, unit="uK_CMB")
Map-based releases also include standard Galactic masks that can be used for analysis, see the documentation.
!ls /project/projectdirs/sobs/v4_sims/mbs/201904_highres_foregrounds_equatorial/512/masks
mask_equatorial_pol_thr1p0_fsky0p28_ns512.fits mask_equatorial_pol_thr2p0_fsky0p50_ns512.fits mask_equatorial_pol_thr3p0_fsky0p62_ns512.fits mask_equatorial_pol_thr3p5_fsky0p67_ns512.fits mask_equatorial_pol_thr4p0_fsky0p71_ns512.fits mask_equatorial_pol_thr4p5_fsky0p75_ns512.fits
mask_filename = "mask_equatorial_pol_thr4p0_fsky0p71_ns512.fits"
mask = hp.read_map("/project/projectdirs/sobs/v4_sims/mbs/201904_highres_foregrounds_equatorial/512/masks/"+mask_filename,
dtype=np.bool)
NSIDE = 512 ORDERING = RING in fits file INDXSCHM = IMPLICIT
mask
array([ True, True, True, ..., False, False, False])
hp.mollview(
mask
, coord="CG")
m_P_masked = m_P.copy()
m_P_masked[np.logical_not(mask)] = hp.UNSEEN
hp.mollview(m_P_masked)
The mapsims.SONoiseSimulator
class, which is used to simulate noise, can also automatically download hitmaps that were simulated by TOAST in time domain.
noise = mapsims.SONoiseSimulator(nside=512)
hitmaps, ave_nhits = noise.get_hitmaps(tube=channel.tube)
/global/common/software/sobs/mbs/lib/python3.7/site-packages/mapsims/utils.py:71: UserWarning: Access data from /global/cfs/cdirs/sobs/www/so_mapsims_data/v0.2/healpix/ST0_UHF1_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz warnings.warn(f"Access data from {full_path}") /global/common/software/sobs/mbs/lib/python3.7/site-packages/mapsims/utils.py:71: UserWarning: Access data from /global/cfs/cdirs/sobs/www/so_mapsims_data/v0.2/healpix/ST0_UHF2_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz warnings.warn(f"Access data from {full_path}")
For this simulation set the hitmaps are the same for both channels but in principle they might be different.
(hitmaps[0] - hitmaps[1]).sum()
0.0
m_P_survey = m_P.copy()
m_P_survey[hitmaps[0] == 0] = hp.UNSEEN
hp.mollview(m_P_survey)