SPAM tutorial

This tutorial was created for the 3GC3 workshop "The Elephants In The Room" in Port Alfred, South Africa, in February 2013.

Introduction

This wiki page contains a tutorial that aims to demonstrate how a direction-dependent (DD) ionospheric phase calibration can be performed following a particular scheme called Source Peeling and Atmospheric Modeling (SPAM). It addresses one of the key issues for high-resolution (HR; ≥ few tens of km baselines, ≤30" resolution), low-frequency (LF; ≤ few hundred MHz) radio interferometric observations. To first order, the ionosphere adds time-variable, DD phase errors to the visibilities. Phase errors can vary significantly over the angular extend of the wide field-of-view (FOV; ≥ few degrees), causing radio sources to appear shifted in position and/or distorted in shape. Since traditional self-calibration (SC) will determine only one phase correction per antenna per time interval for the whole FOV, it fails in the presence of DD ionospheric phase errors. SPAM attempts to measure the ionospheric phase errors towards the apparent brightest sources in the FOV, fit a phase screen model to these measurements, and image the full FOV while applying the DD gain phase corrections from the phase screen model.

This tutorial attempts to take the user through some of the essential steps in SPAM processing: measuring, modeling and applying the ionospheric phases. These steps are generic to any type of ionospheric calibration scheme. SPAM is a python-based set of scripts that accesses AIPS functionality and files through the ParselTongue interface (itself based on ObitTalk). A 2009 paper describing the method can be found here. A link to some SPAM silliness can be found here.

In SPAM, several choices are were on how these steps are implemented. Please keep in mind that building this particular implementation was started in 2005 with some of the better software that was available at the time. With the significant increase in development of algorithms and software in the area of DD calibration and imaging over the last few years, it is now possible to do this more efficiently (and possibly more accurately). Still, SPAM served the purpose of demonstrating that, in principle, improvements of the image quality can be obtained by following this strategy. Over the last few years, SPAM has been successfully applied to many visibility data sets from the Giant Metrewave Radio telescope ( GMRT) at various frequencies (150, 240, 330 & 610 MHz), and several data sets from the Very Large Array ( VLA) at 74 and 330 MHz.

A small example visibility data set was prepared that (in principle) will allow the participant to work through this tutorial within the given time. These data, one hour of GMRT 153 MHz data on the merging galaxy cluster CIZA J2242.8+5301 (containing the 'sausage' relic; van Weeren et al. 2011), were part of a longer observation in November 2009 during moderately disturbed ionospheric conditions.

Fig. 1: Combined 610 MHz radio (red) and X-ray (blue) image of the galaxy cluster CIZA J2242.8+5301.

Accessing the software

For convenience, we will run an AIPS session and a parseltongue (PT) session in parallel, both accessing the same data areas. All necessary software is installed on both JAKE and ELWOOD, so you can work on either machine. Please know that for software reasons the number of (collaborating groups of) participants per machine is LIMITED TO EIGHT.

From your own computer (running an X-server), open a Linux shell:
ssh -X <username>@jake.ru.ac.za

After logging in, check if you're part of the group 'aipsusers' by typing 'groups' in the shell. If not, go see your sysadmin.

Copy the AIPS and PT startup scripts and some other stuff to your local bin directory:
cd /home/<username>
mkdir bin
cp /home/intema/bin/* bin

If the local bin directory did not exist before, to add this new directory to your path, close and re-open your Linux shell.

Create a project directory where we'll be performing the tutorial:
mkdir <your project directory>

Like above, open a second Linux shell on JAKE (the AIPS shell from now on) and execute the following script:
start_aips.sh /home/<username>/<your project directory> 

This will start up AIPS POPS (the ancient AIPS user interface), an AIPS TV window, and an AIPS message window. When asked for a user ID in AIPS POPS, enter your favorite random number between 100 and 9999, but do remember it!

Next, open a third Linux shell (the PT shell from now on) and execute the following script:
start_parseltongue.sh /home/<username>/<your project directory> <AIPS user ID>

This starts up PT in a regular Python shell. PT will ask you for your AIPS user ID again, then return with the Python prompt ">>>". SPAM is a module that is loaded into python. To do this, in the PT shell enter:
from spam import *

If all went well, Python will return with its prompt without complaints. If you go back to the regular Linux shell (the first one), you can browse (ls) the contents of the project directory. The scripts that start up AIPS and PT have created a whole series of sub-directories for convenience. In our case we'll only be using a subset of these directories.

If you want to leave AIPS, type 'kleenex'. If you want to leave PT, type 'exit()' or <CTRL>-D. Better not exit AIPS while PT is still running, since the AIPS exit cleans up files that PT still likes to see.

Loading the visibility data

The visibility data set we'll be working on is pre-processed (calibrated & flagged, and target field field data split off), and has a self-calibration gain table attached to it.
  • Download the data (CIZA.SMALL.UV.FITS) from here and put it in the 'fits' sub-directory.
  • Also download the accompanying AIPS boxfile (CIZA.SMALUV.BOX; explanation will follow) from here and put it in the 'datfil' sub-directory.
  • NOTE: These files are also available locally in /home/intema/data
Next, load the data into the AIPS filesystem by copying the following lines into the PT shell:
aips_disk = 1
target_name = 'CIZA'
uv_file_name = 'fits/CIZA.SMALL.UV.FITS'
uv = get_aips_file( aips_disk, target_name, 'UVDATA', -1, 'UV' )
read_fits_uv( uv_file_name, uv )

The function get_aips_file() creates a new, unique reference to a UV file on AIPS disk 1. The function read_fits_uv() is a wrapper around the AIPS task FITLD, which reads in visibility data in the UVFITS format. The generous output of all AIPS tasks is printed in the same PT shell. This output will become enormous pretty quickly. There are ways of redirecting this, but we'll leave it for now.

The data is now present on AIPS disk 1, which is physically located in sub-directory 'work01'. To see this file in AIPS, go to the AIPS shell and type:
indisk 1
pcat
getname 1
imheader

Like in the PT shell, the AIPS shell now shows the basic structure of the visibility data set, which is a single-source data set (target name = CIZA), has 1 stokes (stokes I), and 1 IF (spectral window) consisting of 14 frequency channels of 0.39 MHz wide. The channels are defined with a negative offset, going down from 155 MHz. Without having to do the math yourself, the central frequency and bandwidth (in MHz) can be obtained in the PT shell:
get_central_frequency( uv ) / 1.e6
get_bandwidth( uv ) / 1.e6

So we have 5.5 MHz centered on 153 MHz. This means we don't have to worry about wide-bandwidth effects. The visibility data we loaded already has a self-calibration gain phase table (SN table in AIPS) attached to it.

For purposes of ionospheric calibration, it is important to know that the (time- and direction-) constant instrumental phase contribution to the visibilities has been estimated and removed. At no point in the pre-processing have any time-variable gain phase solutions been applied, since this would seriously hamper the ionospheric phase modeling we intend to do. The gain phase solutions in any calibration table that we'll determine and look at is therefore assumed to contain a pure ionospheric component only.

When attempting a traditional direction-independent phase calibration in the presence of DD ionospheric phase errors, the self-calibration returns a flux-weighted average phase correction per antenna for the whole FOV. These corrections may be approximately correct for the apparent brightest source in the field, but deteriorate when moving away from this source.

Imaging the data using the self-calibration gain table

The default way to do wide-field imaging in AIPS (task IMAGR) is to tile the FOV up into small patches (facets) for which the w-term is approximately constant. For large arrays and large FOVs, the necessary number of required facets can become quite large (few hundred), generating some bookkeeping challenges. SPAM attempts to keep much of this and other bookkeeping issues away from the user. The available imaging functions control the AIPS IMAGR task by asking for higher-level inputs, and use these and some information derived from the visibility data directly to feed suitable imaging parameters to IMAGR.

The boxfile 'CIZA.SMALUV.BOX' that was put in the 'datfil' directory contains the definition of 73 facets that cover the primary beam (PB) area, as well as 17 additional facets centered on bright sources in and near the PB area. Centering bright point sources on pixel centers helps reduce some of the limitations of CLEAN (e.g., see Cotton & Uson 2008), while including nearby bright outliers reduces sidelobe noise in the main image. The same boxfile also contains the definition of many more CLEAN boxes that are to be used during imaging.

To get a (dirty) preview and noise estimate of our central image, copy the following in the PT shell, either command by command or all at once:
imagr_params = { 'robust' : -1. }
cpb_facets = image_cpb_facets( uv, apply_solutions = True, imagr_params = imagr_params )
cpb_noise = measure_cpb_noise( uv, cpb_facets, keep_image = True )
remove_facets( cpb_facets )

The Python dictionary 'imagr_params' is used to feed some parameters directly to the AIPS task IMAGR. In this case, we specify the weighting scheme to be between robust and uniform. This puts somewhat more weight on longer baselines, since GMRT data has a tendency to be over-represented by the short (central square km) baselines. The function image_cpb_facets() uses AIPS IMAGR to make dirty images of the central 7 facets in the FOV, covering roughly 1/4th of the full field diameter at the most sensitive (and crowded) part of the primary beam. These same 7 central facets are used throughout the processing to measure the background RMS noise. In my opinion, the center of the image is a much more sensible place to measure noise than somewhere far from the field center. The function measure_cpb_noise() stitches the 7 facets together into one image and robustly measures the background RMS noise. Afterwards, remove_facets() deletes the individual facets, but the stitched image is saved.

The measured central noise in the dirty image is around 12.5 mJy/beam for a 28" x 18" beam. To see the image, copy this into the PT shell:
cpb_image = get_aips_file( 1, 'CPB', 'FLATN', 0, 'MA' )
plot_image( cpb_image, z_range = [ -0.01,0.1 ] )

You can also look at this image in AIPS in greyscale by typing (in the AIPS shell):
pcat
getname <number in Cat column of CPB.FLATN.100>
pixrange -0.01, 0.1
tvini
tvlod

The sausage relic is there, really! Obviously, this image needs to be CLEANed; this was just to get a taste. We'll start the imaging and deconvolution of the full FOV. Type in the PT shell:
pb_facets = image_clean_pb_facets( uv,
add_boxes = False,
imagr_params = imagr_params )

The function image_clean_pb_facets() applies the self-calibration solutions while imaging and CLEANing the FOV in an iterative way. The AIPS IMAGR gets called several times, each time CLEANing deeper into the image down to 2 times the latest noise measurement, until the noise improvement between two iterations is less than 5%. CLEANing only gets done in our static, predefined set of CLEAN boxes. However, if we would choose to do so, we can also let the function image_clean_pb_facets() build up its own set.

Now that the imaging has started, this will give us a few minutes to look at the self-calibration gain table (SN table) that is attached to the visibility data. In the AIPS shell type:
task 'SNPLT'
default
indisk 1
getname 1
inext 'SN'
invers 0
dotv 1
nplots 8
factor 0.3
optype 'PHAS'
pixrange -180,180
grchan 0
timerang 1,11,5,0,1,12,5,0
tvini
go

The AIPS TV window will show the gain phase solutions per antenna as a function of time.If you click the window, use 'A' to hold the current 8 antennas, 'B' to go to the next 8 antennas, and 'D' to stop. Type 'go' again to restart. The antennas with numbers 1-14 labelled Cxx are the central square antennas. Our reference antenna 25 is not one of them, but one just outside the central square. We don't see too much ionospheric phase structure on the central antennas, since they're so close to the reference antenna. The antennas with numbers 15-19, 20-24 and 25-30 are the E(ast), S(outh) and W(est) arm antennas, respectively, with increasing numbers meaning further away from the central square. One can see the ionospheric phase structure increase in amplitude when moving down an arm. One can also see a little phase jump occuring on antenna 16. This data could be flagged, but won't make much of an impact when left in.

When the imaging in the PT window is finished, the 90 facets are to be combined into a single image (PT):
sc1_image = combine_facets( uv, pb_facets )
sc1_image.rename( klass = 'SC1' )

The function combine_facets() is mainly a wrapper around the AIPS task FLATN, but also makes each facet image round before combining them. This makes sense because w-term errors (and other errors) increase with radial distance from the facet center. This time, we keep the individual facet images, as they will act as our vehicles to carry around DD calibration tables later. The noise of the CLEANed image has already been measured, as well as the total CLEANed flux (PT):
restore_parameter( sc1_image, 'residual_rms' ) * 1.e3
restore_parameter( sc1_image, 'total_clean_flux' )

So the total CLEAN flux is around 36.5 Jy, and the measured central noise in the CLEANed image is around 4.1 mJy/beam. Like above, you can view the resulting image on the AIPS TV window (AIPS):
pcat
getname <number in Cat column of PB.SC1.100>
pixrange -0.01,0.1
tvini
tvlod

Well, there is the sausage relic! Instead of using the AIPS TV, you may want to look at the resulting image in your favorite FITS viewer (I've provided 'ds9' in the bin directory). For this purpose we export the image in FITS format to the 'fits' sub-directory as 'CIZA.SC1.FITS' (PT):
write_fits_image( sc1_image, 'fits/' + target_name + '.SC1.FITS' )

You can also look at the image using the simple Python-based viewer we've used before, although this is somewhat slow for larger images (~3.5k x 3.5k in this case) (PT):
plot_image( sc1_image, z_range = [ -0.01,0.1 ] )

I encourage you to take a careful look around the whole image. Several of the brighter compact sources have lots of noticeable artificial patterns around them, which are residual sidelobes. These occur mainly due to DD calibration errors, as we will see below.

CIZA SC1.jpg

Fig. 2: GMRT 153 MHz image of the sausage after self-calibration.

Peeling of apparent bright sources

With the self-calibration gain table and local sky model from the CLEANed image, we have a good starting point on which we want to improve. We have now arrived at the measurementstep, where we attempt to 'measure' the DD ionospheric phases by self-calibrating on individual sources in the FOV. We'll follow an iterative approach, where we self-calibrate on the brightest sources in peak flux order. This technique is generally referred to as 'peeling'. While peeling a particular source, the contribution to the visibilities of all other sources will be removed (or at best reduced) by a-priori subtracting their best available model while temporarily applying the best available gain calibration. This creation of residual visibilities is done with (PT):
sub_pb_uv = subtract_model( uv, pb_facets )

For the first peeled source, our best model is the CLEAN component model from the image we just made, and our best calibration is the self-calibration table we started out with. For the second peeled source, the situation is almost the same, except that the first peeled source now has an updated CLEAN model and calibration used for subtraction. There is the logical increase in use of updated (peeled) source models and calibrations while subtracting when moving down the flux order.

We'll start up the peeling now, as it will take some time (15-30 minutes depending on CPU load) to finish. While it is running, you can have your coffee break, or read over some more explanations below. Copy the following in your PT shell:
reference_antenna = 25
calib_params = { 'uvrang' : [ 1., 1.e6 ] }
integration_time = restore_parameter( uv, 'integration_time' )

[ ppb_uv, ppb_facets ] = peel_pbo_facets( sub_pb_uv, pb_facets,
    calib_params = calib_params,
    reference_antenna = reference_antenna,
    phase_interval_min = integration_time / 60.,
    phase_interval_max = 2.,
    amplitude_snr = 300.,
    imagr_params = imagr_params,
    max_peel_count = 10 )

The peeling function peel_pbo_facets() is quite a complex beast. It performs several important steps:
  • Derive a list of candidate sources for peeling from the latest image, and order them
  • Combine nearby sources & match against a astrometry reference catalog (NVSS in our case)
  • Then per source:
    • Identify and add back the source CLEAN model to the residual visibilities
    • Phase self-calibrate this source several rounds, and possibly finish with an amplitude & phase self-calibration
    • Before each round of self-calibration, center source model on reference catalog source position
    • Monitor self-calibration progress, and reject when diverging or lack of image improvement
    • Subtract updated source CLEAN model to create updated residual visibilities
Candidate selection requires at least an image signal-to-noise ratio (SNR) of 15 within a single gain phase solution interval, based on naive sqrt(time) scaling. Sources within a few synthesized beams are grouped together. For use of the NVSS catalog, frequency and resolution differences can cause some errors in the assumed reference positions, especially for resolved sources. Unfortunately there are few other catalogs that are as complete as NVSS. We specify that the gain phase solution interval can vary between 16.1 sec (integration interval) and 2 minutes. Initially, we'll allow for a bit of slack on the flux criteria, since improving the calibration by peeling may lift a source above the threshold. Flux from fainter sources in the near-vicinity (within a facet radius) of a peeled source is also included in the self-calibration to boost the SNR a bit. Gain phase solutions are determined relative to the self-calibration solutions (like a perturbation). Per source, a final amplitude calibration is attempted when the local image SNR (or dynamic range) including all data is larger than 300, since that is roughly the SNR at which residual sidelobes stick clearly above the noise.

We choose to stop peeling when we have 10 peeled sources. This is to prevent 'over-fitting' our data, which can potentially lower the noise in the image artificially by suppressing unCLEANed, real emission near the detection threshold. The peeled sources are kept in a series of images, similar regular facet images, with an updated CLEAN model (AIPS CC table) and peeling gain calibration table (AIPS SN table) stored with each image as a CC and SN table, respectively. To look at the SN tables in AIPS, we have to copy them from the images to a visibility data set. For this we use the returned residual visibility data set (ppb_uv), which has all peeled + other sources subtracted. Therefore, we can dispose of the previous residual data set (PT):
copy_solution_tables( ppb_facets, ppb_uv )
sub_pb_uv.zap()

We'll have a look at the peeling gain calibration results below.

Selection of peeling results to use

We are now at the boundary between measuring and modeling. Before we fit proceed to ionospheric phase models to our peeling gain phases, we're going to manually check the peeling results first (PT):
determine_facet_overlap( ppb_facets )
for i in range( 1, 1 + restore_parameter( ppb_facets, 'facet_count' ) ):
    facet = get_facet( ppb_facets, i )
    print '%2d: %6.4f Jy, %4.2f min, %s, %s' % ( 
        i, get_model_flux( facet ),
        restore_parameter( facet, 'solution_interval' ),
        restore_parameter( facet, 'match_type' ),
        repr( get_facet_overlap( facet ) ) )
    plot_image( facet, z_range = [ -0.01,0.1 ] )
    plot_source( pb_facets, get_radec( facet ), z_range = [ -0.01,0.1 ] )

This bit of script will help assess the peeled source quality. Per peeled source, it will print some info, show the peeled source image, and show the NVSS source location in a (previous) self-calibration image facet (blue dot). Each printed line shows (i) the source ID, (2) the total CLEAN flux, (3) the final gain phase solution interval used, (4) some info on the source structure, (5) some info on the NVSS matching, and (6) a list of other peeled source facets with which this one overlaps.

When going through the images, look for sources that may have been peeled twice (2nd time as a nearby neighbor), or where the NVSS match is shifted or mismatched. Please also notice that for many of the sources the residual sidelobes have decreased noticeably. If life is fully deterministic, you will find that peeled source #9 has a wrong identification, i.e. the matching with NVSS picked the wrong source. So best to leave this one out for the ionospheric phase model fitting, since the mis-identification causes an artificial spatial phase gradient over the array for this source only (PT):
selected_facet_list = [ 1,2,3,4,5,6,7,8,10 ]

We will also have a look at the peeling gain phase tables, using the AIPS TV window. This is a bit awkward to do, but it gets the job done. First, we plot the time-series of the self-calibration gain phase solutions of only the outer antennas (in blue). If any problems exist with the calibration, they tend to show up on the outer antennas first (AIPS):
task 'SNPLT'
default
indisk 1
ucat
getn <number in Cat column of CIZA.PEELPB.100>
inext 'SN'
nplots 8
antennas 17,18,19,23,24,28,29,30,0
timerang 1,11,5,0,1,12,5,0
dotv 1
factor 0.3
optype 'PHAS'
pixrange -180,180
grchan 0
invers 1
go

We have fixed all the axes, so we can easily overplot (in yellow) the solutions for peeled source #1 (AIPS):
grchan 1
invers 1 + 1
go

Purple shows places where self-calibration and peeling solutions overlap. Repeating these lines with 'invers 1 + x' where x is the peeled source # will take you through all peeling gain phase solutions. If all goes well, you'll notice that there is some variation between the phase solutions of different peeled sources. Source #9 should show the largest deviation due to the mis-identification. Based on experience, I judged the quality of all tables (except for source #9) sufficient for ionospheric phase model fitting, which is what we'll do next.

Fitting the ionospheric phase model

The series of selected peeled sources represent a set of phase measurements over time of the GMRT antennas towards multiple directions on the sky. For any peeled source there is a strong correlation between the phase solutions for neighboring antennas, representative of the single physical process (radio wave propagation through the ionosphere) on the basis of this. We define a simple model with which we try to reproduce these observed ionospheric phases, allowing us to predict phase corrections in any direction later on. We'll start fitting the model to the observations, with an explanation of the inputs following below (PT):
fit_ionospheric_pmkl_model( uv, ppb_facets,
    facet_list = selected_facet_list,
    order = 15 )

For each integration time we fit an independent phase model. Our model is a combination of a 2nd order polynomial spatial phase screen to absorb large-scale fluctuations (de-trending), and a dual turbulent-layer ionosphere model, at 250 and 350 km height above the Earth's surface. Each turbulent layer is defined to contribute equally to the observed phase structure. The model name PMKL stands for Polynomial plus Multi-layer Karhunen-Loeve. The total number of free parameters used for describing the model is 20, which is 5 for the polynomial and 15 for the turbulence model.

The fitting routine will try to identify spurious antennas and/or sources, and repeat the fit without them. Successful polynomial fits are propagated as priors to the next fit. The fitting routine also down-weights data points that are very clustered (kind of uniform weighting). The fit results are stored with the visibility data (AIPS NI and OB tables), which can be accessed later to predict phase corrections in any direction. An additional new gain phase solution table (SN table) is stored that contains estimates of systematic (= instrumental) phase residuals per antenna. The latter is especially useful when instrumental phases are varying (phase jumps, phase drifts) over the duration of the observation.

The fitting routine prints per integration time the average residual phase RMS per antenna, for both the polynomial model and the PMKL model. From experience I find that reasonable fits have an RMS around 20 degrees RMS. Really good fits can have an RMS around 10 degrees. Fits become spurious above ~30 degrees.

When the fitting is done, we move the residual gain phase table to ppb_uv for convenience. Also, we plot the antenna average RMS phase residuals as a function of time.
call_aips_task( 'TACOP', indata = uv, outdata = ppb_uv, inext = 'SN', ncount = 1 )
uv.zap_table( 'SN', 0 )
plot_fit_errors( uv )

Fits that failed for some reason (e.g., too few data points) are plotted at phase = 0.

The fitted models can be made visible by depicting the integrated phase screen that a single antenna 'sees' as a function of time. Such an exercise is included as optional at the end of the tutorial.

Generating ionospheric phase corrections

Now the hard work is almost done, as we get to apply the phase model for imaging. With an ionospheric phase model in place, we can predict the phase corrections per antenna for arbitrary viewing directions in our FOV. Since we have quantized our FOV with a grid of facets, we will predict the phase corrections for each center of facet and apply it while imaging the whole facet. Typically, we will introduce a phase error at the edge of any facet that is smaller than the uncertainty of the phase model or the model fit itself. Also remember that the brightest sources have dedicated facets centered on their position, therefore minimizing this potential error.
generate_solution_tables_from_pmkl_fit_table( uv, pb_facets, 
    rms_limit = 25.,
    reference_antenna = reference_antenna )

The phase RMS limit of 25 degrees will reject ~2 time interval, while ~6 more were already rejected earlier. Each of the 90 facet images now has a gain solution table (SN table) attached to it, which will be applied to the visibility data while imaging and CLEANing the facet.

From the ionosphere model we can also derive ionospheric phase delays, which we have included in the generated gain tables by default. This becomes more important for larger fractional bandwidths than ours, but it doesn't hurt to include them even here.

Model gain tables can also be generated for the peeled sources, which enables a direct comparison between measurement (peeling) and model. Such an exercise is included as optional at the end of the tutorial.

Re-imaging the data using ionospheric phase calibration

OK, time to turn the crank again. The imaging routine below re-images and CLEANs the 90 facets covering the FOV, but now using a different gain solution table for each facet. Since this is not directly possible in AIPS, the SPAM code executes the major imaging cycle outside IMAGR. Running this will generate a HUGE amount of AIPS output, so let thy be warned. Running this should take 30-40 minutes (PT):
re_image_clean_pb_facets( uv, pb_facets,
    imagr_params = imagr_params )

sp1_image = combine_facets( uv, pb_facets )
sp1_image.rename( klass = 'SP1' )
write_fits_image( sp1_image, 'fits/' + target_name + '.SP1.FITS' )

restore_parameter( sp1_image, 'residual_rms' ) * 1.e3
restore_parameter( sp1_image, 'total_clean_flux' )

Check the image in AIPS (PB.SP1.100) or in your favorite FITS viewer (fits/CIZA.SP1.FITS). Compared to the self-calibration image, the central noise level has gone down from 4.1 to 3.4 mJy/beam (~17% reduction). When comparing the two images, one can clearly see a reduction of the residual sidelobes around bright sources. There is also a general increase of a few percent in the peak fluxes of sources. These two observations suggest that sources are more 'in focus', and scatter less of their power into sidelobes. The total CLEANed flux is 36.5 Jy, which is exactly the same as before.

Re-imaging while including peeling solutions

The sidelobes around bright sources have reduced, but are still significant. This is probably caused by a combination of issues. Remembering the less affected images coming out of the peeling, it is likely that the remaining sidelobes are caused by:
  • deviations between the peeling phase solutions and model phase solutions, since the model fits were not perfect.
  • amplitude errors, since the brightest peeled sources were also amplitude-calibrated
Since the mechanism for using different gain calibration tables for different facets is already in place, we can choose for the brightest sources to directly use the peeling gain tables. All it requires is to have dedicated facets for the brightest sources (we already have that) and to copy the gain tables from the peeled source facets to the regular facets. As a test, we will do this for the ~6 sources with a peak-to-noise ratio of at least 300 (PT):
peel_facet_list = []
cpb_noise = restore_parameter( uv, 'cpb_noise' )
for i in selected_facet_list:
    facet = get_facet( ppb_facets, i )
    peak_flux = get_image_maximum( facet )[ 0 ]
    if ( peak_flux / cpb_noise > 300. ):
        peel_facet_list.append( i )

print peel_facet_list

To be consistent, we'll include phase delays derived from the ionosphere model in the peeling gain tables (PT):
generate_solution_tables_from_pmkl_fit_table( uv, ppb_facets,
    rms_limit = 25.,
    reference_antenna = reference_antenna,
    include_phases = False,
    facet_list = peel_facet_list )

for i in peel_facet_list:
    facet = get_facet( ppb_facets, i )
    combine_solutions( facet )

The actual replacing of the model solutions with peeling solutions happens here (PT):
replace_model_solutions_with_peel_solutions( uv, pb_facets, ppb_facets, facet_list = peel_facet_list )

And finally we can re-image again for 30-40 minutes. To make it fair, we'll CLEAN as deep as the previous image by manually setting the last measured noise back to 4.1 mJy/beam (PT):
store_parameter( uv, 'cpb_noise', 0.0041 )
re_image_clean_pb_facets( uv, pb_facets,
    imagr_params = imagr_params )

sp1a_image = combine_facets( uv, pb_facets )
sp1a_image.rename( klass = 'SP1A' )
write_fits_image( sp1a_image, 'fits/' + target_name + '.SP1A.FITS' )

restore_parameter( sp1a_image, 'residual_rms' ) * 1.e3
restore_parameter( sp1a_image, 'total_clean_flux' )

Please use AIPS or your favorite FITS viewer to compare this image with the previous ones. The latest noise level is 3.2 mJy/beam, so the noise went down another 6%, which is now 22% compared to the SC image. A total flux of 35.7 Jy got CLEANed, which is a bit (~2%) less than the other maps. This might be a result of introducing gain amplitude calibrations into the imaging. However, doing a rough comparison of peak fluxes shows that there's no obvious systematic decrease in the flux scale in the latest image. This might also be a result of introducing too many calibration parameters, artificially reducing the noise. However, the low-level diffuse emission in the central region does not seem to have changed too much, which would be the first place one would notice. A third option might be that fewer sidelobes got CLEANed in the latest image, which would (of course) be a good thing. It is certain that at least several CLEAN boxes have been (automatically) placed on sidelobes of bright sources. Processing the images with a source extraction tool (e.g., PyBDSM - see talk/tutorial by David) and comparing the results can give more clues where the difference comes from.

Rounding up

The main part of this tutorial is now done. We've witnessed that, for this particular data set, applying DD ionospheric calibration instead of self-calibration gives you:
  • A background RMS noise reduction of order 20%
  • A reduction of residual sidelobe structure near bright sources
  • A moderate increase in source peak fluxes (few %)
There still remains quite some structure in the GMRT image background, but is now less concentrated around bright sources. We may improve on our current result by taking the data through another round (or multiple rounds) of peeling, modeling and imaging. We can safely assume that the current calibration and sky model gives us a better starting point than what self-calibration provided at the start of peeling. Typically, when processing low-frequency GMRT or VLA data, rounds of calibration and imaging need to be complemented with careful editing of the residual visibility data. In our case here, the limited UV-coverage of a single hour observation, and the resulting dirty beam / point spread function with high sidelobes, also provides an extra challenge.

The next two sections are left as optional.

Comparing model phases with peeling phases

To see how well the fitted phase model represents the peeling 'measurements', we will generate model phase tables for the peeled sources (PT):
generate_solution_tables_from_pmkl_fit_table( uv, ppb_facets, rms_limit = 25., reference_antenna = reference_antenna )

To view them in AIPS, we need to copy the gain tables from the peeled source images to a visibility data set. Again, we'll use the residual visibility data after peeling:
copy_solution_tables( ppb_facets, ppb_uv )

First, we plot the time-series of the peeling gain phase solutions for the first source of only the outer antennas (AIPS):
tget snplt
tvini
grchan 0
invers 1 + 1
go

We now overplot (in yellow) the model gain phase solutions for the same source (AIPS):
grchan 1
invers 12 + 1
go

One can repeat this for the other peeled sources by repeating the steps above, replacing the 'invers 1/12 + 1' with 'invers 1/12 + x', where x is the peeled source #. Don't forget that we disqualified source #9. In general, there is a good match between measurements and model. For individual sources there is sometimes a small offset, increasing along an arm. This is a logical signature for a source for which there is a small offset between the true position and the assumed position from NVSS.

Making a phase screen movie

This part requires that the tools 'mplayer and 'mencoder' are available on your system (if not, you can still make the individual movie frames). If we consider the FOV of one antenna, we can make a series of phase screen images derived from our fitted model, one for each integration time. Times that were rejected during the model fitting will get a zero phase screen. In the PT window type:
image_prefix = './prtfil/CIZA_'
make_pmkl_phase_screen_images( uv, reference_antenna,
    facets = ppb_facets,
    facet_list = selected_facet_list,
    print_info = True,
    field_factor=1.5,
    grid_size = 128,
    image_prefix = image_prefix )

The PNG images, one per time frame, are written to the 'prtfil' sub-directory. Each image has the name 'CIZA_xxxxx.png', where xxxxx is a 5-digit sequence number starting at zero. The execution may take a while, as the model is now applied to larger grid than previously. The images use a cyclic color scheme to capture the wrapping nature of phases. In the image, the sky positions of the peeled sources used for model fitting are depicted by '+'-signs.

To combine the images into a movie (10 fps), type in PT:
image_count = len( get_time_list( uv ) )
image_file_name_list = [ image_prefix + '%05d.png' % ( i ) for i in range( image_count ) ]
movie_file_name = image_prefix[ 0 : -1 ]
make_movie_from_images( image_file_name_list, movie_file_name )

The function make_movie_from_images() calls mencoder with the appropriate settings. The resulting movie is 'prtfil/CIZA.avi', which can be played with mplayer.

This topic: Main > TWikiUsers > HuibIntema > HuibIntemaWebHome > HuibIntemaSpamTutorial
Topic revision: 2014-04-28, HuibIntema
This site is powered by FoswikiCopyright © by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding NRAO Public Wiki? Send feedback