First off, we'll be using pylab to plot up the solutions
In [7]:
%pylab inline
The default method for accessing calibrations is from a measurement set format (which is a directory).
In [114]:
from askap.accessors.calibration import Source
In [112]:
ls -ld *.tab/
In [113]:
src = Source('calibdata.tab')
The most recent solution ID is available as a property:
In [15]:
src.most_recent_solution_id
Out[15]:
If you're just intrested in the most recent solution you can get the accessor with
In [103]:
ac = src.most_recent_accessor
Raw access to the bandpass, gain, leakage, and full Jones matrix is available with the get() commands, e.g.
In [18]:
ant = 3
beam= 0
chan = 4
In [19]:
ac.get_gain(ant, beam)
Out[19]:
Aaah, but waht is this freaky JonesJTerm? - it just has 4 useful properties
In [20]:
dir(ac.get_gain(ant, beam))
Out[20]:
In [21]:
g = ac.get_gain(ant, beam)
print 'Gain was', g.g1, g.g2, 'Valid flags:', g.g1IsValid, g.g2IsValid
In [22]:
bp = ac.get_bandpass(ant, beam, chan)
In [23]:
print 'Bandpass was', bp.g1, bp.g2, 'Valid flags:', bp.g1IsValid, bp.g2IsValid
In [26]:
lk = ac.get_leakage(ant, beam)
dir(lk)
Out[26]:
In [28]:
print 'Leakage was', lk.d12, lk.d21, 'Valid flags:', lk.d12IsValid, lk.d21IsValid
In [105]:
ac.get_jones(ant, beam, chan)
Out[105]:
If you want access to the whole data as one big numpy block, try this:
In [31]:
bp_g1, bp_g2 = ac.bandpass[:,:,:]
In [32]:
bp_g1.shape
Out[32]:
The bandpass property is sliceable, and returns g1 and g2 as a tuple. They are each masked numpy arrays of size Nant x Nbeam x Nchan. If you feel the need, you can plot them up too (although this case is not especially interesting):
In [61]:
subplot(2,2,1)
plot(abs(bp_g1[:, 0, :]).T)
ylabel('Amplitude')
title('G1')
subplot(2,2,3)
plot(degrees(angle(bp_g1[:, 0, :]).T))
ylabel('Phase (deg)')
xlabel('channel')
subplot(2,2,2)
plot(abs(bp_g2[:, 0, :]).T)
title('G2')
subplot(2,2,4)
plot(degrees(angle(bp_g2[:, 0, :]).T))
xlabel('channel')
Out[61]:
Maybe un-transpose that puppy up a little
In [106]:
subplot(2,2,1)
plot(abs(bp_g1[:, 0, :]))
ylabel('Amplitude')
title('G1')
subplot(2,2,3)
plot(degrees(angle(bp_g1[:, 0, :])))
ylabel('Phase (deg)')
xlabel('Antenna')
subplot(2,2,2)
plot(abs(bp_g2[:, 0, :]))
title('G2')
subplot(2,2,4)
plot(degrees(angle(bp_g2[:, 0, :])))
xlabel('Antenna')
Out[106]:
Hmm, looks like the same solution for every channel - yeah, we saw that before!
In [41]:
jones = ac.jones[:,:,:]
In [42]:
type(jones)
Out[42]:
In [43]:
jones.shape
Out[43]:
In [45]:
d12, d21= ac.leakage[:,:]
In [46]:
d12.shape
Out[46]:
Currently there's no way to fix this, as the underlying library only tells you you're out of bounds by throwing an exception
In [47]:
g1, g2 = ac.gain[:,:]
In [48]:
g1.shape
Out[48]:
Gain does the same for this dataset
In [49]:
ac2 = Source('caldata1.tab').most_recent_accessor
In [53]:
g1, g2 = ac2.gain[:,:]
In [54]:
g1.shape
Out[54]:
Aah, look, the gain is a sensible size in this example
In [60]:
plot(abs(g1))
plot(abs(g2))
ylabel('Amplitude')
xlabel('Antenna')
figure()
plot(degrees(angle(g1)))
plot(degrees(angle(g2)))
ylabel('Phase (deg)')
xlabel('Antenna')
Out[60]:
In [69]:
d12, d21 = ac2.leakage[:,:]
In [72]:
d12.shape
Out[72]:
But leakage isn't....
Parsets are just text files
In [92]:
ls -l *.in
Here's what a parset-formatted gain solution looks like
In [101]:
with open('rndgains.in', 'rU') as f:
lines = f.readlines()
print '*'*20, 'First 15 lines', '*'*20
print ''.join(lines[0:15])
print '\n\n', '*'*20, 'Last 15 lines', '*'*20
print ''.join(lines[-15:])
If the data is contained in a parset, you can get it by setting type='parset' in the Source constructor:
In [73]:
src = Source('rndgains.in', type='parset')
In [107]:
ac3 = src.most_recent_accessor
Note that the leakage slice is still huge (128x128), even though there are only 36 antennas and 1 beam in the parset:
In [108]:
d12, d21 = ac3.leakage[:,:]
In [109]:
d12.shape
Out[109]:
Note how plotting abs(d12) gives only 36 valid values, but plotting angle(d12) gives 128 antennas - somehow the masking is not working. Oh well, caveat emptor.
In [110]:
subplot(2,2,1)
plot(abs(d12[:, 0]))
ylabel('Amplitude')
title('d12')
subplot(2,2,3)
plot(np.degrees(np.angle(d12[:, 0])))
ylabel('Phase (deg)')
xlabel('Antenna')
subplot(2,2,2)
plot(abs(d21[:, 0]))
title('d21')
subplot(2,2,4)
plot(degrees(angle(d21[:, 0])))
xlabel('Antenna')
Out[110]: