RFI Inspection Daily RTP Notebook¶

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import glob
import os
from astropy import units
from copy import deepcopy
from pyuvdata import UVFlag
from SSINS import INS
from SSINS import version as SSINS_ver
import matplotlib.colors as colors
from matplotlib import cm
from IPython.display import display, HTML

%matplotlib inline
%config InlineBackend.figure_format = 'retina'
display(HTML("<style>.container { width:100% !important; }</style>"))
In [2]:
# Use environment variables to figure out path to data
JD = os.environ['JULIANDATE']
data_path = os.environ['DATA_PATH']
print(f'JD = {JD}')
print(f'data_path = "{data_path}"')
JD = int(JD)
JD = 2459930
data_path = "/mnt/sn1/2459930"
In [3]:
from astropy.time import Time
utc = Time(JD, format='jd').datetime
print(f'Date: {utc.month}-{utc.day}-{utc.year}')
Date: 12-16-2022
In [4]:
uvf = UVFlag(f'{data_path}/zen.{JD}.total_stage_1_threshold_flags.h5')
In [5]:
plt.figure(figsize=(16,12))
plt.imshow(uvf.flag_array[:,:,0], aspect='auto', interpolation='none',
           extent=[uvf.freq_array[0] / 1e6, uvf.freq_array[-1] / 1e6, 
                   uvf.time_array[-1] - JD, uvf.time_array[0] - JD])
plt.xlabel('Frequency (MHz)')
plt.ylabel(f'JD - {JD}')
ax2 = plt.gca().twinx()
ax2.set_ylim([uvf.lst_array[0] * 12 / np.pi, uvf.lst_array[-1] * 12 / np.pi])
ax2.invert_yaxis()
ax2.set_ylabel('LST (hours)')
ax3 = plt.gca().twiny()
ax3.set_xlim([0, uvf.Nfreqs - 1])
ax3.set_xlabel('Channel');

Figure 1(a): Full day of XRFI flags¶

Yellow is flagged. Blue is unflagged.

In [6]:
ssins_dirs_sorted = sorted(glob.glob(f"{data_path}/zen.{JD}*.SSINS"))
init_ssins_path = glob.glob(f"{ssins_dirs_sorted[0]}/*flags.h5")
ssins_uvf = UVFlag(init_ssins_path)
for path in ssins_dirs_sorted[1:]:
    new_path = glob.glob(f"{path}/*flags.h5")[0]
    ssins_uvf += UVFlag(new_path)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In [6], line 3
      1 ssins_dirs_sorted = sorted(glob.glob(f"{data_path}/zen.{JD}*.SSINS"))
      2 init_ssins_path = glob.glob(f"{ssins_dirs_sorted[0]}/*flags.h5")
----> 3 ssins_uvf = UVFlag(init_ssins_path)
      4 for path in ssins_dirs_sorted[1:]:
      5     new_path = glob.glob(f"{path}/*flags.h5")[0]

File ~/mambaforge/envs/RTP/lib/python3.10/site-packages/pyuvdata/uvflag/uvflag.py:521, in UVFlag.__init__(self, indata, mode, copy_flags, waterfall, history, label, use_future_array_shapes, run_check, check_extra, run_check_acceptability)
    518 self.label = ""  # Added to at the end
    519 if isinstance(indata, (list, tuple)):
    520     self.__init__(
--> 521         indata[0],
    522         mode=mode,
    523         copy_flags=copy_flags,
    524         waterfall=waterfall,
    525         history=history,
    526         label=label,
    527         use_future_array_shapes=use_future_array_shapes,
    528         run_check=run_check,
    529         check_extra=check_extra,
    530         run_check_acceptability=run_check_acceptability,
    531     )
    532     if len(indata) > 1:
    533         for i in indata[1:]:

IndexError: list index out of range
In [7]:
plt.figure(figsize=(16,12))
plt.imshow(ssins_uvf.flag_array[:,:,0], aspect='auto', interpolation='none',
           extent=[ssins_uvf.freq_array[0] / 1e6, ssins_uvf.freq_array[-1] / 1e6, 
                   ssins_uvf.time_array[-1] - JD, ssins_uvf.time_array[0] - JD])
plt.xlabel('Frequency (MHz)')
plt.ylabel(f'JD - {JD}')
ax2 = plt.gca().twinx()
ax2.set_ylim([ssins_uvf.lst_array[0] * 12 / np.pi, ssins_uvf.lst_array[-1] * 12 / np.pi])
ax2.invert_yaxis()
ax2.set_ylabel('LST (hours)')
ax3 = plt.gca().twiny()
ax3.set_xlim([0, ssins_uvf.Nfreqs - 1])
ax3.set_xlabel('Channel');
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [7], line 2
      1 plt.figure(figsize=(16,12))
----> 2 plt.imshow(ssins_uvf.flag_array[:,:,0], aspect='auto', interpolation='none',
      3            extent=[ssins_uvf.freq_array[0] / 1e6, ssins_uvf.freq_array[-1] / 1e6, 
      4                    ssins_uvf.time_array[-1] - JD, ssins_uvf.time_array[0] - JD])
      5 plt.xlabel('Frequency (MHz)')
      6 plt.ylabel(f'JD - {JD}')

NameError: name 'ssins_uvf' is not defined
<Figure size 1600x1200 with 0 Axes>

Figure 1(b): Full day of SSINS flags¶

FM is manually flagged at the beginning.

In [8]:
plt.figure(figsize=(16,12))
flag_table = np.zeros_like(uvf.flag_array[:, :, 0]).astype(float)
flag_table[np.logical_and(ssins_uvf.flag_array[:, :, 0], uvf.flag_array[:, :, 0])] = 0.75
flag_table[np.logical_and(ssins_uvf.flag_array[:, :, 0], np.logical_not(uvf.flag_array[:, :, 0]))] = 0.5
flag_table[np.logical_and(np.logical_not(ssins_uvf.flag_array[:, :, 0]), uvf.flag_array[:, :, 0])] = 0.25

# Prepare a colormap.
cmap = plt.cm.colors.ListedColormap(
    ["slategray", "darkturquoise", "plum", "lemonchiffon"]
)
print(cmap.__dict__)


cax = plt.imshow(flag_table, aspect='auto', interpolation='none', 
                 extent=[uvf.freq_array[0] / 1e6, uvf.freq_array[-1] / 1e6, 
                 uvf.time_array[-1] - JD, uvf.time_array[0] - JD],
                 cmap=cmap,
                 vmin=0, vmax=1)
plt.xlabel('Frequency (MHz)')
plt.ylabel(f'JD - {JD}')
ax2 = plt.gca().twinx()
ax2.set_ylim([ssins_uvf.lst_array[0] * 12 / np.pi, ssins_uvf.lst_array[-1] * 12 / np.pi])
ax2.invert_yaxis()
ax2.set_ylabel('LST (hours)')
ax3 = plt.gca().twiny()
ax3.set_xlim([0, ssins_uvf.Nfreqs - 1])
ax3.set_xlabel('Channel');

cbar_ticklabels = ["Flagged in Neither", "Flagged in XRFI", "Flagged in SSINS", "Flagged in Both"]

# Configure the colorbar so that labels are at the center of each section.
cbar = plt.colorbar(cax)
cbar_ticks = np.arange(0.125, 1.125, 0.25)
cbar.set_ticks(cbar_ticks)
cbar.set_ticklabels(cbar_ticklabels);
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [8], line 3
      1 plt.figure(figsize=(16,12))
      2 flag_table = np.zeros_like(uvf.flag_array[:, :, 0]).astype(float)
----> 3 flag_table[np.logical_and(ssins_uvf.flag_array[:, :, 0], uvf.flag_array[:, :, 0])] = 0.75
      4 flag_table[np.logical_and(ssins_uvf.flag_array[:, :, 0], np.logical_not(uvf.flag_array[:, :, 0]))] = 0.5
      5 flag_table[np.logical_and(np.logical_not(ssins_uvf.flag_array[:, :, 0]), uvf.flag_array[:, :, 0])] = 0.25

NameError: name 'ssins_uvf' is not defined
<Figure size 1600x1200 with 0 Axes>

Figure 1(c): Flag Agreement/Disagreement¶

In [9]:
xrfi_dirs = sorted(glob.glob(f'{data_path}/zen.{JD}.?????.stage_1_xrfi'))
print(f'Found {len(xrfi_dirs)} directories containing XRFI intermediate data products.')
files1 = [glob.glob(f'{d}/*combined_metrics1.h5')[0] for d in xrfi_dirs]
print(f'Found {len(files1)} combined round 1 XRFI metrics files.')
files2 = [glob.glob(f'{d}/*combined_metrics2.h5')[0] for d in xrfi_dirs]
print(f'Found {len(files2)} combined round 2 XRFI metrics files.')
uvf1 = UVFlag(files1)
uvf2 = UVFlag(files2)
uvf2.metric_array = np.where(np.isinf(uvf2.metric_array), uvf1.metric_array,
                             uvf2.metric_array)
Found 1849 directories containing XRFI intermediate data products.
Found 1849 combined round 1 XRFI metrics files.
Found 1849 combined round 2 XRFI metrics files.
In [10]:
plt.figure(figsize=(16,12))
max_abs = 100
if np.max(uvf2.metric_array) > max_abs:
    extend = 'max'
    if np.min(uvf2.metric_array) < -max_abs:
        extend = 'both'
elif np.min(uvf2.metric_array) < -max_abs:
    extend = 'min'    
else:
    extend = 'neither'

plt.imshow(uvf2.metric_array[:,:,0], aspect='auto', cmap='RdBu_r', interpolation='none', 
           norm=colors.SymLogNorm(linthresh=1,vmin=-max_abs, vmax=max_abs), 
           extent=[uvf.freq_array[0] / 1e6, uvf.freq_array[-1] / 1e6, 
                   uvf.time_array[-1] - JD, uvf.time_array[0] - JD])
plt.colorbar(pad=.07, extend=extend,
             label='RFI Detection Significance ($\sigma$s)')
plt.title('Combined XRFI Metrics')
plt.xlabel('Frequency (MHz)')
plt.ylabel(f'JD - {JD}')
ax2 = plt.gca().twinx()
ax2.set_ylim([uvf.lst_array[0] * 12 / np.pi, uvf.lst_array[-1] * 12 / np.pi])
ax2.invert_yaxis()
ax2.set_ylabel('LST (hours)')
ax3 = plt.gca().twiny()
ax3.set_xlim([0, uvf.Nfreqs - 1])
ax3.set_xlabel('Channel');

Figure 2(a): Combined XRFI Detection Significance¶

This figure shows round 2 XRFI metrics (mean filter outliers) combined in quadrature. When flagged in round 1 of XRFI, round 1's combined median filter metrics are used instead.

In [11]:
def make_sig_arr(ins, sig_arr, event_arr):
    for event in ins.match_events:
        nomask = np.logical_not(ins.metric_ms.mask[event[:2]])
        if event.sig is not None:
            event_arr[event[:2]] = event.sig
        else:
            event_arr[event[:2]][nomask] = ins.metric_ms[event[:2]][nomask]
        sig_arr[event[:2]][nomask] = ins.metric_ms[event[:2]][nomask]
        ins.metric_array[event[:2]] = np.ma.masked
        ins.metric_ms[:, event[1]] = ins.mean_subtract(freq_slice=event.freq_slice)
    nomask = np.logical_not(ins.metric_ms.mask)
    sig_arr[nomask] = ins.metric_ms[nomask]
    event_arr[nomask] = ins.metric_ms[nomask]
            
        
    return(sig_arr, event_arr)
In [12]:
init_ssins_path_data = glob.glob(f"{ssins_dirs_sorted[0]}/*data.h5")[0]
init_ssins_path_match_events = glob.glob(f"{ssins_dirs_sorted[0]}/*match_events.yml")[0]
init_ssins = INS(init_ssins_path_data, match_events_file=init_ssins_path_match_events)
init_sig_arr = np.ma.copy(init_ssins.metric_ms)
init_event_arr = np.ma.copy(init_sig_arr)
sig_arr, event_arr = make_sig_arr(init_ssins, init_sig_arr, init_event_arr)

for ssins_dir in ssins_dirs_sorted[1:]:
    ssins_path_data = glob.glob(f"{ssins_dir}/*data.h5")[0]
    ssins_path_match_events = glob.glob(f"{ssins_dir}/*match_events.yml")[0]
    ssins = INS(ssins_path_data, match_events_file=ssins_path_match_events)
    init_ssins += ssins
    
    new_sig_arr = np.ma.copy(ssins.metric_ms)
    new_event_arr = np.ma.copy(new_sig_arr)
    
    new_sig_arr, new_event_arr = make_sig_arr(ssins, new_sig_arr, new_event_arr)
    
    sig_arr = np.concatenate([sig_arr, new_sig_arr], axis=0)
    event_arr = np.concatenate([event_arr, new_event_arr], axis=0)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [12], line 12
     10 ssins_path_match_events = glob.glob(f"{ssins_dir}/*match_events.yml")[0]
     11 ssins = INS(ssins_path_data, match_events_file=ssins_path_match_events)
---> 12 init_ssins += ssins
     14 new_sig_arr = np.ma.copy(ssins.metric_ms)
     15 new_event_arr = np.ma.copy(new_sig_arr)

File ~/mambaforge/envs/RTP/lib/python3.10/site-packages/pyuvdata/uvflag/uvflag.py:1922, in UVFlag.__iadd__(self, other, axis, run_check, check_extra, run_check_acceptability)
   1896 def __iadd__(
   1897     self,
   1898     other,
   (...)
   1902     run_check_acceptability=True,
   1903 ):
   1904     """In place add.
   1905 
   1906     Parameters
   (...)
   1920 
   1921     """
-> 1922     self.__add__(
   1923         other,
   1924         inplace=True,
   1925         axis=axis,
   1926         run_check=run_check,
   1927         check_extra=check_extra,
   1928         run_check_acceptability=run_check_acceptability,
   1929     )
   1930     return self

File ~/mambaforge/envs/RTP/lib/python3.10/site-packages/SSINS/incoherent_noise_spectrum.py:523, in INS.__add__(self, other, inplace, axis, run_check, check_extra, run_check_acceptability)
    520 mask_uvf_this = this._make_mask_copy()
    521 mask_uvf_other = other._make_mask_copy()
--> 523 mask_uvf = super(INS, mask_uvf_this).__add__(mask_uvf_other,
    524                                              inplace=False,
    525                                              axis=axis,
    526                                              run_check=run_check,
    527                                              check_extra=check_extra,
    528                                              run_check_acceptability=run_check_acceptability)
    530 if inplace:
    531     super(INS, this).__add__(other, inplace=True, axis=axis,
    532                              run_check=run_check,
    533                              check_extra=check_extra,
    534                              run_check_acceptability=run_check_acceptability)

File ~/mambaforge/envs/RTP/lib/python3.10/site-packages/pyuvdata/uvflag/uvflag.py:1890, in UVFlag.__add__(self, other, inplace, axis, run_check, check_extra, run_check_acceptability)
   1887 this.Ntimes = np.unique(this.time_array).size
   1889 if run_check:
-> 1890     this.check(
   1891         check_extra=check_extra, run_check_acceptability=run_check_acceptability
   1892     )
   1893 if not inplace:
   1894     return this

File ~/mambaforge/envs/RTP/lib/python3.10/site-packages/pyuvdata/uvbase.py:587, in UVBase.check(self, check_extra, run_check_acceptability, ignore_requirements)
    582 else:
    583     # Check the shape of the parameter value. Note that np.shape
    584     # returns an empty tuple for single numbers.
    585     # eshape should do the same.
    586     if not np.shape(param.value) == eshape:
--> 587         raise ValueError(
    588             "UVParameter {param} is not expected shape. "
    589             "Parameter shape is {pshape}, expected shape is "
    590             "{eshape}.".format(
    591                 param=p, pshape=np.shape(param.value), eshape=eshape
    592             )
    593         )
    594     # Quantity objects complicate things slightly
    595     # Do a separate check with warnings until a quantity based
    596     # parameter value is created
    597     if isinstance(param.value, Quantity):
    598         # check if user put expected type as a type of quantity
    599         # not a more generic type of number.

ValueError: UVParameter _flag_array is not expected shape. Parameter shape is (3458, 1536, 4), expected shape is (3457, 1536, 4).
In [13]:
plt.figure(figsize=(16,12))
    
cmap = cm.plasma
cmap.set_bad('white')

plt.imshow(init_ssins.metric_array[:,:,0], aspect='auto', cmap=cmap, interpolation='none', 
           vmax=4e4, 
           extent=[ssins_uvf.freq_array[0] / 1e6, ssins_uvf.freq_array[-1] / 1e6, 
                   ssins_uvf.time_array[-1] - JD, ssins_uvf.time_array[0] - JD])
plt.colorbar(pad=.07, extend='neither',
             label='SSINS (Corr. Units)')
plt.title('SSINS')
plt.xlabel('Frequency (MHz)')
plt.ylabel(f'JD - {JD}')
ax2 = plt.gca().twinx()
ax2.set_ylim([ssins_uvf.lst_array[0] * 12 / np.pi, ssins_uvf.lst_array[-1] * 12 / np.pi])
ax2.invert_yaxis()
ax2.set_ylabel('LST (hours)')
ax3 = plt.gca().twiny()
ax3.set_xlim([0, ssins_uvf.Nfreqs - 1])
ax3.set_xlabel('Channel');
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [13], line 8
      3 cmap = cm.plasma
      4 cmap.set_bad('white')
      6 plt.imshow(init_ssins.metric_array[:,:,0], aspect='auto', cmap=cmap, interpolation='none', 
      7            vmax=4e4, 
----> 8            extent=[ssins_uvf.freq_array[0] / 1e6, ssins_uvf.freq_array[-1] / 1e6, 
      9                    ssins_uvf.time_array[-1] - JD, ssins_uvf.time_array[0] - JD])
     10 plt.colorbar(pad=.07, extend='neither',
     11              label='SSINS (Corr. Units)')
     12 plt.title('SSINS')

NameError: name 'ssins_uvf' is not defined
<Figure size 1600x1200 with 0 Axes>

Figure 2(b): SSINS XX Waterfall¶

In [14]:
fig, ax = plt.subplots(figsize=(16,12), nrows=2)
max_abs = 100
if np.max(sig_arr) > max_abs:
    extend = 'max'
    if np.min(sig_arr) < -max_abs:
        extend = 'both'
elif np.min(sig_arr) < -max_abs:
    extend = 'min'    
else:
    extend = 'neither'
    
cmap = cm.RdBu_r
cmap.set_bad('black')

cax = [None, None]

cax[0] = ax[0].imshow(sig_arr[:,:,0], aspect='auto', cmap=cmap, interpolation='none',
                      norm=colors.SymLogNorm(linthresh=1,vmin=-max_abs, vmax=max_abs), 
                      extent=[ssins_uvf.freq_array[0] / 1e6, ssins_uvf.freq_array[-1] / 1e6, 
                              ssins_uvf.time_array[-1] - JD, ssins_uvf.time_array[0] - JD])

cax[1] = ax[1].imshow(event_arr[:,:,0], aspect='auto', cmap=cmap, interpolation='none',
                      norm=colors.SymLogNorm(linthresh=1,vmin=-max_abs, vmax=max_abs), 
                      extent=[ssins_uvf.freq_array[0] / 1e6, ssins_uvf.freq_array[-1] / 1e6, 
                              ssins_uvf.time_array[-1] - JD, ssins_uvf.time_array[0] - JD])
fig.colorbar(cax[0], pad=.07, extend=extend, ax=ax[0],
             label='SSINS Detection Significance (Sample) [$\sigma$s]')
fig.colorbar(cax[1], pad=.07, extend=extend, ax=ax[1],
             label='SSINS Detection Significance (Event) [$\sigma$s]')
fig.suptitle('SSINS Detection Significance')
ax[0].set_xlabel('Frequency (MHz)')
ax[0].set_ylabel(f'JD - {JD}')

ax[1].set_xlabel('Frequency (MHz)')
ax[1].set_ylabel(f'JD - {JD}')

ax2 = ax[0].twinx()
ax2.set_ylim([ssins_uvf.lst_array[0] * 12 / np.pi, ssins_uvf.lst_array[-1] * 12 / np.pi])
ax2.invert_yaxis()
ax2.set_ylabel('LST (hours)')
ax3 = ax[0].twiny()
ax3.set_xlim([0, ssins_uvf.Nfreqs - 1])
ax3.set_xlabel('Channel');

ax4 = ax[1].twinx()
ax4.set_ylim([ssins_uvf.lst_array[0] * 12 / np.pi, ssins_uvf.lst_array[-1] * 12 / np.pi])
ax4.invert_yaxis()
ax4.set_ylabel('LST (hours)')
ax5 = ax[1].twiny()
ax5.set_xlim([0, ssins_uvf.Nfreqs - 1])
ax5.set_xlabel('Channel');
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [14], line 19
     13 cmap.set_bad('black')
     15 cax = [None, None]
     17 cax[0] = ax[0].imshow(sig_arr[:,:,0], aspect='auto', cmap=cmap, interpolation='none',
     18                       norm=colors.SymLogNorm(linthresh=1,vmin=-max_abs, vmax=max_abs), 
---> 19                       extent=[ssins_uvf.freq_array[0] / 1e6, ssins_uvf.freq_array[-1] / 1e6, 
     20                               ssins_uvf.time_array[-1] - JD, ssins_uvf.time_array[0] - JD])
     22 cax[1] = ax[1].imshow(event_arr[:,:,0], aspect='auto', cmap=cmap, interpolation='none',
     23                       norm=colors.SymLogNorm(linthresh=1,vmin=-max_abs, vmax=max_abs), 
     24                       extent=[ssins_uvf.freq_array[0] / 1e6, ssins_uvf.freq_array[-1] / 1e6, 
     25                               ssins_uvf.time_array[-1] - JD, ssins_uvf.time_array[0] - JD])
     26 fig.colorbar(cax[0], pad=.07, extend=extend, ax=ax[0],
     27              label='SSINS Detection Significance (Sample) [$\sigma$s]')

NameError: name 'ssins_uvf' is not defined

Figure 2(c) SSINS XX detection significance¶

Shamelessly copied/pasted code from 2(a). Could write a function in the future. Slightly wrong time axis due to missing integration around chunk boundary in SSINS.

In [15]:
# Load in the flags from each round of XRFI flagging
round_1_flag_files = [
    glob.glob(os.path.join(xrfi_dir, "*.flags1.h5"))[0]
    for xrfi_dir in xrfi_dirs
]
round_2_flag_files = [f.replace("flags1", "flags2") for f in round_1_flag_files]
thresh_flag_file = f'{data_path}/zen.{JD}.total_stage_1_threshold_flags.h5'
round_1_uvf = UVFlag(round_1_flag_files)
round_2_uvf = UVFlag(round_2_flag_files)
round_3_uvf = UVFlag(thresh_flag_file)
In [16]:
# Load in the data.
round_1_flags = round_1_uvf.flag_array[...,0]
round_2_flags = round_2_uvf.flag_array[...,0]
round_3_flags = round_3_uvf.flag_array[...,0]

# For plotting convenience.
flags = dict(zip(range(1,4), (round_1_flags, round_2_flags, round_3_flags)))
unique_flags = {
    1: round_1_flags,
    2: round_2_flags & ~round_1_flags,
    3: round_3_flags & ~round_1_flags & ~round_2_flags,
}

# Construct an array that can be color-coded by when flags were introduced.
combined_flags = np.zeros(round_1_flags.shape, dtype=np.float)
for round_, flags_ in unique_flags.items():
    combined_flags[flags_] = round_ / len(unique_flags)
    
# Prepare different plot labels.
flag_labels = (
    "Flagged in Median Filter Round",
    "Flagged in Mean Filter Round",
    "Flagged in Thresholding",
)
cbar_ticklabels = ("Unflagged",) + tuple(
    "in\n".join(flag_label.split("in "))
    for flag_label in flag_labels
)

# Prepare a colormap.
cmap = plt.cm.colors.ListedColormap(
    ["slategray", "darkturquoise", "plum", "lemonchiffon"]
)

# Useful plot metadata.
lsts = np.unique(round_3_uvf.lst_array)
times = np.unique(round_3_uvf.time_array)
freqs = np.unique(round_3_uvf.freq_array)
chans = np.arange(freqs.size)
lsts_hr = lsts * units.day.to("hr") * units.rad.to("cycle")
freqs_MHz = freqs / 1e6
plot_times = times - float(JD)
extent = (freqs_MHz[0], freqs_MHz[-1], plot_times[-1], plot_times[0])
`np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
In [17]:
# Make a waterfall showing different flagging products.
fig = plt.figure(figsize=(16,12))
ax = fig.add_subplot(111)
ax.set_xlabel("Frequency (MHz)", fontsize=12)
ax.set_ylabel(f"JD - {JD}", fontsize=12)
ax.set_xlim(*extent[:2])
ax.set_ylim(*extent[2:])
cax = ax.imshow(
    combined_flags,
    aspect="auto",
    extent=extent,
    interpolation='none',
    cmap=cmap,
)

# Add labels on the top and right axes.
twinx = ax.twinx()
twiny = ax.twiny()
twinx.set_ylabel("LST (hours)", fontsize=12)
twinx.set_ylim(lsts_hr.max(), lsts_hr.min())
twiny.set_xlabel("Channel", fontsize=12)
twiny.set_xlim(chans.min(), chans.max())

# Configure the colorbar so that labels are at the center of each section.
cbar = fig.colorbar(cax)
cbar_ticks = np.linspace(combined_flags.min(), combined_flags.max(), len(flags) + 2)
cbar_ticks = 0.5 * (cbar_ticks[1:] + cbar_ticks[:-1])
cbar.set_ticks(cbar_ticks)
cbar.set_ticklabels(cbar_ticklabels);
In [18]:
plt.figure(figsize=(16,12))
max_abs = 100
if np.max(uvf2.metric_array) > max_abs:
    extend = 'max'
    if np.min(uvf2.metric_array) < -max_abs:
        extend = 'both'
elif np.min(uvf2.metric_array) < -max_abs:
    extend = 'min'    
else:
    extend = 'neither'

plt.imshow(uvf2.metric_array[:,:,0], aspect='auto', cmap='RdBu_r',
           norm=colors.SymLogNorm(linthresh=1,vmin=-max_abs, vmax=max_abs), 
           extent=[uvf.freq_array[0] / 1e6, uvf.freq_array[-1] / 1e6, 
                   uvf.time_array[-1] - JD, uvf.time_array[0] - JD])
plt.colorbar(pad=.07, extend=extend,
             label='RFI Detection Significance ($\sigma$s)')
plt.title('Combined XRFI Metrics')
plt.xlabel('Frequency (MHz)')
plt.ylabel(f'JD - {JD}')
ax2 = plt.gca().twinx()
ax2.set_ylim([uvf.lst_array[0] * 12 / np.pi, uvf.lst_array[-1] * 12 / np.pi])
ax2.invert_yaxis()
ax2.set_ylabel('LST (hours)')
ax3 = plt.gca().twiny()
ax3.set_xlim([0, uvf.Nfreqs - 1])
ax3.set_xlabel('Channel');

Figure 3: Flag Evolution¶

This figure shows how the flags are built at each step in the initial XRFI flagging pipeline. Note that the completely flagged sections at the beginning and end of the night are expected. Main thing to look out for is if there are channels which are highly flagged after the second round of flagging but have not been completely flagged after day thresholding.

In [19]:
# Collapse the flags along each axis and plot the result.
fig = plt.figure(figsize=(15,10))
axes = fig.subplots(2, gridspec_kw={"hspace": 0.35})
twin_axes = [ax.twiny() for ax in axes]

# Set the plot labels.
axes[0].set_xlabel(f"JD - {JD}", fontsize=12)
axes[0].set_ylabel("Fraction of Channels Flagged", fontsize=12)
axes[1].set_xlabel("Frequency (MHz)", fontsize=12)
axes[1].set_ylabel(
    "Fraction of Integrations Flagged\nin Total Day Thresholded Flags", fontsize=12
)
twin_axes[0].set_xlabel("LST (hour)", fontsize=12)
twin_axes[0].set_xlim(lsts_hr[0], lsts_hr[-1])
twin_axes[1].set_xlabel("Channel", fontsize=12)
twin_axes[1].set_xlim(chans.min(), chans.max())

# Plot the channel occupancy as a function of time.
for label, flag_array in zip(flag_labels, flags.values()):
    axes[0].plot(plot_times, flag_array.astype(np.float).mean(axis=1), label=label)
axes[0].plot(plot_times, ssins_uvf.flag_array.astype(np.float).mean(axis=(1, -1)), 
             label="Flagged by SSINS")

# Plot the flagging fraction as a function of frequency.
axes[1].plot(
    freqs_MHz, round_3_flags.astype(np.float).mean(axis=0), color="k", ms=1, lw=0, marker="o",
    label="XRFI frequency occupancy",
)
# Do the same with the SSINS flags
axes[1].plot(freqs_MHz, ssins_uvf.flag_array.astype(np.float).mean(axis=(0, -1)), ms=1, lw=0,
             marker="o", label="SSINS frequency occupancy")
    
axes[0].legend();
axes[1].legend()
`np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [19], line 21
     19 for label, flag_array in zip(flag_labels, flags.values()):
     20     axes[0].plot(plot_times, flag_array.astype(np.float).mean(axis=1), label=label)
---> 21 axes[0].plot(plot_times, ssins_uvf.flag_array.astype(np.float).mean(axis=(1, -1)), 
     22              label="Flagged by SSINS")
     24 # Plot the flagging fraction as a function of frequency.
     25 axes[1].plot(
     26     freqs_MHz, round_3_flags.astype(np.float).mean(axis=0), color="k", ms=1, lw=0, marker="o",
     27     label="XRFI frequency occupancy",
     28 )

NameError: name 'ssins_uvf' is not defined

Figure 4: Flagging Occupancies¶

The top plot shows the fraction of channels flagged at each integration for each set of flags. The bottom plot shows the fraction of integrations flagged as a function of frequency for the total thresholded flags.

Metadata¶

In [20]:
from hera_qm import __version__
print(__version__)
2.0.5.dev11+g87299d5
In [21]:
print(f"SSINS version info: {SSINS_ver.construct_version_info()}")
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In [21], line 1
----> 1 print(f"SSINS version info: {SSINS_ver.construct_version_info()}")

AttributeError: module 'SSINS.version' has no attribute 'construct_version_info'