Accessing the Gains used in the SDP pipeline calibrations

It is possible to access the gain tables produced by the pipeline if you use the mvf4/katdal data format.

Unfortunately it is not possible to access the tables if you use the mvf4toms.py script to create a CASA file as we don’t have working code available to convert our solutions to CASA tables yet.

Below is some example code which I hope will illustrate how to access the solutions, there are two different methods you can use so you can chose which suits your purposes best.

import katdal
import numpy as np

fname = '123.rdb?token'
f= katdal.open(fname)

#telstate is where all the metadata about the observation is stored
telstate = f.source.telstate
# polarisation ordering for calibration solutions
pol=telstate['cal_pol_ordering']
antenna_names=telstate['cal_antlist']

# The pipeline solutions are stored in telstate 
# You can get these using e.g below , st = start time, et= end time
# so you can select based on time range
gsol = telstate.get_range('cal_product_G', st=0)

# This produces a list of tuples, where the first element of the 
# tuple is the pipeline solution and the second element is the
# timestamp e.g. 

print("The gain solution at time {} is an array with shape {}"\
      .format(gsol[0][1], gsol[0][0].shape))

# The gain solution has shape (npols, nants),
# the ordering of the polarisation products is given by
# telstate['cal_pol_ordering']
# i.e .if telstate['cal_pol_ordering'] = ['v', 'h'] then g[0,:] is v 
# the ant ordering is given by telstate['cal_antlist']

# To retrieve the gains as an array with time on the first axis
# you can use for example the following function  
def get_cal_product(ts_attrs, cal, st):
    vals, times = np.array([]), []
    product = ts_attrs.get_range(cal, st=st)
    if len(product) > 0:
        vals, times = zip(*product)
        vals = np.array(vals)
    return vals, times
    
gsolns, gtimes = get_cal_product(telstate, 'cal_product_G', 0)
# now gsolns has shape (ntimes, npols, nants)

# Delays can be retrieved in a similar fashion, and have the same shape as gains
# they are measured in seconds 
ksolns, ktimes = get_cal_product(telstate, 'cal_product_K', 0)

# The Bandpass is also similar but it is appended by a number in range(0,4)
# 4K wideband only has 'cal_product_B0' for the whole bandpass
# 32K modes have cal_products_B0 to B3, each of these products
# represents a quarter of the bandpass
# Bandpass has shape (nchans, npols, nants)
bsolns, btimes = get_cal_product(telstate, 'cal_product_B0', 0)

#Alternatively if you use the applycal method you can access
# the corrections (which are 1/solutions) as sensors for each antenna
# and these are already interpolated onto the timestamps of the observtion.
# e.g 
f= katdal.open(fname, applycal='all')

Kcorr=f.sensor['Calibration/Corrections/l1/K/m000h']
Bcorr=f.sensor['Calibration/Corrections/l1/B/m000h']
Gcorr=f.sensor['Calibration/Corrections/l1/G/m000h']

# Here each antenna has a correction sensor with axis (ntimes, nchans), the 
# ntimes will match the timestamps in the observation 
# e.g try
print('The corrections shape is {}'.format(Kcorr.shape))
print('The visibility shape is {}'.format(f.vis.shape))
print('There are {} timestamps for both the \
sensor and the visibilities'.format(Kcorr.shape[0]))

# If you select specific timestamps using the katdal select function
# ie f.select it will alter the selection on the correction sensor as well 
# e.g try
f.select(scans=1)
Kcorr=f.sensor['Calibration/Corrections/l1/K/m000h']
print('The corrections shape is {}'.format(Kcorr.shape))
print('The visibility shape is {}, so the timestamps match'.format(f.vis.shape))