Calibrated, Flagged, and LST-Binned HERA Internal Data Release
On behalf of the HERA Analysis Team
Introduction¶
H1C IDR 2.2 is a re-release of 18 days of calibrated data products from late 2017, part of HERA's first science observing season. It follows up on IDR 2.1, our first internal data release. That release was described in this memo at reionization.org. This memo builds on that memo, and though the overlap is substantial, ID2.2 has a number of improvements:
- File format change: We've moved form using miriad files to uvh5 files. This enables improved partial I/O.
- Python 3 support: All modules (at least for now) support both Python 2 and 3 and the pipeline was run on Python 3.7.
- 4-pol data files, 2-pol calibration files: All data files now contain all four instrumental polarizations, all calibration files now contain both antenna polarizations.
- Bad antennas included in calibration files: Flagged antennas are now included in calibration files, through their flags are all set to
True. Previously, they were simply excluded. - Fewer visibility files saved: With the exception of delay-filtered data, we no longer produce any of the day-by-day calibrated/flagged data products that can be trivially reconstructed from calibration files and raw data. This saves considerable disk space. It is not true of LST-binned data.
- Improved redundant calibration: Redundant calibration is now performed with an generalized firstcal, a new pure-python
omnical, and proper normalization of $\chi^2 / DoF$. The visibility solutions fromomnicalare updated to include absolute calibration in flagging later in the pipeline. - Improved absolute calibration: The externally calibrated set of "model" visibilities has been improved by including more sources and ignoring short baselines (under 40 m). Also, abscal calibration files now contain their own $\chi^2$ and $\chi^2$ per antenna, which compare calibrated data to the abscal model.
- Improved RFI flagging: RFI flagging has been systematized and rationalized to produce a single flag waterfall (per time and frequency, but not per antenna or baseline) using metrics of "outlierness" from both raw data and calibration data products, combined in quadrature. These metrics are also used to pick entire channels or integrations to remove, based on the statistics of a whole day.
- Improved calibration smoothing: Instead of smoothing in time and then in frequency, we now perform a 2D smoothing using an iterative CLEAN-like algorithm.
- Various new ancillary day-by-data data products:
- Multi-frequency synthesis imaging using CASA
- Fits to cable reflections
- Estimates of the noise on calibrated data
- Extracted raw autocorrelations
- Absolutely-calibrated omnical visibility solutions
- UPDATE 12/3/19: This notebook has been updated to run on the most recent version of
hera_cal, which now supports polarizations denoted by cardinal directions ('ee','nn', etc.) in favor of'xx'and'yy'polarizations whenx_orientationis defined in the data files.
This jupyter notebook memo is designed to be run when logged into NRAO. It is available as a .ipynb via github. In addition to the standard python libraries, it requires pyuvdata, hera_cal, and hera_qm. For more information on access, see the HERA wiki. For ongoing discussions, join the #hera-analysis Slack channel or our
telecons on Wednesdays at 10am Pacific in the HERA Zoom room.
from __future__ import print_function, division, absolute_import
import numpy as np
The day-by-day data, along with the pipeline settings file used (makeflow/idr2_2.toml), and the LST-binned data can be found in:
!ls /lustre/aoc/projects/hera/H1C_IDR2/IDR2_2
analysis_folder = '/lustre/aoc/projects/hera/H1C_IDR2/IDR2_2'
Finding and Loading Data¶
H1C IDR 2.2 includes 18 nearly-consecutive nights of data. All data products are sorted by JD into the following folders on the NRAO servers, which also contain softlinks to raw H1C IDR 2 data in the .uvh5 format and antenna metrics released by the commissioning team. Each night has 73 raw visibility files. Each file has 4 polarizations, 1024 frequency channels and (usually) 60 integrations, each 10.7374 seconds.
More information about the data products blessed by the commissioning team team are on the HERA wiki.
import glob, os
for night in sorted(glob.glob(os.path.join(analysis_folder,'245*'))):
print(night, 'has', len(glob.glob(night + '/*HH.uvh5')), 'raw visibility files.')
The HERA collaboration uses pyuvdata (for more, see pyuvdata's documentation) as a standard pythonic interface between data and calibration file formats and in-memory objects. The HERA Analysis team uses subclasses of pyuvdata objects with a number of extra capabilities, namely HERAData and HERACal objects.
For a more complete demonstration of the capabilities of the hera_cal.io module, refer to this example notebook.
Raw HERA visibility files for IDR2 blessed by the commissioning team have the format zen.2458098.43124.HH.uvh5 which can be interpreted as:
zen: zenith pointing2458098: Julian date integer43124: first five decimal digits of the time of the first integration in the fileHH: "HERA Hex", visibilities purely between HERA core antennasuvh5: file extension foruvh5files, a hdf5-based,pyuvdata-compliant data format for visibilities
bl = (25, 51, 'ee')
sample_data_file = os.path.join(analysis_folder, '2458098/zen.2458098.43124.HH.uvh5')
# this is the baseline and file we examine through this notebook
from hera_cal.io import HERAData
hd = HERAData(sample_data_file)
print('This file has', len(hd.times), 'integrations', 'and', len(hd.freqs), 'frequency channels.')
HERAData objects (when initialized with uvh5 files) automatically load relevant metadata and store in the object, including times, freqs, antpos, bls, pols, and lsts.
data, flags, nsamples = hd.read(bls=[bl]) # only loads a single bl; default loads all bls
The data, flags, and nsamples loaded through hera_cal.io are stored in DataContainer objects, which act like dictionaries mapping baseline tuples in the (ant1,ant2,pol) format like (65,71,'ee') to waterfalls of shape (Ntimes, Nfreqs). DataContainers support getting data and flags via [] and normal dictionary functions like .keys(), .values(), .items(), .has_key(). They also know to conjugate data when the reversed baseline key is provided and they abstract away polarization capitalization convetions:
print(np.all(data[25, 51, 'ee'] == np.conj(data[51, 25, 'ee'])))
print(np.all(data[25, 51, 'EE'] == data[25, 51, 'ee']))
This then allows for easy data access and plotting, such as:
from matplotlib import pyplot as plt
%matplotlib inline
plt.figure(figsize=(12,3), dpi=100)
plt.imshow(np.angle(data[bl]) / (~flags[bl]), aspect='auto', cmap='twilight',
extent=[hd.freqs[0]/1e6, hd.freqs[-1]/1e6, hd.times[-1]-2458098, hd.times[0]-2458098])
plt.ylabel('JD - 2458098')
plt.xlabel('Frequency (MHz)')
plt.title('Raw Visibilities: ' + str(bl))
plt.colorbar(label='Phase (Radians)', aspect=8, pad=.025);
The Array During IDR2 and Bad Antennas¶
The data taken for IDR2.2 began on 2458098 and ended on 2458116 with 52 antennas (~40 working well). Bad antennas, are summarized in this folder:
bad_ants_folder = '/users/heramgr/hera_software/hera_opm/pipelines/h1c/idr2/v2/bad_ants/'
bad_ants = np.loadtxt(os.path.join(bad_ants_folder, '2458098.txt')).astype(int)
print('Bad antennas on 2458098:', bad_ants)
These bad antennas be include those identified by the commissioning team as bad as well as several that were found to be particularly non-redundant in an earlier round of analysis. We can visualize the array on 2458098 as follows:
plt.figure(figsize=(5,4), dpi=100)
plt.scatter(np.array(list(hd.antpos.values()))[:,0], np.array(list(hd.antpos.values()))[:,1],
c=[ant in bad_ants for ant in hd.antpos.keys()], s=400, cmap='jet')
for ant,pos in hd.antpos.items():
plt.text(pos[0],pos[1],str(ant), va='center', ha='center', color='w')
plt.xlabel("East-West [meters]")
plt.ylabel("North-South [meters]")
plt.title('Antenna Positions on 2458098 (Red = Flagged)');
plt.axis('equal')
plt.tight_layout()
All data in this IDR come from antennas in the Southwest sector of the split HERA core.
Additionally, we inherit antenna metrics from the commissioning team. These come in json files with the format zen.2458098.43124.HH.uv.ant_metrics.json. These can be read hera_qm.ant_metrics.load_antenna_metrics().
The list of bad antennas generated by ant_metrics is generally a subset of the flagged antennas, but both are ORed together on a file-by-file basis. These antennas are left in all data/calibration files, but their flag arrays are all True.
Pipeline Overview¶
The IDR 2.2 pipeline was run using hera_opm, which interfaces with the Makeflow Workflow System. The configuration files with the specific parameters used, associated bash scripts, and other files related bad antennas can all be found here.
Below is a flowchart (which also lives in hera_opm) detailing all the analysis steps in IDR 2.2, including the interdependencies of the pipeline and the outputs at various steps.
from IPython.display import IFrame
IFrame('https://nbviewer.jupyter.org/github/HERA-Team/hera_opm/blob/' +
'master/pipelines/h1c/idr2/v2/Analysis_Flowchart.pdf', 980, 1650)
Mainline Day-by-Day Data Products¶
The mainline daily analysis pipeline includes several steps of calibration, RFI-flagging, and an optional delay-filtering step.
Redcal: Redundant-Baseline Calibration¶
Redundant-baseline calibration solves for antenna gains and visibility solutions in order to minimize the difference between the observations and the model visibilities, expressed in terms of $\chi^2$:
$\chi^2(t,\nu) = \sum_{i,j} |V_{ij}^\text{obs}(t,\nu) - g_i(t,\nu) g_j^*(t,\nu) V_{i-j}^\text{sol}(t,\nu)|^2 / \sigma^2_{ij}$.
The idea of redundant calibration is to solve for $g_i$ and $V_{ij}^\text{sol}$ simultaneously, taking advantage of the fact that there are far more total visibility measurements than unique baselines. This process takes place in three steps, firstcal, logcal, and lincal which iteratively searches for the minimum value of $\chi^2$ by first restricting degrees of freedom (firstcal finds one delay and one phase offset per antenna), finding an approximate per-frequency calibration solution with logcal, and then refining that solution with lincal.
This procedure operates on ee and nn visibilities independently, cross-polarized visibilities are not used. Each time, frequency, and polarization has four degeneracies that cannot be solved for with redundant calibration (because they do not affect $\chi^2$):
- Overall amplitude
- Overall phase
- North-South phase gradient
- East-West phase gradient
For a more detailed, pedagogical explanation of the degeneracies of redundant-baseline calibration, see Section 2 of this paper.. A thorough overview of the results from redundant calibration of HERA data in IDR 2.2 is forthcoming in Dillon et al. (in prep.).
New in IDR 2.2:
firstcalhas been rewritten to use thelinsolvepackage and to iteratively solve for phase offsets (as opposed to just delays).lincalis no longer uses theomnicalpackage but rather a pure python implementation of it inhera_cal.redcal, which is both faster and more accurate (see HERA Memo #50).- $\chi^2 / DoF$ is now properly normalized (see HERA Memo #61).
- The omnical visibility solutions are subsequently updated to include absolute calibration and flagging. See below.
As the above diagram indicates, the output of redundant calibration includes firstcal calibration solutions, omnical calibration solutions, and omnical visibility solutions for unique baselines. The latter is output as a .uvh5 file, albeit with fewer baselines than the full data (only one per unique separation). The first two are saved as .calfits files. Calibration solution in .calfits files can be read with pyuvdata or using the HERACal subclass in an analogous way to HERAData. HERACal's .read() function returns gains, flags, and $\chi^2$ per antenna, and overall $\chi^2$. In the omnical calibration solutions, that $\chi^2$ is a properly normalized $\chi^2/DOF$.
from hera_cal.io import HERACal
from hera_cal.utils import split_bl
ant = split_bl(bl)[0]
# load firstcal gains and flags
firstcal_file = os.path.join(analysis_folder, '2458098/zen.2458098.43124.HH.first.calfits')
hc = HERACal(firstcal_file)
fc_gains, fc_flags, _, _ = hc.read()
# load redundant calibration gains, flags, chi^2 per antenna, and chi^2 / DoF
omnical_file = os.path.join(analysis_folder, '2458098/zen.2458098.43124.HH.omni.calfits')
hc = HERACal(omnical_file)
oc_gains, oc_flags, oc_quals, oc_total_quals = hc.read()
Note that relevant information about the files produced (git hashes, command line invocation, etc.) by the pipeline is available in their histories:
print(hc.history)
def waterfall_column(waterfalls, flags, titles, clims=None, clabels=None, cmaps=None, ylims=None,
ylabel='JD - 2458098', extents=None, hspace=.1, figsize=(12,6), dpi=100):
'''Useful plotting function for the IDR 2.2 memo.'''
if clims is None:
clims = [None for i in range(len(waterfalls))]
if clabels is None:
clabels = [None for i in range(len(waterfalls))]
if cmaps is None:
cmaps = [None for i in range(len(waterfalls))]
if ylims is None:
ylims = [None for i in range(len(waterfalls))]
if not any(isinstance(ex, list) for ex in extents):
extents = [extents for i in range(len(waterfalls))]
fig, axes = plt.subplots(len(waterfalls), 1, sharex=True, squeeze=True, figsize=figsize, dpi=dpi)
plt.subplots_adjust(hspace=hspace)
for ax, wf, f, t, clim, clabel, cmap, ylim, ex in zip(axes, waterfalls, flags, titles,
clims, clabels, cmaps, ylims, extents):
with np.errstate(divide='ignore', invalid='ignore'):
im = ax.imshow(wf / ~f, aspect='auto', extent=ex, cmap=cmap)
plt.colorbar(im, ax=ax, label=clabel, aspect=8, pad=.025)
if ax == axes[-1]:
ax.set_xlabel('Frequency (MHz)')
im.set_clim(clim)
ax.set_ylabel(ylabel)
ax.set_ylim(ylim)
props = dict(boxstyle='round', facecolor='white', alpha=0.8)
ax.text(0.02, 0.9, t, transform=ax.transAxes, fontsize=14, verticalalignment='top', bbox=props)
waterfall_column([np.angle(fc_gains[ant]), np.angle(oc_gains[ant]), oc_total_quals[ant[1]]],
[fc_flags[ant], oc_flags[ant], oc_flags[ant]],
['Firstcal: {}'.format(ant), 'Omnical: {}'.format(ant), 'Omnical $\chi^2$ / DOF'],
clims=[None, None, [0, 3]],
clabels=['Phase (Radians)', 'Phase (Radians)', '$\chi^2$ / DoF (Unitless)'],
cmaps = ['twilight', 'twilight', 'inferno'],
extents=[hc.freqs[0] / 1e6, hc.freqs[-1] / 1e6, hc.times[-1] - 2458098, hc.times[0] - 2458098])
Here we see redundant baseline calibration solutions, which are dominated by an overall delay. We also see that $\chi^2 / DoF$ is generically greater than 1, indicating non-redundancy that seems to vary both in frequency and time. We also see clear evidence of narrow-band RFI which has not yet been excised. The first and last 50 channels were flagged by hand, a few more channels also got flagged at the high end of the band for producing nans in the calibration solution, which usually occurs when there are 0s in the data. At this stage, we also flag any data where the sun is above the horizon, which explains why the IDR contains many completely flagged files. While some of the integrations appear discrepant, especially after firstcal, this turns out to be discontinuities in the degenerate subspace that are subsequently fixed by abscal.
Next, we compare redundantly calibrated data to the visibility solution for all redundant baselines in the same group. To do that, we use hera_cal.apply_cal.calibrate_in_place() to calibrate the data we've already loaded above and we use hera_cal.redcal.get_reds() to convert antenna positions into a list of redundant baseline groups, which are themselves lists of baselines. The first baseline in the group is used as the key in the omnical visibility solution.
# apply calibration solution to data
from copy import deepcopy
from hera_cal.apply_cal import calibrate_in_place
redcal_data, redcal_flags = deepcopy(data), deepcopy(flags)
calibrate_in_place(redcal_data, oc_gains, data_flags=redcal_flags, cal_flags=oc_flags)
# find the mapping of baselines to the unique baseline key used in the omnical visibility solutions
from hera_cal.redcal import get_reds
reds = get_reds(hd.antpos, pols=['ee'])
red_dict = {bl: bl_group[0] for bl_group in reds for bl in bl_group}
omnical_vis_file = os.path.join(analysis_folder, '2458098/zen.2458098.43124.HH.omni_vis.uvh5')
hd_oc = HERAData(omnical_vis_file)
omnivis, omnivis_flags, _ = hd_oc.read(bls=[red_dict[bl]]) # note, this only has the first baseline in each redundant group
waterfall_column([np.angle(redcal_data[bl]), np.angle(omnivis[red_dict[bl]])],
[redcal_flags[bl], omnivis_flags[red_dict[bl]]],
['Redundantly Calibrated: {}'.format(bl), 'Omnical Visibility Solution: {}'.format(bl)],
clabels=['Phase (Radians)', 'Phase (Radians)'],
cmaps = ['twilight', 'twilight'],
extents=[hd.freqs[0] / 1e6, hd.freqs[-1] / 1e6, hd.times[-1] - 2458098, hd.times[0] - 2458098])
Much of the difference between the individual calibrated visibility and the visibility solution is noise; the visibility solution is approximately equal to the average of all redundant baselines in a group.
Abscal: Fixing Redcal Degeneracies with Externally Calibrated Visibilities¶
The next step in calibration is absolute calibration or abscal. Absolute calibration attaches the omnical gain solutions to a sky reference, allowing us to fix the degeneracies. In fact, absolute calibration is restricted to only modify the degenerate subspace of the omnical solutions; only four degrees of freedom per polarization per time per frequency are allowed to be modified.
Instead of using a point source model, abscal takes as its input a set of externally calibrated visibilities which we call the "abscal model," though simulated visibilities could also work. Individual integrations in the model are phased to match the nearest integration in the data. As in redundant calibration, abscal first calibrates a delay and phase slope across the array (a subset of the degeneracies) before solving for the degeneracies directly, namely a per-frequency and per-time phase slope and amplitude.
The overall phase degeneracy cannot be solved by abscal because it has no impact on calibrated visibilities. To resolve it, we pick a reference antenna and rephase all other gains so that it has zero phase for all frequencies and times. We choose a reference antenna per-file and per-polarization that 1) has no more flags than any other antenna and 2) introduces the least phase noise into the other gains. Later in the calibration smoothing step, we will rephase to a single reference antenna for the whole night using the same criteria.
Because polarizations are calibrated separately without reference to any $en$- or $ne$-visibilities, we have no reason to expect that the cross-polarized visibilties are properly calibrated. This does not affect our analysis of psuedo-Stokes I or Q images or power spectra, but it does mean that our pseudo-Stokes U and V images and power spectra are not reliable at this point.
New in IDR 2.2:
In IDR 2.2, we added a notion of abscal $\chi^2$, defined as:
$\chi^2_\text{abs}(t,\nu) = \sum_{i,j} |V_{ij}^\text{model}(t,\nu) - g^\text{abs}_i(t,\nu) g^{*\text{abs}}_j(t,\nu) V_{i,j}^\text{raw}(t,\nu)|^2 / \sigma^2_{ij}$,
where $\sigma^2$ is the noise on the absolute calibrated visibilities (as inferred from the autocorrelations). Unlike the omnical $\chi^2$, this is not properly normalized per DoF. However, it is still a useful metric of discrepancy between the data after redcal and abscal and the externally calibrated visibilities.
For IDR2.2, the data was calibrated by looking at three days that span the LST range of what will eventually be H1C IDR3, namely 2458042, 2458116, and 2458161. Each day was calibrated in CASA using a different calibrator field from GLEAM and only baselines > 40 m in length:
- 2458042 was calibrated on a field centered at 2.0167 hours of LST (including GLEAM J0200-3053)
- 2458116 was calibrated on a field centered at 5.2367 hours of LST (including GLEAM J0513-3028)
- 2458161 was calibrated on a field centered at 14.4428 hours of LST (including GLEAM J1425-2959)
The calibration solutions were assumed to be stable over the course of each calibration day. Then, to form one synthetic day externally calibrated visibilities, these days are "cross-faded" into each other, effectively a weighted average of the three days with weights that fall linearly away from the calibration field and are 0 when the sun is up.
Finally, the abscal model was low-pass delay filtered (using a Tukey window with $\alpha=0.2$) at either the baseline's horizon + 50 ns, or 150 ns, whichever is longer. This lowered the noise on the model, removed any high-frequency spectral structure we don't expect to be associated with the sky, and allowed interpolation over RFI-contaminated channels. Unfortunately, this also risks spreading out any isolated but unflagged RFI, which we see some evidence of in the the abscal $\chi^2$ below.
For more on absolute calibration, see HERA Memo #42 for the IDR 2.1 implementation or Kern et al. (in prep.).
# Load abscal calibration solutions
abscal_file = os.path.join(analysis_folder, '2458098/zen.2458098.43124.HH.abs.calfits')
hc = HERACal(abscal_file)
ac_gains, ac_flags, ac_quals, ac_total_quals = hc.read()
waterfall_column([np.abs(ac_gains[(ant)]), np.angle(ac_gains[ant]), ac_total_quals[ant[1]]],
[ac_flags[ant], ac_flags[ant], ac_flags[ant]],
['Abscal Gain Amplitude: {}'.format(ant), 'Abscal Gain Phase: {}'.format(ant), 'Abscal $\chi^2$'],
clims=[[0, .04], None, [0, 5]],
clabels=['Amplitude (Unitless)', 'Phase (Radians)', '$\chi^2$ (Unnormalized)'],
cmaps = ['inferno', 'twilight', 'inferno'],
extents=[hc.freqs[0] / 1e6, hc.freqs[-1] / 1e6, hc.times[-1] - 2458098, hc.times[0] - 2458098])
We can also see that the absolutely calibrated data now closely match the reference at the same LSTs.
# Apply abscal calibration to data
abscal_data, abscal_flags = deepcopy(data), deepcopy(flags)
calibrate_in_place(abscal_data, ac_gains, data_flags=abscal_flags, cal_flags=ac_flags)
# Find corresponding LSTs in the externally calibrated Abscal Model
from hera_cal.abscal import match_times
model_files = glob.glob('/lustre/aoc/projects/hera/nkern/idr3_abscal_models/full_model/' +
'zen.2458042.5*.HH.uvRXLS.uvh5') # cheatingly slightly by only checking the .5* files
matched_model_files = sorted(set(match_times(sample_data_file, model_files, filetype='uvh5')))
hdm = HERAData(matched_model_files)
model, model_flags, _ = hdm.read(bls=[bl])
# compare redundantly and then absolutely calibrated data to externally calibrated reference
waterfall_column([np.angle(abscal_data[bl]), np.angle(model[bl])],
[abscal_flags[bl], model_flags[bl]],
['Absolutely Calibrated: {}'.format(bl), 'Externally Calibrated Reference: {}'.format(bl)],
clabels=['Phase (Radians)', 'Phase (Radians)'],
ylabel='LST (Radians)',
cmaps = ['twilight', 'twilight'],
ylims = [[model.lsts[-1], model.lsts[0]], None],
extents=[[hd.freqs[0] / 1e6, hd.freqs[-1] / 1e6, hd.lsts[-1], hd.lsts[0]],
[model.freqs[0] / 1e6, model.freqs[-1] / 1e6, model.lsts[-1], model.lsts[0]]])
RFI Flagging¶
Our Radio Frequency Interference (RFI) flagging algorithm has largely been overhauled for IDR2.2, though it still relies fundamentally on the basic idea of looking for highly-significant outliers in detrended data and then also flagging less significant outliers that neighbor them in time or frequency, the so-called "watershed" algorithm. This relies on the idea that RFI events are compact in frequency or time or both.
New in IDR 2.2:
The following algorithm is used for RFI flagging on a per-observation basis.
- Compute metric waterfalls of calibration data products using median filters. We begin by flagging on the
omnicalgains,omnical$\chi^2/DoF$, theomnicalvisibility solutions, theabscalgains, and theabscal$\chi^2$. Each baseline, antenna, or polarization's waterfall is median filtered (i.e. the local median subtracted off in a 8 integration by 8 channel box) and each each data point is compared to the local median and median absolute deviation to compute a metric of "outlierness" in number of $\sigma$s. These "metrics" are averaged in quadrature to produce a single waterfall for each data product.
- Flag and watershed on metrics waterfalls individually and collectively to produce a set of initial flags. $5\sigma$ outliers on each waterfall are flagged, and then any $2\sigma$ outliers that neighbor a flag are also flagged. This process is performed both on individual calibration data products to look for specific issues with the associated step, and then on a equally-weighted average in quadrature of all 5 waterfalls to look for low level RFI. These are all ORed together to produce an initial set of flags.
- Use the initial flags to compute metric waterfalls of both calibration products and raw data using mean filters. The process is repeated on the five calibration data products above, but now also on the raw data itself (these were excluded from the first step because they are slow to work with). However, because the data now already have flags, we use a mean filter instead of a median filter, which can better handle setting the weight of flagged data to zero. Since the data are already pre-flagged, the worst outliers that would throw off this calculation are already gone. This is effectively computing a number of $\sigma$s for each time and frequency in terms of a weighted local standard deviation. These are then averaged in quadrature to form metrics waterfalls.
- Flag and watershed on new metrics waterfalls individually and collectively and OR those flags with the initial flags. Once again, we look for significant outliers ($> 5\sigma$) and minor neighboring outliers (growing to encompass all $> 2\sigma$ outliers that touch $> 5\sigma$ outliers or each other) and flag both on individual metrics waterfalls and on the quadrature average of all six waterfalls. All the resulting flags are ORed together with the initial flags.
- Perform day-long thresholding and broadcasting of flags. The metrics waterfalls used for flagging individual frequencies and times are subsequently analyzed for an entire day. We examine median metrics across frequency or across time, renormalize, and look for significant outliers (in this case $> 7\sigma$) or neighbors ($ >3\sigma$) to flag for all times or frequencies. This is performed on all the different metrics used above and also ORed together with the previous flagging.
- Update and save flags in updated
.flagged_abs.calfitsfile. Also save various intermediate metrics and flagging products as HDF5 files (readable withpyuvdata.UVFlag).*
# load flagged calibration solutions
flagged_abscal_file = os.path.join(analysis_folder, '2458098/zen.2458098.43124.HH.flagged_abs.calfits')
hc = HERACal(flagged_abscal_file)
fac_gains, fac_flags, fac_quals, fac_total_quals = hc.read()
waterfall_column([np.abs(fac_gains[ant]), np.angle(fac_gains[ant]), fac_total_quals[ant[1]]],
[fac_flags[ant], fac_flags[ant], fac_flags[ant]],
['Flagged Abscal Amplitude: {}'.format(ant),
'Flagged Abscal Phase: {}'.format(ant),
'Flagged Abscal $\chi^2$'],
clims=[[0, .04], None, [0, 5]],
clabels=['Amplitude (Unitless)', 'Phase (Radians)', '$\chi^2$ (Unnormalized)'],
cmaps = ['inferno', 'twilight', 'inferno'],
extents=[hc.freqs[0] / 1e6, hc.freqs[-1] / 1e6, hc.times[-1] - 2458098, hc.times[0] - 2458098])
This algorithm seems to catch all obvious RFI and also flags the impact of RFI in the abscal model (as seen in the abscal $\chi^2$).
Calibration Smoothing¶
To reduce the number of degrees of freedom in the calibration solution, we rely on the assumption that the instrument does not have spectral structure on small frequency scales and that it does not evolve quickly in time. Our smoothing algorithm loads in a whole day of calibration solutions simultaneously and put them on a time grid, leaving flagged holes for any missing files or missing time.
New in IDR 2.2:
Most notably, we have implemented 2D smoothing. In IDR 2.1, we smoothed in time and then in frequency. Here we do both simultaneously using a 2D CLEANing algorithm. We restrict our CLEANed solution in delay/fringe-rate space, filtering on 10 MHz (i.e. 100 ns) and 30 minute (0.556 mHz) scales. This filtering is done with a Tukey window with $\alpha = .5$ in frequency and a top-hat window in time.
The calibration smoothing step also checks for antennas that are often flagged as bad while others are not. We count the number of unflagged visibilities on each antenna and any antenna with less than 50% of the maximum among all antennas is flagged before smoothing is performed.
Here is an example of a calibration smoothing in action:
smoothed_abscal_file = os.path.join(analysis_folder, '2458098/zen.2458098.43124.HH.smooth_abs.calfits')
hc = HERACal(smoothed_abscal_file)
sac_gains, sac_flags, _, _ = hc.read()
# find antenna that was picked as the reference by looking for unflagged antennas with 0 phase
sc_refant = {pol: sorted([(np.mean(np.abs(np.angle(sac_gains[ant]))), ant)
for ant in sac_gains if ant[1]==pol and not np.all(sac_flags[ant])])[0][1]
for pol in hc.pols}
for pol in hc.pols:
print('Antenna {} chosen as the reference for {}.'.format(sc_refant[pol][0], pol))
# rephase abscal to to have the same reference antenna as smooth_cal
rephased_abscal = fac_gains[ant] * np.abs(fac_gains[sc_refant[ant[1]]]) / fac_gains[sc_refant[ant[1]]]
smoothcal_gains = deepcopy(sac_gains[ant])
waterfall_column([np.abs(rephased_abscal), np.abs(smoothcal_gains)],
[fac_flags[ant], sac_flags[ant]],
['Abscal Amplitude: {}'.format(ant), 'Smoothcal Amplitude: {}'.format(ant)],
clabels=['Amplitude (Unitless)', 'Amplitude (Unitless)'],
cmaps = ['inferno', 'inferno'],
figsize=(12,4),
clims=[[0, .04], [0, .04]],
extents=[hc.freqs[0] / 1e6, hc.freqs[-1] / 1e6, hc.times[-1] - 2458098, hc.times[0] - 2458098])
# replace flags with nans in order to take a nanmedian over time
rephased_abscal[fac_flags[ant]] = np.nan
smoothcal_gains[sac_flags[ant]] = np.nan
# plot median smooth calibration
plt.figure(figsize=(12,3), dpi=100)
plt.plot(hc.freqs / 1e6, np.nanmedian(rephased_abscal, axis=0).real, '.', label='Flagged Abscal')
plt.plot(hc.freqs / 1e6, np.nanmedian(smoothcal_gains, axis=0).real, label='Smoothcal')
plt.title('Median Abscal and Smoothcal Gains')
plt.ylabel('Real Part of Gains (Unitless)')
plt.xlabel('Frequency (MHz)')
plt.ylim(-.05, .05)
plt.legend();
Delay Filter¶
At this point in the pipeline, we are ready to LST-bin our calibrated data. However, we also produce calibrated and delay-filtered data. These files, the .OCRSD.uvh5 files are the only full-sized data files produced before LST-binning. OCRSD is a legacy nomenclature meaning Omicaled, absCaled, xRfi, Smoothcaled, Delay-filtered data.
Delay filtering is performed independently for each integration and each baseline in the calibrated and flagged data. Delay filtering removes Fourier content outside the baseline's horizon delay, plus a 15 ns buffer. To account for RFI gaps, it is done using the CLEAN algorithm with a $10^{-9}$ tolerance and a Tukey window with $\alpha=0.5$.
This procedure is largely the same as in IDR2.1 (the only significant update was improved memory management using partial I/O enabled by .uvh5).
delay_filtered_file = os.path.join(analysis_folder, '2458098/zen.2458098.43124.HH.OCRSD.uvh5')
hd_dfil = HERAData(delay_filtered_file)
dfil_data, dfil_flags, _ = hd_dfil.read(bls=[bl])
# Apply smooth_abs calibration to data
sac_data, sac_data_flags = deepcopy(data), deepcopy(flags)
calibrate_in_place(sac_data, sac_gains, data_flags=sac_data_flags, cal_flags=sac_flags)
waterfall_column([np.abs(sac_data[bl]), np.abs(dfil_data[bl]), np.angle(sac_data[bl]), np.angle(dfil_data[bl])],
[sac_data_flags[bl], dfil_flags[bl], sac_data_flags[bl], dfil_flags[bl]],
['Calibrated Amplitude: {}'.format(bl), 'Delay-Filtered Amplitude: {}'.format(bl),
'Calibrated Phase: {}'.format(bl), 'Delay-Filtered Phase: {}'.format(bl)],
clabels=['Amplitude (Jy)', 'Amplitude (Jy)', 'Phase (Radians)', 'Phase (Radians)'],
clims=[[0, 150], [0, 15], None, None],
cmaps = ['inferno', 'inferno', 'twilight', 'twilight'],
extents=[hd_dfil.freqs[0] / 1e6, hd_dfil.freqs[-1] / 1e6,
hd_dfil.times[-1] - 2458098, hd_dfil.times[0] - 2458098],
figsize=(12,7))
(Note the different color scale on the delay-filtered amplitude.)
While the delay filter vastly reduces the amplitude of the visibility to much closer to the noise level, the result is not purely noise-like, as can be seen in the spectral structure in the delay-filtered phase. This may be, in part, due to some of the non-fringing structure and the cable reflections, both examined in the two Kern et al. (submitted) papers examining systematics.
Ancillary Day-by-Day Data Products¶
The main goal of analysis is to calibrate and RFI-flag to as to enable LST-binning and the resultant reduction in data volume. However, we perform several additional steps that produce useful data products for quality assessment of the data and various open research questions.
Update omnical Visibility Solutions with Absolute Calibration and Flagging¶
New in IDR 2.2:
Previously, omnical visibility solutions we only used to look for RFI. Now, we update the omnical visibilities by filling in their degenerate components using either the results of abscal or the results of smooth_cal (we do both and leave it as an open question whether there's any advantage or one or the other). This step is performed by taking the average of the ratio of omnical gain products to abscal/smooth_cal gain products in a redundant baseline group. This isolates the degenerate subspace, which we then use to update the visibility solutions. Like the omnical visibility solutions, these .uvh5 data products contain one visibility per unique baseline group and an updated nsamples.
While we do no subsequent processing of these data, they are meant to help enable research questions about whether this data reductions step (which is substantial for HERA) can help us get to the power spectrum quicker without losing information.
# load representative baseline from updated omnical visibility solutions
fac_vis_file = os.path.join(analysis_folder, '2458098/zen.2458098.43124.HH.flagged_abs_vis.uvh5')
hd_fac = HERAData(fac_vis_file)
facvis, facvis_flags, _ = hd_fac.read(bls=[red_dict[bl]])
sac_vis_file = os.path.join(analysis_folder, '2458098/zen.2458098.43124.HH.smooth_abs_vis.uvh5')
hd_sac = HERAData(sac_vis_file)
sacvis, sacvis_flags, _ = hd_sac.read(bls=[red_dict[bl]])
waterfall_column([np.angle(abscal_data[bl]), np.angle(facvis[red_dict[bl]]), np.angle(sacvis[red_dict[bl]])],
[sacvis_flags[red_dict[bl]], sacvis_flags[red_dict[bl]], sacvis_flags[red_dict[bl]]],
['Absolutely Calibrated and Flagged: {}'.format(bl),
'Omnical Solution Updated with Flagged Abscal: {}'.format(bl),
'Omnical Solution Updated with Smoothed Abscal: {}'.format(bl)],
clabels=['Phase (Radians)', 'Phase (Radians)', 'Phase (Radians)'],
cmaps = ['twilight', 'twilight', 'twilight'],
extents=[hd.freqs[0] / 1e6, hd.freqs[-1] / 1e6, hd.times[-1] - 2458098, hd.times[0] - 2458098])
waterfall_column([np.abs(abscal_data[bl]), np.abs(facvis[red_dict[bl]]), np.abs(sacvis[red_dict[bl]])],
[sacvis_flags[red_dict[bl]], sacvis_flags[red_dict[bl]], sacvis_flags[red_dict[bl]]],
['Absolutely Calibrated and Flagged: {}'.format(bl),
'Omnical Solution Updated with Flagged Abscal: {}'.format(bl),
'Omnical Solution Updated with Smoothed Abscal: {}'.format(bl)],
clabels=['Amplitude (Jy)', 'Amplitude (Jy)', 'Amplitude (Jy)'],
cmaps = ['inferno', 'inferno', 'inferno'],
extents=[hd.freqs[0] / 1e6, hd.freqs[-1] / 1e6, hd.times[-1] - 2458098, hd.times[0] - 2458098])
Raw Autocorrelations and Noise¶
New in IDR 2.2:
To enable future calibration tests, the pipeline extracts autocorrelations. While in principle this is duplicating data and is otherwise trivial with partial I/O, this step anticipates IDR 3.1 and beyond where we cannot keep the full data set staged on lustre. For future tests of calirabtion and noise, it was requested to store these raw data products alongside calibration solutions. These autocorrelations are stored in .autos.uvh5 and only have autocorrelation baseline keys.
Using calibrated autocorrelations, we can predict the noise variance on visibilities, $\sigma_{ij}^2$. Namely,
$\sigma_{ij}^2 = V_{ii} V_{jj}$ $ / $ $B t$
where $B$ is the bandwidth of a channel and $t$ is the integration time. Instead of computing this quantity for all baselines, we instead compute and save $\sigma_{ii}$ where
$\sigma_{ij} \equiv \sqrt{\sigma_{ii} \sigma_{jj}} = \left(V_{ii} / \sqrt{Bt}\right) \left( V_{jj} / \sqrt{Bt} \right)$.
These quantities, $\sigma_{ii}$, are stored in .noise_std.uvh5 files. Though they are technically per-antenna, we felt it more sensible to store them as visibility data files (since the units are Jy) with autocorrelation keys instead of storing them in .calfits files.
from hera_cal.utils import join_bl
ant1, ant2 = split_bl(bl)
auto_bl1 = join_bl(ant1, ant1)
auto_bl2 = join_bl(ant2, ant2)
# Load autocorrelation
autos_file = os.path.join(analysis_folder, '2458098/zen.2458098.43124.HH.autos.uvh5')
hd_autos = HERAData(autos_file)
autos, auto_flags, _ = hd_autos.read(bls=[auto_bl1, auto_bl2])
# Calibrate autocorrelation
cal_autos, cal_auto_flags = deepcopy(autos), deepcopy(auto_flags)
calibrate_in_place(cal_autos, sac_gains, data_flags=cal_auto_flags, cal_flags=sac_flags)
# Load inferred noise on calibrated data
noise_file = os.path.join(analysis_folder, '2458098/zen.2458098.43124.HH.noise_std.uvh5')
hd_noise = HERAData(noise_file)
noise, noise_flags, _ = hd_noise.read(bls=[auto_bl1, auto_bl2])
bl_noise = np.sqrt(noise[auto_bl1] * noise[auto_bl2])
bl_noise_flags = noise_flags[auto_bl1] | noise_flags[auto_bl2]
waterfall_column([np.abs(autos[auto_bl1]), np.abs(cal_autos[auto_bl1]),
np.abs(noise[auto_bl1]), np.abs(bl_noise)],
[auto_flags[auto_bl1], cal_auto_flags[auto_bl1],
noise_flags[auto_bl1], bl_noise_flags],
['Raw Autocorrelations: {}'.format(auto_bl1),
'Smooth-Calibrated Autocorrelations: {}'.format(auto_bl1),
'Inferred Antenna Noise: {}.'.format(ant1),
'Inferred Visibility Noise: {}.'.format(bl)],
clabels=['Amplitude\n(Uncalibrated)', 'Amplitude (Jy)', 'Amplitude (Jy)', 'Amplitude (Jy)'],
clims=[[0, 15], None, [0,15], [0,15]],
cmaps = ['inferno', 'inferno', 'inferno', 'inferno'],
extents=[hd_autos.freqs[0] / 1e6, hd_autos.freqs[-1] / 1e6,
hd_autos.times[-1] - 2458098, hd_autos.times[0] - 2458098],)
We can also check our inferred value for the noise on visibilities by checking them against a sequential difference of the data. In this case, we use hera_cal.noise.interleaved_noise_variance_estimate() to estimate the noise on the data by subtracting 0.5 times the next and previous channels from the data. Averaging in time over the file, we see that these two estimates of the noise agree.
# Estimate noise from visibility data using interleaved frequencies
from hera_cal.noise import interleaved_noise_variance_estimate
data_with_nans = deepcopy(sac_data[bl])
data_with_nans[sac_data_flags[bl]] = np.nan
noise_var_est = interleaved_noise_variance_estimate(data_with_nans, kernel=[[-.5, 1, -.5]])
with np.errstate(divide='ignore', invalid='ignore'):
interleaved_noise = np.sqrt(np.nanmean(noise_var_est, axis=0))
# Estimate noise on baseline using autocorrelations
var_with_nans = noise[auto_bl1] * noise[auto_bl2]
var_with_nans[sac_data_flags[bl]] = np.nan
autocorrelation_noise = np.sqrt(np.nanmean(var_with_nans, axis=0))
# Plot Results
plt.figure(figsize=(12,3), dpi=100)
plt.plot(hd.freqs / 1e6, interleaved_noise, label='Interleaved Noise Estimate')
plt.plot(hd.freqs / 1e6, autocorrelation_noise, lw=3, label='Noise Inferred from Autocorrelations')
plt.xlim(100,200)
plt.xlabel('Frequency (MHz)')
plt.ylabel('Amplitude (Jy)')
plt.legend();
Reflection Fitting¶
New in IDR 2.2:
Calibration smoothing erases small scale spectral structure in calibration solutions, some of which we believe is real. For example, we see clear evidence from the autocorrelations of ripples in the calibration solutions attributable to reflection systematics. An extra trip down and back up 150 m cables produce delays of 1250 ns, assuming a speed of light in the cables of $0.8c$. Such reflections were seen in every antenna at different delays and amplitudes in HERA Memo #64.
In the two forthcoming Kern et al. papers on HERA systematics, we develop and test a method for mitigating these systematics by fitting the reflection in delay space and then running an iterative, non-linear optimizer to refine the estimate of the reflection delay, amplitude, and phase.
In theory, the calibration solutions found for the cable reflections could be incorporated into the final, smoothed calibration solution by multiplying (an example of which we show below). In practice, we feel that this move is risky until the stability of the fits and of the cable reflections themselves are more thoroughly explored. Therefore, we provide calibration solutions of just cable reflections in .reflections.calfits files in order to support open research questions about these reflections.
reflections_file = os.path.join(analysis_folder, '2458098/zen.2458098.43124.HH.reflections.calfits')
hc = HERACal(reflections_file)
rfl_gains, rfl_flags, _, _ = hc.read()
waterfall_column([np.abs(rfl_gains[ant]), np.abs(rfl_gains[ant] * sac_gains[ant])],
[rfl_flags[ant], sac_flags[ant]],
['Cable Reflection Fit: {}'.format(ant), 'Smoothcal x Reflections: {}'.format(ant)],
clabels=['Gain Amplitude\n(Unitless)', 'Gain Amplitude\n(Unitless)'],
cmaps = ['inferno', 'inferno'],
extents=[hd.freqs[0] / 1e6, hd.freqs[-1] / 1e6, hd.times[-1] - 2458098, hd.times[0] - 2458098],
figsize=(12,4))
Multi-Frequency Synthesis Imaging with CASA¶
New in IDR 2.2:
For this IDR, all files are also passed through a multi-frequency synthesis imagine routine in CASA. These images are largely for diagnostic purposes, they are not meant for direct comparison to external catalogs or for any sort of power spectrum analysis. We image both in instrumental polarizations, i.e. ee and nn, and in pseudo-Stokes I, Q, U, and V. The various imaging products are stored as .fits in the *.calibrated.uvh5_image/ folder corresponding to each file.
The image products are dirty images (no CLEANing has been done, i.e. niter=0) and are not primary-beam corrected. We use Briggs weighting with robust=0. We synthesize the image using channels 150 to 900, corresponding to 114.6 - 187.9 MHz. The images are 512 by 512 pixels centered at zenith with 500 arcsecond pixels.
from astropy.io import fits
from astropy import wcs
from hera_cal.utils import polnum2str
image_fits = os.path.join(analysis_folder, '2458098/zen.2458098.41632.HH.calibrated.uvh5_image/' +
'zen.2458098.41632.HH.calibrated.uvh5.image.image.fits')
hdulist = fits.open(image_fits)
coords = wcs.WCS(hdulist[0].header, naxis=[wcs.WCSSUB_CELESTIAL])
fig, axes = plt.subplots(2, 2, figsize=(10, 10), dpi=100, subplot_kw={'projection': coords})
plt.subplots_adjust(hspace=0, wspace=.3)
for i, (ax, image_data) in enumerate(zip(np.ravel(axes), hdulist[0].data[:, 0, ...])):
im = ax.imshow(image_data, cmap='inferno')
ax.grid()
ax.set_xlim([150, image_data.shape[0] - 150])
ax.set_ylim([150, image_data.shape[1] - 150])
pol = polnum2str(int(i * hdulist[0].header['CDELT4'] + hdulist[0].header['CRVAL4'])).replace('p', 'Pseudo-')
plt.colorbar(im, ax=ax, label='Jy/Beam', fraction=0.046, pad=0.04)
props = dict(boxstyle='round', facecolor='white', alpha=0.8)
ax.text(0.05, 0.95, 'Stokes ' + pol, transform=ax.transAxes, fontsize=14, verticalalignment='top', bbox=props)
This field can be compared to the one used in HERA Memo #42 to calibrate IDR 2.1, thought a different field was used to calibrate this data (see above). Overall, the agreement in terms of source positions and flux scale is pretty good, though a more detailed analysis of the image data products is warranted.
LST-Binned Data Products¶
LST-Binning¶
Both before and after delay filtering, the data is LST-binned into a single group and into two groups of interleaved nights. During LST binning, each integration is assigned to the the nearest LST bin, but is rephased to account for the slight difference in LST between its bin center and the LST bin's center. LST bins were chosen to be twice as long as raw data bins (21.475 seconds = 0.00156597 radians of a sidreal day), so every LST bin gets two data points for every input night.
Next, the list of rephased data points that would go into any given LST and frequency bin are examined for outliers. For relatively robust outlier detection, we compute the median absoute deviation (MAD) of all the points that are to be binned together and use that to compute a modified z-score. As long as there are at least 5 data points to bin together, we then flag all data that with a modified z-score greater than 5, which corresponds to $5\sigma$ outliers if the distribution were Gaussian.
Finally, we compute the standard deviation of both the real and imaginary parts of the unflagged data and store those in the .STD. file. We also save the number of samples that went into each binned data point in both the .LST. and .STD. files. We do not further perform a cut on the number of samples or any other form of post-LST binning RFI flagging. The procedure is repeated with the delay-filtered data as well.
For the nights split into two groups, we have picked:
- grp1: 2458098, 2458101, 2458103, 2458105, 2458107, 2458109, 2458111, 2458113, 2458115
- grp2: 2458099, 2458102, 2458104, 2458106, 2458108, 2458110, 2458112, 2458114, 2458116
The data can be found in /lustre/aoc/projects/hera/H1C_IDR2/IDR2_1/LSTBIN/{one/two}_group/grp{N}(_dly)/ and have the format like zen.grp1.of1.LST.1.31552.HH.OCRSL.uvh5. This means:
grp1: 1-indexed group of nights that the LST binned data products are divided intoof2: number of groups of nights that the data was divided into before separately LST-binning each group togetherLST: LST-binned data. The other alternative isSTD, which encodes the standard deviation of the real and imaginary parts (separately) of the binned data in its real and imaginary data array1.31552: starting LST of the data in radiansOCRSL: this denotes calibrated, RFI-flagged, LST-binned data, but not delay-filtered. The LST-binned results also includeOCRSDLfiles, which have been delay-filtered before LST-binning. Those are in folders that look likegrp1_dlyinstead of justgrp1. L is for LST-binning.
Here we show a comparison of LST-binned data with data from a single night:
lstbinned_file = os.path.join(analysis_folder,'LSTBIN/one_group/grp1/zen.grp1.of1.LST.1.31552.HH.OCRSL.uvh5')
hd_lst = HERAData(lstbinned_file)
lst_data, lst_flags, lst_nsamples = hd_lst.read(bls=[bl, auto_bl1, auto_bl2])
# compare redundantly and then absolutely calibrated data to externally calibrated reference
waterfall_column([np.angle(sac_data[bl]), np.angle(lst_data[bl]), lst_nsamples[bl]],
[sac_data_flags[bl], lst_flags[bl], np.zeros_like(lst_flags[bl])],
['2458098: {}'.format(bl), 'LST-binned: {}'.format(bl), 'LST-binning samples: {}'.format(bl)],
clabels=['Phase (Radians)', 'Phase (Radians)', 'Number of Samples'],
ylabel='LST (Radians)',
cmaps = ['twilight', 'twilight', 'inferno'],
ylims = [[hd_lst.lsts[-1], hd_lst.lsts[0]], None, None],
extents=[[hd.freqs[0] / 1e6, hd.freqs[-1] / 1e6, hd.lsts[-1], hd.lsts[0]],
[hd_lst.freqs[0] / 1e6, hd_lst.freqs[-1] / 1e6, hd_lst.lsts[-1], hd_lst.lsts[0]],
[hd_lst.freqs[0] / 1e6, hd_lst.freqs[-1] / 1e6, hd_lst.lsts[-1], hd_lst.lsts[0]]])
In general, the these results look reasonable. The flagging is substantially more aggressive than in IDR2.1, which is reflected in the lower and more variable number of samples above, through there are still a couple of "clean" windows around 125 MHz and around 160 MHz.
We can also compare the observed noise level to what we excepted from the autocorrelations (both from individual days and from the LST-binned autocorrelations). This match is pretty good, though there's still some evidence for temporal structure in the STD files, especially in the real/imaginary components.
lststd_file = os.path.join(analysis_folder,'LSTBIN/one_group/grp1/zen.grp1.of1.STD.1.31552.HH.OCRSL.uvh5')
hd_std = HERAData(lststd_file)
std_data, std_flags, _ = hd_std.read(bls=[bl])
# predict noise from lst-binned autocorrelations
from hera_cal.noise import predict_noise_variance_from_autos
tint = np.median(np.diff(hd.times)) * 24 * 3600
lst_bl_noise = np.sqrt(predict_noise_variance_from_autos(bl, lst_data, dt=tint))
lst_bl_flags = lst_flags[bl] | lst_flags[auto_bl1] | lst_flags[auto_bl2]
# compare redundantly and then absolutely calibrated data to externally calibrated reference
waterfall_column([np.abs(bl_noise), np.abs(lst_bl_noise), np.abs(std_data[bl])],
[bl_noise_flags, lst_bl_flags, std_flags[bl]],
['Predicted Noise from Autos: ' + str(bl),
'Predicted Noise from LST-Binned Autos: ' + str(bl),
'|STDEV| over Nights: {}'.format(bl)],
clabels=['$\sigma$ (Jy)', '$\sigma$ (Jy)', '$\sigma$ (Jy)'],
ylabel='LST (Radians)',
clims=[[0,15], [0,15], [0,15]],
ylims = [[hd_lst.lsts[-1], hd_lst.lsts[0]], None, None],
cmaps = ['inferno', 'inferno', 'inferno'],
extents=[[hd.freqs[0] / 1e6, hd.freqs[-1] / 1e6, hd.lsts[-1], hd.lsts[0]],
[hd_std.freqs[0] / 1e6, hd_std.freqs[-1] / 1e6, hd_std.lsts[-1], hd_std.lsts[0]],
[hd_std.freqs[0] / 1e6, hd_std.freqs[-1] / 1e6, hd_std.lsts[-1], hd_std.lsts[0]]])
The Path Forward¶
The goal of IDR 2.2 is to provide calibration, flagging, and LST-binning to support an astrophysically interesting power spectrum upper limit. It is not meant to be the last word in HERA analysis. Hopefully, this release is "good enough" to move forward, but it is important to document the issues as they arise and continue to cross-check our results and better understand the systematics in the data.
Known Issues and Future Improvements¶
- Abscal model flagging: The abscal model shows evidence for low-level RFI going unflagged, contaminating nearly channels due to the low-pass delay filter. While our subsequent flagging algorithm appears to catch these outliers, better flagging of the externally calibrated data set would remove the the need for this. Also, we expect this issue to lead to RFI flags that repeat at the same LST day after day. We have an open Github issue to address this.
- Antennas flagged in abscal model are flagged throughout the data: As it stands, any antennas that are completely flagged in the abscal model are flagged in the data as well. This is probably too conservative. We have an open Github issue to address this.
- Abscal model improvements: A number possible improvements to the externally calibrated abscal model have been suggested. These include:
- Fringe-rate filtering the abscal model (along with low-pass delay filtering it) to improve signal to noise.
- Including
enandnevisibilities. As of now, our the externally calibrated data used as our abscal model has no polarization information, so we cannot calibrate the absolute phase of x relative to y. This makes pseudo-Stokes U and V difficult to interpret. - Calibrating off a limited set of trustworthy, redundantly-averaged (or possibly simulated) baselines. This might be the only practicable way to make abscal work for RTP in the next season.
- Abscal $\chi^2$ is not normalized nor particularly well-understood. We could use a systematic investigation of the discrepancies between the abscal model and the absolutely calibrated data.
- Cross-talk. Even after the delay filter, we see clear evidence of persistent structure on short baselines. This was explored in the two forthcoming Kern et al. papers. At this point, at this point, our plan is to apply these techniques to LST-binned data before forming power spectra, though it remains an open question whether this is best done before LST-binning---which is to say, does the cross-talk repeat from day to day?
- Fringe structure in STD files. The real and imaginary components of the STD files appear to have fringe structure in them. This is likely attributable to un/mis-calibrated gain variations from day to day, though we don't know that for sure. Does this mean that data is not integrating down properly? Is it also a problem in the delay-filtered, LST-binned, data?
Open Questions¶
Cable Reflections. We solve for cable reflections at every integration, but our current plan is to ignore these results and only solve for and then calibrate out cable reflections in LST-binned data. How accurate and reliable are our cable reflection fits in this IDR? And do the cable reflections exhibit any temporal evolution or additional spectral structure (beyond the three parameter model used to constrain them?)
Temporal structure. How stable are the calibration solutions in time (per feed)? What kind of temporal structure do they exhibit and how repeatable is that structure from day to day? What is the optimal temporal smoothing scale for calibration solutions? (This is under investigation by UCB undergraduate Max Lee and a memo on this question is forthcoming.)
Spectral structure. What kind of spectral structure do our calibration solutions exhibit? How does it vary from antenna to antenna? From integration to integration? From day to day? Is the spectral structure in our calibration solutions consistent with that expected from beam modeling of spectral structure? What is the optimal spectral smoothing scale for calibration solutions? Is there some spectral form other than delay modes on which we should be filtering (e.g. Hermite polynomials? PCA/Weiner filter?).
RFI Report. What are the persistent sources of RFI? How much of the band are we losing? Are there any major transient broadband events and can they be tracked to any source (e.g. MeerKAT, as we saw in 2458109)?
Post-LST RFI. Is there any low level RFI that we are not detecting? Are we overflagging? Does the LST-binned data product show any evidence for RFI? Should we be removing RFI after LST binning?
Versus Sky-Based Calibration. How do the calibration solutions compare to those generated purely from a sky-based calibration?
Relative e-to-n Phase Calibration. Right now, our abscal model has no en or ne visibilities, which means that the relative phase of x to y as a function of frequency is uncalibrated. Is this a problem? How is it affecting Stokes U and V?
Imaging. We have produced MFS images, but we have not examined them systematically. Do our maps look reasonable? Do source fluxes and spectra look right? What about psuedo-Stokes I, Q, U, and V images?
Noise integration. Does the noise integrate down from day to day?
LST results. Does the difference between a given day's data and the corresponding LST-binned data look noise-like? Are there some days that are worse than others?
Redundancy. How redundant are redundant baselines? How do we interpret $\chi^2$? Does the noise integrate down inside redundant baseline groups?
Repeatability of calibrated data. How repeatable are redundant baselines from day to day after calibration? Can or should we be LST-binning updated omnical visibility solutions?
Beams. Does our beam model make sense? Can we constrain it with foregrounds moving through the beam?
Delay filtering. We use a CLEAN-based filtering in calibration smoothing and in the delay transform steps. Is that the best algorithm for it? Are we using the best tapering function?
Delay spectrum. How does the delay spectrum look? Is it noise-like at k~.2 and beyond after foreground filtering on the EoR baselines? How does it depend on LST?
If you want to work on any of these questions (or already are!) or are curious about this IDR, please join our HERA analysis telecon on Tuesdays at 10am Pacific on Zoom or the #hera-analysis Slack channel.