<aside> <img src="/icons/reorder_gray.svg" alt="/icons/reorder_gray.svg" width="40px" />

Menu

</aside>

Overview

Logistics

Program

Reaching the venue

Before the workshop

After the workshop

Hands-on tutorials

Tutorial 1: Visualization with CARTA

Tutorial 2: Continuum Imaging

Tutorial 3: Line Imaging

Tutorial 4: Self-Calibration

Tutorial 5: ALMA Archive

Lectures

Lecture 1: ALMA Overview

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 tutorial3_script_line_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_Line_Imaging_CASA_6.6.1>
Adapted for ALMA Data Reduction Workshop at Universidad de Chile 
Pietro Curone ([email protected]) 29 October 2025

</aside>

Dataset

TW Hya self-calibrated N2H+ line+cont measurement set (435 MB, already contained in the twhya_firstlook.tar you can download from the Before the workshop page)

twhya_selfcal.ms.tar (Download)

Your working directory should look like:

Tutorial3_line_imaging >> ls
tutorial3_script_line_imaging.py
twhya_selfcal.ms.tar

Untar with:

# In bash
tar -xvzf twhya_selfcal.ms.tar

Tutorial3_line_imaging >> ls
tutorial3_script_line_imaging.py
twhya_selfcal.ms
twhya_selfcal.ms.tar

Open CASA

Open a CASA terminal

Tutorial3_line_imaging >> casa

A Python prompt will appear like this:

CASA <1>:

Import useful libraries

# In CASA
import os
import numpy as np

Getting information on the dataset

listobs(vis='twhya_selfcal.ms',
        listfile='twhya_selfcal.ms.txt',
        overwrite=True)

Prepare to subtract the continuum

View amplitude vs. channel:

plotms(vis='twhya_selfcal.ms',
       xaxis='channel',
       yaxis='amp',
       field='TW Hya',
       avgspw=False,
       avgtime='1e9',
       avgscan=True,
       avgbaseline=True,
       showgui=False,
       plotfile='twhya_selfcal.ms_amp-v-channel.png')

Note that the amplitudes away from the spectral line are well above zero, showing the continuum emission.

View amplitude vs. frequency (LSRK, corrected for Earth’s movement):

plotms(vis='twhya_selfcal.ms',
       xaxis='frequency',
       yaxis='amp',
       field='TW Hya',
       avgspw=False,
       avgtime='1e9',
       avgscan=True,
       avgbaseline=True,
       transform=True,            
       freqframe='LSRK',          # LSRK frame
       showgui=False,
       plotfile='twhya_selfcal.ms_amp-v-frequency.png')

We know that N2H+ line is at 372.67249 GHz

Perform continuum subtraction

os.system('rm -rf twhya_selfcal.ms.contsub')

uvcontsub(vis = 'twhya_selfcal.ms',
          field = 'TW Hya',
          fitspec = '0:0~209;271~383',
          fitorder = 0,
          outputvis='twhya_selfcal.ms.contsub')

Inspect the continuum-subtracted data

listobs(vis='twhya_selfcal.ms.contsub',
        listfile='twhya_selfcal.ms.contsub.txt',
        overwrite=True)

View amplitude vs. channel again:

plotms(vis='twhya_selfcal.ms.contsub',
       xaxis='channel',
       yaxis='amp',
       field='TW Hya',
       avgspw=False,
       avgtime='1e9',
       avgscan=True,
       avgbaseline=True,
       showgui=False,
       plotfile='twhya_selfcal.ms.contsub_amp-v-channel.png')

Note that because the amplitude is positive definite, you don't expect the line-free channels to go entirely down to zero. They should just average to a very small value that reflects the noise.

Image the science target

Line imaging is like continuum imaging with the additional dimension of spectral information.

Set the rest frequency of the N2H+ molecular transition we are imaging:

restfreq = '372.67249GHz' # N2H+ rest frequency <https://splatalogue.online/#/basic>

To auto-mask and image non-interactively:

imagename = 'twhya_n2hp'

os.system(f'rm -rf {imagename}.*') # Remove old versions in case you've run this before

tclean(vis='twhya_selfcal.ms.contsub',
       imagename=imagename,
       field='TW Hya',
       spw='0',
       specmode='cube',
       nchan=15,
       start='0.0km/s',
       width='0.5km/s',
       outframe='LSRK',
       restfreq=restfreq,
       deconvolver= 'hogbom',
       gridder='standard',
       imsize=250,
       cell='0.1arcsec',
       weighting='briggsbwtaper',
       robust=0.5,
       restoringbeam='common',   # Ensure a constant beam size for each plane in the cube
       niter=1000000,    # Just choose a high number
       nsigma = 3, # CLEAN down to a threshold of 3 times the rms noise
       usemask='auto-multithresh', # Automasking
       noisethreshold=4.0,  # Automasking parameter, lower than default (5.0) because emission here is faint
       fastnoise=False,   # Better noise calculation
       interactive=False)

To export to a fits file

exportfits(imagename = f'twhya_n2hp.image',
           fitsimage = f'twhya_n2hp.fits')

Perform primary beam correction

os.system('rm -rf twhya_n2hp.pbcor.image')
impbcor(imagename='twhya_n2hp.image',
        pbimage='twhya_n2hp.pb',
        outfile='twhya_n2hp.pbcor.image')