VLBI Data Reduction and Imaging Tutorial
Advanced Tutorial 3: Introduction to VLBI data reduction and imaging
Originally developed for the 10th European Radio Interferometry School in Granada (IAA-CSIC, September 30 – October 4, 2024)
Overview
This is a hands-on introduction to VLBI data reduction using CASA. The tutorial dataset will be updated with data from the e-EVN observation conducted during the school.
Note: This is a simplified introduction. For more detailed guides, please refer to:
Data Download
https://archive.jive.eu/scripts/arch.php?exp=RSM07
You can go to the fits idi page where you will find the information to download the FITS-IDI files containing the data.
The credentials for accessing the data are:
username: rsm07
password: m7PQOTdy2Mgl
Download the a-priori flag table from here.
Download the a-priori calibration table from here.
Dataset Information
Observation: RSM07 continuum dataset
Type: EVN network monitoring test with phase-referencing
Duration: 3 hours
Method: Telescopes nodding between calibrator (C) and target (T) in C-T-C-T pattern
| Source Type | Source Name |
|---|---|
| Calibrator (C) | J1848+3219 |
| Target (T) |
3C395 |
| Fringe-finder (FF) |
3C345 |
Data Reduction in CASA
For extra details on what these steps do, check the slides for the calibration and flagging.
Setup
Navigate to your working directory and start CASA:
cd ~/ERIS/eris-evn/
casaStep 1: Loading and Inspecting Data
Load the data:
default(importfitsidi)
fitsidifile= ['rsm07_1_1.IDI1','rsm07_1_1.IDI2']
vis='rsm07.ms'
constobsid=True
scanreindexgap_s=15
importfitsidi()
Load flag table:
default(flagdata)
vis='rsm07.ms'
mode='list'
inpfile='rsm07_casa.flags'
flagdata()Note that the rsm07_casa.flags file that we provide here is just the ANTAB file provided in the EVN Archive (under Pipeline) and converted to CASA format, which can be done via scripts following the steps in the EVN data reduction guide. Here we have already done it for you.
Understand your data:
default(listobs)
vis='rsm07.ms'
listfile='rsm07.listobs'
overwrite=True
listobs()The output will provide to you the information about the sources, antennas that were participating in the observation, and the different scans times, and frequency setup.
Inspect uncalibrated data:
default(plotms)
vis='rsm07.ms'
xaxis='frequency'
yaxis='amp'
field='3C345'
avgtime='3600'
antenna='EF&*'
gridrows=3
gridcols=3
iteraxis='antenna'
coloraxis='scan'
plotms()Step 2: A Priori Amplitude Calibration
Generate system temperature information:
default(gencal)
vis='rsm07.ms'
caltable='rsm07.tsys'
caltype='tsys'
uniform=False
gencal()
Generate gain curves (SKIP THIS STEP!):
default(gencal)
vis="rsm07.ms"
caltable= "rsm07.gcal"
caltype="gc"
gencal()WARNING: instead of creating the gain curve table, we will be using the one linked above (download it from here and uncompress). There is a small bug in the GC information that we appended in this observation.
Inspect Tsys corrections:
default(plotms)
vis='rsm07.tsys'
xaxis='frequency'
yaxis='tsys'
gridrows=3
gridcols=3
coloraxis='corr'
iteraxis='antenna'
plotms()
xaxis='time'
plotms()As we are analyzing the project before the observing run actually finished, you will see that all Tsys values are the same along time, but different for each telescope. This is because the information is only generated after the run finishes. The Tsys values there are proportional to the sensitivity of each antenna.
Step 3: Data Inspection and Flagging
Plot amplitude vs frequency on sensitive baseline:
default(plotms)
vis='rsm07.ms'
xaxis='channel'
yaxis='amp'
field = '3C345'
avgtime='3600'
antenna='EF&JB'
iteraxis='spw'
gridrows=1
gridcols=4
plotms()
Flag band edges:
default(flagdata)
vis='rsm07.ms'
mode='manual'
spw='*:0~3;60~63'
flagdata()We remove the channels that are at the edge of each subbands, as the sensitivity there quickly drops to zero.
Flag autocorrelations:
default(flagdata)
vis='rsm07.ms'
mode='manual'
autocorr=True
flagdata()
Flag Onsala (O8) as did not observe due to winds:
default(flagdata)
vis='rsm07.ms'
mode='manual'
antenna='O8'flagdata()
Step 4: Fringe Fitting
Identify good timerange and check data:
default(plotms)
vis='rsm07.ms'
xaxis='frequency'
yaxis='phase'
ydatacolumn='data'
antenna='EF&*'
coloraxis='corr'
timerange='15:48:00~15:50:00'
averagedata=True
avgtime='120'
iteraxis='antenna'
gridrows=3
gridcols=3
plotms()You can check different scans on the fringe finder and the time range to pick would be a few minutes that contain data that look stable and without bad data.
Calculate delays (SBD - Single Band Delays):
default(fringefit)
vis='rsm07.ms'
caltable='rsm07.sbd'
timerange='15:48:00~15:50:00'
solint='inf'
zerorates=True
refant='EF'
minsnr=10
gaintable=['rsm07.gcal','rsm07.tsys']
interp=['nearest','nearest,nearest']
parang=True
fringefit()
Apply SBD solutions:
default(applycal)
vis='rsm07.ms'
gaintable=['rsm07.gcal','rsm07.tsys','rsm07.sbd']
interp=['nearest','nearest,nearest','nearest']
parang=True
applycal()
Check corrected data:
tget(plotms)
ydatacolumn='corrected'
plotms()Phases should not look around zero for the time when we corrected via the SBD.
Multi-band delays and rates (MBD):
default(fringefit)
vis='rsm07.ms'
caltable='rsm07.mbd'
field='J1848+3219,3C345'
solint='120s'
zerorates=False
refant='EF,MC'
combine='spw'
minsnr=7
gaintable=['rsm07.gcal','rsm07.tsys','rsm07.sbd']
interp=['nearest','nearest,nearest','nearest']
parang=True
fringefit()
Inspect MBD solutions:
default(plotms)
vis='rsm07.mbd'
xaxis='time'
yaxis='delay'
gridrows=4
gridcols=4
coloraxis='corr'
iteraxis='antenna'
plotms()
Apply phase referencing:
default(applycal)
vis='rsm07.ms'
field='J1848+3219,3C395,3C345'
gaintable=['rsm07.gcal','rsm07.tsys','rsm07.sbd','rsm07.mbd']
interp=['nearest', 'nearest,nearest','nearest','linear']
spwmap=[[],[],[],[0,0,0,0]]
parang=True
applycal()Step 5: Bandpass Calibration
Calculate bandpass:
default(bandpass)
vis='rsm07.ms'
caltable='rsm07.bpass'
field='3C345'
gaintable=['rsm07.gcal','rsm07.tsys','rsm07.sbd','rsm07.mbd']
interp=['nearest','nearest,nearest','nearest','linear']
solnorm=True
solint='inf'
refant='EF'
bandtype='B'
spwmap=[[],[],[],[0,0,0,0]]
parang=True
bandpass()
Inspect bandpass solutions:
default(plotms)
vis='rsm07.bpass'
xaxis='frequency'
yaxis='amp'
coloraxis='corr'
iteraxis='antenna'
plotms()
Apply all calibrations:
default(applycal)
vis='rsm07.ms'
field='J1848+3219, 3C395'
gaintable=['rsm07.gcal','rsm07.tsys','rsm07.sbd','rsm07.mbd','rsm07.bpass']
interp=['nearest','nearest,nearest','nearest','linear','nearest,nearest']
spwmap=[[],[],[],[0,0,0,0],[]]
parang=True
applycal()Step 6: Export Data
Create calibrated datasets:
# Calibrator
default(mstransform)
vis='rsm07.ms'
outputvis='rsm07_J1848+3219.ms'
field='J1848+3219'
datacolumn='corrected'
keepflags=True
chanaverage=True
chanbin=64
mstransform()
# Target
default(mstransform)
vis='rsm07.ms'
outputvis='rsm07_3C395.ms'
field='3C395'
datacolumn='corrected'
keepflags=True
chanaverage=True
chanbin=64
mstransform()Step 7: Imaging the data
For more details on this process, look at the slides on imaging and self-calibration.
tclean(vis='***',
imagename='***',
specmode='mfs',
niter=1500,
cycleniter=300,
threshold=0,
imsize=['**','**'],
cell='**',
weighting='**',
deconvolver='**',
savemodel='modelcolumn',
interactive=True)In MacOS:
from casagui.apps import run_iclean
run_iclean(vis=target+'.ms',
imagename='***',
specmode='mfs',
niter=1500,
cycleniter=300,
threshold=0,
imsize=['**','**'],
cell='**',
weighting='**',
deconvolver='**',
savemodel='modelcolumn')
ft(vis=target+'.ms',
model='***.model',
usescratch=True)
Step 6: Self Calibration
After imaging, the MS will now contain the model that you computed for your source.
You can look at the model in plots:
plotms(vis=target+'.ms', xaxis='uvdist', yaxis='amp', ydatacolumn='model', correlation='RR,LL', avgchannel='16', coloraxis='spw', plotfile='', overwrite=True)
Then you can derive the calibration table. First start with a long solution interval and to correct only for the phases.
gaincal(vis=target+'.ms', calmode='p',caltable='target-sc-p1', field=target, solint='10min',refant='EF', minblperant=3, minsnr=5)
The calibration table should contain corrections that are on a small level, and you can check with:
plotms(vis='target-sc-p1', gridcols=2, gridrows=3, xaxis='time', yaxis='phase', iteraxis='antenna', coloraxis='spw', plotrange=[-1,-1,-180,180])
Apply the calibration:
applycal(vis=target+'.ms', gaintable=['target-sc-p1'])
Then you can repeat the cycle of imaging (it should look slightly better now), and self calibration in phase with a shorted interval (2min, as you do not expect faster changes due to the atmosphere).
After another imaging round, if the solutions (in plots) seem to converge to just small improvements, one can run the self-calibration in amplitude and phase. In this case, running gaincal with calmode='ap' . For amplitude, you do not expect fast changes, but only on the order of ~30 min or more. Therefore a longer solution interval of half an hour or one hour is enough.
The final image with this corrections should look much better!
