Skip to content

Commit 69ca3c8

Browse files
timcornwelltimcornwell
committed
Add configuration to Gaintable
1 parent 6af7089 commit 69ca3c8

5 files changed

Lines changed: 109 additions & 29 deletions

File tree

‎.gitignore‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,3 +24,4 @@ dask-worker-space
2424
*.egg
2525
/data_models/.idea
2626
*.fits
27+
/sourcetree.txt

‎data_models/memory_data_models.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -137,7 +137,7 @@ class GainTable:
137137

138138
def __init__(self, data=None, gain: numpy.array = None, time: numpy.array = None, interval=None,
139139
weight: numpy.array = None, residual: numpy.array = None, frequency: numpy.array = None,
140-
receptor_frame: ReceptorFrame = ReceptorFrame("linear"), phasecentre=None):
140+
receptor_frame: ReceptorFrame = ReceptorFrame("linear"), phasecentre=None, configuration=None):
141141
""" Create a gaintable from arrays
142142
143143
The definition of gain is:
@@ -177,6 +177,7 @@ def __init__(self, data=None, gain: numpy.array = None, time: numpy.array = None
177177
self.frequency = frequency
178178
self.receptor_frame = receptor_frame
179179
self.phasecentre = phasecentre
180+
self.configuration = configuration
180181

181182
def size(self):
182183
""" Return size in GB

‎processing_components/calibration/operations.py‎

Lines changed: 22 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -67,8 +67,8 @@ def create_gaintable_from_blockvisibility(vis: BlockVisibility, timeslice=None,
6767
gain_residual = numpy.zeros([ntimes, nfrequency, nrec, nrec])
6868

6969
gt = GainTable(gain=gain, time=gain_time, interval=gain_interval, weight=gain_weight, residual=gain_residual,
70-
frequency=gain_frequency,
71-
receptor_frame=receptor_frame, phasecentre=vis.phasecentre)
70+
frequency=gain_frequency, receptor_frame=receptor_frame, phasecentre=vis.phasecentre,
71+
configuration=vis.configuration)
7272

7373
assert isinstance(gt, GainTable), "gt is not a GainTable: %r" % gt
7474
assert_vis_gt_compatible(vis, gt)
@@ -242,24 +242,37 @@ def qa_gaintable(gt: GainTable, context=None) -> QA:
242242
return QA(origin='qa_gaintable', data=data, context=context)
243243

244244

245-
def gaintable_plot(gt: GainTable, ax, title='', value='amp', **kwargs):
245+
def gaintable_plot(gt: GainTable, ax, title='', value='amp', ants=None, channels=None,
246+
label_max=10, **kwargs):
246247
""" Standard plot of gain table
247248
248249
:param gt: Gaintable
249250
:param ax: matplotlib axes
250251
:param value: 'amp' or 'phase'
252+
:param ants: Antennas to plot
253+
:param channels: Channels to plot
251254
:param kwargs:
252255
:return:
253256
"""
254-
nants = gt.nants
255-
for ant in range(nants):
256-
amp = numpy.abs(gt.gain[:, ant, 0, 0, 0])
257+
if ants is None:
258+
ants = range(gt.nants)
259+
if channels is None:
260+
channels = range(gt.nchan)
261+
for ant in ants:
262+
if gt.configuration is not None:
263+
label = gt.configuration.names[ant]
264+
else:
265+
label = ''
266+
amp = numpy.abs(gt.gain[:, ant, channels, 0, 0])
257267
if value == 'amp':
258-
ax.plot(gt.time[gt.time > 0.0][amp > 0.0], 1.0 / amp[amp > 0.0], '.')
268+
ax.plot(gt.time, amp, '.', label=label)
259269
else:
260-
angle = numpy.angle(gt.gain[:, ant, 0, 0, 0])
261-
ax.plot(gt.time[gt.time > 0.0][amp > 0.0], angle[amp > 0.0], '.')
270+
angle = numpy.angle(gt.gain[:, ant, channels, 0, 0])
271+
ax.plot(gt.time, angle, '.', label=label)
262272

273+
if gt.configuration is not None:
274+
if len(gt.configuration.names) < label_max:
275+
ax.legend()
263276
ax.set_title(title)
264277
ax.set_xlabel('Time (s)')
265278
ax.set_ylabel(value)

‎processing_components/image/operations.py‎

Lines changed: 15 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -124,7 +124,7 @@ def qa_image(im, context="") -> QA:
124124

125125

126126
def show_image(im: Image, fig=None, title: str = '', pol=0, chan=0, cm='Greys', components=None,
127-
vmin=None, vmax=None):
127+
vmin=None, vmax=None, vscale=1.0):
128128
""" Show an Image with coordinates using matplotlib, optionally with components
129129
130130
:param im: Image
@@ -135,38 +135,39 @@ def show_image(im: Image, fig=None, title: str = '', pol=0, chan=0, cm='Greys',
135135
:param components: Optional components
136136
:param vmin: Clip to this minimum
137137
:param vmax: Clip to this maximum
138+
:param vscale: scale max, min by this amount
138139
:return:
139140
"""
140141
import matplotlib.pyplot as plt
141142

142143
assert isinstance(im, Image), im
143-
if not fig:
144-
fig = plt.figure()
145-
plt.clf()
146-
147-
fig.add_subplot(111, projection=im.wcs.sub([1,2]))
148144

145+
fig = plt.figure()
146+
ax = fig.add_subplot(1, 1, 1, projection=im.wcs.sub([1,2]))
147+
149148
if len(im.data.shape) == 4:
150149
data_array = numpy.real(im.data[chan, pol, :, :])
151150
else:
152151
data_array = numpy.real(im.data)
153152

154153
if vmax is None:
155-
vmax = numpy.max(data_array)
154+
vmax = vscale * numpy.max(data_array)
156155
if vmin is None:
157-
vmin = numpy.min(data_array)
156+
vmin = vscale * numpy.min(data_array)
158157

159-
plt.imshow(data_array, origin = 'lower', cmap = cm, vmax=vmax, vmin=vmin)
158+
cm = ax.imshow(data_array, origin = 'lower', cmap = cm, vmax=vmax, vmin=vmin)
160159

161-
plt.xlabel(im.wcs.wcs.ctype[0])
162-
plt.ylabel(im.wcs.wcs.ctype[1])
163-
plt.title(title)
164-
plt.colorbar()
160+
ax.set_xlabel(im.wcs.wcs.ctype[0])
161+
ax.set_ylabel(im.wcs.wcs.ctype[1])
162+
ax.set_title(title)
163+
164+
fig.colorbar(cm, orientation='vertical', shrink=0.7)
165165

166166
if components is not None:
167167
for sc in components:
168168
x, y = skycoord_to_pixel(sc.direction, im.wcs, 0, 'wcs')
169-
plt.plot(x, y, marker='+', color='red')
169+
ax.plot(x, y, marker='+', color='red')
170+
170171
return fig
171172

172173

@@ -213,11 +214,6 @@ def smooth_image(model: Image, width=1.0):
213214
:param width: Kernel in pixels
214215
215216
"""
216-
# TODO: Remove filter when astropy fixes convolve
217-
import warnings
218-
warnings.simplefilter(action='ignore', category=FutureWarning)
219-
import astropy.convolution
220-
221217
assert isinstance(model, Image), model
222218

223219
kernel = astropy.convolution.kernels.Gaussian2DKernel(width)

‎tests/processing_components/test_calibration_plot.py‎

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
""" Unit tests for calibration solution
2+
3+
4+
"""
5+
import logging
6+
import unittest
7+
8+
import astropy.units as u
9+
import matplotlib.pyplot as plt
10+
import numpy
11+
from astropy.coordinates import SkyCoord
12+
13+
from data_models.memory_data_models import Skycomponent
14+
from data_models.polarisation import PolarisationFrame
15+
from processing_components.calibration.operations import create_gaintable_from_blockvisibility, \
16+
gaintable_summary, gaintable_plot
17+
from processing_components.imaging.base import predict_skycomponent_visibility
18+
from processing_components.simulation.configurations import create_named_configuration
19+
from processing_components.simulation.testing_support import simulate_gaintable
20+
from processing_components.visibility.base import create_blockvisibility
21+
22+
log = logging.getLogger(__name__)
23+
24+
25+
class TestCalibrationPlot(unittest.TestCase):
26+
def setUp(self):
27+
numpy.random.seed(180555)
28+
29+
def actualSetup(self, sky_pol_frame='stokesIQUV', data_pol_frame='linear', f=None, vnchan=3, ntimes=30):
30+
self.lowcore = create_named_configuration('LOWBD2', rmax=50.0)
31+
self.times = (numpy.pi / 43200.0) * numpy.linspace(0.0, 30.0, ntimes)
32+
self.frequency = numpy.linspace(1.0e8, 1.1e8, vnchan)
33+
self.channel_bandwidth = numpy.array(vnchan * [self.frequency[1] - self.frequency[0]])
34+
35+
if f is None:
36+
f = [100.0, 50.0, -10.0, 40.0]
37+
38+
if sky_pol_frame == 'stokesI':
39+
f = [100.0]
40+
41+
self.flux = numpy.outer(numpy.array([numpy.power(freq / 1e8, -0.7) for freq in self.frequency]), f)
42+
43+
# The phase centre is absolute and the component is specified relative (for now).
44+
# This means that the component should end up at the position phasecentre+compredirection
45+
self.phasecentre = SkyCoord(ra=+180.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000')
46+
self.compabsdirection = SkyCoord(ra=+181.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000')
47+
self.comp = Skycomponent(direction=self.compabsdirection, frequency=self.frequency, flux=self.flux,
48+
polarisation_frame=PolarisationFrame(sky_pol_frame))
49+
self.vis = create_blockvisibility(self.lowcore, self.times, self.frequency, phasecentre=self.phasecentre,
50+
channel_bandwidth=self.channel_bandwidth, weight=1.0,
51+
polarisation_frame=PolarisationFrame(data_pol_frame))
52+
self.vis = predict_skycomponent_visibility(self.vis, self.comp)
53+
54+
def test_plot_gaintable_scalar(self):
55+
self.actualSetup('stokesI', 'stokesI', f=[100.0])
56+
gt = create_gaintable_from_blockvisibility(self.vis)
57+
log.info("Created gain table: %s" % (gaintable_summary(gt)))
58+
gt = simulate_gaintable(gt, phase_error=0.1, amplitude_error=0.1)
59+
plt.clf()
60+
fig, ax = plt.subplots(1, 1)
61+
gaintable_plot(gt, ax, value='amp')
62+
plt.show(block=False)
63+
fig, ax = plt.subplots(1, 1)
64+
gaintable_plot(gt, ax, value='phase')
65+
plt.show(block=False)
66+
67+
68+
if __name__ == '__main__':
69+
unittest.main()

0 commit comments

Comments
 (0)