VLA 5 GHz continuum survey of Seyfert galaxies 5.7.2: Difference between revisions

From CASA Guides
Jump to navigationJump to search
 
(21 intermediate revisions by 6 users not shown)
Line 1: Line 1:
[[Category:VLA]] [[Category:Calibration]] [[Category:Imaging]] [[Category:Self-Calibration]] [[Category:Continuum]]
[[Category:VLA]] [[Category:Calibration]] [[Category:Imaging]] [[Category:Self-Calibration]] [[Category:Continuum]]


{{VLA Intro}}
{{Checked 5.7.2}}


This particular tutorial demonstrates how to calibrate a VLA continuum survey, consisting of a flux calibrator, target sources, and their respective phase calibrators. What distinguishes this example is that it includes multiple sources and multiple phase calibrators.
 
Disclaimer: Due to continuous CASA software updates, GUI images may look different on more recent versions than those shown here.
 
 
This particular tutorial demonstrates how to calibrate a VLA continuum survey, consisting of a flux denisty scale calibrator, target sources, and their respective complex gain calibrators. What distinguishes this example is that it includes multiple target sources and multiple complex gain calibrators.




== Importing VLA Archive Data ==
== Importing VLA Archive Data ==


Download your data from the [http://archive.cv.nrao.edu: VLA Archive]; in this example we'll use the publicly available survey [https://archive.nrao.edu/archive/ArchiveQuery?QUERY_ID=9999&QUERY_MODE=Prepare+Download&PROTOCOL=HTML&SORT_PARM=Starttime&SORT_ORDER=Asc&LOCKMODE=PROJECT&SITE_CODE=AOC&DBHOST=CHEWBACCA&PROJECT_CODE=AG733&OBSERVER=&TELESCOPE=ALL&TIMERANGE1=&TIMERANGE2=&PASSWD=&QUERYTYPE=ARCHIVE&SUBMIT=Submit+Query: AG733], which consists of VLA C-array, C-band (5 GHz) continuum observations of a handful of Seyfert galaxies.  
Download your data from the [http://archive.nrao.edu: VLA Archive]; in this example we'll use the publicly available survey [https://archive.nrao.edu/archive/ArchiveQuery?QUERY_ID=9999&QUERY_MODE=Prepare+Download&PROTOCOL=HTML&SORT_PARM=Starttime&SORT_ORDER=Asc&LOCKMODE=PROJECT&SITE_CODE=AOC&DBHOST=CHEWBACCA&PROJECT_CODE=AG733&OBSERVER=&TELESCOPE=ALL&TIMERANGE1=&TIMERANGE2=&PASSWD=&QUERYTYPE=ARCHIVE&SUBMIT=Submit+Query: AG733], which consists of VLA C-configuration, C-band (5 GHz) continuum observations of a handful of Seyfert galaxies.  


With some sensible fiddling of the web interface, you should be able to download two archive formatted files, one from 06-Dec-09 and one from 06-Dec-13, with names like '''AG733_A061209.xp1''' and '''AG733_B061213.xp1''' or '''2006-12/vla2006-12-09.dat''' and '''2006-12/vla2006-12-13.dat'''.
With some sensible fiddling of the web interface, you should be able to download two archive formatted files, one from 06-Dec-09 and one from 06-Dec-13, with names like '''AG733_A061209.xp1''' and '''AG733_B061213.xp1''' or '''2006-12/vla2006-12-09.dat''' and '''2006-12/vla2006-12-13.dat'''.


[[Image:Screenshot-Log_Messages_(ashlesha--users-jgallimo-AG733-casapy.log).png|thumb|Example snippet of listobs output for program AG733.]]
[[Image:PreEVLA_Seyfert_Example_ListObs.png|thumb|Example snippet of listobs output for program AG733.]]


It's helpful to make a subdirectory to store the VLA archive data and the CASA measurement sets before starting CASA. Linux-level commands are given assuming you are using bash as your default shell.
It's helpful to make a subdirectory to store the VLA archive data and the CASA measurement sets before starting CASA. We call this subdirectory AG733 and place the downloaded files in it. Linux-level commands are given assuming you are using bash as your default shell.


<source lang="bash">
<source lang="bash">
# In bash
# In bash
mkdir ~/AG733
cd AG733
mv AG733_* ~/AG733
casa
cd ~/AG733
casapy
</source>
</source>


To load the data in CASA, we employ the '''{{importvla}}''' command.[[Image:casa_vla_antenna_plot.png|thumb|Graphics window produced by {{plotants}}. Shown are the antennas of the VLA for the program AG733.]]
To load the data in CASA, we employ the '''{{importvla}}''' command.[[Image:PreEVLA_Seyfert_PlotAnts.png|thumb|Graphics window produced by {{plotants}}. Shown are the antennas of the VLA for the program AG733.]]


<source lang="python">
<source lang="python">
Line 38: Line 40:
</source>
</source>


The output measurement set, i.e., the data to be worked on, is ''ag733.ms.'' Note that the (interesting) output of {{listobs}} and {{vishead}} go to the CASA message log window.
The output measurement set, i.e., the data to be worked on, is ''ag733.ms.'' Note that the output of {{listobs}} and {{vishead}} go to the CASA message log window.


{{listobs}} is similar to LISTR in AIPS and provides a list of scans, sources and their respective indices, and antennas. You can also get a more useful plot of antennas, akin to the output of the AIPS task PRTAN, using {{plotants}}. The antennas plot aids the selection of a reference antenna; the antenna ''VA05'' looks like a reasonable choice.
{{listobs}} provides a list of scans, sources and their respective indices, and antennas. You can also get a plot of antennas using {{plotants}}. The antennas plot aids the selection of a reference antenna; the antenna ''VA05'' looks like a reasonable choice.


<source lang="python">
<source lang="python">
Line 48: Line 50:
</source>
</source>


<!-- {{plotants bug}}
You may want to '''close''' {{plotants}} before moving on in the script to avoid locked tables.
 
''' <pre style="color: red"> Note: Some older VLA data may give errors in more recent CASA versions when attempting to generate a plot of antennas. If this occurs, use the listobs task and choose an antenna that was located on a pad close to the center of the array (Pads numbers increase with distance from the center of the array, ie. "VLA:_N1" and "EVLA:W2" are pads close to the center of the array, where as, "VLA:_N18" is not).
</pre>'''


It is also currently important to '''close''' {{plotants}} before moving on in the script to avoid locked tables. -->


== Editing the Measurement Set ==
== Editing the Measurement Set ==
Line 61: Line 65:
</source>
</source>


[[plotms|Plotms]] brings up the following graphical display.
[[plotms|Plotms]] brings up a graphical display that allows you to plot various aspects of your data.


[[Image:Screenshot-PlotMS.png|300px|The plotms GUI.]]
[[Image:PreEVLA_Seyfert_PlotMS.png|300px|The plotms GUI with proper field and timerange options selected.]]


In this example, the measurement set was selected by navigating to and highlighting the measurement set ag733.ms (which will look like a directory to your file browser, because measurement sets are actually directory structures). A specific subset of data (the flux calibrator 3C48 = 0137+331 on the first observing day) was selected by specifying in the appropriate "Data -> Selection" boxes,
In this example, the measurement set was selected by navigating to and highlighting the measurement set ag733.ms (which will look like a directory to your file browser, because measurement sets are actually directory structures). A specific subset of data (the flux density scale calibrator 3C48 = 0137+331 on the first observing day) was selected by specifying in the appropriate "Data -> Selection" boxes:


* field = 0137+331, and
* field = 0137+331, and
* timerange = 2006/12/09/0:0:0~23:0:0.
* timerange = 2006/12/09/0:0:0~23:0:0.


For this specific field, we see no obvious bad data. Note that the low points are actually cross-polarization data; if you want to only look at the right and left polarizations, enter
For this specific field, we see no obvious bad data. Note that the low points are actually cross-polarization data; you can see this by entering:


* corr = RR, LL.
* corr = RR, LL


For the sake of this example, we'll clip out the early part of the observations:  
For the sake of this example, we'll clip out the early part of the observations, say, all the points before 04:00:45:


Firstly, box the region to be flagged by using the "Mark Regions" button [[Image:MarkRegionsButton.png]]. Next, go the the "Flag" tab at top, and select "Extend flags" and "Correlation", so that we will delete the RL and LR data as well (despite the fact that we're not plotting it). Then, press the "Go" button, or "Flag" button [[Image:FlagThoseData.png‎]] to zap those data.  
Firstly, box the region to be flagged by using the "Mark Regions" button [[Image:MarkRegionsButton.png]]. Next, go to the "Flag" tab at top, and select "Extend flags" and "Correlation", so that we will delete the RL and LR data as well (despite the fact that we're not plotting it). Then, press the "Go" or "Flag" button [[Image:FlagThoseData.png‎]] to zap those data.  


The {{viewer}} tool provides utilities to flag data in a manner similar to the AIPS task TVFLG; see the tutorial [[Data flagging with viewer]].  
The {{viewer}} tool provides utilities to flag data in a manner similar to the AIPS task TVFLG; see the tutorial [[Data flagging with viewer]].  


If we leave the field and timerange filter boxes blank and plot again, we see a number of visibilities with extremely large amplitudes. By highlighting and locating the bad points (described in detail in the [[data flagging with plotms]] tutorial), we find most of the bad data can be attributed to one antenna, 'VA15', at specific times, and one time range where all antennas contribute. In the text below we flag these points using {{flagcmd}} and {{tflagdata}}; you can also flag data interactively in the {{plotms}} display (though this is generally not encouraged since it does not save backup versions of the flags, and it can be easy to flag data by mistake or miss data that's not plotted).  
If we leave the field and timerange filter boxes blank and plot again, we see a number of visibilities with extremely large amplitudes. If you have trouble seeing the points, adjust their size by going to "Display", setting "Unflagged Points Symbol" to "Custom", and selecting a larger pixel size for the "circle", "square", or "diamond" options and replotting. Highlight these few extremely high amplitude points and flag them as before. Now we'll reload the plot by selecting the "Reload" checkbox and clicking "Plot". When the plot reloads with its new adjusted scale, we notice that there are even more bad points. By highlighting and locating these new bad points (described in detail in the [[data flagging with plotms]] tutorial), we find most of the bad data can be attributed to one antenna, 'VA15', at specific times, and one time range where all antennas contribute. Close inspection will show there is also some bad data on VA05 and VA22 as well. You may need to use the Zoom Tool [[File:casaplotms-zoom.png]] to see all of the bad points clearly. In the text below we flag these points using {{flagdata}}; you can also flag data interactively in the {{plotms}} display (though this is generally not encouraged since it does not save backup versions of the flags, and it can be easy to flag data by mistake or miss data that's not plotted).  


<source lang="python">
<source lang="python">
# In CASA
# In CASA
flagcmd(vis='ag733.ms', inpmode='cmd',
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/09/04:45:10.0~04:45:20.0')
      command=["antenna='VA15' timerange='2006/12/09/04:45:10.0~04:45:20.0'",
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/09/06:12:30.0~06:14:10.0')
                "antenna='VA15' timerange='2006/12/09/06:12:30.0~06:14:10.0'",
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/13/04:32:40.0~04:35:30.0')
                "antenna='VA15' timerange='2006/12/13/04:32:40.0~04:35:30.0'",
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/09/06:12:30.0~06:12:40.0')
                "antenna='VA15' timerange='2006/12/09/06:12:30.0~06:12:40.0'",
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/09/06:13:00.0~06:13:40.0')
                "antenna='VA15' timerange='2006/12/09/06:13:00.0~06:13:40.0'",
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/09/06:13:50.0~06:14:00.0')
                "antenna='VA15' timerange='2006/12/09/06:13:50.0~06:14:00.0'",
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/09/06:14:00.0~06:14:10.0')
                "antenna='VA15' timerange='2006/12/09/06:14:00.0~06:14:10.0'"],  
flagdata(vis='ag733.ms', antenna='VA05', timerange='2006/12/09/04:30:00.0~04:30:10.0')
      action='apply')
flagdata(vis='ag733.ms', antenna='VA22', timerange='2006/12/09/06:28:10.0~06:28:20.0')
tflagdata(vis='ag733.ms', timerange='2006/12/09/06:24:50.0~06:25:00.0')
flagdata(vis='ag733.ms', timerange='2006/12/09/06:24:50.0~06:25:00.0')
</source>
</source>


Notice that you can flag in one line multiple time ranges for a given antenna. Plot the data again to make sure the amplitudes look reasonable (if {{plotms}} is already open, just hit Shift+Plot to reload with the current flags):  
Plot the data again to make sure the amplitudes look reasonable (if {{plotms}} is already open, just hit Shift+Plot to reload with the current flags):  


<source lang="python">
<source lang="python">
Line 107: Line 111:
By default, {{plotms}} plots all of the data, including cross-polarization. All those low-amplitude points in your data are probably not junk. They are real cross-polarization data. Flagging those data will cause grief later on. Take a look at the tutorial [[Data flagging with plotms#Cross-Pol Data of Bright Sources | data flagging with plotms]].
By default, {{plotms}} plots all of the data, including cross-polarization. All those low-amplitude points in your data are probably not junk. They are real cross-polarization data. Flagging those data will cause grief later on. Take a look at the tutorial [[Data flagging with plotms#Cross-Pol Data of Bright Sources | data flagging with plotms]].


Tread lightly!  And, whenever possible, flag with {{tflagdata}} and keep careful track of what you flag in a separate file.
Tread lightly!  And, whenever possible, flag with {{flagdata}} and keep careful track of what you flag in a separate file.


==A Python Interlude==
==A Python Interlude==
Line 114: Line 118:


<source lang="python">
<source lang="python">
ampcallist = ['0137+331']
# In CASA
phasecallist = ['2250+143', '0119+321', '0237+288', '0239-025',
FDS_callist = ['0137+331']
CG_callist = ['2250+143', '0119+321', '0237+288', '0239-025',
                 '0323+055', '0339-017', '0423-013']  
                 '0323+055', '0339-017', '0423-013']  
sourcelist = ['NGC7469', 'MRK0993', 'MRK1040', 'NGC1056', 'NGC1068',
sourcelist = ['NGC7469', 'MRK0993', 'MRK1040', 'NGC1056', 'NGC1068',
               'NGC1194', 'NGC1241', 'NGC1320', 'F04385-0828',
               'NGC1194', 'NGC1241', 'NGC1320', 'F04385-0828',
               'NGC1667']  
               'NGC1667']  
allcals = ampcallist + phasecallist
allcals = FDS_callist + CG_callist
print allcals
print allcals
</source>
</source>


It will also be helpful later to have a way of referencing which phase calibrator should be assigned to which source. We can define the matching explicitly by using a ''dictionary'' variable.
It will also be helpful later to have a way of referencing which complex gain calibrator should be assigned to which source. We can define the matching explicitly by using a ''dictionary'' variable.


<source lang="python">
<source lang="python">
# In CASA
calDict = {'NGC7469':'2250+143',
calDict = {'NGC7469':'2250+143',
           'MRK0993':'0119+321',  
           'MRK0993':'0119+321',  
Line 143: Line 149:
== Calibration ==
== Calibration ==


=== Set the Flux Scale ===
=== Set the Flux Denisty Scale ===
As in AIPS, CASA sets the flux scale using an observation of a flux calibrator (here, 0137+331) and the task {{setjy}}.  For this example, we'll determine the calibration in reference to a model image of the calibrator source. Look at the inputs of {{setjy}} and get it ready to go!
As in AIPS, CASA sets the flux density scale using an observation of a flux density scale calibrator (here, 0137+331) and the task {{setjy}}.  For this example, we'll determine the calibration in reference to a model image of the calibrator source. Look at the inputs of {{setjy}} and get it ready to go!


<source lang="python">  
<source lang="python">  
# in CASA
# In CASA
# perhaps restore setjy back to its defaults
# restore setjy back to its defaults
default setjy
default setjy
# shows several of the current and possible input parameters for the current task
inp
inp
vis = 'ag733.ms'
vis = 'ag733.ms'
Line 159: Line 166:
casapath = os.environ.get('CASAPATH').split()[0]
casapath = os.environ.get('CASAPATH').split()[0]
# generate full path to the image
# generate full path to the image
modimage = casapath + '/data/nrao/VLA/CalModels/3C48_C.im'
model = casapath + '/data/nrao/VLA/CalModels/3C48_C.im'
go
go
</source>
</source>


That example was set up to demonstrate some AIPS heritage. You could equally well run {{setjy}} the following way.
That example was set up to demonstrate some AIPS heritage. You could equally well run {{setjy}} the following way:
 
<source lang="python">  
<source lang="python">  
# in CASA
# In CASA
setjy('ag733.ms', field = '0137+331', modimage = casapath + '/data/nrao/VLA/CalModels/3C48_C.im')
import os
casapath = os.environ.get('CASAPATH').split()[0]
setjy(vis = 'ag733.ms', field = '0137+331', model = casapath + '/data/nrao/VLA/CalModels/3C48_C.im')
</source>
</source>


=== Determine calibration solutions ===
=== Determine calibration solutions ===


At this stage the data have an overall flux scaling determined, but full gain solutions aren't there yet. The relevant task is {{gaincal}} (analogous to the AIPS task CALIB).
At this stage the data have an overall flux density scaling determined, but full gain solutions aren't there yet. The relevant task is {{gaincal}}, which will produce the table ''cal.G.''
 
Rather than generate solution tables (SN tables in AIPS) that are attached to the measurement set, {{gaincal}} will produce a separate table, which we'll call ''cal.G.''


<source lang="python">
<source lang="python">
# In CASA
rmtables('cal.G') # get rid of any old, unmeritorious versions of the calibration
rmtables('cal.G') # get rid of any old, unmeritorious versions of the calibration
rmtables('gc.cal')# will give a warning if the tables don't already exist, which can be safely ignored
# First we'll generate an antenna zenith-angle dependent gain curve calibration table
gencal(vis='ag733.ms', caltype='gc', caltable='gc.cal')
#
#
# Some caution on the following step. The solution interval is set to 'inf,' meaning to average over entire
# Some caution on the following step. The solution interval is set to 'inf,' meaning to average over entire
# scans of each calibrator in turn. If a given scan is however longer than the coherence time (esp. on long baselines),  
# scans of each calibrator in turn. If a given scan is however longer than the coherence time (esp. on long baselines),  
# some decorrelation will occur and the flux calibration will be thrown off. Suggest, e.g., solint='1min' or solint='2min'  
# some decorrelation will occur and the flux density scale calibration will be thrown off. Suggest, e.g., solint='1min' or solint='2min'  
# as reasonable alternatives to solint='inf'. Trial and error may be your friend here.
# as reasonable alternatives to solint='inf'. Trial and error may be your friend here.
#
#
gaincal(vis='ag733.ms',caltable='cal.G', field=','.join(allcals), solint='inf',
gaincal(vis='ag733.ms',caltable='cal.G', field=','.join(allcals), solint='inf',
         refant='VA05', append=False, gaincurve=False)
         refant='VA05', gaintable=['gc.cal'], append=False)
plotcal('cal.G',yaxis='amp',timerange='2006/12/09/04:0:0~06:30:0')
plotms(vis='cal.G', coloraxis='Antenna1', yaxis='amp', xaxis='time', timerange='2006/12/09/04:0:0~06:30:0')
</source>
</source>


There's some python trickery buried in this example. If you are unfamiliar with the string method ''join,'' try the following example. Recall that we had defined the variable allcals above in [[#A_Python_Interlude|A Python Interlude]].
There's some python trickery buried in this example. If you are unfamiliar with the string method ''join,'' try the following example. Recall that we had defined the variable allcals above in [[#A_Python_Interlude|A Python Interlude]].
<source lang="python">
<source lang="python">
# In CASA
print allcals  
print allcals  
print ','.join(allcals)
print ','.join(allcals)
</source>
</source>


It's worth playing with the timerange, etc. in {{plotcal}} to explore the solutions. You may notice that the antenna 'EA26' looks a little junky on the second day. Here's how to flag it manually and then repeat the calibration (be sure to close the {{plotcal}} window displaying the previous solutions first).
It's worth playing with the timerange, etc. in {{plotms}} to explore the solutions. You may notice that antenna 'EA26' does not appear in the solutions on the second day. For completeness, we'll flag it manually from the measurement set and then repeat the calibration.


<source lang="python">
<source lang="python">
tflagdata(vis='ag733.ms',antenna='EA26', timerange='2006/12/13/0:0:0~24:0:0')
# In CASA
flagdata(vis='ag733.ms',antenna='EA26', timerange='2006/12/13/0:0:0~24:0:0')
rmtables('cal.G')  
rmtables('cal.G')  
gaincal(vis='ag733.ms',caltable='cal.G', field=','.join(allcals), solint='inf',
gaincal(vis='ag733.ms',caltable='cal.G', field=','.join(allcals), solint='inf',
         refant='VA05', append=False, gaincurve=False)
         refant='VA05', gaintable=['gc.cal'], append=False )
</source>
</source>


Line 207: Line 221:


<source lang="python">
<source lang="python">
tget plotcal
# In CASA
plotcal()
tget plotms
plotms()
</source>
</source>


Finally, bootstrap the flux scale of the flux calibrator onto the phase calibrators. AIPS called it GETJY, but CASA calls it {{fluxscale}}.
Using the task {{fluxscale}}, bootstrap the flux density scale onto the complex gain calibrators.


<source lang="python">
<source lang="python">
# In CASA
rmtables('cal.Gflx') # remove any old versions of the calibration table, if necessary
rmtables('cal.Gflx') # remove any old versions of the calibration table, if necessary
myFluxes = fluxscale(vis='ag733.ms', caltable='cal.G', reference=','.join(ampcallist),  
myFluxes = fluxscale(vis='ag733.ms', caltable='cal.G', reference=','.join(FDS_callist),  
                     transfer=','.join(phasecallist), fluxtable='cal.Gflx', append=False)
                     transfer=','.join(CG_callist), fluxtable='cal.Gflx', append=False)
</source>
</source>


Here, the calibration table ''cal.G'' is modified and stored as 'cal.Gflx.' The Python dictionary "myFluxes" will also contain information about the flux scaling and other useful tidbits (type "myFluxes" to see its contents).  So far, solutions have been generated only for the calibrators, and they have not yet been transfered to the sources.
Here, the calibration table ''cal.G'' is modified and stored as 'cal.Gflx.' The Python dictionary "myFluxes" will also contain information about the flux density scaling and other useful tidbits (type "myFluxes" to see its contents).  So far, solutions have been generated only for the calibrators, and they have not yet been transfered to the sources.


We could equally well have used reference='0137+337', but use of the ''join'' method allows for the more general possibility of multiple flux calibrators.
We could equally well have used reference='0137+337', but use of the ''join'' method allows for the more general possibility of multiple flux density scale calibrators.


=== Apply the Calibrations ===
=== Apply the Calibrations ===


In AIPS, the task CLCAL applies calibration solutions to a calibration table. In CASA the equivalent task is {{applycal}}.  We loop over sources and calibrators to properly match them:
The task in CASA is {{applycal}}.  We loop over sources and calibrators to properly match them:
<!-- for now, just use applycal, since accum appears to be broken in 3.4 (see CAS-4240)
 
However, we need to apply the calibration ''incrementally,'' so that calibration solutions for sources are properly matched with their respective phase calibrator.  The task {{accum}} performs incremental calibration.
 
In this example, the solutions in ''cal.Gflx'' will be incrementally compiled into a new calibration table called ''multisource.gcal.''
 
<source lang="python">
# first, generate a "blank" table.
accum(vis='ag733.ms',tablein='',accumtime=10.0,incrtable='cal.Gflx',
      interp='linear', caltable='multisource.gcal',
      field=','.join(allcals), calfield=','.join(allcals))
 
# now accumulate solutions onto that table
# loop over the source and calibrator pairs in the calDict dictionary variable
for source, calibrator in calDict.iteritems():
    accum(vis='ag733.ms',tablein='multisource.gcal',accumtime=10.0,
          incrtable='cal.Gflx',
          interp='linear', field=source,
          calfield=calibrator)
 
# Finally, apply the calibrations to the measurement set
default(applycal)
applycal(vis='ag733.ms', gaintable='multisource.gcal',
        gaincurve=False)
</source>
 
-->


<source lang="python">
<source lang="python">
# In CASA
for source, calibrator in calDict.iteritems():
for source, calibrator in calDict.iteritems():
     applycal(vis='ag733.ms', field=source,
     applycal(vis='ag733.ms', field=source,
             gaintable='cal.Gflx', gainfield=calibrator,
             gaintable='cal.Gflx', gainfield=calibrator)
            gaincurve=False)
</source>
</source>


=== Baseline-based corrections ===
=== Baseline-based corrections ===


[[File:NGC1068_without_blcal.png|thumb|Initial clean image of NGC1068 ''without'' blcal.]]  
[[File:PreEVLA_Seyfert_IntClean_Noblcal.png|thumb|Initial tclean image of NGC1068 ''without'' blcal.]]  
[[File:NGC1068_with_blcal.png|thumb|Initial clean image of NGC1068 ''with'' blcal applied.]]
[[File:PreELA_Seyfert_IntClean_blcal.png|thumb|Initial tclean image of NGC1068 ''with'' blcal applied.]]


Although it makes only modest improvements for the present example, observations including a mix of EVLA and VLA telescopes may benefit from performing baseline-based corrections on a strong calibrator. Now that we have initial antenna-based solutions, we can use {{blcal}}.
Although it makes only modest improvements for the present example, observations including a mix of EVLA and VLA telescopes may benefit from performing baseline-based corrections on a strong calibrator. Now that we have initial antenna-based solutions, we can use {{blcal}}.


<source lang="python">
<source lang="python">
# In CASA
default('blcal')
default('blcal')
vis = 'ag733.ms'
vis = 'ag733.ms'
# output baseline-based calibration solutions
# output baseline-based calibration solutions
caltable = 'cal.BL'  
caltable = 'cal.BL'  
# use the strong flux calibrator to determine baseline-based solutions
# use the strong flux density scale calibrator to determine baseline-based solutions
field = '0137+331'  
field = '0137+331'  
# generate a solution for each calibrator scan
# generate a solution for each calibrator scan
Line 282: Line 272:
# calibrator for the BL calibrator; in this example, gain cal and bl cal are the same
# calibrator for the BL calibrator; in this example, gain cal and bl cal are the same
interp = 'nearest'
interp = 'nearest'
gaincurve = False
# do the first day's data
# do the first day's data
timerange = '2006/12/09/0:0:0~24:0:0'
timerange = '2006/12/09/0:0:0~24:0:0'
Line 291: Line 280:
</source>
</source>


This operation generates the calibration table ''cal.BL''. Solutions can be inspected with {{plotcal}}. Now we have to use {{applycal}} again to apply both the antenna-based and baseline-based calibration solutions to the data.
This operation generates the calibration table ''cal.BL''. Solutions can be inspected with {{plotcal}} (this functionality may transfer to {{plotms}} in the near future). Now we have to use {{applycal}} again to apply both the antenna-based and baseline-based calibration solutions to the data.


<source lang="python">
<source lang="python">
# In CASA
for source, calibrator in calDict.iteritems():
for source, calibrator in calDict.iteritems():
     applycal(vis='ag733.ms', field=source,
     applycal(vis='ag733.ms', field=source,
             gaintable=['cal.Gflx','cal.BL'], gainfield=calibrator,
             gaintable=['cal.Gflx','cal.BL'], gainfield=calibrator)
            gaincurve=False)
</source>
</source>


Did the {{blcal}} operation make a difference here? Getting ahead of ourselves a little, the plots at right show the initial {{clean}} images of NGC1068, one with {{blcal}} applied, and one without, and both images are displayed with the same stretch. There does appear a slight reduction of artifacts with {{blcal}} applied, but of course your mileage will vary.
Did the {{blcal}} operation make a difference here? Getting ahead of ourselves a little, the plots at right show the initial {{tclean}} images of NGC1068, one with {{blcal}} applied, and one without, and both images are displayed with the same stretch. There does appear a slight reduction of artifacts with {{blcal}} applied, but of course your mileage will vary.


=== Splitting the calibrated source data from the multisource measurement set ===
=== Splitting the calibrated source data from the multisource measurement set ===
Line 307: Line 296:


<source lang="python">
<source lang="python">
# In CASA
splitfile = 'NGC1667.split.ms'
splitfile = 'NGC1667.split.ms'
rmtables(splitfile) # get rid of any old versions before splitting
rmtables(splitfile) # get rid of any old versions before splitting
Line 312: Line 302:
</source>
</source>


Alternatively, you can take advantage of scripting to loop over the sources.
Instead of having to manually repeat the split task for each source, you can take advantage of scripting to loop over the sources:


<source lang="python">
<source lang="python">
# In CASA
for source in sourcelist:
for source in sourcelist:
     splitfile = source + '.split.ms'
     splitfile = source + '.split.ms'
Line 324: Line 315:
== Imaging ==
== Imaging ==


[[Image:NGC1667_example.png|thumb|An example desktop display of ''{{viewer}},'' showing the image display panel, the data display options window, and the data loader.]]
[[Image:PreEVLA_Seyfert_ExampleViewer.png|thumb|An example desktop display of ''{{viewer}},'' showing the image display panel, the data display options window, and the data loader.]]


Use {{clean}} for [http://www.aoc.nrao.edu/events/synthesis/2008/lectures/wilner_synthesis08.pdf imaging] and [http://adsabs.harvard.edu/abs/1974A%26AS...15..417H clean deconvolution]. Like IMAGR in AIPS, it's worth some time mulling over all of the options. A complete list of keywords is given in the example, below. Most of the options should be relatively familiar to the AIPS user, but one keyword bears mentioning in this example: '''calready = True'''. By selecting this option, the CLEAN model is transformed back to (''u'',''v'') space and stored as as a new column in the visibility data. The transformed CLEAN model will be necessary for [[#Self-calibration|self calibration]], described below.
Use {{tclean}} for [http://www.aoc.nrao.edu/events/synthesis/2008/lectures/wilner_synthesis08.pdf imaging] and [http://adsabs.harvard.edu/abs/1974A%26AS...15..417H clean deconvolution]. Like IMAGR in AIPS, it's worth some time mulling over all of the options. A list of some of the various keywords is given in the example, below. Most of the options should be relatively familiar to the AIPS user.  <!--, but one keyword bears mentioning in this example: '''calready = True'''. By selecting this option, the CLEAN model is transformed back to (''u'',''v'') space and stored as as a new column in the visibility data. The transformed CLEAN model will be necessary for [[#Self-calibration|self calibration]], described below. -->


For this snapshot survey, the example chooses a fairly conservative, slow clean to accommodate the less-than-ideal beam. You could probably get away with more efficient cleaning, but it's better to play it safe in scripted self-calibration.
For this snapshot survey, the example chooses a fairly conservative, slow clean to accommodate the less-than-ideal beam. You could probably get away with more efficient cleaning, but it's better to play it safe in scripted self-calibration.


<source lang="python">
<source lang="python">
default(clean)
# In CASA
mode                =      'mfs'         
default(tclean)
niter              =       2000       
savemodel          =      'modelcolumn'
gain                =       0.1       
specmode            =      'mfs'         
threshold          = '4.5e-5Jy'       
niter              =     2000       
psfmode            =   'clark'      
gain                =     0.1       
imagermode          =  'csclean'   
threshold          =     '4.5e-5Jy'       
cyclefactor   =         3      
deconvolver        =     'clark_exp'          
cyclespeedup  =        -1     
cyclefactor         =     3            
multiscale          =        []     
interactive        =      False       
interactive        =      False       
mask                =         []       
mask                =     []       
imsize              = [1024, 1024]       
imsize              =     [1024, 1024]       
cell                = ['0.75arcsec', '0.75arcsec']
cell                =     ['0.75arcsec', '0.75arcsec']
stokes              =       'I'         
stokes              =     'I'         
weighting          = 'natural'         
weighting          =     'natural'         
uvtaper            =      False      
uvtaper            =      []      
pbcor              =      False       
pbcor              =      False       
minpb              =       0.1    
pblimit            =     0.1  
usescratch         =      False      
selectdata         =      False                
async              =      False     
vis                =      'NGC1068.split.ms'
for source in sourcelist:
imagename           =     'NGC1068_img'
    splitfile = source + '.split.ms'
rmtables(imagename + '*')
    vis = splitfile
tclean()
    imagename = source + '_img'
    rmtables(imagename + '*') # get rid of earlier versions of the image
    clean()
</source>
</source>


Take a look at the image using {{viewer}}.
Take a look at the image using {{viewer}}. From the GUI, select the .image file you want to view, select ''Raster Image'', and then dismiss and close your way through to the display.


<source lang="python">
<source lang="python">
# In CASA
viewer()
viewer()
</source>
</source>


From the GUI, select the .image file you want to view, select ''Raster Image'', and then dismiss and close your way through to the display.
In the interest of time for this tutorial, we've chosen to only image NGC1068; however, should one be interested, all of the sources can be imaged by looping over a similar {{tclean}} command.
 
<source lang="python">
# In CASA
default(tclean)
savemodel          =      'modelcolumn'
specmode            =      'mfs'       
niter              =      2000     
gain                =      0.1     
threshold          =      '4.5e-5Jy'     
deconvolver        =      'clark_exp'           
cyclefactor        =      3             
interactive        =      False     
mask                =      []     
imsize              =      [1024, 1024]     
cell                =      ['0.75arcsec', '0.75arcsec']
stokes              =      'I'       
weighting          =      'natural'       
uvtaper            =      []     
pbcor              =      False     
pblimit            =      0.1
selectdata          =      False                 
for source in sourcelist:
    splitfile = source + '.split.ms'
    vis = splitfile
    imagename = source + '_img'
    rmtables(imagename + '*') # get ride of earlier versions of the image
    tclean()
</source>


==Self-calibration==
==Self-calibration==
There's not much more that can be done for the fainter sources, but the brighter sources will almost certainly show some residual calibration artifacts that can be cleaned up with self-calibration. The strategy is to repeat using the primary calibration application {{gaincal}}, only this time we use the source itself as a calibrator.


At this point we have a set of calibrated visibilities (the "splits") and a set of images. There's not much more that can be done for the fainter sources, but the brighter sources will almost certainly show some residual calibration artifacts that can be cleaned up with self-calibration. The strategy is to repeat using the primary calibration application {{gaincal}}, only this time we use the source itself as a calibrator.
NGC 1068 is a good example of a difficult source. Firstly, it's very bright, so residual complex gain calibration errors should be obvious. Secondly, like the larger fraction of the sky, it's located painfully close to the celestial equator, meaning that the (''u'',''v'') coverage is particularly poor.
 
NGC 1068 is a good example of a difficult source. Firstly, it's very bright, so residual phase calibration errors should be obvious. Secondly, like the larger fraction of the sky, it's located painfully close to the celestial equator, meaning that the (''u'',''v'') coverage is particularly poor.  


To view the calibration artifacts, use {{viewer}} to display the image of NGC 1068, but adjust the display levels so that the displayed image saturates at 10 mJy. To do this, select the leftmost button bearing a wrench icon. The "Data Display Options" GUI opens, and you can enter the desired display stretch in the '''Data Range''' text box. To cutoff the stretch at 30 mJy,
To view the calibration artifacts, use {{viewer}} to display the image of NGC 1068, but adjust the display levels so that the displayed image saturates at 30 mJy. To do this, select the leftmost button bearing a wrench icon. The "Data Display Options" GUI opens, and you can enter the desired display stretch in the '''Data Range''' text box. To cutoff the stretch at 30 mJy,


* Data Range: [-0.005,0.03]
* Data Range: [-0.005,0.03]


You can see some residual systematics in the background and sidelobe artifacts that suggest the phase calibration isn't as good as it might be.
You can see some residual systematics in the background and sidelobe artifacts that suggest the complex gain calibration isn't as good as it might be.


The process of self-calibration will involve iterating between (''i'') improved {{clean}} models of the source (which are updated in the "split" visibility file so long as calready = True) and (''i''+1) improved calibrations generated by {{gaincal}}.  For simplicity, in this example we will repeatedly generate new calibration tables so that we can avoid using ''accum.''
The process of self-calibration will involve iterating between (''i'') improved {{tclean}} models of the source <!--(which are updated in the "split" visibility file so long as calready = True)--> and (''i''+1) improved calibrations generated by {{gaincal}}.  For simplicity, in this example we will repeatedly generate new calibration tables so that we can avoid using {{accum}}.


=== Step ''i'' ===
=== Step ''i'' ===


In principle, we could start with the initial {{clean}} model that we generated above. For this example, we'll take a more cautious approach and re-run {{clean}} interactively.
In principle, we could start with the initial {{tclean}} model that we generated above. For this example, we'll take a more cautious approach and re-run {{tclean}} interactively.


<source lang="python">
<source lang="python">
# In CASA
source = 'NGC1068' # set up a global variable to store the source we're working on
source = 'NGC1068' # set up a global variable to store the source we're working on
tget clean
tget tclean
vis = source + '.split.ms'
vis = source + '.split.ms'
imagename = source + '_img0'
imagename = source + '_img0'
calready = True # make sure we store the model in the measurement set
interactive = True
interactive = True
cycleniter = 10
go
go
</source>
</source>
Line 405: Line 421:


<source lang="python">
<source lang="python">
default('gaincal') # reset inputs for gaincal (akin to RESTORE in AIPS)
# In CASA
default('gaincal') # reset inputs for gaincal
vis = source + '.split.ms'
vis = source + '.split.ms'
caltable = source + '.gcal0'
caltable = source + '.gcal0'
Line 415: Line 432:
field = source
field = source
gainfield = source # are these necessary?
gainfield = source # are these necessary?
gaincurve = False
gaincal() # or go
gaincal() # or go!
</source>
</source>


=== Pause for Reflection and Troubleshooting ===
=== Pause for Reflection and Troubleshooting ===


Self-calibration is a method to fine-tune your original calibration, and, as such, should produce relatively small phase corrections assuming the primary calibration was applied correctly. Now is a good time to run {{plotcal}} on your new calibration table (''NGC1068.gcal0'') to make sure the phase corrections range from a few to perhaps 10 or 20 degrees. If you see phase corrections greatly exceeding 20 degrees, then there's a chance of either
Self-calibration is a method to fine-tune your original calibration, and, as such, should produce relatively small complex gain corrections assuming the primary calibration was applied correctly. Now is a good time to run {{plotms}} on your new calibration table (''NGC1068.gcal0'') to make sure the phase corrections range from a few to perhaps 10 or 20 degrees. If you see phase corrections greatly exceeding 20 degrees, then there's a chance the primary calibration wasn't applied to the split database correctly (you could check this by searching your CASA Log file for errors during the applycal tasks). To check the phase corrections, run the following and step through each antenna: 
 
* the primary calibration wasn't applied to the split database correctly, and should be checked, or
* the model wasn't updated in the last iteration of {{clean}}, meaning you should verify that ''calready = True.''


<!--* the primary calibration wasn't applied to the split database correctly, and should be checked, or -->
<!--* the model wasn't updated in the last iteration of {{clean}}, meaning you should verify that ''calready = True.''-->
<source lang="python">
<source lang="python">
plotcal(caltable, iteration='antenna', xaxis='time', yaxis='phase')
# In CASA
plotms(caltable, iteraxis='antenna', xaxis='time', yaxis='phase', plotrange=[-1,-1,-180,180])
</source>
</source>


=== Step ''i'' + 2 ===
=== Step ''i'' + 2 ===


Assuming {{plotcal}} gave reasonable results, next apply the calibration.
Assuming {{plotms}} gave reasonable results, next apply the calibration.


<source lang="python">
<source lang="python">
# In CASA
default('applycal')
default('applycal')
vis = source + '.split.ms'
vis = source + '.split.ms'
Line 442: Line 459:
</source>
</source>


And now we can generate a new {{clean}} model.
And now we can generate a new {{tclean}} model.


<source lang="python">
<source lang="python">
tget clean
# In CASA
tget tclean
imagename=source + '_img1' # choose an image root name that matches the future cal table to facilitate bookkeeping
imagename=source + '_img1' # choose an image root name that matches the future cal table to facilitate bookkeeping
go
go
Line 454: Line 472:
Assuming things are going well, we can now continue to repeat the process until we are satisfied with the image.  
Assuming things are going well, we can now continue to repeat the process until we are satisfied with the image.  


[[Image:N1068final.png|thumb|Final, self-calibrated image of NGC1068.]]
[[Image:PreEVLA_Seyfert_Final.png|thumb|Final, self-calibrated image of NGC1068.]]


<source lang="python">
<source lang="python">
# In CASA
tget gaincal
tget gaincal
caltable = source +'.gcal1'
caltable = source +'.gcal1'
Line 462: Line 481:
go
go


# good idea to plotcal() again here!
# good idea to plotms() again here!


tget applycal
tget applycal
Line 468: Line 487:
go
go


tget clean
tget tclean
imagename = source + '_img2'
imagename = source + '_img2'
go
go
Line 482: Line 501:


<source lang="python">
<source lang="python">
default('gaincal') # reset inputs for gaincal (akin to RESTORE in AIPS)
# In CASA
default('gaincal') # reset inputs for gaincal
vis = source + '.split.ms'
vis = source + '.split.ms'
caltable = source + '.gcalap'
caltable = source + '.gcalap'
Line 490: Line 510:
combine = ''
combine = ''
minsnr = 3.0
minsnr = 3.0
gaintable = source + '.gcal1' # apply your best phase calibration on the fly.  
gaintable = source + '.gcal1' # apply your best phase calibration on the fly.
gaincal() # or go!
gaincal() # or go
</source>
</source>


Here, we are applying the best phase calibration solution on the fly by setting the ''gaintable'' parameter, so that this A&P self calibration iteration will only make small tweaks on top of the phase solutions we just worked so hard to find. Check your solutions in {{plotcal}}---amplitudes should not vary much and should be close to unity. Phase tweaks should be minimal and close to zero. If the solutions look reasonable, go ahead and apply the calibration:
Here, we are applying the best phase calibration solution on the fly by setting the ''gaintable'' parameter, so that this 'ap' self calibration iteration will only make small tweaks on top of the phase solutions we just worked so hard to find. Check your solutions in {{plotms}}---amplitudes should not vary much and should be close to unity. Complex gain tweaks should be minimal and close to zero. If the solutions look reasonable, go ahead and apply the calibration:


<source lang="python">
<source lang="python">
# In CASA
default('applycal')
default('applycal')
vis = source + '.split.ms'
vis = source + '.split.ms'
Line 504: Line 525:
</source>
</source>


The critical difference here from phase-only self calibration iterations is that we are applying ''two'' gain tables---the best phase solutions (source + '.gcal1') and the A&P solutions (source + '.gcalap'). Now you should be ready for one final imaging run with {{clean}}!  For consistency, we will run it with the same command we used the first time:
The critical difference here from phase-only self calibration iterations is that we are applying ''two'' gain tables---the best phase solutions (source + '.gcal1') and the A&P solutions (source + '.gcalap'). Now you should be ready for one final imaging run with {{tclean}}!  For consistency, we will run it with the same command we used the first time:


<source lang="python">
<source lang="python">
default(clean)
# In CASA
mode                =      'mfs'         
default(tclean)
niter              =       2000       
savemodel          =      'modelcolumn'
gain                =       0.1       
specmode            =      'mfs'         
threshold          = '4.5e-5Jy'       
niter              =     2000       
psfmode            =   'clark'      
gain                =     0.1       
imagermode          =  'csclean'   
threshold          =     '4.5e-5Jy'       
cyclefactor   =         3      
deconvolver        =     'clark_exp'          
cyclespeedup  =        -1     
cyclefactor         =     3                  
multiscale          =        []     
interactive        =      False       
interactive        =      False       
mask                =         []       
mask                =     []       
imsize              = [1024, 1024]       
imsize              =     [1024, 1024]       
cell                = ['0.75arcsec', '0.75arcsec']
cell                =     ['0.75arcsec', '0.75arcsec']
stokes              =       'I'         
stokes              =     'I'         
weighting          = 'natural'         
weighting          =     'natural'         
uvtaper            =      False      
uvtaper            =      []      
pbcor              =      False       
pbcor              =      False       
minpb              =       0.1       
pblimit            =     0.1       
calready            =      True     
selectdata          =      False  
async              =      False    
vis                 =     'NGC1068.split.ms'
vis = 'NGC1068.split.ms'
imagename           =     'NGC1068_img4'
imagename = 'NGC1068_img4'
rmtables(imagename + '*') # get rid of earlier versions of the image
rmtables(imagename + '*') # get rid of earlier versions of the image
clean()
tclean()
</source>
</source>


You can use the {{viewer}} to inspect your final image:


{{Checked 3.4.0}}
<source lang="python">
# In CASA
viewer('NGC1068_img4.image')
</source>


[[VLA Tutorials | &#8629; '''VLA Tutorials''']]
{{Checked 5.7.2}}


Updated: [[User:Mkrauss|Miriam Krauss]], 15 June 2012
[[Pre-upgrade VLA Tutorials | &#8629; '''Pre-upgrade VLA Tutorials''']]

Latest revision as of 17:16, 2 February 2023


Last checked on CASA Version 5.7.2.


Disclaimer: Due to continuous CASA software updates, GUI images may look different on more recent versions than those shown here.


This particular tutorial demonstrates how to calibrate a VLA continuum survey, consisting of a flux denisty scale calibrator, target sources, and their respective complex gain calibrators. What distinguishes this example is that it includes multiple target sources and multiple complex gain calibrators.


Importing VLA Archive Data

Download your data from the VLA Archive; in this example we'll use the publicly available survey AG733, which consists of VLA C-configuration, C-band (5 GHz) continuum observations of a handful of Seyfert galaxies.

With some sensible fiddling of the web interface, you should be able to download two archive formatted files, one from 06-Dec-09 and one from 06-Dec-13, with names like AG733_A061209.xp1 and AG733_B061213.xp1 or 2006-12/vla2006-12-09.dat and 2006-12/vla2006-12-13.dat.

Example snippet of listobs output for program AG733.

It's helpful to make a subdirectory to store the VLA archive data and the CASA measurement sets before starting CASA. We call this subdirectory AG733 and place the downloaded files in it. Linux-level commands are given assuming you are using bash as your default shell.

# In bash
cd AG733
casa

To load the data in CASA, we employ the importvla command.

Graphics window produced by plotants. Shown are the antennas of the VLA for the program AG733.
# In CASA
# Before importing, get rid of any original version of the measurement set
if os.path.exists('ag733.ms'): rmtables('ag733.ms')
# let glob find the archive files, store in a python global "myFiles"
from glob import glob as filesearch 
myFiles = filesearch("AG733*.xp?")
importvla(archivefiles=myFiles,vis='ag733.ms')
listobs('ag733.ms')
vishead('ag733.ms')

The output measurement set, i.e., the data to be worked on, is ag733.ms. Note that the output of listobs and vishead go to the CASA message log window.

listobs provides a list of scans, sources and their respective indices, and antennas. You can also get a plot of antennas using plotants. The antennas plot aids the selection of a reference antenna; the antenna VA05 looks like a reasonable choice.

# In CASA
vis = 'ag733.ms'
plotants()

You may want to close plotants before moving on in the script to avoid locked tables.

 Note: Some older VLA data may give errors in more recent CASA versions when attempting to generate a plot of antennas. If this occurs, use the listobs task and choose an antenna that was located on a pad close to the center of the array (Pads numbers increase with distance from the center of the array, ie. "VLA:_N1" and "EVLA:W2" are pads close to the center of the array, where as, "VLA:_N18" is not).


Editing the Measurement Set

This technique is detailed in the tutorial Data flagging with plotms. Here's a shorter version.

# In CASA
plotms()

Plotms brings up a graphical display that allows you to plot various aspects of your data.

The plotms GUI with proper field and timerange options selected.

In this example, the measurement set was selected by navigating to and highlighting the measurement set ag733.ms (which will look like a directory to your file browser, because measurement sets are actually directory structures). A specific subset of data (the flux density scale calibrator 3C48 = 0137+331 on the first observing day) was selected by specifying in the appropriate "Data -> Selection" boxes:

  • field = 0137+331, and
  • timerange = 2006/12/09/0:0:0~23:0:0.

For this specific field, we see no obvious bad data. Note that the low points are actually cross-polarization data; you can see this by entering:

  • corr = RR, LL

For the sake of this example, we'll clip out the early part of the observations, say, all the points before 04:00:45:

Firstly, box the region to be flagged by using the "Mark Regions" button MarkRegionsButton.png. Next, go to the "Flag" tab at top, and select "Extend flags" and "Correlation", so that we will delete the RL and LR data as well (despite the fact that we're not plotting it). Then, press the "Go" or "Flag" button FlagThoseData.png to zap those data.

The viewer tool provides utilities to flag data in a manner similar to the AIPS task TVFLG; see the tutorial Data flagging with viewer.

If we leave the field and timerange filter boxes blank and plot again, we see a number of visibilities with extremely large amplitudes. If you have trouble seeing the points, adjust their size by going to "Display", setting "Unflagged Points Symbol" to "Custom", and selecting a larger pixel size for the "circle", "square", or "diamond" options and replotting. Highlight these few extremely high amplitude points and flag them as before. Now we'll reload the plot by selecting the "Reload" checkbox and clicking "Plot". When the plot reloads with its new adjusted scale, we notice that there are even more bad points. By highlighting and locating these new bad points (described in detail in the data flagging with plotms tutorial), we find most of the bad data can be attributed to one antenna, 'VA15', at specific times, and one time range where all antennas contribute. Close inspection will show there is also some bad data on VA05 and VA22 as well. You may need to use the Zoom Tool Casaplotms-zoom.png to see all of the bad points clearly. In the text below we flag these points using flagdata; you can also flag data interactively in the plotms display (though this is generally not encouraged since it does not save backup versions of the flags, and it can be easy to flag data by mistake or miss data that's not plotted).

# In CASA
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/09/04:45:10.0~04:45:20.0')
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/09/06:12:30.0~06:14:10.0')
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/13/04:32:40.0~04:35:30.0')
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/09/06:12:30.0~06:12:40.0')
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/09/06:13:00.0~06:13:40.0')
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/09/06:13:50.0~06:14:00.0')
flagdata(vis='ag733.ms', antenna='VA15', timerange='2006/12/09/06:14:00.0~06:14:10.0')
flagdata(vis='ag733.ms', antenna='VA05', timerange='2006/12/09/04:30:00.0~04:30:10.0')
flagdata(vis='ag733.ms', antenna='VA22', timerange='2006/12/09/06:28:10.0~06:28:20.0')
flagdata(vis='ag733.ms', timerange='2006/12/09/06:24:50.0~06:25:00.0')

Plot the data again to make sure the amplitudes look reasonable (if plotms is already open, just hit Shift+Plot to reload with the current flags):

# In CASA
plotms()

Careful with Cross-Pols

By default, plotms plots all of the data, including cross-polarization. All those low-amplitude points in your data are probably not junk. They are real cross-polarization data. Flagging those data will cause grief later on. Take a look at the tutorial data flagging with plotms.

Tread lightly! And, whenever possible, flag with flagdata and keep careful track of what you flag in a separate file.

A Python Interlude

CASA uses python as an interactive wrapper and scripting language. Python offers the ability to define global variables on the fly, and it's worth using such variables to take care of bookkeeping. First we'll define some list variables that store calibrator and source information.

# In CASA
FDS_callist = ['0137+331']
CG_callist = ['2250+143', '0119+321', '0237+288', '0239-025',
                '0323+055', '0339-017', '0423-013'] 
sourcelist = ['NGC7469', 'MRK0993', 'MRK1040', 'NGC1056', 'NGC1068',
              'NGC1194', 'NGC1241', 'NGC1320', 'F04385-0828',
              'NGC1667'] 
allcals = FDS_callist + CG_callist
print allcals

It will also be helpful later to have a way of referencing which complex gain calibrator should be assigned to which source. We can define the matching explicitly by using a dictionary variable.

# In CASA
calDict = {'NGC7469':'2250+143',
           'MRK0993':'0119+321', 
           'MRK1040':'0237+288', 
           'NGC1056':'0237+288',
           'NGC1068':'0239-025',
           'NGC1194':'0323+055',
           'NGC1241':'0323+055',
           'NGC1320':'0339-017',
           'F04385-0828':'0423-013',
           'NGC1667':'0423-013'}

Those variables being defined, we can proceed with the calibration (and not have to type source and calibrator names over and over).

Calibration

Set the Flux Denisty Scale

As in AIPS, CASA sets the flux density scale using an observation of a flux density scale calibrator (here, 0137+331) and the task setjy. For this example, we'll determine the calibration in reference to a model image of the calibrator source. Look at the inputs of setjy and get it ready to go!

 
# In CASA
# restore setjy back to its defaults
default setjy
# shows several of the current and possible input parameters for the current task
inp
vis = 'ag733.ms'
field = '0137+331'
# For the model image, need to specify the correct directory, something like the following
#
# load the CASA path from environment variable
import os
casapath = os.environ.get('CASAPATH').split()[0]
# generate full path to the image
model = casapath + '/data/nrao/VLA/CalModels/3C48_C.im'
go

That example was set up to demonstrate some AIPS heritage. You could equally well run setjy the following way:

 
# In CASA
import os
casapath = os.environ.get('CASAPATH').split()[0]
setjy(vis = 'ag733.ms', field = '0137+331', model = casapath + '/data/nrao/VLA/CalModels/3C48_C.im')

Determine calibration solutions

At this stage the data have an overall flux density scaling determined, but full gain solutions aren't there yet. The relevant task is gaincal, which will produce the table cal.G.

# In CASA
rmtables('cal.G') # get rid of any old, unmeritorious versions of the calibration
rmtables('gc.cal')# will give a warning if the tables don't already exist, which can be safely ignored

# First we'll generate an antenna zenith-angle dependent gain curve calibration table
gencal(vis='ag733.ms', caltype='gc', caltable='gc.cal') 

#
# Some caution on the following step. The solution interval is set to 'inf,' meaning to average over entire
# scans of each calibrator in turn. If a given scan is however longer than the coherence time (esp. on long baselines), 
# some decorrelation will occur and the flux density scale calibration will be thrown off. Suggest, e.g., solint='1min' or solint='2min' 
# as reasonable alternatives to solint='inf'. Trial and error may be your friend here.
#
gaincal(vis='ag733.ms',caltable='cal.G', field=','.join(allcals), solint='inf',
        refant='VA05', gaintable=['gc.cal'], append=False)
plotms(vis='cal.G', coloraxis='Antenna1', yaxis='amp', xaxis='time', timerange='2006/12/09/04:0:0~06:30:0')

There's some python trickery buried in this example. If you are unfamiliar with the string method join, try the following example. Recall that we had defined the variable allcals above in A Python Interlude.

# In CASA
print allcals 
print ','.join(allcals)

It's worth playing with the timerange, etc. in plotms to explore the solutions. You may notice that antenna 'EA26' does not appear in the solutions on the second day. For completeness, we'll flag it manually from the measurement set and then repeat the calibration.

# In CASA
flagdata(vis='ag733.ms',antenna='EA26', timerange='2006/12/13/0:0:0~24:0:0')
rmtables('cal.G') 
gaincal(vis='ag733.ms',caltable='cal.G', field=','.join(allcals), solint='inf',
        refant='VA05', gaintable=['gc.cal'], append=False )

Go back and re-plot the calibration; this example shows even more AIPS heritage.

# In CASA
tget plotms
plotms()

Using the task fluxscale, bootstrap the flux density scale onto the complex gain calibrators.

# In CASA
rmtables('cal.Gflx') # remove any old versions of the calibration table, if necessary
myFluxes = fluxscale(vis='ag733.ms', caltable='cal.G', reference=','.join(FDS_callist), 
                     transfer=','.join(CG_callist), fluxtable='cal.Gflx', append=False)

Here, the calibration table cal.G is modified and stored as 'cal.Gflx.' The Python dictionary "myFluxes" will also contain information about the flux density scaling and other useful tidbits (type "myFluxes" to see its contents). So far, solutions have been generated only for the calibrators, and they have not yet been transfered to the sources.

We could equally well have used reference='0137+337', but use of the join method allows for the more general possibility of multiple flux density scale calibrators.

Apply the Calibrations

The task in CASA is applycal. We loop over sources and calibrators to properly match them:

# In CASA
for source, calibrator in calDict.iteritems():
    applycal(vis='ag733.ms', field=source,
             gaintable='cal.Gflx', gainfield=calibrator)

Baseline-based corrections

Initial tclean image of NGC1068 without blcal.
Initial tclean image of NGC1068 with blcal applied.

Although it makes only modest improvements for the present example, observations including a mix of EVLA and VLA telescopes may benefit from performing baseline-based corrections on a strong calibrator. Now that we have initial antenna-based solutions, we can use blcal.

# In CASA
default('blcal')
vis = 'ag733.ms'
# output baseline-based calibration solutions
caltable = 'cal.BL' 
# use the strong flux density scale calibrator to determine baseline-based solutions
field = '0137+331' 
# generate a solution for each calibrator scan
solint = 'inf' 
gaintable = 'cal.Gflx'
# calibration table with the best antenna-based solutions, so far
gainfield = '0137+331' 
# calibrator for the BL calibrator; in this example, gain cal and bl cal are the same
interp = 'nearest'
# do the first day's data
timerange = '2006/12/09/0:0:0~24:0:0'
blcal()
# now the second day
timerange = '2006/12/13/0:0:0~24:0:0'
blcal()

This operation generates the calibration table cal.BL. Solutions can be inspected with plotcal (this functionality may transfer to plotms in the near future). Now we have to use applycal again to apply both the antenna-based and baseline-based calibration solutions to the data.

# In CASA
for source, calibrator in calDict.iteritems():
    applycal(vis='ag733.ms', field=source,
             gaintable=['cal.Gflx','cal.BL'], gainfield=calibrator)

Did the blcal operation make a difference here? Getting ahead of ourselves a little, the plots at right show the initial tclean images of NGC1068, one with blcal applied, and one without, and both images are displayed with the same stretch. There does appear a slight reduction of artifacts with blcal applied, but of course your mileage will vary.

Splitting the calibrated source data from the multisource measurement set

To split the calibrated source data out of the measurement set, use the CASA task split.

# In CASA
splitfile = 'NGC1667.split.ms'
rmtables(splitfile) # get rid of any old versions before splitting
split(vis='ag733.ms',outputvis=splitfile, datacolumn='corrected', field='NGC1667')

Instead of having to manually repeat the split task for each source, you can take advantage of scripting to loop over the sources:

# In CASA
for source in sourcelist:
    splitfile = source + '.split.ms'
    rmtables(splitfile)
    split(vis='ag733.ms',outputvis=splitfile, 
          datacolumn='corrected', field=source)

Imaging

An example desktop display of viewer, showing the image display panel, the data display options window, and the data loader.

Use tclean for imaging and clean deconvolution. Like IMAGR in AIPS, it's worth some time mulling over all of the options. A list of some of the various keywords is given in the example, below. Most of the options should be relatively familiar to the AIPS user.

For this snapshot survey, the example chooses a fairly conservative, slow clean to accommodate the less-than-ideal beam. You could probably get away with more efficient cleaning, but it's better to play it safe in scripted self-calibration.

# In CASA
default(tclean)
savemodel           =      'modelcolumn'
specmode            =      'mfs'        
niter               =      2000      
gain                =      0.1       
threshold           =      '4.5e-5Jy'      
deconvolver         =      'clark_exp'            
cyclefactor         =      3              
interactive         =      False      
mask                =      []       
imsize              =      [1024, 1024]      
cell                =      ['0.75arcsec', '0.75arcsec']
stokes              =      'I'        
weighting           =      'natural'        
uvtaper             =      []      
pbcor               =      False      
pblimit             =      0.1 
selectdata          =      False                  
vis                 =      'NGC1068.split.ms'
imagename           =      'NGC1068_img'
rmtables(imagename + '*')
tclean()

Take a look at the image using viewer. From the GUI, select the .image file you want to view, select Raster Image, and then dismiss and close your way through to the display.

# In CASA
viewer()

In the interest of time for this tutorial, we've chosen to only image NGC1068; however, should one be interested, all of the sources can be imaged by looping over a similar tclean command.

# In CASA
default(tclean)
savemodel           =      'modelcolumn'
specmode            =      'mfs'        
niter               =      2000      
gain                =      0.1       
threshold           =      '4.5e-5Jy'      
deconvolver         =      'clark_exp'            
cyclefactor         =      3              
interactive         =      False      
mask                =      []       
imsize              =      [1024, 1024]      
cell                =      ['0.75arcsec', '0.75arcsec']
stokes              =      'I'        
weighting           =      'natural'        
uvtaper             =      []      
pbcor               =      False      
pblimit             =      0.1 
selectdata          =      False                  
for source in sourcelist:
    splitfile = source + '.split.ms'
    vis = splitfile
    imagename = source + '_img'
    rmtables(imagename + '*') # get ride of earlier versions of the image
    tclean()

Self-calibration

There's not much more that can be done for the fainter sources, but the brighter sources will almost certainly show some residual calibration artifacts that can be cleaned up with self-calibration. The strategy is to repeat using the primary calibration application gaincal, only this time we use the source itself as a calibrator.

NGC 1068 is a good example of a difficult source. Firstly, it's very bright, so residual complex gain calibration errors should be obvious. Secondly, like the larger fraction of the sky, it's located painfully close to the celestial equator, meaning that the (u,v) coverage is particularly poor.

To view the calibration artifacts, use viewer to display the image of NGC 1068, but adjust the display levels so that the displayed image saturates at 30 mJy. To do this, select the leftmost button bearing a wrench icon. The "Data Display Options" GUI opens, and you can enter the desired display stretch in the Data Range text box. To cutoff the stretch at 30 mJy,

  • Data Range: [-0.005,0.03]

You can see some residual systematics in the background and sidelobe artifacts that suggest the complex gain calibration isn't as good as it might be.

The process of self-calibration will involve iterating between (i) improved tclean models of the source and (i+1) improved calibrations generated by gaincal. For simplicity, in this example we will repeatedly generate new calibration tables so that we can avoid using accum.

Step i

In principle, we could start with the initial tclean model that we generated above. For this example, we'll take a more cautious approach and re-run tclean interactively.

# In CASA
source = 'NGC1068' # set up a global variable to store the source we're working on
tget tclean
vis = source + '.split.ms'
imagename = source + '_img0'
interactive = True
cycleniter = 10
go

The trick to interactive cleaning is to draw clean boxes around the sources as they appear in the residual image. By default, clean boxes are set by right-clicking on the image and dragging. The box can be edited until you double right-click inside the box, which then sets that box for cleaning. Press the green icon to continue cleaning the image; on the next major cycle, a new residual map will be displayed, affording the option to add clean boxes as new sources appear where they were once swamped by sidelobes.

Of course, for self-calibration, it pays to stop the cleaning early when artifacts start to appear in the residuals!

Step i + 1

To start, we should already have a clean model for NGC 1068 stored in the visibility database. Now we can set up gaincal.

# In CASA
default('gaincal') # reset inputs for gaincal
vis = source + '.split.ms'
caltable = source + '.gcal0'
gaintype = 'G'
calmode = 'p' # phase-only self-calibration at this stage
solint = '10 min' # set up a fairly coarse correction interval in this early stage
minsnr = 3.0
gaintable = ''
field = source
gainfield = source # are these necessary?
gaincal() # or go

Pause for Reflection and Troubleshooting

Self-calibration is a method to fine-tune your original calibration, and, as such, should produce relatively small complex gain corrections assuming the primary calibration was applied correctly. Now is a good time to run plotms on your new calibration table (NGC1068.gcal0) to make sure the phase corrections range from a few to perhaps 10 or 20 degrees. If you see phase corrections greatly exceeding 20 degrees, then there's a chance the primary calibration wasn't applied to the split database correctly (you could check this by searching your CASA Log file for errors during the applycal tasks). To check the phase corrections, run the following and step through each antenna:

# In CASA
plotms(caltable, iteraxis='antenna', xaxis='time', yaxis='phase', plotrange=[-1,-1,-180,180])

Step i + 2

Assuming plotms gave reasonable results, next apply the calibration.

# In CASA
default('applycal')
vis = source + '.split.ms'
gaintable = source + '.gcal0'
gainfield = source
go

And now we can generate a new tclean model.

# In CASA
tget tclean
imagename=source + '_img1' # choose an image root name that matches the future cal table to facilitate bookkeeping
go

Further Steps

Assuming things are going well, we can now continue to repeat the process until we are satisfied with the image.

Final, self-calibrated image of NGC1068.
# In CASA
tget gaincal
caltable = source +'.gcal1'
solint = "5 min" # halve the solution interval. 
go

# good idea to plotms() again here!

tget applycal
gaintable = source + '.gcal1'
go

tget tclean
imagename = source + '_img2'
go

And so on, decreasing the solution interval as much as the S/N and sampling allow.

Note that gaincal always uses the 'data' column (not the recently self-calibrated 'corrected' column) to find gain solutions. Therefore, these iterations simply serve to improve the data model, and gain solutions at each iteration are independent (e.g., the solution from the second self calibration iteration does not accumulate on top of the solution from the first interval. Use the CASA task accum if you wish to accumulate gain solutions).

One Last Iteration: Amplitude & Phase Self Calibration

After you have improved the phase calibration as much as possible, one last iteration where both phase and amplitude are calibrated can be beneficial. A&P self calibration should be performed with long solution intervals---typically the length of a scan.

# In CASA
default('gaincal') # reset inputs for gaincal
vis = source + '.split.ms'
caltable = source + '.gcalap'
gaintype = 'G'
calmode = 'ap' # it's finally time to calibrate amplitudes, too
solint = 'inf' # one solution interval per scan (in conjunction with combine='') 
combine = ''
minsnr = 3.0
gaintable = source + '.gcal1' # apply your best phase calibration on the fly.
gaincal() # or go

Here, we are applying the best phase calibration solution on the fly by setting the gaintable parameter, so that this 'ap' self calibration iteration will only make small tweaks on top of the phase solutions we just worked so hard to find. Check your solutions in plotms---amplitudes should not vary much and should be close to unity. Complex gain tweaks should be minimal and close to zero. If the solutions look reasonable, go ahead and apply the calibration:

# In CASA
default('applycal')
vis = source + '.split.ms'
gaintable = [source + '.gcal1', source + '.gcalap']
gainfield = source
go

The critical difference here from phase-only self calibration iterations is that we are applying two gain tables---the best phase solutions (source + '.gcal1') and the A&P solutions (source + '.gcalap'). Now you should be ready for one final imaging run with tclean! For consistency, we will run it with the same command we used the first time:

# In CASA
default(tclean)
savemodel           =      'modelcolumn'
specmode            =      'mfs'        
niter               =      2000      
gain                =      0.1       
threshold           =      '4.5e-5Jy'      
deconvolver         =      'clark_exp'           
cyclefactor         =      3                    
interactive         =      False      
mask                =      []       
imsize              =      [1024, 1024]      
cell                =      ['0.75arcsec', '0.75arcsec']
stokes              =      'I'        
weighting           =      'natural'        
uvtaper             =      []      
pbcor               =      False      
pblimit             =      0.1      
selectdata          =      False    
vis                 =      'NGC1068.split.ms'
imagename           =      'NGC1068_img4'
rmtables(imagename + '*') # get rid of earlier versions of the image
tclean()

You can use the viewer to inspect your final image:

# In CASA
viewer('NGC1068_img4.image')

Last checked on CASA Version 5.7.2.

Pre-upgrade VLA Tutorials