This script makes a fake data set, fills it with a number of point components, and then images it using a variety of algorithms. See imaging-fits for a similar notebook that checks for errors in the recovered properties of the images.
The measurement equation for a wide field of view interferometer is:
$$V(u,v,w) =\int \frac{I(l,m)}{\sqrt{1-l^2-m^2}} e^{-2 \pi j (ul+um + w(\sqrt{1-l^2-m^2}-1))} dl dm$$We will show various algorithms for computing approximations to this integral. Calculation of the visibility V from the sky brightness I is called predict, and the inverese is called invert.
In [ ]:
%matplotlib inline
import os
import sys
sys.path.append(os.path.join('..', '..'))
from data_models.parameters import arl_path
results_dir = arl_path('test_results')
from matplotlib import pylab
import numpy
from astropy.coordinates import SkyCoord
from astropy import units as u
from astropy.wcs.utils import pixel_to_skycoord
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from data_models.polarisation import PolarisationFrame
from wrappers.serial.image.iterators import image_raster_iter
from processing_library.image.operations import create_w_term_like
# Use serial wrappers by default
from wrappers.serial.visibility.base import create_visibility, create_visibility, create_visibility_from_rows
from wrappers.serial.skycomponent.operations import create_skycomponent
from wrappers.serial.image.operations import show_image, export_image_to_fits
from wrappers.serial.visibility.iterators import vis_timeslice_iter
from wrappers.serial.simulation.configurations import create_named_configuration
from wrappers.serial.imaging.base import invert_2d, create_image_from_visibility, \
predict_skycomponent_visibility, advise_wide_field
from wrappers.serial.visibility.iterators import vis_timeslice_iter
from wrappers.serial.imaging.weighting import weight_visibility
from wrappers.serial.visibility.iterators import vis_timeslices
from wrappers.arlexecute.griddata.kernels import create_awterm_convolutionfunction
from wrappers.arlexecute.griddata.convolution_functions import apply_bounding_box_convolutionfunction
# Use arlexecute for imaging
from wrappers.arlexecute.execution_support.arlexecute import arlexecute
from workflows.arlexecute.imaging.imaging_arlexecute import invert_list_arlexecute_workflow
import logging
log = logging.getLogger()
log.setLevel(logging.DEBUG)
log.addHandler(logging.StreamHandler(sys.stdout))
doplot = True
In [ ]:
pylab.rcParams['figure.figsize'] = (12.0, 12.0)
pylab.rcParams['image.cmap'] = 'rainbow'
Construct the SKA1-LOW core configuration
In [ ]:
lowcore = create_named_configuration('LOWBD2-CORE')
We create the visibility.
This just makes the uvw, time, antenna1, antenna2, weight columns in a table
In [ ]:
times = numpy.array([-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0]) * (numpy.pi / 12.0)
frequency = numpy.array([1e8])
channel_bandwidth = numpy.array([1e7])
reffrequency = numpy.max(frequency)
phasecentre = SkyCoord(ra=+15.0 * u.deg, dec=-45.0 * u.deg, frame='icrs', equinox='J2000')
vt = create_visibility(lowcore, times, frequency, channel_bandwidth=channel_bandwidth,
weight=1.0, phasecentre=phasecentre, polarisation_frame=PolarisationFrame("stokesI"))
Advise on wide field parameters. This returns a dictionary with all the input and calculated variables.
In [ ]:
advice = advise_wide_field(vt, wprojection_planes=1)
Plot the synthesized UV coverage.
In [ ]:
if doplot:
plt.clf()
plt.plot(vt.data['uvw'][:, 0], vt.data['uvw'][:, 1], '.', color='b')
plt.plot(-vt.data['uvw'][:, 0], -vt.data['uvw'][:, 1], '.', color='r')
plt.xlabel('U (wavelengths)')
plt.ylabel('V (wavelengths)')
plt.show()
plt.clf()
plt.plot(vt.data['uvw'][:, 0], vt.data['uvw'][:, 2], '.', color='b')
plt.xlabel('U (wavelengths)')
plt.ylabel('W (wavelengths)')
plt.show()
plt.clf()
plt.plot(vt.data['time'][vt.u>0.0], vt.data['uvw'][:, 2][vt.u>0.0], '.', color='b')
plt.plot(vt.data['time'][vt.u<=0.0], vt.data['uvw'][:, 2][vt.u<=0.0], '.', color='r')
plt.xlabel('U (wavelengths)')
plt.ylabel('W (wavelengths)')
plt.show()
plt.clf()
n, bins, patches = plt.hist(vt.w, 50, normed=1, facecolor='green', alpha=0.75)
plt.xlabel('W (wavelengths)')
plt.ylabel('Count')
plt.show()
Show the planar nature of the uvw sampling, rotating with hour angle
Create a grid of components and predict each in turn, using the full phase term including w.
In [ ]:
npixel = 512
cellsize=0.001
facets = 4
flux = numpy.array([[100.0]])
vt.data['vis'] *= 0.0
model = create_image_from_visibility(vt, npixel=512, cellsize=0.001, npol=1)
spacing_pixels = npixel // facets
log.info('Spacing in pixels = %s' % spacing_pixels)
spacing = 180.0 * cellsize * spacing_pixels / numpy.pi
centers = -1.5, -0.5, +0.5, +1.5
comps=list()
for iy in centers:
for ix in centers:
pra = int(round(npixel // 2 + ix * spacing_pixels - 1))
pdec = int(round(npixel // 2 + iy * spacing_pixels - 1))
sc = pixel_to_skycoord(pra, pdec, model.wcs)
log.info("Component at (%f, %f) %s" % (pra, pdec, str(sc)))
comp = create_skycomponent(flux=flux, frequency=frequency, direction=sc,
polarisation_frame=PolarisationFrame("stokesI"))
comps.append(comp)
predict_skycomponent_visibility(vt, comps)
Make the dirty image and point spread function using the two-dimensional approximation:
$$V(u,v,w) =\int I(l,m) e^{2 \pi j (ul+um)} dl dm$$Note that the shape of the sources vary with position in the image. This space-variant property of the PSF arises from the w-term neglected in the two-dimensional invert.
In [ ]:
arlexecute.set_client(use_dask=True)
In [ ]:
dirty = create_image_from_visibility(vt, npixel=512, cellsize=0.001,
polarisation_frame=PolarisationFrame("stokesI"))
vt = weight_visibility(vt, dirty)
future = invert_list_arlexecute_workflow([vt], [dirty], context='2d')
dirty, sumwt = arlexecute.compute(future, sync=True)[0]
if doplot:
show_image(dirty)
print("Max, min in dirty image = %.6f, %.6f, sumwt = %f" % (dirty.data.max(), dirty.data.min(), sumwt))
export_image_to_fits(dirty, '%s/imaging-wterm_dirty.fits' % (results_dir))
This occurs because the Fourier transform relationship between sky brightness and visibility is only accurate over small fields of view.
Hence we can make an accurate image by partitioning the image plane into small regions, treating each separately and then glueing the resulting partitions into one image. We call this image plane partitioning image plane faceting.
$$V(u,v,w) = \sum_{i,j} \frac{1}{\sqrt{1- l_{i,j}^2- m_{i,j}^2}} e^{-2 \pi j (ul_{i,j}+um_{i,j} + w(\sqrt{1-l_{i,j}^2-m_{i,j}^2}-1))} \int I(\Delta l, \Delta m) e^{-2 \pi j (u\Delta l_{i,j}+u \Delta m_{i,j})} dl dm$$
In [ ]:
dirtyFacet = create_image_from_visibility(vt, npixel=512, cellsize=0.001, npol=1)
future = invert_list_arlexecute_workflow([vt], [dirtyFacet], facets=4, context='facets')
dirtyFacet, sumwt = arlexecute.compute(future, sync=True)[0]
if doplot:
show_image(dirtyFacet)
print("Max, min in dirty image = %.6f, %.6f, sumwt = %f" % (dirtyFacet.data.max(), dirtyFacet.data.min(), sumwt))
export_image_to_fits(dirtyFacet, '%s/imaging-wterm_dirtyFacet.fits' % (results_dir))
That was the best case. This time, we will not arrange for the partitions to be centred on the sources.
In [ ]:
dirtyFacet2 = create_image_from_visibility(vt, npixel=512, cellsize=0.001, npol=1)
future = invert_list_arlexecute_workflow([vt], [dirtyFacet2], facets=2, context='facets')
dirtyFacet2, sumwt = arlexecute.compute(future, sync=True)[0]
if doplot:
show_image(dirtyFacet2)
print("Max, min in dirty image = %.6f, %.6f, sumwt = %f" % (dirtyFacet2.data.max(), dirtyFacet2.data.min(), sumwt))
export_image_to_fits(dirtyFacet2, '%s/imaging-wterm_dirtyFacet2.fits' % (results_dir))
Another approach is to partition the visibility data by slices in w. The measurement equation is approximated as:
$$V(u,v,w) =\sum_i \int \frac{ I(l,m) e^{-2 \pi j (w_i(\sqrt{1-l^2-m^2}-1))})}{\sqrt{1-l^2-m^2}} e^{-2 \pi j (ul+um)} dl dm$$If images constructed from slices in w are added after applying a w-dependent image plane correction, the w term will be corrected.
The w-dependent w-beam is:
In [ ]:
if doplot:
wterm = create_w_term_like(model, phasecentre=vt.phasecentre, w=numpy.max(vt.w))
show_image(wterm)
plt.show()
In [ ]:
dirtywstack = create_image_from_visibility(vt, npixel=512, cellsize=0.001, npol=1)
future = invert_list_arlexecute_workflow([vt], [dirtywstack], vis_slices=101, context='wstack')
dirtywstack, sumwt = arlexecute.compute(future, sync=True)[0]
show_image(dirtywstack)
plt.show()
print("Max, min in dirty image = %.6f, %.6f, sumwt = %f" %
(dirtywstack.data.max(), dirtywstack.data.min(), sumwt))
export_image_to_fits(dirtywstack, '%s/imaging-wterm_dirty_wstack.fits' % (results_dir))
The w-term can also be viewed as a time-variable distortion. Approximating the array as instantaneously co-planar, we have that w can be expressed in terms of $u,v$
$$w = a u + b v$$Transforming to a new coordinate system:
$$ l' = l + a (\sqrt{1-l^2-m^2}-1))$$$$ m' = m + b (\sqrt{1-l^2-m^2}-1))$$Ignoring changes in the normalisation term, we have:
$$V(u,v,w) =\int \frac{I(l',m')}{\sqrt{1-l'^2-m'^2}} e^{-2 \pi j (ul'+um')} dl' dm'$$To illustrate this, we will construct images as a function of time. For comparison, we show difference of each time slice from the best facet image. Instantaneously the sources are un-distorted but do lie in the wrong location.
In [ ]:
for rows in vis_timeslice_iter(vt):
visslice = create_visibility_from_rows(vt, rows)
dirtySnapshot = create_image_from_visibility(visslice, npixel=512, cellsize=0.001, npol=1, compress_factor=0.0)
future = invert_list_arlexecute_workflow([visslice], [dirtySnapshot], context='2d')
dirtySnapshot, sumwt = arlexecute.compute(future, sync=True)[0]
print("Max, min in dirty image = %.6f, %.6f, sumwt = %f" %
(dirtySnapshot.data.max(), dirtySnapshot.data.min(), sumwt))
if doplot:
dirtySnapshot.data -= dirtyFacet.data
show_image(dirtySnapshot)
plt.title("Hour angle %.2f hours" % (numpy.average(visslice.time) * 12.0 / 43200.0))
plt.show()
This timeslice imaging leads to a straightforward algorithm in which we correct each time slice and then sum the resulting timeslices.
In [ ]:
dirtyTimeslice = create_image_from_visibility(vt, npixel=512, cellsize=0.001, npol=1)
future = invert_list_arlexecute_workflow([vt], [dirtyTimeslice], vis_slices=vis_timeslices(vt, 'auto'),
padding=2, context='timeslice')
dirtyTimeslice, sumwt = arlexecute.compute(future, sync=True)[0]
show_image(dirtyTimeslice)
plt.show()
print("Max, min in dirty image = %.6f, %.6f, sumwt = %f" %
(dirtyTimeslice.data.max(), dirtyTimeslice.data.min(), sumwt))
export_image_to_fits(dirtyTimeslice, '%s/imaging-wterm_dirty_Timeslice.fits' % (results_dir))
Finally we try w-projection. For a fixed w, the measurement equation can be stated as as a convolution in Fourier space.
$$V(u,v,w) =G_w(u,v) \ast \int \frac{I(l,m)}{\sqrt{1-l^2-m^2}} e^{-2 \pi j (ul+um)} dl dm$$where the convolution function is:
$$G_w(u,v) = \int \frac{1}{\sqrt{1-l^2-m^2}} e^{-2 \pi j (ul+um + w(\sqrt{1-l^2-m^2}-1))} dl dm$$Hence when gridding, we can use the transform of the w beam to correct this effect while gridding.
In [ ]:
dirtyWProjection = create_image_from_visibility(vt, npixel=512, cellsize=0.001, npol=1)
gcfcf = create_awterm_convolutionfunction(model, nw=101, wstep=800.0/101, oversampling=8,
support=60,
use_aaf=True)
future = invert_list_arlexecute_workflow([vt], [dirtyWProjection], context='2d', gcfcf=[gcfcf])
dirtyWProjection, sumwt = arlexecute.compute(future, sync=True)[0]
if doplot:
show_image(dirtyWProjection)
print("Max, min in dirty image = %.6f, %.6f, sumwt = %f" % (dirtyWProjection.data.max(),
dirtyWProjection.data.min(), sumwt))
export_image_to_fits(dirtyWProjection, '%s/imaging-wterm_dirty_WProjection.fits' % (results_dir))
In [ ]: