Skip to content

Commit 2076fcb

Browse files
timcornwelltimcornwell
committed
First version pointing_simulation notebook
1 parent a19d514 commit 2076fcb

9 files changed

Lines changed: 336 additions & 8 deletions

File tree

‎data_models/memory_data_models.py‎

Lines changed: 17 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -264,8 +264,8 @@ def __init__(self, data=None, pointing: numpy.array = None, time: numpy.array =
264264
nchan = pointing.shape[2]
265265
assert len(frequency) == nchan, "Discrepancy in frequency channels"
266266
desc = [('pointing', '>f16', (nants, nchan, nrec, 2)),
267-
('weight', '>f8', (nants, nchan, nrec, nrec)),
268-
('residual', '>f8', (nchan, nrec, nrec)),
267+
('weight', '>f8', (nants, nchan, nrec, 2)),
268+
('residual', '>f8', (nchan, nrec, 2)),
269269
('time', '>f8'),
270270
('interval', '>f8')]
271271
data = numpy.zeros(shape=[nrows], dtype=desc)
@@ -691,10 +691,25 @@ def __init__(self, image=None, components=None, gaintable=None, mask=None, fixed
691691
"""
692692
if components is None:
693693
components = []
694+
695+
if image is not None:
696+
assert isinstance(image, Image), image
694697
self.image = image
698+
699+
if components is not None:
700+
assert isinstance(components, list)
701+
for comp in components:
702+
assert isinstance(comp, Skycomponent), comp
695703
self.components = [sc for sc in components]
704+
705+
if gaintable is not None:
706+
assert isinstance(gaintable, GainTable), gaintable
696707
self.gaintable = gaintable
708+
709+
if mask is not None:
710+
assert isinstance(mask, Image), mask
697711
self.mask = mask
712+
698713
self.fixed = fixed
699714

700715
def __str__(self):

‎processing_components/image/operations.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -113,6 +113,7 @@ def qa_image(im, context="") -> QA:
113113
data = {'shape': str(im.data.shape),
114114
'max': numpy.max(im.data),
115115
'min': numpy.min(im.data),
116+
'maxabs': numpy.max(numpy.abs(im.data)),
116117
'rms': numpy.std(im.data),
117118
'sum': numpy.sum(im.data),
118119
'medianabs': numpy.median(numpy.abs(im.data)),

‎processing_components/simulation/pointing.py‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -48,10 +48,10 @@ def create_gaintable_from_pointingtable(vis, sc, pt, vp, vis_slices=None, scale=
4848
for ant in range(nant):
4949
worldloc = [float((comp.direction.ra.rad + pointing_ha[0, ant, 0, 0, 0])*r2d) ,
5050
float((comp.direction.dec.rad + pointing_ha[0, ant, 0, 0, 1])*r2d),
51-
1.0, 1e8]
51+
1.0, 1e9]
5252
try:
5353
pixloc = vp.wcs.wcs_world2pix([worldloc], 0)[0].astype('int')
54-
antgain[ant] = scale * vp.data[pixloc[3], pixloc[2], pixloc[1], pixloc[0]]
54+
antgain[ant] = 1.0 / (scale * vp.data[pixloc[3], pixloc[2], pixloc[1], pixloc[0]])
5555
number_good += 1
5656
except:
5757
number_bad += 1
@@ -61,7 +61,7 @@ def create_gaintable_from_pointingtable(vis, sc, pt, vp, vis_slices=None, scale=
6161
gaintables[icomp].phasecentre = comp.direction
6262

6363
if number_bad > 0:
64-
log.warning("create_gaintable_from_screen: %d points are inside the voltage pattern image" % (number_good))
65-
log.warning("create_gaintable_from_screen: %d points are outside the voltage pattern image" % (number_bad))
64+
log.warning("create_gaintable_from_pointingtable: %d points are inside the voltage pattern image" % (number_good))
65+
log.warning("create_gaintable_from_pointingtable: %d points are outside the voltage pattern image" % (number_bad))
6666

6767
return gaintables

‎processing_components/simulation/testing_support.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -924,8 +924,8 @@ def simulate_pointingtable(pt: PointingTable, pointing_error, static_pointing_er
924924
if seed is not None:
925925
numpy.random.seed(seed)
926926

927-
log.debug("simulate_pointingtable: Simulating pointing error = %.4f (rad)"
928-
% (pointing_error))
927+
log.debug("simulate_pointingtable: Simulating pointing error = %.6f (rad), static pointing error = %.6f (rad)"
928+
% (pointing_error, static_pointing_error))
929929

930930
ntimes, nant, nchan, nrec, _ = pt.data['pointing'].shape
931931
if pointing_error > 0.0:

‎workflows/notebooks/pointing_simulation.ipynb‎

Lines changed: 302 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,302 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"metadata": {
7+
"ExecuteTime": {
8+
"end_time": "2019-04-06T14:48:31.808478Z",
9+
"start_time": "2019-04-06T14:48:30.637288Z"
10+
},
11+
"pycharm": {},
12+
"scrolled": false
13+
},
14+
"outputs": [],
15+
"source": [
16+
"\n",
17+
"import sys\n",
18+
"%matplotlib inline\n",
19+
"\n",
20+
"from data_models.parameters import arl_path\n",
21+
"\n",
22+
"results_dir = arl_path('test_results')\n",
23+
"\n",
24+
"import numpy\n",
25+
"\n",
26+
"from astropy.coordinates import SkyCoord\n",
27+
"from astropy import units as u\n",
28+
"\n",
29+
"from matplotlib import pyplot as plt\n",
30+
"\n",
31+
"from data_models.polarisation import PolarisationFrame\n",
32+
"from data_models.memory_data_models import Skycomponent, SkyModel\n",
33+
"\n",
34+
"from wrappers.serial.visibility.base import create_blockvisibility\n",
35+
"from wrappers.serial.image.operations import show_image, qa_image\n",
36+
"from wrappers.serial.simulation.testing_support import create_named_configuration, simulate_pointingtable\n",
37+
"from wrappers.serial.imaging.primary_beams import create_vp, create_pb\n",
38+
"from wrappers.serial.skycomponent.operations import apply_beam_to_skycomponent\n",
39+
"from wrappers.serial.imaging.base import create_image_from_visibility, advise_wide_field\n",
40+
"from processing_components.calibration.pointing import create_pointingtable_from_blockvisibility\n",
41+
"from processing_components.simulation.pointing import create_gaintable_from_pointingtable\n",
42+
"from wrappers.arlexecute.visibility.base import copy_visibility\n",
43+
"\n",
44+
"\n",
45+
"from wrappers.arlexecute.visibility.coalesce import convert_blockvisibility_to_visibility\n",
46+
"\n",
47+
"from workflows.arlexecute.skymodel.skymodel_arlexecute import predict_skymodel_list_arlexecute_workflow\n",
48+
"from workflows.arlexecute.imaging.imaging_arlexecute import invert_list_arlexecute_workflow\n",
49+
"\n",
50+
"from wrappers.arlexecute.execution_support.arlexecute import arlexecute\n",
51+
"\n",
52+
"import logging\n",
53+
"log = logging.getLogger()\n",
54+
"log.setLevel(logging.DEBUG)\n",
55+
"log.addHandler(logging.StreamHandler(sys.stdout))\n"
56+
]
57+
},
58+
{
59+
"cell_type": "code",
60+
"execution_count": null,
61+
"metadata": {
62+
"ExecuteTime": {
63+
"end_time": "2019-04-06T14:48:34.011538Z",
64+
"start_time": "2019-04-06T14:48:31.811472Z"
65+
}
66+
},
67+
"outputs": [],
68+
"source": [
69+
"arlexecute.set_client(use_dask=True)"
70+
]
71+
},
72+
{
73+
"cell_type": "code",
74+
"execution_count": null,
75+
"metadata": {
76+
"ExecuteTime": {
77+
"end_time": "2019-04-06T14:48:34.122482Z",
78+
"start_time": "2019-04-06T14:48:34.014916Z"
79+
},
80+
"pycharm": {}
81+
},
82+
"outputs": [],
83+
"source": [
84+
"nfreqwin = 1\n",
85+
"ntimes = 1\n",
86+
"rmax = 1000.0\n",
87+
"if nfreqwin > 1:\n",
88+
" frequency = numpy.linspace(0.8e9, 1.2e9, nfreqwin)\n",
89+
" channel_bandwidth = numpy.array(nfreqwin * [frequency[1] - frequency[0]])\n",
90+
"else:\n",
91+
" frequency = [1e9]\n",
92+
" channel_bandwidth = [0.4e9]\n",
93+
"times = numpy.linspace(-numpy.pi / 4.0, numpy.pi / 4.0, ntimes)\n",
94+
"\n",
95+
"phasecentre = SkyCoord(ra=+15.0 * u.deg, dec=-45.0 * u.deg, frame='icrs', equinox='J2000')\n",
96+
"midcore = create_named_configuration('MIDR5', rmax=rmax)\n",
97+
"\n",
98+
"block_vis = create_blockvisibility(midcore, times, frequency=frequency,\n",
99+
" channel_bandwidth=channel_bandwidth, weight=1.0, phasecentre=phasecentre,\n",
100+
" polarisation_frame=PolarisationFrame(\"stokesI\"), zerow=True)\n"
101+
]
102+
},
103+
{
104+
"cell_type": "code",
105+
"execution_count": null,
106+
"metadata": {
107+
"ExecuteTime": {
108+
"end_time": "2019-04-06T14:48:34.218959Z",
109+
"start_time": "2019-04-06T14:48:34.124530Z"
110+
},
111+
"pycharm": {}
112+
},
113+
"outputs": [],
114+
"source": [
115+
"wprojection_planes=1\n",
116+
"vis = convert_blockvisibility_to_visibility(block_vis)\n",
117+
"advice=advise_wide_field(vis, guard_band_image=5.0, delA=0.02, \n",
118+
" wprojection_planes=wprojection_planes)\n",
119+
"\n",
120+
"vis_slices = advice['vis_slices']\n",
121+
"npixel=16384\n",
122+
"cellsize=advice['cellsize']\n",
123+
"pb_cellsize=advice['primary_beam_fov']/npixel"
124+
]
125+
},
126+
{
127+
"cell_type": "code",
128+
"execution_count": null,
129+
"metadata": {
130+
"ExecuteTime": {
131+
"end_time": "2019-04-06T14:48:34.231069Z",
132+
"start_time": "2019-04-06T14:48:34.222902Z"
133+
}
134+
},
135+
"outputs": [],
136+
"source": [
137+
"from matplotlib import pylab\n",
138+
"\n",
139+
"pylab.rcParams['figure.figsize'] = (12.0, 12.0)\n",
140+
"pylab.rcParams['image.cmap'] = 'rainbow'\n",
141+
"doplot=False"
142+
]
143+
},
144+
{
145+
"cell_type": "code",
146+
"execution_count": null,
147+
"metadata": {
148+
"ExecuteTime": {
149+
"end_time": "2019-04-06T14:51:03.403059Z",
150+
"start_time": "2019-04-06T14:48:34.234002Z"
151+
}
152+
},
153+
"outputs": [],
154+
"source": [
155+
"original_component = [Skycomponent(flux=[[1.0]], direction=phasecentre, \n",
156+
" frequency=frequency, polarisation_frame=PolarisationFrame('stokesI'))]\n",
157+
"vp = create_image_from_visibility(block_vis, npixel=npixel, frequency=frequency,\n",
158+
" nchan=nfreqwin, cellsize=pb_cellsize, phasecentre=phasecentre)\n",
159+
"model = create_image_from_visibility(block_vis, npixel=512, frequency=frequency,\n",
160+
" nchan=nfreqwin, cellsize=cellsize, phasecentre=phasecentre)\n",
161+
"if doplot:\n",
162+
" pb = create_pb(vp, 'MID', pointingcentre=pb_direction)\n",
163+
" show_image(pb)\n",
164+
" \n",
165+
"offset = [0.595, 0.0]\n",
166+
"pb_direction = SkyCoord(ra=(+15.0 + offset[0]) * u.deg, \n",
167+
" dec=(-45.0 + offset[1]) * u.deg, frame='icrs', equinox='J2000')\n",
168+
"\n",
169+
"vp = create_vp(vp, 'MID', pointingcentre=pb_direction)\n",
170+
"pt = create_pointingtable_from_blockvisibility(block_vis, vp)\n",
171+
"\n",
172+
"no_error_pt = simulate_pointingtable(pt, 0.0, 0.0)\n",
173+
"no_error_gt = create_gaintable_from_pointingtable(block_vis, original_component, no_error_pt, vp)\n",
174+
"\n",
175+
"no_error_sm=SkyModel(components=original_component, gaintable=no_error_gt[0])\n",
176+
"no_error_vis = copy_visibility(vis)\n",
177+
"no_error_vis = predict_skymodel_list_arlexecute_workflow(no_error_vis, [no_error_sm], context='2d', docal=True)\n",
178+
"no_error_vis=arlexecute.compute(no_error_vis, sync=True)[0]\n"
179+
]
180+
},
181+
{
182+
"cell_type": "code",
183+
"execution_count": null,
184+
"metadata": {
185+
"ExecuteTime": {
186+
"end_time": "2019-04-06T14:51:03.415162Z",
187+
"start_time": "2019-04-06T14:51:03.406050Z"
188+
}
189+
},
190+
"outputs": [],
191+
"source": [
192+
"print(numpy.average(no_error_vis.vis))"
193+
]
194+
},
195+
{
196+
"cell_type": "markdown",
197+
"metadata": {
198+
"pycharm": {}
199+
},
200+
"source": [
201+
"Generate the model from the GLEAM catalog, including application of the primary beam."
202+
]
203+
},
204+
{
205+
"cell_type": "code",
206+
"execution_count": null,
207+
"metadata": {
208+
"ExecuteTime": {
209+
"end_time": "2019-04-06T14:51:11.148096Z",
210+
"start_time": "2019-04-06T14:51:03.418008Z"
211+
},
212+
"pycharm": {},
213+
"scrolled": false
214+
},
215+
"outputs": [],
216+
"source": [
217+
"\n",
218+
"fwhm = 1.0\n",
219+
"pes = [1.0, 3.0, 10.0, 30.0, 100.0, 300.0, 1000.0]\n",
220+
"qas = []\n",
221+
"for pe in pes:\n",
222+
" pointing_error = pe * numpy.pi/(3600.0*180.0)\n",
223+
" static_pointing_error = 0.0*numpy.pi/(3600.0*180.0)\n",
224+
" error_pt = simulate_pointingtable(pt, pointing_error=pointing_error, static_pointing_error=static_pointing_error)\n",
225+
" error_gt = create_gaintable_from_pointingtable(block_vis, original_component, error_pt, vp)\n",
226+
"\n",
227+
" error_sm=SkyModel(components=original_component, gaintable=error_gt[0])\n",
228+
"\n",
229+
" error_vis = copy_visibility(vis)\n",
230+
" error_vis = predict_skymodel_list_arlexecute_workflow(error_vis, [error_sm], context='2d', docal=True)\n",
231+
" error_vis=arlexecute.compute(error_vis, sync=True)[0]\n",
232+
"\n",
233+
" error_vis.data['vis']-=no_error_vis.data['vis']\n",
234+
" \n",
235+
" dirty = invert_list_arlexecute_workflow([error_vis], [model], '2d', padding=4)\n",
236+
" dirty, sumwt = arlexecute.compute(dirty, sync=True)[0]\n",
237+
" bpix=129\n",
238+
" epix=384\n",
239+
" print(bpix,epix)\n",
240+
" print(dirty.data.shape)\n",
241+
"\n",
242+
" dirty.data = dirty.data[:,:,bpix:epix,bpix:epix]\n",
243+
" qa=qa_image(dirty)\n",
244+
" print(qa)\n",
245+
" qas.append(qa)\n",
246+
" if True:\n",
247+
" show_image(dirty)\n",
248+
" plt.show()\n"
249+
]
250+
},
251+
{
252+
"cell_type": "code",
253+
"execution_count": null,
254+
"metadata": {
255+
"ExecuteTime": {
256+
"end_time": "2019-04-06T14:52:40.289627Z",
257+
"start_time": "2019-04-06T14:52:39.448339Z"
258+
},
259+
"scrolled": false
260+
},
261+
"outputs": [],
262+
"source": [
263+
"plt.clf()\n",
264+
"for field in ['maxabs', 'rms', 'medianabs']:\n",
265+
" plt.loglog(pes, [q.data[field] for q in qas], '.', label=field)\n",
266+
"plt.xlabel('Pointing error (arcsec)')\n",
267+
"plt.ylabel('Error (Jy)')\n",
268+
"plt.title('Error for 1Jy point source at half power point of MID primary beam')\n",
269+
"plt.legend()\n",
270+
"plt.show()"
271+
]
272+
},
273+
{
274+
"cell_type": "code",
275+
"execution_count": null,
276+
"metadata": {},
277+
"outputs": [],
278+
"source": []
279+
}
280+
],
281+
"metadata": {
282+
"kernelspec": {
283+
"display_name": "Python 3",
284+
"language": "python",
285+
"name": "python3"
286+
},
287+
"language_info": {
288+
"codemirror_mode": {
289+
"name": "ipython",
290+
"version": 3
291+
},
292+
"file_extension": ".py",
293+
"mimetype": "text/x-python",
294+
"name": "python",
295+
"nbconvert_exporter": "python",
296+
"pygments_lexer": "ipython3",
297+
"version": "3.6.1"
298+
}
299+
},
300+
"nbformat": 4,
301+
"nbformat_minor": 2
302+
}

‎wrappers/arlexecute/imaging/primary_beams.py‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,9 @@
33
"""
44
from processing_components.imaging.primary_beams import ft_disk
55
from processing_components.imaging.primary_beams import create_pb
6+
from processing_components.imaging.primary_beams import create_vp
67
from processing_components.imaging.primary_beams import mosaic_pb
78
from processing_components.imaging.primary_beams import create_pb_generic
9+
from processing_components.imaging.primary_beams import create_vp_generic
810
from processing_components.imaging.primary_beams import create_low_test_beam
11+
from processing_components.imaging.primary_beams import create_low_test_vp

0 commit comments

Comments
 (0)