Image Continuum CASA 6.5.4

This guide should be used after completing Prepare the data for Imaging. You should have created calibrated_final.ms prior to proceeding. Commands for this guide can be found in scriptForImaging_template.py available on github.

This guide features CARTA, the “Cube Analysis and Rendering Tool for Astronomy,” which is the new NRAO visualization tool for images and cubes. The CASA viewer (imview) has not been maintained for a few years and will be removed from future versions of CASA. We strongly recommend using CARTA, as it provides a much more efficient, stable, and feature rich user experience. A comparison of the CASA viewer and CARTA, as well as instructions on how to use CARTA at NRAO, is provided in the CARTA section of the CASA docs.

Check CASA version

This template is for use in CASA versions 6.5.4 and greater. The following code checks the version of CASA and exits if it is less than 6.5.4. You can download the appropriate version of CASA from Obtaining CASA . See Prepare for Imaging section on the main page of this guide for more information.

# in CASA
import re
try:
    import casalith
except:
    print("Script requires CASA 6.0 or greater")

if casalith.compare_version("<",[6,5,4]):
    print("Please use CASA version greater than or equal to 6.5.4 with this script")

Create an Averaged Continuum MS

Figure 1: spw 0

Figure 2: Plots showing each spectral window for TW Hydra with channel vs amplitude. Spw 0 has a strong spectral line and spw 3 is a continuum only window.

Appropriate averaging of a measurement set can reduce overall data volume, making imaging faster with tclean. Since the continuum image is formed by essentially summing the entire bandwidth, we can spectrally average the input measurement set prior to imaging without significantly affecting the final imaging results. Below, we outline a procedure to create a spectrally averaged measurement set for continuum imaging.

The first step is to identify which spectral windows (spws) you would like to include in the continuum. In general, you should include all spws with bandwidths greater than 250MHz to maximize the available continuum bandwidth. Once you have determined which spws you would like to use to form the continuum ms, set the contspw variable:

# in CASA
finalvis='calibrated_final.ms' # This is your output ms from the data
                               # preparation script.
# Set spws to be used to form continuum
contspws = '0,1,2,3'

# Set the weighting parameters for the upcoming tclean commands. These are discussed in detail in the Image Parameters section below
weighting = 'briggs'
robust = 0.5

The next step is to identify and flag the spectral lines in all spectral windows that you will use to make the continuum image. You will split and average the data with the spectral lines flagged so that your final spectrally averaged continuum ms only contains continuum emission. Two methods commonly used to identical spectral line emission are: 1) a channel vs. amplitude plot of the visibilities and 2) a dirty image of the cube.

To use the first method, create a channel vs. amplitude plot using plotms. An example command is given below. You may wish to change the averaging and the uvrange to identify extended emission.

Figures 1 and 2 show example plots of a spw with and without lines to flag. Figure 1 shows a clear line at channel 1700 whereas Figure 2 shows no line emission and just continuum.

Warning: If you apply channel averaging greater than 1, the numbers displayed on the channel axis will represent the channel bin number, rather than the averaged channel.

# in CASA
plotms(vis=finalvis, xaxis='channel', yaxis='amplitude',
   	ydatacolumn='data',
   	avgtime='1e8', avgscan=True, avgchannel='1',
   	iteraxis='spw' )

The second, and more accurate, method is to use tclean to make a quick dirty image cube of channels with niter set to zero and mode=’channel’. Inspecting this channel cube for line emission gives a better defined channel range to flag. The basic command to create dirty image cubes is given below. If you are going to use this method, you will need to set the imaging parameters before running the tclean command. For now, we will set the parameters specifically, but the link (which takes you to another section of this guide) will explain how we arrive at their values. This command should be repeated for each spw and science field.

#In CASA 
testimagename='testImage'

for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt']:
    rmtables(testimagename+ext)

field=['2'] #list all target fields
spw=['0','1','2','3'] #list all target spw’s
veltype='radio'
cell='0.25arcsec'
imsize=[128,128]
gridder='standard'

for i in field:
      for j in spw:  
            tclean(vis=finalvis,
            	  imagename=testimagename+'Field_'+str(i)+'_spw_'+str(j), 
            	  field=str(i),
      	          spw=str(j),
                  # phasecenter=phasecenter, # uncomment if mosaic and set to appropriate field number
                  # phasecenter='TRACKFIELD' # uncomment if imaging an ephemeris object, the phasecenter needs to be TRACKFIELD, not a field number as above.
      	          specmode='cube',
                  veltype=veltype,
                  nchan=-1,
      	          outframe='lsrk', # velocity reference frame. See science goals
      	          niter=0,
      	          interactive=True,
      	          cell=cell,
      	          imsize=imsize, 
      	          weighting=weighting, 
      	          robust=robust,
                  pbcor=True,
                  restoringbeam='common',
      	          gridder=gridder)

Figure 3: The weight spectrum vs Frequency for the calibrated_final.ms for each science spectral window. Notice how the weights in each individual spectral window do not vary wildly and have no outliers but the overall weight measure is different for each spectral window due to the different observing frequency.

The first step to flagging the spectral line channels in your data is to use the flagmanager task to save the state of the data before any flagging is applied. You will need to revert back to the non spectral line flagged dataset when line imaging is done later on. In addition, if you accidentally flag any data you can easily correct this by restoring the before_cont_flags file using the flagmanager.

# in CASA
flagmanager(vis=finalvis,mode='save',
versionname='before_cont_flags')

Initialize the per-channel (or spectral) weights in the ms using initweights. This step ensures that when the flagged and unflagged channels are combined, the appropriate weighting is given to the final set of averaged channels. Figure 3 shows a representation of the weightspectrum vs. frequency plot you should see after initializing the spectral weights.

# in CASA
initweights(vis=finalvis,wtmode='weight',dowtsp=True)

The exact spectral window and channel ranges to flag in the flagchannels variable needs to be specified. For example, if spws 2&3 have a line between channels 1201 and 2199 and spectral windows 0 and 1 are line-free the command would be:

# in CASA
flagchannels='0:1660~1760,2:2100~2200' # modify the channel range for your dataset

Figure 4: Amp vs channel for spectral window 0 for TW Hydra with the spectral line flagged out.

Next, use the task flagdata to apply these flags. After this is done, check that the flags were applied correctly by using plotms to inspect the flagged ms. Figure 4 shows an example of what you should expect to see.

# in CASA
flagdata(vis=finalvis,mode='manual',
      	spw=flagchannels,flagbackup=False)

# check that flags are as expected, NOTE must check reload on plotms
# gui if its still open.
plotms(vis=finalvis,yaxis='amp',xaxis='channel',
   	avgchannel='1',avgtime='1e8',avgscan=True,iteraxis='spw')

Now you can spectrally average the channels together to reduce the size of the continuum ms. The number of channels that you can average together is limited by what is referred to as "bandwidth smearing". When creating an interferometer image, we assume that the emission is essentially monochromatic. If you average large numbers of channels together, this is no longer an appropriate assumption and results in radial smearing in the image that increases in magnitude from the delay tracking center. For a detailed discussion, consult Bridle and Schwab, 1999, Synthesis Imaging in Radio Astronomy II, 180, 371. For a short derivation of the relevant quantities, see How to Calculate Bandwidth Smearing. Conservative values for averaging channels in various ALMA bands are suggested below. As a rough guide, the number of channels you can average together is:

• Bands 3, 4, 5, and 6: 125MHz/ (Channel Width (MHz))
• Bands 7, 8, 9 and 10: 250MHz/ (Channel Width (MHz))

In general, you should make sure that the number of channels you are averaging together is an integer multiple of the original total number of channels. For example, if you have 128 channels and a channel width of 15.625 MHz in band 6, you can average together 8 channels at a time. The resulting ms will contain 16 final channels each with a channel width of 125MHz.

If you want to be more careful, the following table will tell you the maximum bandwidth each averaged channel can be to avoid bandwidth smearing worse than about 1%. For example, at Band 7 (373 GHz) in a compact configuration (Bmax=500m) bandwidth smearing is relatively unimportant even for wide bandwidths (2.1 GHz), whereas for extended configurations (Bmax=10km) the maximum bandwidth at the same Band 7 should be of order 100 MHz or less.

This table lists the maximum bandwidth allowed for a reduction in peak response to a point source over the field of view of 1% for a a square and Gaussian bandpass for various observing frequencies and baselines for different two bandpass types. See How to Calculate Bandwidth Smearing for more information.

Finally, the task split is used to average the channels together to produce the spectrally averaged continuum data set.

# in CASA
contvis='calibrated_final_cont.ms'
rmtables(contvis)
os.system('rm -rf ' + contvis + '.flagversions')
split(vis=finalvis,
     spw=contspws,      
     outputvis=contvis,
     width=[256,8,8,8], # number of channels to average together. The final channel width should be less than 125MHz in Bands 3, 4, and 6 
     # and 250MHz in Bands 7, 8, 9 and 10.
     datacolumn='data')

Now you should check the weights of the new continuum measurement set. The ratio of the weights value between Time Domain Mode (TDM) and Frequency Domain Mode (FDM) windows should be approximately equal to the ratio of the channel widths. For more information on the correlator modes TDM and FDM, see section 5.1 of the ALMA Technical Handbook. In addition, any heavily flagged channels should have their weights scaled by the ratio of the unflagged bandwidth to the bandwidth per output channel. For more information about data weights, see the Data Weights and Combination guide.

# in CASA
# update the antenna and field parameters for your dataset
plotms(vis=contvis, yaxis='wtsp',xaxis='freq',spw='',antenna='DA42',field='2')

Finally, we need to use the flagmanager tasks to restore the ms file to its original unflagged state, so that later we can do continuum subtraction and line imaging.

Figure 5: Amplitude vs UVDistance plot of the continuum, colored by spw.

# in CASA
# If you flagged any line channels, restore the previous flags
flagmanager(vis=finalvis,mode='restore',
        	versionname='before_cont_flags')

It is a good practice to inspect the final continuum ms to make sure that it is correct. In general, the distribution of amplitude vs. uv distance should be smooth. Recall this plot will be a horizontal line for a point source, while the short uv-distances will have higher amplitudes for a more extended object.

# in CASA
plotms(vis=contvis,xaxis='uvdist',yaxis='amp',coloraxis='spw')

Image Parameters

Before imaging, you should set a few key parameters that you will use throughout the rest of the script.

Set the field id for the science target you are interested in imaging. You can use the listobs file generated during the imaging prep step to find this information. If you are imaging a mosaic then you will have to select all the fields. For example, field = ‘3~25’. You can also set the field variable to the name of the source, as long as it matches the name in the listobs file. Do not leave the field parameter blank. If you leave the field parameter blank, tclean will attempt to image all sources in the measurement set.

# in CASA
# update for your ms
field='2'

Uncomment the gridder that is relevant to your dataset. Set the phase center by field id or coordinates if you are imaging a mosaic. Check the spatial setup in the weblog (for pipeline calibrations) or qa2 report (for manual calibrations) to find the phase center. You should choose the central field for the phase center in order to get the best results.

If you are unsure which field to use for the phase center after looking at the weblog then you may use the following Analysis Utilities command:

au.pickCellSize('calibrated_final.ms',imsize=True)

This will return a four element array with that contains the calculated cell size, the X axis number of pixels, the Y axis number of pixels, and the field number that is most centered in the mosaic. You may use this as the phase center field id in your script. If you haven't installed Analysis Utilities, see Obtaining Analysis Utilities for instructions.

Another method of finding the central field number uses the Analysis Utilities plotmosaic au.plotmosaic(‘calibrated_final.ms’). This method will produce a plot of all fields with their corresponding field numbers plotted on the sky. You can also set the phase center with the coordinates at the center of this plot.

For ephemeris objects such as planets, the outframe should be set to phasecenter = 'TRACKFIELD', not a field number as above.

# in CASA
# gridder='standard' # uncomment if single field 
# gridder='mosaic' # uncomment if mosaic or if combining one 7m and one 12m pointing.
# phasecenter=3 # uncomment and set to field number for phase
                # center. Note lack of ''.  Use the weblog to
                # determine which pointing to use. Remember that the
                # field ids for each pointing will be re-numbered
                # after your initial split. You can also specify the
                # phase center using coordinates, e.g.,
                # phasecenter='J2000 19h30m00 -40d00m00'.
# phasecenter = 'TRACKFIELD' # If imaging an ephemeris object (planet, etc), the phasecenter needs to be TRACKFIELD, not a field number as above.

Figure 6: Amplitude versus UV wave for the calibrated_final.ms showing the maximum uv wave point for determining cell size.

Next, determine the cell (or pixel) size. To do this, you need to determine the approximate resolution of your observations. Then you divide the resolution in arcsec by 5 to 8 to adequately sample the PSF. The resolution of a particular interferometer observations can be estimated from the length of the longest baseline:

resolution (radian) ~ (observed wavelength) / (length of longest baseline)

If the baseline is measured in wavelengths, this becomes

resolution (radian) ~ 1.0 / (length of longest baseline in wavelengths)

To convert from radians to arcsec, we multiple by 206265.0 to obtain:

resolution(arcsec) ~ 206265.0/(longest baseline in wavelengths) 

To determine the longest baseline in wavelengths use plotms with xaxis=’uvwave’ and yaxis=’amp’ on your ms file. Figure 6 shows an example of this plot. It is generally better to oversample your beam than to undersample, particularly for observations with poor uv-coverage.

cellsize(arc seconds) = resolution(arcsec)/7 

The next step is to determine the image size in pixels. There are two methods of doing this depending on if you are imaging one field or a mosaic.

If this is a single field, the image size can usually be approximated to be the same size as the primary beam of the telescope. The ALMA 12m primary beam in arcsec scales as 6300 / nu[GHz] and the ALMA 7m primary beam in arcsec scales as 10608 / nu[GHz], where nu[GHz] is the sky frequency. However, if there is significant point source and/or extended emission beyond the edges of your initial images, you should increase the imsize to incorporate more emission.

For mosaics, make the image substantially larger than the mosaic footprint. Use au.plotmosaic(finalvis) to find the mosaic footprint in arcseconds. Padding the imsize substantially avoids artifacts in the image. You should be able to see the edges of the outside fields being cut off in an appropriately padded image. In Figure 7 you will see a mosaic field under padded with the edges not visible. Figure 8 shows an image made with a larger imsize of the same mosaic with the edges clearly visible.

Figure 7

Figure 8: A mosaic of Antennae showing under padding and correctly padded images. For the correctly padded image you can see the edges of the primary beam field whereas the under padded does not cut off.

Note that the imsize parameter is in PIXELS, not arcsec, so you will need to divide the image size in arcsec by the pixel size to determine a value for imsize.

If you would like to check any of these calculations you may use the following command, au.pickCellSize('calibrated_final.ms', imsize=True), in CASA. If you haven't installed Analysis Utilities, see Obtaining Analysis Utilities for instructions. This will return a four element array of the calculated cell size, the x axis imsize, the y axis imsize, and the central field id number.

# in CASA
cell='0.25arcsec' # cell size for imaging.
imsize = [128,128] # size of image in pixels.

When imaging lines, you will need to set two specific velocity parameters called outframe and veltype. Outframe is the coordinate system used for the observation. If you have access to the original proposal, this can be found in the Observing Tool (OT) under field setup. A list of acceptable outframes that can be used in CASA can be found at https://help.almascience.org/kb/articles/what-are-the-frequency-reference-frames-in-casa. Note: heliocentric(hel) is deprecated in CASA. Use barycentric(bary) in this case. The most common choices are 'bary' and 'lsrk'. Usually 'bary' is used for sources where z>0.2 ('extragalactic") and 'lsrk is used for 'galactic' sources. For ephemeris objects, the outframe should be set to a blank string, for example outframe = '', as you have likely already regridded to the source velocity in cvel() or can allow tclean to do it on the fly.

You will also have to set the veltype for the tclean command. This variable has only two options available, radio and optical. Due to an interaction between the ALMA Observing Tool and CASA, set the veltype to radio. Even if the object has an optically defined velocity, the sensitivity calculation uses the radio definition. This will avoid confusion in comparing the achieved sensitivity to the expected sensitivity.

# in CASA
outframe='lsrk' # velocity reference frame. See science goals.
veltype='radio' # velocity type.

The last four parameters that must be set are associated with how tclean will weight and clean the data. During the imaging process, the individual visibilities are placed on a uv grid and then combined. Weighting determines how tclean will combine the uv gridded to produce an image. There are several weighting schemes that can be used:

• Natural: Each cell in the uv plane is weighted by the number of points in the cell. This weighting scheme is the default in tclean. It results in a lower noise image at the expense of worse angular resolution.

• Uniform: Each cell in the uv plane has the same weight. In this case, the sidelobes will be reduced because the imaging plane is filled in more uniformly. This weighting scheme also gives more weight to longer baselines resulting in higher resolution at the expense of higher noise. It can emphasize visibilities with calibration errors.

• Briggs: This weighting scheme is a combination of natural and uniform weighting. The weighting is controlled by the robust parameter. A robust parameter of 2 gives natural weighting, -2 uniform weighting, and a number in between a combination of natural and uniform. Refer to Brigg's Thesis for more information. Currently, robust is set to 0.5, which gives a good compromise between natural and uniform weighting. You may find, after imaging, that you have to decrease the noise or angular resolution based on the science goals. Playing with the robust parameter can affect your final noise in the image and also the angular resolution. If you are making a mosaic image do not use a robust value smaller than 0 as this may introduce major artifacts in the images including uneven noise across the image. If you choose to do any form of uvtapering to the data in tclean(), set robust to 2 (Natural weighting) to avoid upweighting points that are going to be downweighted by uv-taper.

The parameters niter and threshold provide two ways to stop the tclean process. The niter parameter is the maximum number of iterations allowed. After this limit has been reached, tclean will terminate. The threshold parameter sets a flux density threshold for the tclean. When the peak residual is below this value, tclean terminates. When cleaning interactively, we recommend setting the niter parameter to a large(say 1000), but not too large value (say 10000 or 100000) so that tclean terminates eventually in the case of human error. The threshold can be left at the default of 0.0mJy for interactive tclean. For non-interactive tclean, the recommendation is to set the threshold to several times the noise in the final image. To determine the threshold for cube images, run tclean with niter set to zero and inspect the resulting image in CARTA. In a line-free channel, select a region and look at the statistics panel to determine the noise level. For more detailed instructions, see the “Image the Spectral Line Data” section of EVLA Spectral Line Imaging Analysis IRC+10216.

# in CASA
weighting = 'briggs'
robust=0.5
niter=1000
threshold = '0.0mJy'

For more information on the various options available in tclean, refer to David Wilner’s Imaging and Deconvolution presentation from the 2016 Synthesis Imaging Workshop. The table below provides helpful pointers to relevant slides.

Imaging Parameters

Pixel and Image Size	Slides 40-41
Weighting	Slides 42-48, 60-61
Deconvolution Algorithms	Slides 50-53, 62

Imaging the Continuum

Figure 9: During interactive tclean, the GUI will show an initial dirty image. From this GUI, create the tclean mask.

Figure 10: After a mask is created, the green arrow will be illuminated. Press this to begin the first 100 iterations. Once cleaning is complete, press the red "X" to finish.

Figure 11: The final residual image.

Now that you have set all of the imaging parameters you will need in tclean, you can proceed to imaging the continuum. The First Look at Imaging CASAGuide gives an introduction to cleaning and imaging. You may also find the Automasking Guide useful in allowing tclean to generate the mask used for cleaning automatically.

We provide an abbreviated set of commands here. If you are imaging a mosaic, the phasecenter parameter should be set. Refer to the Image Parameters section of this guide for instructions on how to determine this for your project. Type “help tclean()” in CASA if you would like to explore the possible parameters of tclean. Specode=’mfs’ sets the spectral gridding type to multi-frequency synthesis and creates a continuum image. If you are imaging the aggregate continuum in Band 3 or 4 and have a fractional bandwidth larger than 10%, you should consider using multi-term multi-frequency synthesis (deconvolver='mtmfs' and nterms=2). This tclean mode accounts for spatial spectral index variations and especially important to include if you have a spatially resolved, high S/N source. For more information on multi-term multi-frequency synthesis, see Rau, U., & Cornwell, T.J. 2011, A&A, 532, AA71. Clean interactively as the threshold is set at 0 mJy. The mask parameter may also be added if you have an existing file. You can create a mask from the dirty image using the instructions at Create_a_Clean_Mask_from_Continuum_Image_or_Moment_Cube.

In CASA 5.4 and later, tclean calls with gridder = 'mosaic' have an additional parameter mosweight with a default of True. When mosweight = True, the gridder weights each field in the mosaic independently. The mosweight parameter is particularly important for mosaics with non-uniform sensitivity, with rectangular shapes, or when using more uniform values of robust Briggs weighting. For more information on mosweight, please see the tclean documentation.

# in CASA
contvis = 'calibrated_final_cont.ms'         
contimagename = 'calibrated_final_cont'

# in CASA
tclean(vis=contvis,
       imagename=contimagename,
       field=field,
       #  phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
       #  mosweight=True, # uncomment if mosaic     
       specmode='mfs',
       deconvolver='hogbom', 
       # Uncomment the below to image with nterms>1. Use if fractional bandwidth is >10%.
       #deconvolver='mtmfs',
       #nterms=2,
       imsize = imsize, 
       cell= cell, 
       weighting = weighting,
       robust = robust,
       niter = niter, 
       threshold = threshold,
       interactive = True,
       gridder = gridder,
       pbcor = True,
       usepointing=False)

Figure 9 shows the tclean GUI that will appear. If no source is apparent, no cleaning should be done. Press the red “X” to complete the task. If a source is apparent, create a mask around it using the mouse. Once a mask is created, the green arrow will be illuminated and you can begin the first round of cleaning. Figure 10 shows what a mask would look like for our example source. The logger window gives you vital information on the progression of tclean. Once the cycle is complete, the residuals will appear in the GUI. You should continue to iterate until the region inside the mask matches the noise outside the mask. You may need multiple cycles depending on the complexity of the source. Figure 11 shows an example of when we chose to stop cleaning.

The red "X" will stop tclean where you are, the blue arrow will stop the interactive part of tclean, but continue to clean non-interactively until reaching the number of iterations requested (niter) or the flux density threshold (whichever comes first), and the green circle arrow will clean until it reaches the "iterations" parameter on the left side of the green area. These are the only safe exit buttons to use during tclean. DO NOT CTRL-C OR KILL TCLEAN WHILE IT IS RUNNING. If you do this, it is very likely that your ms will be corrupted.

Figure 12: The final continuum image.

After creating the continuum image, you should now open the image in CARTA. If using NRAO machines, you can open a new terminal tab, cd to the working directory, then type:

# in bash
carta --no_browser

Copy the output URL into a browser to view your CARTA session. Select and load <name.image>, and check the noise and resolution to make sure the results match the expected values. Use the ALMA Sensitivity Calculator to estimate the expected sensitivity.

For each image it creates, tclean generates several images with the name imagename+extension. If you re-run tclean with the same imagename, tclean will use the existing files as a starting point, continuing the tclean where you left off. To start completely from scratch, either change the imagename or delete all the files from the previous tclean run. Note that CASA retains some image information in memory, so to truly delete the images, you need to run the rmtables command. You will also want to remove datacolumns that were added by tclean with clearcal and delmod. See below for an example.

# in CASA
clearcal(contvis)
delmod(contvis)
for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt']:
    rmtables(contimagename+ext)

Once you are happy with your continuum image(s), you can continue to Self-Calibration Template or Spectral Line Imaging Template. If you do not wish to self-calibrate or create line cubes, continue with this guide to create primary beam corrected images and fits files.

Export the images

Use exportfits to create fits files of the *.flux and *.pbcor files.

# in CASA
import glob

myimages = glob.glob("*.pbcor")
for image in myimages:
    exportfits(imagename=image, fitsimage=image+'.fits',overwrite=True)

myimages = glob.glob("*.pb")
for image in myimages:
    exportfits(imagename=image, fitsimage=image+'.fits',overwrite=True)

Create Diagnostic PNGs

The First Look at Image Analysis guide gives an introduction to a variety of options to begin image analysis, including using immoments to create moment maps. The following commands create png files of the continuum image.

# in CASA
os.system("rm -rf *.png")
mycontimages = glob.glob("*mfs*manual.image")
for cimage in mycontimages:
    mymax=imstat(cimage)['max'][0]
    mymin=-0.1*mymax
    outimage = cimage+'.png'
    os.system('rm -rf '+outimage)
    imview(raster={'file':cimage,'range':[mymin,mymax]},out=outimage)

Return to the Main Page