Ok, now that all of the calibration is out of the way we can start with imaging the data. Finally!
This tutorial is designed to be worked through using this scriptFromImaging.py python script. The syntax is the same as the "scriptForCalibration" which was used yesterday but for some steps you will have to calculated and fill in the blank values. These will be noted in red in this document. To proceed you will need the fully calibrated measurement we created yesterday:
'uid___A002_X7fb89e_X6e1.ms.split.cal'. If you didn't get finished or want to start from a clean copy it can be downloaded as a tar fail here. (but be aware it is 45GB).
STEP 0: CONTINUUM SUBTRACTION
The first thing to consider is the level of continuum emission in the source. To perform accurate science with detected spectral lines we need to work with continuum subtracted data; conversely to do good science with the source continuum emission we need to use only line free channels!. To guage the level of continuum and find the line free channels we need to make data cubes of the each of the science spectral windows1. This is time consuming for this dataset so we've done that bit for you. (You can download the full cubes here: SPW0, SPW1, SPW2 and SPW3).
In practice once you had full cubes of the data you can inspect the spectra of your target and note down the line free channels. Figure 1 gives you an idea of the regions of the target spectrum which we need to use in continuum subtraction.
Fig.1 - Simple spectrum plotted as Channel vs Flux Density. The line free channel ranges are noted by the black lines and channel values.
For this simple case the line free channel are 0 to 90, 150 to 280, 317 to 370 and 430 to 510. In CASA syntax this would be written as SPW No.:0~90;150~280;317~370;430~510 (where 'SPW No.' would be replaced by the SPW number and remember these are not the values to use with our data).
Unfortunately things are not this simple for IRAS16293 as sources A and B have very different spectra (see Figures 2 and 3 for sources B and A in SPW3) and source B is incredibly line rich. So for the purposes of this tutorial we will use the below values in the script as line free. These values were created using CASA task Adam is developing (LUMBERJACK on Github) to find line free channels algorithmically.
Fig. 2 & 3 - Spectrum at IRAS16293 sources B (left) and A (right) in SPW3.
The values from the Lumberjack are here: Line free channels
Using these values we will use the CASA task uvcontsub to fit a continuum level and subtract it directly from the visibility data. To do this enter the line free chan values in to the 'fitSPW' parameter of uvcontsub at mystep = 0 in the script, which looks like this:
uvcontsub(vis = thevis,
field = 'IRAS16293-2422',
fitSPW = [values from within text],
want_cont=False)
You can then set mysteps=[0] and then execute the script.This will create the measurement set 'uid___A002_X7fb89e_X6e1.ms.split.cal.contsub'. We will use this in the imaging of spectral lines in step 2.
STEP 1: IMAGE CONTINUUM
Using the channels we identified as line free we can create an image of the continuum emission in IRAS 16293-2422 with all the avaliable line free bandwidth. To do this we will use the task tclean which is CASA's clean task2 and the data in the non continuum subtracted measurement set (so thevis = 'uid___A002_X7fb89e_X6e1.ms.split.cal'). From this mornings presentation we will also need to calculate the following parameters:
imsize, the image size (in pixels). This should be at least the primary beam size of an ALMA antenna (Diameter = 12m). Preferably twice the PB.
cell, the size of the pixels in arcsec. This should be a fraction (1/3 to ~1/5) the estimated synthesised beam size.
threshold, the limit at which to stop cleaning (in uJy, mJy or Jy).
To calculate the first two values you will need to know (a) the central frequency of the observations, (b) the diameter of an ALMA dish and (c) the
the maximum baseline in the array during the observations. These can be found using listobs (for central freq and the number of dishes) and using the command x=aU.getBaselineLengths('uid___A002_X7fb89e_X6e1.ms.split.cal') for the maximum baseline.
For the threshold you will need to know (1) the diameter of an ALMA dish (again), (2) the number of dishes used, (3) the time on source, (4) the bandwidth of line free channels and (5) the system temperature (Tsys). (1) to (3) can be read directly from the listobs output, (4) can be calculated from the number of line free channels multiplied by the channel width (also in listobs) and (5) can be estimated from the Tsys plots made during calibration yesterday.
Knowing these values and using a couple of the equations in the Fundamentals of Interferometry talk yesterday you can calculate the values needed for clean.
Alternately, if you're pushed for time you can highlight the next line.
Spoilers:Central Freq= 231.25 GHz, Maximum baseline= 558.2 m, Number of dishes = 35, Tsys ~ 80K, Line free bandwidth = 3538 channels at a channel width of 122.070kHz ~ 0.431GHz Once you have calculated these values you can enter them into the relevant lines of code in mystep = 1, which looks like this:
tclean(vis = thevis,
imagename = 'IRAS_16293-2422_sci.spw0_1_2_3.mfs.I.manual',#- Or whatever you would like to call the output image
field = 'IRAS_16293-2422',
stokes = 'I',
spw =[values from within text at step 0],
outframe = 'LSRK',
specmode = 'mfs',
nterms = 1,
imsize = [X, Y],
cell = 'Z.Aarcsec',
deconvolver = 'hogbom',
niter = 10000,
weighting = 'briggs',
robust = 0.5,
gridder = 'standard',
pbcor = True,
threshold = 'Q.VmJy',
interactive = True)
You can then set mysteps=[1] and execute the script. This will start an interactive clean during which you can draw clean boxes around the sources of emission and clean down until you are happy with the residuals. As a rule of thumb we want the emission within our clean boxes to look like the emission outside out clean boxes. Figures 4, 5 and 6 show examples of this.
Fig. 4, 5 & 6 - Continuum image residual after 0 clean iterations (Fig. 4 left), two major cycle (Fig. 5 centre) and when clean box residuals are smaller than outside (Fig. 6 left).
When you have finished cleaning your continuum image open it in the CASA viewer with the command viewer and explore the image.
ASIDE: A bit of explanation for the parameters used:
specmode = 'mfs', MFS stands for multi-frequency synthesis. This tells tclean that you are making a continuum image and to produce an image with a single channel which uses all the data in the specified bandwidth (SPW parameter) to denfine the sky brightness distribution.
deconvolver = 'hogbom', This is the algorithm used in deconvolution during clean minor cycles. 'Hogbom' is the 'Classic' clean algorithm, it tends to be the fastest but is best suited to point like/ minimally extended sources. For IRAS16293 that is fine, but for more extended sources there are other algorithms e.g. Multiscale which will do a better job. See the CASA help documentation for more details on this.
weighting = 'briggs', This parameter dictates how the visibilties are weighted during Clean. Typically options are:
'natural', where visibilties are weighted based on noise on each visibility and summed when gridded in uv-space. This produces the best SNR ratio but with a larger synthesised beam that other options.
'uniform', as natural but the weights are recalculated after gridding to give equal weights across uv-space. This increase resolution at the cost of SNR.
'briggs', a sliding scale weighting somewhere between 'natural' and 'uniform'.
robust = 0.5, This parameter is used in conjunction with weighting='briggs' to set the sliding scale between 'natural'-like (robust = 2.0) and 'uniform'-like weighting (robust = -2.0). A value of 0.5 is somwhere between the two, but closer to natural weighting.
STEP 2: IMAGE A LINE
Here we switch to using the continuum subtracted measurement set made in step 0 (thevis = 'uid___A002_X7fb89e_X6e1.ms.split.cal.contsub'). For this example we will use a line of SO2 found in spectral window 3 (see Figure 7) 3.
Fig.7 - Zoom in of the SO2 line in SPW3. Click on this to see a larger version.
The line is detected to peak at nu=221.963GHz in source A and the rest frequecy of the line is nu0=221.9652GHz. This equatest to a Vlsr for IRAS16293 of ~2.9km/s (which is compatible with published values Vlsr source A = 1.5 to 3.2 km/s [Pineda et al 2012, and reference therein]).
For this step we need to change specmode from 'mfs' to 'cube' so that tclean will create a spectral cube. This change in specmode opens some additional parameters we need to set. Those are:
nchan, the number of channels you wish to have in your cube.
start, the start channel,frequency or velocity you wish to start you cube at.
width, the width of the individual channels within your cube. This must be in compatible units to the value you entered start start. Here you can average up some channels. I would recommend using twice the channel width for this tutorial.
Based on the spectra in Figure 7 select your nchan, width and start values and enter them into the code at step = 2.
The final thing we need to do before running this step is to recalculate the threshold value. Now that we are channelising the data we should replace the bandwidth of line free channels used in Step 1 for the bandwidth of a single channel in your image cube (whatever you set width to).
When you have decided upon these values you can enter them into the relevant lines of code in mystep = 2, which looks like this:
tclean(vis = thevis,
imagename = 'CONTSUB_IRAS_16293-2422_sci.spw3_line_SO2.cube.I.manual',#- Or whatever you would like to call the output image
field = 'IRAS_16293-2422',
stokes = 'I',
spw = '3',
outframe = 'LSRK',
restfreq = '221.963GHz',
specmode = 'cube',
imsize = [X, Y],#- You can use the previously calculated values here or re-calculate for this specific frequency.
cell = 'Z.Aarcsec',#- You can use the previously calculated values here or re-calculate for this specific frequency.
deconvolver = 'hogbom',
niter = 10000,
weighting = 'briggs',
robust = 0.5,
gridder = 'standard',
pbcor = True,
threshold = 'Q.VmJy',
width = '000.00kHz',
start = '000.00GHz',
nchan = W,
interactive = True)
You can then set mysteps=[2] and then execute the script. Again this will start an interactive clean, where you can draw a clean box around emission in each channel of the cube.
ASIDE: Auto-masking
You'll quickly see that this can be very time consuming! But thankfully CASA recently introduced an Auto-masking option. Automasking (called by setting usemask='auto-multithresh') uses various user defined thresholds to 'mimic' what an experienced user would do when drawing clean boxes. To give this method a go by setting mysteps = [3] which has the recommended thresholds set to clean the same line as above, you just need to input your nchan, width, start and threshold values and execute this.
The definitions of the threholds and the effects of changing them are described in detail in the Automasking CASA guide which can be found here: https://casaguides.nrao.edu/index.php/Automasking_Guide
When you have finished cleaning your cube open the CASA viewer with the command viewer and explore the image.
STEP 3: IMAGE MOMENTS
Using the cube generated in Step 2 we can start to preform some image analysis. A useful task for this is immoments. This task allows you generate images of the -1th moment (mean value in spectrum), 0th moment (integrated intensity) up to the 11th moment (coordinate of the miminum pixel in the spectrum, apparently). (Note some of these moments have the typical mathematical meaning, where as some are more like convenience functions).
The most commonly used moments for image analysis are 0, 1 and 2 which are 'integrated value of the spectrum', 'intensity weighted coordinate' and 'intensity weighted dispersion of the coordinate'. The latter two are commonly used to infer 'velocity fields' and 'velocity dispersion' in the source.
To generate such images for IRAS16293 we call immoments as follows:
immoments(imagename=yourImage,#- The name of the file you created in Step 2
moments=[0,1,2],
outfile='CONTSUB_IRAS_16293-2422_sci.spw3_line_SO2.cube.I.manual.image.mom',#- Or whatever you would like to call the output image
excludepix=[X,Y],
chans='z~w')
For this task you need to set the channels over which you want to collapse the image to generate the moments using the chans parameter. It is also useful to set values in either includepix or excludepix, either of these parameters can be used (but not both) to limit the pixel values used in the generation of the moment images. includepix will allow immoments to use only pixel from a low value threshold (e.g. a few times the RMS noise) to a high value and conversely excludepix will only allow use of pixels from a very low value to a high value cut off (again e.g. a few times the RMS noise).
STEP 4: EXPORT TO FITS FILE
Possibly the least exciting step, but as we have seen CASA images are in the CASA image format, which is essentially a directory structure. This is fine if we want to work solely in CASA for all our analysis and will be sharing our results with other who will exclusively use CASA. However, it is more that likely that your colleage/supervisor/tame theorist will ask for an image as a FITS file, (the standard image format for most astronomical images) to use in another package. CASA has the task exportfits which allows us to do this.
In its simplest form you can export a FITS file with the following code:
There are however a few things to keep in mind which might be useful when trying to get your CASA exported FITS file to be loaded properly in external packages. This is due to the fact that by default CASA images have 4 axis, which are listed in the following order:
Axis 1 = x spatial axis,
Axis 2 = y spatial axis,
Axis 3 = Stokes axis (I Q U V),
Axis 4 = frequency axis,
All four typically exist even if they have a value of 1 (i.e. Stokes 'I' only and 1 frequency plane).
Many other packages expect either the Stokes and Frequency axes to be reversed or for the Stokes axis not to exist at all. To over come this issue exportfits has a few of parameters to deal with this. dropstokes = True will drop the Stokes axis entirely, dropdeg = True will drop any axes which have a length of 1 and stokeslast = True will put ensure that the Stokes axis is the last (rather than 3rd) axis in the exported FITS.
mystep=[5] in the code has an example of this.
Finally, if you have any problems then here is a version of the script with the all the answers in: aa_scriptForImaging.py
Optional Extra Step:
STEP 6: SELFCAL THE CONTINUUM
In this step we introduce the technique of self calibration. For bright target sources it is possible to use the model generated in Clean to calibrate (initially) the phases of the target visibilities to improve upon the corrections generated from phase calibration on the phase calibrator. This ultimately leads to an improvement to image quality/SNR.
In this example, we will self calibrate only the solutions in SPW 0 (for time reasons) but the technique can be extended to all SPWs. A complete script for preforming selfcalibration is available here (scriptForSelfCal.py). Note: The first steps in the script deletes any existing model columns in the calibrated measurement set, so don't start this until you're done with the previous 5 steps.
Self cal is an iterative process and in this example we preform two cycles of phase only calibration on the target. The first leads to a large SNR improvement and the second only a marignal one, which tells us we've probably achieved as much as we can. We then finish with an amplitude and phase calibration.
Footnotes:
1) Mysteps 18, 19, 20, 21 in the code will run full cube imaging for the non-continuum subtracted data, run this later... it takes ages!
2) The `t' notes that it has superceded an older clean task, simply known as clean. At some point in the future clean will be deprecated and tclean will replace it.
3) To inspect for other lines we can use the full cubes without continuum subtraction or make some new ones with the continuum subtracted data. Mysteps 22, 23, 24, 25 will run 'full' cubes from the continuum subtracted data... as per footnote 1, run this later as it will take a while!
This workshop has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 730562 [RadioNet]
