<aside> <img src="/icons/reorder_gray.svg" alt="/icons/reorder_gray.svg" width="40px" />
Menu
</aside>
Tutorial 1: Visualization with CARTA
Lecture 2: Basics of Interferometry and Imaging
Lecture 3: Introduction to CASA
Lecture 4: Introduction to Calibration and Self-Calibration
Lecture 5: Science Ready Data Products and Weblog
Lecture 6: NRAO/NAASC services
<aside>
Download tutorial2_script_continuum_imaging.py
This script was written for CASA 6.6.6 (ALMA pipeline version)
Source: First Look at Imaging (CASA Guide)
<https://casaguides.nrao.edu/index.php/First_Look_at_Imaging_CASA_6.6.1>
Adapted for ALMA Data Reduction Workshop at Universidad de Chile
Pietro Curone ([email protected]) 29 October 2025
</aside>
TW Hya calibrated continuum measurement set (435 MB, already contained in the twhya_firstlook.tar you can download from the Before the workshop page)
twhya_calibrated.ms.tar (Download)
Your working directory should look like:
Tutorial2_continuum_imaging >> ls
tutorial2_script_continuum_imaging.py
twhya_calibrated.ms.tar
Untar with:
# In bash
tar -xvzf twhya_calibrated.ms.tar
Tutorial2_continuum_imaging >> ls
tutorial2_script_continuum_imaging.py
twhya_calibrated.ms
twhya_calibrated.ms.tar
Open a CASA terminal
Tutorial2_continuum_imaging >> casa
A Python prompt will appear like this:
CASA <1>:
CASA is running on a python interpreter, so python syntax is valid here.
A complete description of CASA and all its tasks is available in the documentation.
Import useful libraries
# In CASA
import os
import numpy as np
listobs(vis='twhya_calibrated.ms',
listfile='twhya_calibrated.ms.txt',
overwrite=True)
Important sections: Scans, Fields, and Spectral Windows
plotms(vis='twhya_calibrated.ms',
xaxis='u',
yaxis='v',
avgchannel='10000',
avgspw=False,
avgtime='1e9',
avgscan=False,
coloraxis='field',
showgui=True)
Save resulting plot as png file directly:
plotms(vis='twhya_calibrated.ms',
xaxis='u',
yaxis='v',
avgchannel='10000',
avgspw=False,
avgtime='1e9',
avgscan=False,
coloraxis='field',
showgui=False,
plotfile='twhya_calibrated_uvcoverage.png')
Let’s plot the uv-coverage of the science target only
plotms(vis='twhya_calibrated.ms',
xaxis='u',
yaxis='v',
field='TW Hya',
avgchannel='10000',
avgspw=False,
avgtime='1e9',
avgscan=False,
showgui=False,
plotfile='twhya_calibrated_uvcoverage_science_target.png')
Given our averaging, each pair of points represent a single baseline, per spw (we have only one), per scan.
Try setting avgscan=True …
Check the Amplitude vs Baseline (= uv-distance) plot:
plotms(vis='twhya_calibrated.ms',
xaxis='UVdist',
yaxis='Amp',
avgchannel='10000',
avgspw=False,
avgtime='1e9',
avgscan=False,
field='TW Hya',
showgui=False,
plotfile='twhya_calibrated_amp-v-uvdist_science_target.png')
First, in case you've run this task before, let's remove old versions of the dataset that use the same name (CASA tends to not overwrite and to raise a SEVERE error instead).
os.system('rm -rf twhya_science_target_avg.ms')
We average in frequency (8 channels) and time (30s) as we are really only interested in continuum imaging at this point. This reduces the data volume without losing much information.
split(vis='twhya_calibrated.ms', # Old measurement set (435 MB)
outputvis='twhya_science_target_avg.ms', # New measurement set (13 MB)
datacolumn='data',
field='TW Hya',
width='8',
timebin='30s')
Inspect the new MS:
listobs(vis='twhya_science_target_avg.ms',
listfile='twhya_science_target_avg.ms.txt',
overwrite=True)
Let’s estimate the proper cell size and image size:
cell=['0.1arcsec']imsize=250 ) just in caseimagename = 'twhya_cont_dirtyimage'
os.system(f'rm -rf {imagename}.*') # Remove old versions in case you've run this before
tclean(vis='twhya_science_target_avg.ms',
imagename=imagename,
specmode='mfs',
gridder='standard',
deconvolver='hogbom',
cell=['0.1arcsec'],
imsize=250,
weighting='briggs',
robust=0.5,
niter=0, # Dirty image
interactive=False)
The beam size is 0.56’’ x 0.42’’, so our calculation was correct!
Let’s define a circular mask around the source. The center of the source can be established by eye or, more precisely, with a 2D Gaussian fit in CARTA.
Set the mask major axis in CARTA using the ruler.
mask_ra = '11h01m51.837s'
mask_dec = '-34.42.17.216'
mask_pa = 0 # mask position angle in degrees (here 0 as we want a circular mask)
mask_inc = 0 # mask inclination in degrees
mask_maj = 1.3 # semimajor axis of mask in arcsec
mask_min = mask_maj*np.cos(mask_inc * np.pi/180) #semiminor axis of mask in arcsec
mask = f'ellipse[[{mask_ra}, {mask_dec}], [{mask_maj}arcsec, {mask_min}arcsec], {mask_pa}deg]'
imagename = 'twhya_cont_CLEANimage'
os.system(f'rm -rf {imagename}.*') # Remove old versions in case you've run this before
tclean(vis='twhya_science_target_avg.ms',
imagename=imagename,
specmode='mfs',
gridder='standard',
deconvolver='hogbom',
cell=['0.1arcsec'],
imsize=250,
weighting='briggs',
robust=0.5,
mask=mask, # Use the mask we just defined
niter=1000000, # Just choose a high number
nsigma=3, # CLEAN down to a threshold of 3 times the rms noise
fastnoise=False, # Better noise calculation
interactive=False)
To export to a fits file
exportfits(imagename = f'{imagename}.image',
fitsimage = f'{imagename}.fits')
os.system(f'rm -rf {imagename}.pbcor.image')
Perform primary beam correction and generate primary beam corrected image:
impbcor(imagename=f'{imagename}.image',
pbimage=f'{imagename}.pb',
outfile=f'{imagename}.pbcor.image')
usemask='auto-multithresh' instead of mask=mask