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 = 2459966
data_path = "/mnt/sn1/2459966"
In [3]:
from astropy.time import Time
utc = Time(JD, format='jd').datetime
print(f'Date: {utc.month}-{utc.day}-{utc.year}')
Date: 1-21-2023
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:613, in UVFlag.__init__(self, indata, mode, copy_flags, waterfall, history, label, use_future_array_shapes, run_check, check_extra, run_check_acceptability)
    610 self.label = ""  # Added to at the end
    611 if isinstance(indata, (list, tuple)):
    612     self.__init__(
--> 613         indata[0],
    614         mode=mode,
    615         copy_flags=copy_flags,
    616         waterfall=waterfall,
    617         history=history,
    618         label=label,
    619         use_future_array_shapes=use_future_array_shapes,
    620         run_check=run_check,
    621         check_extra=check_extra,
    622         run_check_acceptability=run_check_acceptability,
    623     )
    624     if len(indata) > 1:
    625         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)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In [12], line 3
      1 init_ssins_path_data = glob.glob(f"{ssins_dirs_sorted[0]}/*data.h5")[0]
      2 init_ssins_path_match_events = glob.glob(f"{ssins_dirs_sorted[0]}/*match_events.yml")[0]
----> 3 init_ssins = INS(init_ssins_path_data, match_events_file=init_ssins_path_match_events)
      4 init_sig_arr = np.ma.copy(init_ssins.metric_ms)
      5 init_event_arr = np.ma.copy(init_sig_arr)

File ~/mambaforge/envs/RTP/lib/python3.10/site-packages/SSINS/incoherent_noise_spectrum.py:44, in INS.__init__(self, input, history, label, order, mask_file, match_events_file, spectrum_type, use_integration_weights, nsample_default)
     21 def __init__(self, input, history='', label='', order=0, mask_file=None,
     22              match_events_file=None, spectrum_type="cross",
     23              use_integration_weights=False, nsample_default=1):
     25     """
     26     init function for the INS class.
     27 
   (...)
     41             the uvfits file.
     42     """
---> 44     super().__init__(input, mode='metric', copy_flags=False,
     45                      waterfall=False, history='', label='')
     47     # Used in _data_params to determine when not to return None
     48     self._super_complete = True

File ~/mambaforge/envs/RTP/lib/python3.10/site-packages/pyuvdata/uvflag/uvflag.py:648, in UVFlag.__init__(self, indata, mode, copy_flags, waterfall, history, label, use_future_array_shapes, run_check, check_extra, run_check_acceptability)
    644         del fobj
    646 elif issubclass(indata.__class__, (str, pathlib.Path)):
    647     # Given a path, read indata
--> 648     self.read(
    649         indata,
    650         history,
    651         use_future_array_shapes=use_future_array_shapes,
    652         run_check=run_check,
    653         check_extra=check_extra,
    654         run_check_acceptability=run_check_acceptability,
    655     )
    656 elif issubclass(indata.__class__, UVData):
    657     self.from_uvdata(
    658         indata,
    659         mode=mode,
   (...)
    667         run_check_acceptability=run_check_acceptability,
    668     )

File ~/mambaforge/envs/RTP/lib/python3.10/site-packages/pyuvdata/uvflag/uvflag.py:3314, in UVFlag.read(self, filename, history, use_future_array_shapes, run_check, check_extra, run_check_acceptability)
   3292 """Read in flag/metric data from a HDF5 file.
   3293 
   3294 Parameters
   (...)
   3311 
   3312 """
   3313 # make sure we have an empty object.
-> 3314 self.__init__()
   3315 if isinstance(filename, (tuple, list)):
   3316     self.read(filename[0])

TypeError: INS.__init__() missing 1 required positional argument: 'input'
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 6
      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 'init_ssins' 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 3
      1 fig, ax = plt.subplots(figsize=(16,12), nrows=2)
      2 max_abs = 100
----> 3 if np.max(sig_arr) > max_abs:
      4     extend = 'max'
      5     if np.min(sig_arr) < -max_abs:

NameError: name 'sig_arr' 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.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'