Codes for making plots for O3 anisotropic stochastic paper

Importing packages and setting things up;

Before proceeding make sure that all the packages listed below are already installed

In [1]:
import numpy as np
import healpy as hp
import h5py
import scipy.io as sio
import matplotlib # requires matplotlib version >= 3.2.2
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from scipy.io import loadmat
import matplotlib.patches as mpatches
from matplotlib.lines import Line2D
import warnings
warnings.filterwarnings("ignore", category=UserWarning) 

Parameters for plotting

In [2]:
fontsize = 24
matplotlib.rcParams.update({
    "font.size": fontsize,
    "axes.titlesize": fontsize,
    "axes.labelsize": fontsize,
    "xtick.labelsize": fontsize,
    "ytick.labelsize": fontsize,
    "xtick.major.size": fontsize * 0.8,
    "ytick.major.size": fontsize * 0.8,
    "legend.fontsize": fontsize * 0.9,
    "font.family": "Times new Roman",
    "figure.dpi": 100,
    "savefig.dpi": 300,
    "text.usetex": True,
    "path.simplify": True,
    "figure.figsize": (8, 6)
})

# Class to deal in an easy way with number of digits in colorbars (used for ULs and sigma maps)
# It allows to choose the number of digits for tickslabels and keeping 10^N at the right of the colorbar
# instead of having 1eN for each tickslabel.

class FormatScalarFormatter(matplotlib.ticker.ScalarFormatter):
            def __init__(self, fformat="%1.1f", offset=True, mathText=True):
                self.fformat = fformat
                matplotlib.ticker.ScalarFormatter.__init__(self,useOffset=offset,
                                                        useMathText=mathText)
            def _set_format(self):
                self.format = self.fformat
                if self._useMathText:
                    #self.format = '$%s$' % matplotlib.ticker._mathdefault(self.format)
                    self.format = '$%s$' % ('\\mathdefault{%s}' % self.format)

Variable definitions

In [3]:
# variable corresponding to spectral shapes of signal models 
alphas = ["0", "2over3", "3"]

Figure 2: Top row (BBR SNR maps)

In [4]:
# SNR sky-maps for BBR analysis

# Create figure
fig = plt.figure(figsize=(24, 8))

# Limits for the colorbars
vlims = [2.55, 3.72, 3.63]

# Loop over the *combined_maps.hdf5 files for the three spectral indices 0, 2/3 and 3
for alpha, vlim, index in zip (alphas, vlims, range(0,len(alphas))):
    
    # SNR *.hdf5 file
    snr_file = "./datafiles/bbr_a" + alpha + "_o1_o2_o3_combined_maps.hdf5"
    
    # Read the file and create a dictionary-like variable; the file has three variables SNR (map_snr_total), signal estimate (ptEst_total) and its sigma (sig_total)
    combined_snr = h5py.File(snr_file, "r")

    # Select the map of interest (SNR map)
    snr_map = np.array(combined_snr["map_snr_total"])
    print("alpha:%s, max snr:%.1f" % (alpha, np.max(snr_map)))

    
    # Close the file
    combined_snr.close()
    
    # Plot the sky-map for a given alpha as subplot using healpy Mollweide projection
    hp.mollview(np.real(np.squeeze(snr_map)),rot=(180,0,0),flip='astro', title="",nest=False, format = "%.2f", cbar=False, sub=(1, 3, index+1), margins=[0.0125,0.0125, 0.0125, 0.0125], min=-vlim, max=vlim)
    
    # Add graticule
    hp.graticule(45,90)
    
    # Add useful annotations on the map for right ascension and declination
    plt.annotate('12h',(0.46,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
    plt.annotate('6h',(0.725,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
    plt.annotate('18h',(0.21,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
    plt.annotate(r'$45^{\circ}$',(0.03,0.8),xycoords='axes fraction', fontsize=fontsize)
    plt.annotate(r'$0^{\circ}$',(-0.045,0.47),xycoords='axes fraction', fontsize=fontsize)
    plt.annotate(r'$-45^{\circ}$',(-0.002,0.14),xycoords='axes fraction', fontsize=fontsize) 
    
    # Trick to make colorbars all the same (healpy colorbar does not leave much freedom of customization and it is not used)
    fig = plt.gcf()
    ax = plt.gca()  
    image = ax.get_images()[0]
    cmap = fig.colorbar(image, ax=ax, shrink=0.7, orientation="horizontal", pad=0.15, aspect = 40, ticks=[-2,0,2])
plt.show()
alpha:0, max snr:2.6
alpha:2over3, max snr:2.7
alpha:3, max snr:3.6

Figure 2: Bottom row (BBR Upper limit maps)

In [5]:
# Upper limit sky-maps for BBR analysis

fig = plt.figure(figsize=(24, 8))

# Loop over the upper limit *.mat files for the three spectral indices 0, 2/3 and 3
for alpha, index in zip (alphas, range(0,len(alphas))):

    # upper limit *.mat file
    ul_dir = "./datafiles/bbr_a"+ alpha +"_zl_limits.mat"
    
    # Read the file and create a dictionary-like variable; the file has variables correponding to upplimits from O1+O2, just O3 (for three baselines) and O1+O2+O3
    combined_ul = sio.loadmat(ul_dir)
    ul_map = np.squeeze(np.array(combined_ul["flux_ul_combined"]))
    print("alpha:%s, min UL:%.1e, max UL:%.1e" % (alpha, np.min(ul_map), np.max(ul_map)))
    
    # Plot the sky-map for a given alpha as subplot using healpy Mollweide projection
    hp.mollview(np.real(np.squeeze(ul_map)),rot=(180,0,0),flip='astro', title="",nest=False, format = "%.2g", cbar=False, sub=(1, 3, index+1), margins=[0.0125,0.0125, 0.0125, 0.0125])

    # Add graticule
    hp.graticule(45,90)
    
    # Add useful annotations on the map for right ascension and declination
    plt.annotate('12h',(0.46,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
    plt.annotate('6h',(0.725,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
    plt.annotate('18h',(0.21,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
    plt.annotate(r'$45^{\circ}$',(0.03,0.8),xycoords='axes fraction', fontsize=fontsize)
    plt.annotate(r'$0^{\circ}$',(-0.045,0.47),xycoords='axes fraction', fontsize=fontsize)
    plt.annotate(r'$-45^{\circ}$',(-0.002,0.14),xycoords='axes fraction', fontsize=fontsize) 
    
    
    fig = plt.gcf()
    ax = plt.gca()  
    image = ax.get_images()[0]
    
    # Further colorbar customization
    nticks = 4
    interval = 1./(nticks-1)
    diff= np.max(ul_map) -np.min(ul_map)
    ticks_no_round = [np.min(ul_map)+(diff)*interval*x for x in np.arange(0,nticks)]
    
    cmap = fig.colorbar(image, ax=ax, ticks=ticks_no_round, format=FormatScalarFormatter("%.1f"), shrink=0.7, orientation="horizontal", pad=0.15, aspect = 30)
    
    # Adding units as colorbar label
    cmap.set_label(r'[erg cm$^{-2}$ Hz$^{-1}$ s$^{-1}$]', fontsize = 1.2*fontsize, labelpad = 20)
    cmap.ax.tick_params(labelsize=fontsize)
plt.show()
alpha:0, min UL:1.7e-08, max UL:7.6e-08
alpha:2over3, min UL:8.5e-09, max UL:4.1e-08
alpha:3, min UL:1.3e-10, max UL:1.1e-09