Jupyter Notebook Interface for CASA

This notebook demonstrates a Jupyter notebook interface to CASA, the Common Astronomy Software Applications package developed by NRAO, used for reduction of radio astronomy data such as that of the VLA and ALMA. A standalone version of CASA can be obtained from here.

This notebook runs in a sandbox on a webserver at ASTRON. There is some demonstration data available in this sandbox. The purpose of this notebook is to introduce the notebook interface. To actually work with the notebook interface, you can install it on your own computer through the instructions at https://github.com/aardk/jupyter-casa.

The notebook interface to CASA was developed by Aard Keimpema at the Joint Institute for VLBI ERIC (JIVE). Installation in the sandbox interface (based on tmpnb) was done by Tammo Jan Dijkema at ASTRON. This work is a part of the ASTERICS project, Astronomy ESFRI & Research Infrastructure Cluster ASTERICS - 653477.

This notebook interface was built on top of CASA 5.4.0. Adaptions were made to run it on a newer version of Python than the one distributed with CASA. Future development of this notebook, or upgrading the underlying CASA version, is currently being planned.

The current notebook is taken from the VLA Continuum Tutorial. The dataset inside this sandbox has been averaged down considerably in size, the original was 3GB data, so the results will be worse than in the tutorial.

NB: This container is for demonstration purposes only, it will be deleted automatically after 5 minutes of inactivity

Tutorial, inspect the data

To get a sense of the array, as well as identify an antenna for later use in calibration, we use the task plotants.

We now inspect some basic information about the data. The task listobs can be used to get a listing of the individual scans comprising the observation, the frequency setup, source list, and antenna locations.

The log of this task can be shown by clicking the blue button.

Examining and editing the data

It is always a good idea to examine the data before jumping straight into calibration. Moreover, from the observer's log, we already know that one antenna will need to be flagged because it does not have a C-band receiver. Start by flagging data known to be bad, then examine the data.

From the listobs output above, one may have noticed that the first scan is less than 1 minute long. This first scan can safely be flagged. This can be done using the task flagdata.

The logs (visible by clicking the blue button), indicate that, among other things, the flags that existed in the data set prior to this run will be saved to another file called flagdata_1. Should you ever desire to revert to the data prior to this run, the task flagmanager could be used.

From the observer's log, we know that antenna ea13 does not have a C-band receiver and antenna ea15 had some corrupted data, so they should be flagged as well.

It is common for the array to require a small amount of time to settle down at the start of a scan. Consequently, it has become standard practice to flag the initial samples from the start of each scan. This is known as 'quack' flagging.

Having now done some basic editing of the data, based in part on a priori information, it is time to look at the data to determine if there are any other obvious problems. One task to examine the data themselves is plotms:

The task plotms allows one to select and view the data in many ways. Another useful plot we will make is a datastream plot of the antenna2 in a baseline for the data versus ea01. This shows, assuming that ea01 is in the entire observation, when various antennas drop out.

From this display, you see immediately that the flagging we did earlier of antennas 10 and 12 (ea13 and ea15) has taken effect. For the remaining antennas, you see that 1, 6, and 13 (ea02, ea08, and ea16) are missing some blocks toward the beginning and also toward the end of the run. Antenna 3 (ea04) is missing the last scan (on the polarization calibrator, 3C84) and antenna 23 (ea26) is missing scans near the end. None of these antennas should be chosen as the reference calibrator during the calibration process, below.

Calibration

It is now time to begin calibrating the data. The general data reduction strategy is to derive a series of scaling factors or corrections from the calibrators, which are then collectively applied to the science data. For more discussion of the philosophy, strategy, and implementation of calibration of synthesis data within CASA, see Synthesis Calibration in the CASA Reference Manual.

Recall that the observed visibility $V^{\prime}$ between two antennas $(i,j)$ is related to the true visibility $V$ by:

$V^{\prime}_{i,j}(u,v,f) = b_{ij}(t)\,[B_i(f,t) B^{*}_j(f,t)]\,g_i(t) g_j(t)\,V_{i,j}(u,v,f)\,e^{i [\theta_i(t) - \theta_j(t)]}$

Here, for generality, we show the visibility as a function of frequency $f$ and spatial wave numbers $u$ and $v$. The other terms are:

• $g_i$ and $\theta_i$ are the amplitude and phase portions of what is commonly termed the complex gain. They are shown separately here because they are usually determined separately. For completeness, these are shown as a function of time t to indicate that they can change with temperature, atmospheric conditions, etc.
• $B_i$ is the complex bandpass, the instrumental response as a function of frequency $f$. As shown here, the bandpass may also vary as a function of time.
• $b(t)$ is the often-neglected baseline term. It can be important to include for the highest dynamic range images or shortly after a configuration change at the VLA, when antenna positions may not be known well.

Strictly, the equation above is a simplification of a more general measurement equation formalism, but it is a useful simplification in many cases.

A priori Antenna Position Corrections

NRAO monitors the positions of the VLA antennas on a regular basis. The corrections are then placed into an NRAO database. If updated positions were entered into the database after your observation date, the corrections to the newly measured positions can still be applied during your data reduction process in this step. Any updated positions that were entered into the database before your observations will already be accounted for in your data.

The calculations are inserted via gencal which allows automated lookup of the corrections. To see how to calculate corrections manually, go to the VLA Baseline Corrections site.

In the logger you can see the corrections reported.

Initial Flux Density Scaling

Let's find the available calibrator models with setjy and setting the parameter listmodels=True:

Since any image could be a potential calibrator model, setjy will list all *.im and *.mod images in the working directory. In addition, it will list all models that are provided by NRAO with the CASA package, and they will be picked by their names. We will be using the C-band VLA standard model for 3C286 which is aptly named '3C286_C.im':

Initial phase calibration

Before solving for the bandpass, we will do an initial phase calibration. The reason for this step is to average over the (typically small) variations of phase with time in the bandpass, before solving for the bandpass solution itself. Depending upon frequency and configuration, there could be significant gain variations between different scans of the bandpass calibrator, particularly if the scans happen at much different elevations. One can solve for an initial set of antenna-based gains, which will later be discarded, in order to moderate the effects of variations from integration to integration and from scan to scan on the bandpass calibrator. While amplitude variations with time will have little effect on the bandpass solutions, it is important to solve for phase variations with time to prevent de-correlation when vector averaging the data for computing the final bandpass solution.

We use the CASA task gaincal to solve for phase versus time for the central channels on our three calibrators:

To really see what is going on, we use plotcal to inspect the solutions from gaincal for a single antenna. Note below, that there are phase jumps for antenna ea05 where the phase appears to be oscillating between two states.

We will not be able to transfer calibration for antenna ea05 so we flag it from the data:

For the following bandpass solution we need only solve for our bandpass calibrator, and we will do so now after flagging. The following call to gaincal is similar to the one above, but selects only the bandpass calibrator (using the field parameter). This is the calibration table we will use when solving for the bandpass solution, below.

Delay calibration

The first stage of bandpass calibration involves solving for the antenna-based delays which put a phase ramp versus frequency channel in each spectral window (Figure 3C). The K gain type in gaincal solves for the relative delays of each antenna relative to the reference antenna (parameter refant), so be sure you pick one that is there for this entire scan and good. This is not a full global delay, but gives one value per spw per polarization.

We can plot these solutions (in nanoseconds) as a function of antenna:

These are within about 4 nanoseconds, as expected for the early science observations with the newly upgraded VLA.

Bandpass calibration

This step solves for the complex bandpass, $B_i$. All data with the VLA are taken in spectral line mode, even if the science that one is conducting is continuum, and therefore requires a bandpass solution to account for gain variations with frequency. Solving for the bandpass won't hurt for continuum data, and, for moderate or high dynamic range image, it is essential.

We now form the bandpass, using the phase solutions just derived:

Once again, one can use plotcal to display the bandpass solutions. Note that in the plotcal inputs below, the amplitudes are being displayed as a function of frequency channel.

Gain calibration

The next step is to derive corrections for the complex antenna gains, $g_i$ and $\theta_i$. As discussed above, the absolute magnitude of the gain amplitudes $g_i$ are determined by reference to a standard flux density calibrator.

In principle, one could determine the complex antenna gains for all sources with a single invocation of gaincal; for clarity here, two separate invocations will be used.

In the first step, we derive the appropriate complex gains $g_i$ and $\theta_i$ for the flux density calibrator 3C 286.

In the second step, the appropriate complex gains for a direction on the sky close to the target source will be determined from the phase calibrator J1822-0938. We also determine the complex gains for the polarization calibrator source J0319+4130. These will be solved separately, but in practice could be solved together as there are no gaintables that are time dependent at this point (and thus would risk having cross-source interpolation issues), nor are we doing different solution intervals per source.

If one checks the gain phase solutions using plotcal, one should see smooth solutions for each antenna as a function of time.

This is also a good time to check that our chosen reference antenna (ea21) has good phase stability (i.e., the phase difference between the right and left polarizations is stable with time). This is a prerequisite for accurate polarization calibration. To do this, we plot the complex polarization ratio by selecting poln='/':

As can be seen above, there is a bit of drift (a few degrees here and there), but no phase jumps. This means that ea21 is, indeed, a good choice for reference antenna.

Scaling the Amplitude Gains

While we know the flux density of our primary calibrator (J1331+3030$\equiv$3C 286), the model assumed for the secondary calibrator (J1822-0938) was a point source of 1 Jy located at the phase center. While the secondary calibrator was chosen to be a point source (at least, over some limited range of uv-distance; see the VLA calibrator manual for any u-v restrictions on your calibrator of choice at the observing frequency), its absolute flux density is unknown. Being point-like, secondary calibrators typically vary on timescales of months to years, in some cases by up to 50-100%.

We use the primary (flux) calibrator to determine the system response to a source of known flux density and assume that the mean gain amplitudes for the primary calibrator are the same as those for the secondary calibrator. This allows us to find the true flux density of the secondary calibrator. To do this, we use the task fluxscale, which produces a new calibration table containing properly-scaled amplitude gains for the secondary calibrator.

The task fluxscale will print to the CASA logger the derived flux densities of all calibrator sources specified with the transfer parameter. These are also captured in the return variable from the task. You should examine the output to ensure that it looks sensible. If the data set has more than one spectral window, depending upon where they are spaced and the spectrum of the source, it is possible to find quite different flux densities and spectral indexes for the secondary calibrators.

We plot the rescaled amplitudes from this table:

You can see in these figures that the amplitude gain factors are now similar across sources, compared to the raw factors in the G1 table.

Applying the calibration

Now that we have derived all the calibration solutions, we need to apply them to the actual data, using the task applycal. The measurement set DATA column contains the original data. To apply the calibration we have derived, we specify the appropriate calibration tables, which are then applied to the DATA column, with the results being written in the CORRECTED_DATA column. If the dataset does not already have a CORRECTED_DATA scratch column, then one will be created in the first applycal run.

First, we apply the calibration to each individual calibrator, using the gain solutions derived on that calibrator alone to compute the CORRECTED_DATA. To do this, we iterate over the different calibrators, in each case specifying the source to be calibrated (using the field parameter).

Finally, we apply the calibration to the target fields in the mosaic, linearly interpolating the gain solutions from the secondary calibrator, J1822-0938. In this case however, we want to apply the amplitude and phase gains derived from the secondary calibrator, J1822-0938, since that is close to the target source on the sky and we assume that the gains applicable to the target source are very similar to those derived in the direction of the secondary calibrator. Of course, this is not strictly true, since the gains on J1822-0938 were derived at a different time and in a different position on the sky from the target. However, assuming that the calibrator was sufficiently close to the target, and the weather was sufficiently well-behaved, then this is a reasonable approximation and should get us a sufficiently good calibration that we can later use self-calibration to correct for the small inaccuracies thus introduced.

We should now have fully-calibrated visibilities in the CORRECTED_DATA column of the measurement set, and it is worthwhile pausing to inspect them to ensure that the calibration did what we expected it to. We make some standard plots:

The phase solutions do not look very good, which may be due to over-averaging the data or not doing polarization calibration. In the original tutorial, the phases do look good.

Another nice display is to use plotms to plot the amplitude and phase of the CORRECTED_DATA column against one another for one of the parallel-hand correlations (RR or LL); the signal in the cross-hands, (RL and LR) is much smaller and will be noise-like for an unpolarized calibrator. This should then show a nice oval of visibilities, with some scatter, centered at zero phase and at the amplitude found for that source in fluxscale.

Inspecting the data at this stage may well show up previously-unnoticed bad data. Plotting the corrected amplitude against UV distance or against time is a good way to find such issues. If you find bad data, you can remove them via interactive flagging in plotms (not in the notebook interface yet) or via manual flagging in flagdata once you have identified the offending antennas/baselines/channels/times. When you are happy that all data (particularly on your target source) look good, you may proceed.

Now that the calibration has been applied to the target data we can split off the science targets, creating a new, calibrated measurement set containing all the target fields.

Initial imaging

Now that we have split off the target data into a separate measurement set with all the calibration applied, it's time to make an image. Recall that the visibility data and the sky brightness distribution (a.k.a. image) are Fourier transform pairs.

$I(l,m) = \int V(u,v) e^{[2\pi i(ul + vm)]} dudv$

The $u$ and $v$ coordinates are the baselines measured in units of the observing wavelength, while the $l$ and $m$ coordinates are the direction cosines on the sky. For generality, the sky coordinates are written in terms of direction cosines, but for most VLA (and ALMA) observations they can be related simply to the right ascension ($l$) and declination ($m$). Also recall that this equation is valid only if the $w$ coordinate of the baselines can be neglected. This assumption is almost always true at high frequencies and smaller VLA configurations (such as the 4.6 GHz D-configuration observations here). The $w$ coordinate cannot be neglected at lower frequencies and larger configurations (e.g., 0.33 GHz, A-configuration observations). This expression also neglects other factors, such as the shape of the primary beam. For more information on imaging, see Synthesis Imaging within the CASA Reference Manual.

CASA has a single task clean which both Fourier transforms the data and deconvolves the resulting image. For the purposes of this tutorial, we will make a simple mosaic clean image in Stokes I only.

Setting the appropriate pixel size for imaging depends upon basic optics aspects of interferometry. Using plotms to look at the newly-calibrated, target-only data set:

The maximum baseline is about 16,000 wavelengths, i.e., an angular scale of 12 arcseconds ($\lambda/D=1/16000$). The most effective cleaning occurs with 3-5 pixels across the synthesized beam. For example, a cell size of 2.5 arcseconds will give just under 5 pixels per beam.

The supernova remnant itself is known to have a diameter of order 9 arcminutes, corresponding to about 216 pixels for the chosen cell size. The mosaic was set up with 7 fields, 1 centered on the remnant with 6 flanking fields; the spacing of the fields was chosen based on the size of the antenna primary beam. For the choice of ftmachine='mosaic' (our main mosaicking algorithm), you do not have to fit the mosaic inside the inner quarter of the total image in order to prevent image artifacts arising from aliasing, we just want to have a bit of padding around the outside. Although CASA has the feature that its Fourier transform engine (FFTW) does not require a strict power of 2 for the number of linear pixels in a given image axis, it is somewhat more efficient if the number of pixels on a side is a composite number divisible by any pair of 2 and 3 and/or 5. Because clean internally applies a padding of 1.2 (=3x2/5) choose 480, which is 25 × 3 × 5 (so 480 × 1.2 = 576 = 26 × 32). We therefore set imsize=[480,480] and our mosaic fits comfortably inside the image.

Running clean will take about two minutes.

Interactive cleaning is not (yet) supported in this notebook interface. We can open the viewer to see the result of the imaging.

Note that the image is cut off in a circular fashion at the edges, corresponding to the default minimum primary beam response (minpb) within clean of 0.2.