import mapsims
import numpy as np
import healpy as hp
import pysm3.units as u
from pathlib import Path
%matplotlib inline
import matplotlib.pyplot as plt
The mapsims
package allows to generate maps on-the-fly, so they can be generated as needed inside other pipelines.
The interface is the SONoiseSimulator
class, specify nside
to work with HEALPix maps, shape
and wcs
for CAR (see lesson 6).
The class constructor just initializes the configuration of the class, doesn't load neither spectra nor hitmaps.
mapsims.SONoiseSimulator?
Init signature: mapsims.SONoiseSimulator( nside=None, shape=None, wcs=None, ell_max=None, return_uK_CMB=True, sensitivity_mode='baseline', apply_beam_correction=False, apply_kludge_correction=True, homogeneous=False, no_power_below_ell=None, rolloff_ell=50, survey_efficiency=0.2, full_covariance=True, LA_years=5, LA_noise_model='SOLatV3point1', elevation=50, SA_years=5, SA_one_over_f_mode='pessimistic', sky_fraction=None, cache_hitmaps=True, boolean_sky_fraction=False, instrument_parameters='simonsobs_instrument_parameters_2020.06', ) Docstring: <no docstring> Init docstring: Simulate noise maps for Simons Observatory Simulate the noise power spectrum in spherical harmonics domain and then generate a map in microK_CMB or microK_RJ (based on return_uK_CMB) In the constructor, this object calls the published 20180822 noise simulator and generates the expected noise power spectra for all channels. Then you need to call the `simulate` method with a channel identifier to create a simulated map. Parameters ---------- nside : int nside of HEALPix map. If None, uses rectangular pixel geometry specified through shape and wcs. shape : tuple of ints shape of ndmap array (see pixell.pixell.enmap). Must also specify wcs. wcs : astropy.wcs.wcs.WCS instance World Coordinate System for geometry of map (see pixell.pixell.enmap). Must also specify shape. ell_max : int Maximum ell for the angular power spectrum, if not provided set to 3 * nside when using healpix or 10000 * (1.0 / pixel_height_arcmin) when using CAR, corresponding roughly to the Nyquist frequency. return_uK_CMB : bool True, output is in microK_CMB, False output is in microK_RJ sensitivity_mode : str Value should be threshold, baseline or goal to use predefined sensitivities apply_beam_correction : bool Include the effect of the beam in the noise angular power spectrum apply_kludge_correction : bool If True, reduce the hitcount by a factor of 0.85 to account for not-uniformity in the scanning homogeneous : bool Set to True to generate full-sky maps with no hit-count variation, with noise curves corresponding to a survey that covers a sky fraction of sky_fraction (defaults to 1). no_power_below_ell : int The input spectra have significant power at low :math:`\ell`, we can zero that power specifying an integer :math:`\ell` value here. The power spectra at :math:`\ell < \ell_0` are set to zero. rolloff_ell : int Low ell power damping, see the docstring of `so_noise_models.so_models_v3.SO_Noise_Calculator_Public_v3_1_1.rolloff` survey_efficiency : float Fraction of calendar time that may be used to compute map depth. full_covariance : bool Whether or not to include the intra-tube covariance between bands. If white noise (atmosphere=False) sims are requested, no covariance is included regardless of the value of full_covariance. LA_years : int Total number of years for the Large Aperture telescopes survey LA_noise_model : str Noise model among the ones available in `so_noise_model`, "SOLatV3point1" is default, "SOLatV3" is the model released in 2018 which had a bug in the atmosphere contribution elevation : float Elevation of the scans in degrees, the V3.1.1 noise model includes elevation dependence for the LAT. This should reproduced original V3 results at the reference elevation of 50 degrees. SA_years : int Total number of years for the Small Aperture telescopes survey SA_one_over_f_mode : {"pessimistic", "optimistic", "none"} Correlated noise performance of the detectors on the Small Aperture telescopes sky_fraction : optional,float If homogeneous is True, this sky_fraction is used for the noise curves. cache_hitmaps : bool If True, caches hitmaps. boolean_sky_fraction: bool If True, determines sky fraction based on fraction of hitmap that is zero. If False, determines sky_fraction from <Nhits>. instrument_parameters : Path or str See the help of MapSims File: /global/common/software/sobs/mbs/lib/python3.7/site-packages/mapsims/noise.py Type: type Subclasses:
noise_sim = mapsims.SONoiseSimulator(nside=256)
The simulate
method does all the processing, it loads the HEALPix or CAR hitmaps (from a path when run at NERSC or downloading and caching them locally) and calls the so_noise_models
package to generate the expected spectra given the instrument configuration.
The hitmaps were generated in time domain executing parallel runs with TOAST and then saved to disk. They are only used as a relative weighting, the global noise properties are driven by the spectra from so_noise_models
.
The simulate
method gets the expected spectra, weights them by the sky fraction and then generates isotropic noise over the whole sky with synfast
for HEALPix and sym_expand
for CAR and then scales it by the relative hitmap.
It always simulates a dichroic tube at a time so that also the cross-correlation between the channels is simulated and captured in the output maps.
noise_sim.simulate?
Signature: noise_sim.simulate( tube, output_units='uK_CMB', seed=None, nsplits=1, mask_value=None, atmosphere=True, hitmap=None, white_noise_rms=None, ) Docstring: Create a random realization of the noise power spectrum Parameters ---------- tube : str Specify a tube (for SO: ST0-ST3, LT0-LT6) see the `tubes` attribute output_units : str Output unit supported by PySM.units, e.g. uK_CMB or K_RJ seed : integer or tuple of integers, optional Specify a seed. The seed is converted to a tuple if not already one and appended to (0,0,6,tube_id) to avoid collisions between tubes, with the signal sims and with ACT noise sims, where tube_id is the integer ID of the tube. nsplits : integer, optional Number of splits to generate. The splits will have independent noise realizations, with noise power scaled by a factor of nsplits, i.e. atmospheric noise is assumed to average down with observing time the same way the white noise does. By default, only one split (the coadd) is generated. mask_value : float, optional The value to set in masked (unobserved) regions. By default, it uses the value in default_mask_value, which for healpix is healpy.UNSEEN and for CAR is numpy.nan. atmosphere : bool, optional Whether to include the correlated 1/f from the noise model. This is True by default. If it is set to False, then a pure white noise map is generated from the white noise power in the noise model, and the covariance between arrays is ignored. hitmap : string or map, optional Provide the path to a hitmap to override the default used for the tube. You could also provide the hitmap as an array directly. white_noise_rms : float or tuple of floats, optional Optionally scale the simulation so that the small-scale limit white noise level is white_noise_rms in uK-arcmin (either a single number or a pair for the dichroic array). Returns ------- output_map : ndarray or ndmap Numpy array with the HEALPix or CAR map realization of noise. The shape of the returned array is (2,3,nsplits,)+oshape, where oshape is (npix,) for HEALPix and (Ny,Nx) for CAR. The first dimension of size 2 corresponds to the two different bands within a dichroic tube. See the `band_id` attribute of the Channel class to identify which is the index of a Channel in the array. The second dimension corresponds to independent split realizations of the noise, e.g. it is 1 for full mission. The third dimension corresponds to the three polarization Stokes components I,Q,U The last dimension is the number of pixels File: /global/common/software/sobs/mbs/lib/python3.7/site-packages/mapsims/noise.py Type: method
tube = "LT1"
noise_maps = noise_sim.simulate(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/LT1_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/LT1_UHF2_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/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 "
noise_maps.shape
(2, 1, 3, 786432)
channels = mapsims.parse_channels("tube:" + tube)[0]
channels
(Channel LT1_UHF1, Channel LT1_UHF2)
for ch, m in zip(channels, noise_maps):
hp.mollview(m[0][1], title="EE map " + ch.tag, unit="uK", min=-10, max=10)
We can take a cross-spectrum between the maps of the 2 tubes and compare with the expected spectrum from noise.get_fullsky_noise_spectra
.
ell, ps_T, ps_P = noise_sim.get_fullsky_noise_spectra(tube=tube)
hitmaps, ave_nhits = noise_sim.get_hitmaps(tube=tube)
sky_fraction = (hitmaps[0] > 0).sum() / len(hitmaps[0])
noise_maps.shape
(2, 1, 3, 786432)
x_C_ell = hp.anafast(hp.ma(noise_maps[0][0])*hitmaps[0], hp.ma(noise_maps[1][0])*hitmaps[1]) / sky_fraction / np.mean(hitmaps[0]*hitmaps[1])
x_C_ell.shape
(6, 768)
plt.figure(figsize=(8,5))
plt.loglog(ell, ps_P[2] * sky_fraction, lw=2, label="Input polarization Cross-spectrum")
plt.loglog(x_C_ell[1], label="EE Cross-spectrum", alpha=.5)
plt.loglog(x_C_ell[2], label="BB Cross-spectrum", alpha=.5)
plt.xlabel("$\ell$")
plt.ylabel("$C_\ell [\mu K^2]$")
plt.grid()
plt.ylim(1e-5, 1e-1)
plt.legend();
A single SONoiseSimulator
object is capable of simulating noise for different tubes. See for example a loop to simulate all SAT channels.
tubes = ["ST{}".format(i) for i in range(4)]
for tube in tubes:
noise_maps = noise_sim.simulate(tube)
for ch, m in zip(mapsims.parse_channels("tube:" + tube)[0], noise_maps):
hp.mollview(m[0][1], title="EE map " + ch.tag, unit="$uK$", min=-10, max=10)
/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}") /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 " /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/ST1_MFF1_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/ST1_MFF2_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/ST2_MFS1_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/ST2_MFS2_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/ST3_LF1_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/ST3_LF2_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz warnings.warn(f"Access data from {full_path}")
Override nsplits
to generate multiple splits, the output maps will be a 3 dimensional array where for each split we have the 3 IQU components.
noise_maps_splits = noise_sim.simulate("ST0", nsplits=4)
noise_maps_splits.shape
(2, 4, 3, 786432)