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 = 2459965
data_path = "/mnt/sn1/2459965"
In [3]:
from astropy.time import Time
utc = Time(JD, format='jd').datetime
print(f'Date: {utc.month}-{utc.day}-{utc.year}')

Date: 1-20-2023
In [4]:
uvf = UVFlag(f'{data_path}/zen.{JD}.total_stage_1_threshold_flags.h5')

---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
Cell In [4], line 1
----> 1 uvf = UVFlag(f'{data_path}/zen.{JD}.total_stage_1_threshold_flags.h5')

File ~/mambaforge/envs/RTP/lib/python3.10/site-packages/pyuvdata/uvflag/uvflag.py:556, in UVFlag.__init__(self, indata, mode, copy_flags, waterfall, history, label, use_future_array_shapes, run_check, check_extra, run_check_acceptability)
    552         del fobj
    554 elif issubclass(indata.__class__, (str, pathlib.Path)):
    555     # Given a path, read indata
--> 556     self.read(
    557         indata,
    558         history,
    559         use_future_array_shapes=use_future_array_shapes,
    560         run_check=run_check,
    561         check_extra=check_extra,
    562         run_check_acceptability=run_check_acceptability,
    563     )
    564 elif issubclass(indata.__class__, UVData):
    565     self.from_uvdata(
    566         indata,
    567         mode=mode,
   (...)
    575         run_check_acceptability=run_check_acceptability,
    576     )

File ~/mambaforge/envs/RTP/lib/python3.10/site-packages/pyuvdata/uvflag/uvflag.py:2783, in UVFlag.read(self, filename, history, use_future_array_shapes, run_check, check_extra, run_check_acceptability)
   2781 else:
   2782     if not os.path.exists(filename):
-> 2783         raise IOError(filename + " not found.")
   2785     # update filename attribute
   2786     basename = os.path.basename(filename)

OSError: /mnt/sn1/2459965/zen.2459965.total_stage_1_threshold_flags.h5 not found.
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');

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [5], line 2
      1 plt.figure(figsize=(16,12))
----> 2 plt.imshow(uvf.flag_array[:,:,0], aspect='auto', interpolation='none',
      3            extent=[uvf.freq_array[0] / 1e6, uvf.freq_array[-1] / 1e6, 
      4                    uvf.time_array[-1] - JD, uvf.time_array[0] - JD])
      5 plt.xlabel('Frequency (MHz)')
      6 plt.ylabel(f'JD - {JD}')

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

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 2
      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

NameError: name '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 1850 directories containing XRFI intermediate data products.
Found 1850 combined round 1 XRFI metrics files.
Found 1850 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');

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [10], line 14
      9 else:
     10     extend = 'neither'
     12 plt.imshow(uvf2.metric_array[:,:,0], aspect='auto', cmap='RdBu_r', interpolation='none', 
     13            norm=colors.SymLogNorm(linthresh=1,vmin=-max_abs, vmax=max_abs), 
---> 14            extent=[uvf.freq_array[0] / 1e6, uvf.freq_array[-1] / 1e6, 
     15                    uvf.time_array[-1] - JD, uvf.time_array[0] - JD])
     16 plt.colorbar(pad=.07, extend=extend,
     17              label='RFI Detection Significance ($\sigma$s)')
     18 plt.title('Combined XRFI Metrics')

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

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)
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)

---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
Cell In [15], line 10
      8 round_1_uvf = UVFlag(round_1_flag_files)
      9 round_2_uvf = UVFlag(round_2_flag_files)
---> 10 round_3_uvf = UVFlag(thresh_flag_file)

File ~/mambaforge/envs/RTP/lib/python3.10/site-packages/pyuvdata/uvflag/uvflag.py:556, in UVFlag.__init__(self, indata, mode, copy_flags, waterfall, history, label, use_future_array_shapes, run_check, check_extra, run_check_acceptability)
    552         del fobj
    554 elif issubclass(indata.__class__, (str, pathlib.Path)):
    555     # Given a path, read indata
--> 556     self.read(
    557         indata,
    558         history,
    559         use_future_array_shapes=use_future_array_shapes,
    560         run_check=run_check,
    561         check_extra=check_extra,
    562         run_check_acceptability=run_check_acceptability,
    563     )
    564 elif issubclass(indata.__class__, UVData):
    565     self.from_uvdata(
    566         indata,
    567         mode=mode,
   (...)
    575         run_check_acceptability=run_check_acceptability,
    576     )

File ~/mambaforge/envs/RTP/lib/python3.10/site-packages/pyuvdata/uvflag/uvflag.py:2783, in UVFlag.read(self, filename, history, use_future_array_shapes, run_check, check_extra, run_check_acceptability)
   2781 else:
   2782     if not os.path.exists(filename):
-> 2783         raise IOError(filename + " not found.")
   2785     # update filename attribute
   2786     basename = os.path.basename(filename)

OSError: /mnt/sn1/2459965/zen.2459965.total_stage_1_threshold_flags.h5 not found.
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])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [16], line 4
      2 round_1_flags = round_1_uvf.flag_array[...,0]
      3 round_2_flags = round_2_uvf.flag_array[...,0]
----> 4 round_3_flags = round_3_uvf.flag_array[...,0]
      6 # For plotting convenience.
      7 flags = dict(zip(range(1,4), (round_1_flags, round_2_flags, round_3_flags)))

NameError: name 'round_3_uvf' is not defined
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);

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [17], line 6
      4 ax.set_xlabel("Frequency (MHz)", fontsize=12)
      5 ax.set_ylabel(f"JD - {JD}", fontsize=12)
----> 6 ax.set_xlim(*extent[:2])
      7 ax.set_ylim(*extent[2:])
      8 cax = ax.imshow(
      9     combined_flags,
     10     aspect="auto",
   (...)
     13     cmap=cmap,
     14 )

NameError: name 'extent' is not defined
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');

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [18], line 14
      9 else:
     10     extend = 'neither'
     12 plt.imshow(uvf2.metric_array[:,:,0], aspect='auto', cmap='RdBu_r',
     13            norm=colors.SymLogNorm(linthresh=1,vmin=-max_abs, vmax=max_abs), 
---> 14            extent=[uvf.freq_array[0] / 1e6, uvf.freq_array[-1] / 1e6, 
     15                    uvf.time_array[-1] - JD, uvf.time_array[0] - JD])
     16 plt.colorbar(pad=.07, extend=extend,
     17              label='RFI Detection Significance ($\sigma$s)')
     18 plt.title('Combined XRFI Metrics')

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

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()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [19], line 14
     10 axes[1].set_ylabel(
     11     "Fraction of Integrations Flagged\nin Total Day Thresholded Flags", fontsize=12
     12 )
     13 twin_axes[0].set_xlabel("LST (hour)", fontsize=12)
---> 14 twin_axes[0].set_xlim(lsts_hr[0], lsts_hr[-1])
     15 twin_axes[1].set_xlabel("Channel", fontsize=12)
     16 twin_axes[1].set_xlim(chans.min(), chans.max())

NameError: name 'lsts_hr' 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.dev13+gd6c757c
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'