VLA CASA Imaging-CASA4.5.2: Difference between revisions

From CASA Guides
Jump to navigationJump to search
(Created page with "* '''This CASA guide was designed for CASA 4.5.3''' https://casaguides.nrao.edu/index.php?title=VLA-CASA-Flagging-CASA4.5.3 https://casaguides.nrao.edu/index.php?title=VLA-C...")
 
No edit summary
 
(312 intermediate revisions by 5 users not shown)
Line 1: Line 1:
* '''This CASA guide was designed for CASA 4.5.3'''


https://casaguides.nrao.edu/index.php?title=VLA-CASA-Flagging-CASA4.5.3
== Imaging ==
https://casaguides.nrao.edu/index.php?title=VLA-CASA-Imaging-CASA4.5.3
== Overview ==
 
This CASA guide will cover priori data flagging, including online flagging, shadowing, zero clipping, and quacking. It will also cover auto-flagging RFI (Radio Frequency Interference) via TFcrop (Time-Frequency crop) and rflag.


The guide will then move on to imaging the data, and cover topics on the CLEAN algorithm, MS (Multi-Scale) deconvolution, MS-MFS (Multi-Scale, Multi-Frequency Synthesis), outlier fields, spectral indices, image weights and tapering, small scale bias, primary beam correction, image headers, and imaging in widefield mode, which will make use of w-projection.
This tutorial provides guidance on imaging procedures in CASA. The tutorial covers basic continuum cleaning and the influence of image weights, as well as wide-band and wide-field imaging techniques, multi-scale clean, an outlier field setup, and primary beam correction. Spectral line imaging procedures are explained but not covered in detail. For a more thorough example of spectral line imaging, refer to the [https://casaguides.nrao.edu/index.php/VLA_high_frequency_Spectral_Line_tutorial_-_IRC%2B10216 VLA high frequency Spectral Line tutorial - IRC+10216]. This imaging tutorial concludes with basic image header calls and the conversion of Jy/beam to K surface brightness units through image arithmetic and header manipulation.  


We will be utilizing data taken with the Karl G. Jansky, Very Large Array, of a supernova remnant [http://simbad.u-strasbg.fr/simbad/sim-id?Ident=SNR+G055.7%2B03.4&NbIdent=1&Radius=2&Radius.unit=arcmin&submit=submit+id G055.7+3.4.]. The data were taken on August 23, 2010, in the first D-configuration for which the new wide-band capabilities of the WIDAR (Wideband Interferometric Digital ARchitecture) correlator were available. The 8-hour-long observation includes all available 1 GHz of bandwidth in L-band, from 1-2 GHz in frequency. The observations were made in 8-bit mode, which has the best dynamic range and is thus best suited for rfi excision.
We will be utilizing data, taken with the Karl G. Jansky Very Large Array, of G055.7+3.4, which is a supernova remnant. The data were taken on August 23, 2010, in the first D-configuration for which the wideband capabilities of the WIDAR (Wideband Interferometric Digital ARchitecture) correlator were available. The 8 hour long observation includes all available 1 GHz of bandwidth in L-band, from 1–2 GHz in frequency.  


The guide will often reference the CASA cookbook which is available [http://casa.nrao.edu/docs/cookbook/index.html online] and as a [http://casa.nrao.edu/casa_cookbook.pdf pdf download].  
We will skip the calibration process as examples of calibration can be found in several other guides, including the
[https://casaguides.nrao.edu/index.php?title=VLA_Continuum_Tutorial_3C391 VLA Continuum Tutorial 3C391] and  
[https://casaguides.nrao.edu/index.php?title=VLA_high_frequency_Spectral_Line_tutorial_-_IRC%2B10216 VLA high frequency Spectral Line tutorial - IRC+10216] tutorials.


== Obtaining the Data ==
A copy of the calibrated data <font color=green>(1.2GB)</font> can be downloaded from [http://casa.nrao.edu/Data/EVLA/SNRG55/SNR_G55_10s.calib.tar.gz http://casa.nrao.edu/Data/EVLA/SNRG55/SNR_G55_10s.calib.tar.gz]


A copy of the data <font color=green>(5.1GB)</font> can be downloaded from [http://casa.nrao.edu/Data/EVLA/SNRG55/SNR_G55_10s.tar.gz http://casa.nrao.edu/Data/EVLA/SNRG55/SNR_G55_10s.tar.gz]
Your first step will be to unzip and untar the file in a terminal window (before you start CASA):
These data are already in CASA format and were prepared specifically for this guide.


Once you've downloaded the file, let's untar it:
<source lang="bash">
<source lang="bash">
tar -xzvf SNR_G55_10s.tar.gz
tar -xzvf SNR_G55_10s.calib.tar.gz
</source>
</source>


That will create a directory called SNR_G55_10s.ms <font color=green>(6.1GB)</font>, which is the MS (Measurement Set).
Then start casa as usual via the '''casa''' command, which will bring up the ipython interface and launches the logger.
 


== Starting CASA ==
== The CLEAN Algorithm ==


Start CASA by typing '''casa''' on a terminal command line.  If you have not used CASA before, some helpful tips are available on the [[Getting Started in CASA]] page.
[[Image:CLEAN_Cycle.png|500px|thumb|right|'''Figure 1''' <br /> The CLEAN major and minor cycles, indicating the steps undertaken during gridding, projection algorithms, and creation of images.]]


This guide has been written for CASA version 4.5.2. Please confirm your version before proceeding by checking the message in the command line interface window or the CASA logger after startup.
The CLEAN algorithm, developed by J. Högbom (1974) enabled the synthesis of complex objects, even if they have relatively poor Fourier uv-plane coverage. Poor coverage occurs with partial earth rotation synthesis or with arrays composed of few antennas. The dirty image is formed by a simple Fourier inversion of the sampled visibility data, with each point on the sky being represented by a suitably scaled and centered PSF (Point Spread Function), or dirty beam, which itself is the Fourier inversion of the visibility (u,v) coverage.  


For this tutorial, we will be running tasks using the ''task (parameter = value)'' syntax. When called in this manner, all parameters not explicitly set will use their default values.  
The convolution with the dirty beam creates artifacts in the image and limits the dynamic range. The CLEAN algorithm attempts to remove the dirty beam pattern from the image via deconvolution. This implies that it interpolates from the measured (u,v) points across gaps in the (u,v) coverage. In short, CLEAN provides solutions to the convolution equation by representing radio sources by a number of point sources in an empty field. The brightest points are found by performing a cross-correlation between the dirty image and the PSF. The brightest parts are then subtracted, and the process is repeated again for the next brighter sources. Variants of CLEAN, such as multi-scale CLEAN, take into account extended kernels which may be better suited for extended objects.


== Preliminary Data Evaluation ==
For single pointings, CASA uses the Cotton-Schwab cleaning algorithm in the task {{clean}} (parameter ''imagermode='csclean'''), which breaks the process into major and minor cycles (see Figure 1). To start with, the visibilities are gridded, weighted, and Fourier transformed to create a dirty image. The minor cycle then operates in the image domain to find the clean components that are added to the clean model. The model image is Fourier transformed back to the visibility domain, degridded, and subtracted from the visibilities. This creates a new residual that is then gridded, weighted, and FFT'ed again to the image domain for the next iteration. The gridding, FFT, degridding, and subtraction processes form the major cycle.


As a first step, use {{listobs}} to have a look at the MS:
This iterative process is continued until a stopping criterion is reached, such as a maximum number of clean components or a flux threshold in the residual image.


<source lang="python">
In CASA {{clean}}, two versions of the PSF can be used: parameter ''psfmode = 'hogbom' ''uses the full-sized PSF for subtraction and is a thorough but slow method. Parameter ''psfmode='clark' ''uses a smaller beam patch, which increases the speed. The patch size and length of the minor cycle are internally chosen such that clean converges well without giving up the speed improvement. Parameter ''psfmode='clark' ''is the default option in {{clean}}.  
# In CASA
listobs(vis='SNR_G55_10s.ms', listfile='SNR_G55_10s.listobs')
</source>


*listfile: Creates a text document with details for the observation.
<!--
<pre>
''' 1. Högbom Algorithm '''<br />
#CASA log when listfile is not defined
This algorithm will initially find the strength and position of a peak in a dirty image, subtract it from the dirty image, record this position and maginitude, and repeat for further peaks. The remainder of the dirty image is known as the residuals.  
2016-03-07 16:54:23 INFO listobs  ##########################################
2016-03-07 16:54:23 INFO listobs  ##### Begin Task: listobs            #####
2016-03-07 16:54:23 INFO listobs  listobs(vis="SNR_G55_10s.ms",selectdata=True,spw="",field="",antenna="",
2016-03-07 16:54:23 INFO listobs   uvrange="",timerange="",correlation="",scan="",intent="",
2016-03-07 16:54:23 INFO listobs   feed="",array="",observation="",verbose=True,listfile="",
2016-03-07 16:54:23 INFO listobs   listunfl=False,cachesize=50,overwrite=False)
2016-03-07 16:54:23 INFO listobs  ================================================================================
2016-03-07 16:54:23 INFO listobs     MeasurementSet Name:  /lustre/aoc/sciops/CASA_Tutorials/SNR_G55_10s.ms      MS Version 2
2016-03-07 16:54:23 INFO listobs  ================================================================================
2016-03-07 16:54:23 INFO listobs    Observer: Dr. Sanjay Sanjay Bhatnagar    Project: uid://evla/pdb/1072564 
2016-03-07 16:54:23 INFO listobs  Observation: EVLA
2016-03-07 16:54:25 INFO listobs  Data records: 3780000      Total elapsed time = 26926 seconds
2016-03-07 16:54:25 INFO listobs    Observed from  23-Aug-2010/00:56:36.0  to  23-Aug-2010/08:25:22.0 (UTC)
2016-03-07 16:54:33 INFO listobs  
2016-03-07 16:54:33 INFO listobs    ObservationID = 0        ArrayID = 0
2016-03-07 16:54:33 INFO listobs  Date        Timerange (UTC)          Scan  FldId FieldName    nRows    SpwIds      Average Interval(s)    ScanIntent
2016-03-07 16:54:33 INFO listobs  23-Aug-2010/00:56:36.0 - 00:58:06.0    14      0 J1925+2106    12600  [0,1,2,3]  [10, 10, 10, 10, 10] [CALIBRATE_PHASE#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       00:58:06.0 - 00:59:36.0    15      0 J1925+2106    12600  [0,1,2,3]  [10, 10, 10, 10, 10] [CALIBRATE_PHASE#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       00:59:36.0 - 01:01:05.0    16      0 J1925+2106    12600  [0,1,2,3]  [9.89, 9.89, 9.89, 9.89, 9.89] [CALIBRATE_PHASE#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:01:05.0 - 01:02:35.0    17      0 J1925+2106    12600  [0,1,2,3]  [10, 10, 10, 10, 10] [CALIBRATE_PHASE#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:02:35.0 - 01:04:05.0    18      0 J1925+2106    12600  [0,1,2,3]  [10, 10, 10, 10, 10] [CALIBRATE_PHASE#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:04:05.0 - 01:05:34.0    19      0 J1925+2106    12600  [0,1,2,3]  [9.89, 9.89, 9.89, 9.89, 9.89] [CALIBRATE_PHASE#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:05:34.0 - 01:07:04.0    20      0 J1925+2106    12600  [0,1,2,3]  [10, 10, 10, 10, 10] [CALIBRATE_PHASE#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:07:04.0 - 01:08:34.0    21      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:08:34.0 - 01:10:04.0    22      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:10:04.0 - 01:11:34.0    23      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:11:34.0 - 01:13:03.0    24      1 G55.7+3.4      12600  [0,1,2,3]  [9.89, 9.89, 9.89, 9.89, 9.89] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:13:03.0 - 01:14:33.0    25      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:14:33.0 - 01:16:03.0    26      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:16:03.0 - 01:17:33.0    27      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:17:33.0 - 01:19:02.0    28      1 G55.7+3.4      12600  [0,1,2,3]  [9.89, 9.89, 9.89, 9.89, 9.89] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:19:02.0 - 01:20:32.0    29      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:20:32.0 - 01:22:02.0    30      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:22:02.0 - 01:23:32.0    31      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:23:32.0 - 01:25:01.0    32      1 G55.7+3.4      12600  [0,1,2,3]  [9.89, 9.89, 9.89, 9.89, 9.89] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:25:01.0 - 01:26:31.0    33      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:26:31.0 - 01:28:01.0    34      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:28:01.0 - 01:29:31.0    35      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:29:31.0 - 01:31:00.0    36      1 G55.7+3.4      12600  [0,1,2,3]  [9.89, 9.89, 9.89, 9.89, 9.89] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:31:00.0 - 01:32:30.0    37      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:32:30.0 - 01:34:00.0    38      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:34:00.0 - 01:35:30.0    39      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:35:30.0 - 01:36:59.0    40      1 G55.7+3.4      12600  [0,1,2,3]  [9.89, 9.89, 9.89, 9.89, 9.89] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:36:59.0 - 01:38:29.0    41      0 J1925+2106    12600  [0,1,2,3]  [10, 10, 10, 10, 10] [CALIBRATE_PHASE#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:38:29.0 - 01:39:59.0    42      0 J1925+2106    12600  [0,1,2,3]  [10, 10, 10, 10, 10] [CALIBRATE_PHASE#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:39:59.0 - 01:41:29.0    43      0 J1925+2106    12600  [0,1,2,3]  [10, 10, 10, 10, 10] [CALIBRATE_PHASE#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       01:41:29.0 - 01:42:58.0    44      1 G55.7+3.4      12600  [0,1,2,3]  [9.89, 9.89, 9.89, 9.89, 9.89] [OBSERVE_TARGET#UNSPECIFIED]
<snip>
2016-03-07 16:54:33 INFO listobs       08:11:54.0 - 08:13:24.0  305      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       08:13:24.0 - 08:14:54.0  306      1 G55.7+3.4      12600  [0,1,2,3]  [10, 10, 10, 10, 10] [OBSERVE_TARGET#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       08:14:54.0 - 08:16:24.0  307      2 0542+498=3C147 12600  [0,1,2,3]  [10, 10, 10, 10, 10] [CALIBRATE_AMPLI#UNSPECIFIED,CALIBRATE_BANDPASS#UNSPECIFIED,UNSPECIFIED#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       08:16:24.0 - 08:17:54.0  308      2 0542+498=3C147 12600  [0,1,2,3]  [10, 10, 10, 10, 10] [CALIBRATE_AMPLI#UNSPECIFIED,CALIBRATE_BANDPASS#UNSPECIFIED,UNSPECIFIED#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       08:17:54.0 - 08:19:23.0  309      2 0542+498=3C147 12600  [0,1,2,3]  [9.89, 9.89, 9.89, 9.89, 9.89] [CALIBRATE_AMPLI#UNSPECIFIED,CALIBRATE_BANDPASS#UNSPECIFIED,UNSPECIFIED#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       08:19:23.0 - 08:20:53.0  310      2 0542+498=3C147 12600  [0,1,2,3]  [10, 10, 10, 10, 10] [CALIBRATE_AMPLI#UNSPECIFIED,CALIBRATE_BANDPASS#UNSPECIFIED,UNSPECIFIED#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       08:20:53.0 - 08:22:23.0  311      2 0542+498=3C147 12600  [0,1,2,3]  [10, 10, 10, 10, 10] [CALIBRATE_AMPLI#UNSPECIFIED,CALIBRATE_BANDPASS#UNSPECIFIED,UNSPECIFIED#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       08:22:23.0 - 08:23:52.0  312      2 0542+498=3C147 12600  [0,1,2,3]  [9.89, 9.89, 9.89, 9.89, 9.89] [CALIBRATE_AMPLI#UNSPECIFIED,CALIBRATE_BANDPASS#UNSPECIFIED,UNSPECIFIED#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs       08:23:52.0 - 08:25:22.0  313      2 0542+498=3C147 12600  [0,1,2,3]  [10, 10, 10, 10, 10] [CALIBRATE_AMPLI#UNSPECIFIED,CALIBRATE_BANDPASS#UNSPECIFIED,UNSPECIFIED#UNSPECIFIED]
2016-03-07 16:54:33 INFO listobs           (nRows = Total number of rows per scan)
2016-03-07 16:54:33 INFO listobs  Fields: 3
2016-03-07 16:54:33 INFO listobs  ID  Code Name                RA              Decl          Epoch  SrcId      nRows
2016-03-07 16:54:33 INFO listobs  0    D    J1925+2106          19:25:59.605371 +21.06.26.16218 J2000  0        541800
2016-03-07 16:54:33 INFO listobs  1    NONE G55.7+3.4          19:21:40.000000 +21.45.00.00000 J2000  1        3150000
2016-03-07 16:54:33 INFO listobs  2    N    0542+498=3C147      05:42:36.137916 +49.51.07.23356 J2000  2          88200
2016-03-07 16:54:33 INFO listobs  Spectral Windows:  (5 unique spectral windows and 1 unique polarization setups)
2016-03-07 16:54:33 INFO listobs  SpwID  Name      #Chans  Frame  Ch0(MHz)  ChanWid(kHz)  TotBW(kHz) CtrFreq(MHz) BBC Num  Corrs         
2016-03-07 16:54:33 INFO listobs  0      Subband:0    64  TOPO    1256.000      2000.000    128000.0  1319.0000        4  RR  RL  LR  LL
2016-03-07 16:54:33 INFO listobs  1      Subband:2    64  TOPO    1384.000      2000.000    128000.0  1447.0000        4  RR  RL  LR  LL
2016-03-07 16:54:33 INFO listobs  2      Subband:1    64  TOPO    1648.000      2000.000    128000.0  1711.0000        8  RR  RL  LR  LL
2016-03-07 16:54:33 INFO listobs  3      Subband:0    64  TOPO    1776.000      2000.000    128000.0  1839.0000        8  RR  RL  LR  LL
2016-03-07 16:54:33 INFO listobs  Sources: 15
2016-03-07 16:54:33 INFO listobs  ID  Name                SpwId RestFreq(MHz)  SysVel(km/s)
2016-03-07 16:54:33 INFO listobs  0    J1925+2106          0    -              -           
2016-03-07 16:54:33 INFO listobs  0    J1925+2106          1    -              -           
2016-03-07 16:54:33 INFO listobs  0    J1925+2106          2    -              -           
2016-03-07 16:54:33 INFO listobs  0    J1925+2106          3    -              -           
2016-03-07 16:54:33 INFO listobs  1    G55.7+3.4          0    -              -           
2016-03-07 16:54:33 INFO listobs  1    G55.7+3.4          1    -              -           
2016-03-07 16:54:33 INFO listobs  1    G55.7+3.4          2    -              -           
2016-03-07 16:54:33 INFO listobs  1    G55.7+3.4          3    -              -           
2016-03-07 16:54:33 INFO listobs  2    0542+498=3C147      0    -              -           
2016-03-07 16:54:33 INFO listobs  2    0542+498=3C147      1    -              -           
2016-03-07 16:54:33 INFO listobs  2    0542+498=3C147      2    -              -           
2016-03-07 16:54:33 INFO listobs  2    0542+498=3C147      3    -              -           
2016-03-07 16:54:34 INFO listobs  Antennas: 27:
2016-03-07 16:54:34 INFO listobs  ID  Name  Station  Diam.    Long.        Lat.                Offset from array center (m)                ITRF Geocentric coordinates (m)       
2016-03-07 16:54:34 INFO listobs                                                                East      North    Elevation        x              y              z
2016-03-07 16:54:34 INFO listobs  0    ea01  W09      25.0 m  -107.37.25.2  +33.53.51.0  -521.9416  -332.7766    -1.2001 -1601710.017000 -5042006.925200  3554602.355600
2016-03-07 16:54:34 INFO listobs  1    ea02  E02      25.0 m  -107.37.04.4  +33.54.01.1      9.8240    -20.4293    -2.7806 -1601150.060300 -5042000.619800  3554860.729400
2016-03-07 16:54:34 INFO listobs  2    ea03  E09      25.0 m  -107.36.45.1  +33.53.53.6    506.056    -251.8670    -3.5825 -1600715.950800 -5042273.187000  3554668.184500
2016-03-07 16:54:34 INFO listobs  3    ea04  W01      25.0 m  -107.37.05.9  +33.54.00.5    -27.3562    -41.3030    -2.7418 -1601189.030140 -5042000.493300  3554843.425700
2016-03-07 16:54:34 INFO listobs  4    ea05  W08      25.0 m  -107.37.21.6  +33.53.53.0  -432.1167  -272.1478    -1.5054 -1601614.091000 -5042001.652900  3554652.509300
2016-03-07 16:54:34 INFO listobs  5    ea06  N06      25.0 m  -107.37.06.9  +33.54.10.3    -54.0649    263.8778    -4.2273 -1601162.591000 -5041828.999000  3555095.896400
2016-03-07 16:54:34 INFO listobs  6    ea07  E05      25.0 m  -107.36.58.4  +33.53.58.8    164.9788    -92.8032    -2.5268 -1601014.462000 -5042086.252000  3554800.799800
2016-03-07 16:54:34 INFO listobs  7    ea08  N01      25.0 m  -107.37.06.0  +33.54.01.8    -30.8810    -1.4664    -2.8597 -1601185.634945 -5041978.156586  3554876.424700
2016-03-07 16:54:34 INFO listobs  8    ea09  E06      25.0 m  -107.36.55.6  +33.53.57.7    236.9058  -126.3369    -2.4443 -1600951.588000 -5042125.911000  3554773.012300
2016-03-07 16:54:34 INFO listobs  9    ea10  N03      25.0 m  -107.37.06.3  +33.54.04.8    -39.0773    93.0192    -3.3330 -1601177.376760 -5041925.073200  3554954.584100
2016-03-07 16:54:34 INFO listobs  10  ea11  E04      25.0 m  -107.37.00.8  +33.53.59.7    102.8054    -63.7682    -2.6414 -1601068.790300 -5042051.910200  3554824.835300
2016-03-07 16:54:34 INFO listobs  11  ea12  E08      25.0 m  -107.36.48.9  +33.53.55.1    407.8285  -206.0065    -3.2272 -1600801.926000 -5042219.366500  3554706.448200
2016-03-07 16:54:34 INFO listobs  12  ea13  N07      25.0 m  -107.37.07.2  +33.54.12.9    -61.1037    344.2331    -4.6138 -1601155.635800 -5041783.843800  3555162.374100
2016-03-07 16:54:34 INFO listobs  13  ea15  W06      25.0 m  -107.37.15.6  +33.53.56.4  -275.8288  -166.7451    -2.0590 -1601447.198000 -5041992.502500  3554739.687600
2016-03-07 16:54:34 INFO listobs  14  ea16  W02      25.0 m  -107.37.07.5  +33.54.00.9    -67.9687    -26.5614    -2.7175 -1601225.255200 -5041980.383590  3554855.675000
2016-03-07 16:54:34 INFO listobs  15  ea17  W07      25.0 m  -107.37.18.4  +33.53.54.8  -349.9877  -216.7509    -1.7975 -1601526.387300 -5041996.840100  3554698.327400
2016-03-07 16:54:34 INFO listobs  16  ea18  N09      25.0 m  -107.37.07.8  +33.54.19.0    -77.4346    530.6273    -5.5859 -1601139.485100 -5041679.036800  3555316.533200
2016-03-07 16:54:34 INFO listobs  17  ea19  W04      25.0 m  -107.37.10.8  +33.53.59.1  -152.8599    -83.8054    -2.4614 -1601315.893000 -5041985.320170  3554808.304600
2016-03-07 16:54:34 INFO listobs  18  ea20  N05      25.0 m  -107.37.06.7  +33.54.08.0    -47.8454    192.6015    -3.8723 -1601168.786100 -5041869.054000  3555036.936000
2016-03-07 16:54:34 INFO listobs  19  ea21  E01      25.0 m  -107.37.05.7  +33.53.59.2    -23.8638    -81.1510    -2.5851 -1601192.467800 -5042022.856800  3554810.438800
2016-03-07 16:54:34 INFO listobs  20  ea22  N04      25.0 m  -107.37.06.5  +33.54.06.1    -42.6239    132.8436    -3.5494 -1601173.979400 -5041902.657700  3554987.517500
2016-03-07 16:54:34 INFO listobs  21  ea23  E07      25.0 m  -107.36.52.4  +33.53.56.5    318.0509  -164.1850    -2.6957 -1600880.571400 -5042170.388000  3554741.457400
2016-03-07 16:54:34 INFO listobs  22  ea24  W05      25.0 m  -107.37.13.0  +33.53.57.8  -210.0959  -122.3887    -2.2577 -1601377.009500 -5041988.665500  3554776.393400
2016-03-07 16:54:34 INFO listobs  23  ea25  N02      25.0 m  -107.37.06.2  +33.54.03.5    -35.6245    53.1806    -3.1345 -1601180.861480 -5041947.453400  3554921.628700
2016-03-07 16:54:34 INFO listobs  24  ea26  W03      25.0 m  -107.37.08.9  +33.54.00.1  -105.3447    -51.7177    -2.6037 -1601265.153600 -5041982.533050  3554834.858400
2016-03-07 16:54:34 INFO listobs  25  ea27  E03      25.0 m  -107.37.02.8  +33.54.00.5    50.6641    -39.4835    -2.7273 -1601114.365500 -5042023.151800  3554844.944000
2016-03-07 16:54:34 INFO listobs  26  ea28  N08      25.0 m  -107.37.07.5  +33.54.15.8    -68.9057    433.1889    -5.0602 -1601147.940400 -5041733.837000  3555235.956000
2016-03-07 16:54:34 INFO listobs ##### End Task: listobs              #####
2016-03-07 16:54:34 INFO listobs ##########################################
</pre>


* J1925+2106,     field ID 0: Phase Calibrator;
The accumulated point sources, now residing in a model, is convolved with an idealized CLEAN beam (usually a Gaussian fitted to the central lobe of the dirty beam), creating a CLEAN image. As the final step, the residuals of the dirty image are then added to the CLEAN image.  
* G55.7+3.4,      field ID 1: The Supernova Remnant;
* 0542+498=3C147, field ID 2: Amplitude/Bandpass Calibrator


We can also see that these sources have associated "scan intents", which indicate their function in the observation. Note that you can select sources based on their intents in certain CASA tasks. The various scan intents in this data set are:
''' 2. Clark Algorithm '''<br />
Clark (1980), developed a FFT-based CLEAN algorithm, which more efficiently shifts and scales the dirty beam by approximating the position and strength of components using a small patch of the dirty beam. This algorithm is the default within the {{clean}} task, which involves major and minor cycles.  


* CALIBRATE_PHASE indicates that this is a scan to be used for gain calibration;
The algorithm will first select a beam patch, which will include the highest exterior sidelobes. Points are then selected from the dirty image, which are up to a fraction of the image peak, and are greater than the highest exterior sidelobe of the beam. It will then conduct a list-based Högbom CLEAN, creating a model and convolution with an idealized CLEAN beam. This process is the minor cycle.  
* OBSERVE_TARGET indicates that this is the science target;
* CALIBRATE_AMPLI indicates that this is to be used for flux calibration; and  
* CALIBRATE_BANDPASS indicates that these scans are to be used for bandpass calibration.


Note that 3C147 is to be used for both flux and bandpass calibration.
The major cycle involves transforming the point source model via a FFT (Fast-Fourier Transform), mutiplying this by the weight sampling function (more on this below), and transformed back. This is then subtracted from the dirty image, creating your CLEAN image. The process is then repeated with subsequent minor cycles.


It's important to also note that the antennas have a name and ID associated with them. For example antenna ID 15 is named ea17 ( The "ea" stemming from the Expanded VLA project). When specifying an antenna within a task parameter, we will mainly reference them by name.  
''' 3. Cotton-Schwab Algorithm '''<br />
This is the default imager mode (''csclean''), and is a variant of the Clark algorith in which the major cycle involves the subtraction of CLEAN components of ungridded visibility data. This allows the removal of gridding errors, as well as noise. One advantage is its ability to image and clean many seperate fields simultaneously. Fields are cleaned independently in the minor cycle, and components from all fields cleaned together in the major cycles.


We can see the antenna configuration for this observation by using {{plotants}}:
This algorithm is faster than the Clark algorithm, except when dealing with a large number of visibility samples, due to the re-gridding process it undergoes. It is most useful in cleaning sensitive high-resolution images at lower frequencies where a number of confusing sources are within the primary beam.
-->
In a final step, {{clean}} derives a Gaussian fit to the inner part of the dirty beam, which defines the clean beam. The clean model is then convolved with the clean beam and added to the last residual image to create the final image.


[[Image:plotAnts.png|200px|thumb|right|plotants image showing the antenna configuration during the observation.]]
Note that the CASA team is currently developing a refactored clean task, called  [http://casa.nrao.edu/docs/TaskRef/tclean-task.html tclean]. Task tclean has a better interface, provides new algorithms, and more combinations between imaging algorithms. Task tclean also includes software to parallelize the computations in a multi-processor environment. Eventually, tclean will replace the current {{clean}} task. For this guide, however, we will stick with the original {{clean}} task, as tclean is still in the development and testing phase. Nevertheless, the reader is encouraged to try tclean and send us feedback through the [https://help.nrao.edu/ NRAO helpdesk].
[[Image:amp_v_freq_rawdata.png|200px|thumb|right|plotms image showing strong RFI spikes in amplitude.]]


<source lang="python">
For more details on imaging and deconvolution, we refer to the Astronomical Society of the Pacific Conference Series book entitled [http://www.aspbooks.org/a/volumes/table_of_contents/?book_id=292 Synthesis Imaging in Radio Astronomy II]; the chapter on [http://www.aspbooks.org/a/volumes/article_details/?paper_id=17942 Deconvolution] may prove helpful. Additionally, imaging presentations are available on the [https://science.nrao.edu/science/meetings/2014/14th-synthesis-imaging-workshop/lectures Synthesis Imaging Workshop] and [https://science.nrao.edu/science/meetings/2016/vla-data-reduction/program VLA Data Reduction Workshop] webpages. The [https://casa.nrao.edu/docs/cookbook/index.html  CASA cookbook] chapter on [https://casa.nrao.edu/docs/cookbook/casa_cookbook006.html Synthesis Imaging] provides a wealth of information on the CASA implementation of {{clean}} and related tasks.
# In CASA
plotants(vis='SNR_G55_10s.ms', figfile='SNR_G55_10s.plotants.png')
</source>


This shows that antennas ea01, ea03, and ea18 were on the extreme ends of the west, east, and north arms, respectively.  The antenna position diagram is particularly useful as a guide to help determine which antenna to use as the reference antenna later during calibration. Note that antennas on stations 8 of each arm (N08, E08, W08) do not get moved during array reconfigurations, they can therefore at times be good choices as reference antennas. In this case, we'll probably want to choose something closer to the center of the array.  
Finally, we refer users to the [https://science.nrao.edu/facilities/vla/docs/manuals/oss VLA Observational Status Summary] and the [https://science.nrao.edu/facilities/vla/docs/manuals/obsguide Guide to Observing with the VLA] for information on the VLA capabilities and observing strategies.


We may also inspect the raw data using {{plotms}}. To start with, lets look at a subset of scans on the supernova remnant:


<source lang="python">
== Weights and Tapering ==
# In CASA
plotms(vis='SNR_G55_10s.ms', scan='30,75,120,165,190,235,303', antenna='ea24', xaxis='freq',
      yaxis='amp', coloraxis='spw', iteraxis='scan', correlation='RR,LL')
</source>


* coloraxis='spw': Parameter indicates that a different color will be assigned to each spectral window.
[[Image:SNR-G55-uvcoverage-CASA4.5.2.png|400px|thumb|right|'''Figure 2''' <br /> u,v coverage for the 8 hour observation of the supernova remnant G055.7+3.4]]
* antenna='ea24' : We chose only information for antenna ea24.
* iteraxis='scan': Parameter tells plotms to display a new plot for each scan.
* correlation='RR,LL': We just want to display the right and left circular polarizations, without the cross-hand terms.  


Flipping through to scan 190, we can see that there is significant time and frequency variable RFI present in the observation, as seen by the large spikes in amplitude. In particular, we can see that several spectral windows are quite badly affected. To determine which spectral windows they are, click on the "Mark Regions" tool at the bottom of the plotms GUI (the open box with a green "plus" sign). Use the mouse to select a few of the highest-amplitude points in each of the spectral windows. Click on the "Locate" button (magnifying glass). Information about the selected areas should now display in the logger window:
When imaging data, a map is created associating the visibilities with the image.  
<br />
The sampling function is modified by a weight function that defines the shape and size of the PSF. Weighting therefore provides some control over the spatial resolution and the surface brightness sensitivity of the map, where either direction can be emphasized.
<br />
<!-- <math> S(u,v) \to S(u,v)W(u,v) </math>.  
<pre>
Due to the fact that each VLA antenna performs slightly differently, different weights should be applied to each antenna. Therefore, the weight column in the data table reflects how much weight each corrected data sample should receive.  
Frequency in [1.30662 1.31407] or [1.32377 1.33569] or [1.67866 1.69581], Amp in [0.16871 0.217097] or [0.153226 0.186613] or [0.172581 0.235968]:
-->
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:19:58.0 BL=ea16@W02 & ea24@W05[14&22] Spw=0 Chan=27 Freq=1.31  Corr=RR X=1.31  Y=0.17187  (1718/144/1718)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:19:58.0 BL=ea19@W04 & ea24@W05[17&22] Spw=0 Chan=27 Freq=1.31  Corr=LL X=1.31  Y=0.178607  (1975/144/1975)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:19:58.0 BL=ea22@N04 & ea24@W05[20&22] Spw=0 Chan=37 Freq=1.33  Corr=RR X=1.33  Y=0.18105  (2250/144/2250)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:08.0 BL=ea19@W04 & ea24@W05[17&22] Spw=0 Chan=27 Freq=1.31  Corr=LL X=1.31  Y=0.21228  (1975/145/1975)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:08.0 BL=ea21@E01 & ea24@W05[19&22] Spw=0 Chan=27 Freq=1.31  Corr=LL X=1.31  Y=0.181632  (2103/145/2103)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:08.0 BL=ea22@N04 & ea24@W05[20&22] Spw=0 Chan=37 Freq=1.33  Corr=RR X=1.33  Y=0.167068  (2250/145/2250)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:28.0 BL=ea19@W04 & ea24@W05[17&22] Spw=0 Chan=27 Freq=1.31  Corr=LL X=1.31  Y=0.216252  (1818323722/147/1975)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:28.0 BL=ea22@N04 & ea24@W05[20&22] Spw=0 Chan=37 Freq=1.33  Corr=RR X=1.33  Y=0.182341  (1818323997/147/2250)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:38.0 BL=ea07@E05 & ea24@W05[6&22]  Spw=0 Chan=27 Freq=1.31  Corr=LL X=1.31  Y=0.175032  (1866691864/148/695)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:38.0 BL=ea09@E06 & ea24@W05[8&22]  Spw=0 Chan=27 Freq=1.31  Corr=RR X=1.31  Y=0.176728  (1866692119/148/950)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:38.0 BL=ea16@W02 & ea24@W05[14&22] Spw=0 Chan=27 Freq=1.31  Corr=LL X=1.31  Y=0.170435  (1866692888/148/1719)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:38.0 BL=ea19@W04 & ea24@W05[17&22] Spw=0 Chan=27 Freq=1.31  Corr=LL X=1.31  Y=0.210123  (1866693144/148/1975)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:38.0 BL=ea21@E01 & ea24@W05[19&22] Spw=0 Chan=27 Freq=1.31  Corr=LL X=1.31  Y=0.185065  (1866693272/148/2103)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:38.0 BL=ea22@N04 & ea24@W05[20&22] Spw=0 Chan=37 Freq=1.33  Corr=RR X=1.33  Y=0.158359  (1866693419/148/2250)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:38.0 BL=ea24@W05 & ea27@E03[22&25] Spw=0 Chan=27 Freq=1.31  Corr=RR X=1.31  Y=0.190659  (1866693783/148/2614)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:48.0 BL=ea21@E01 & ea24@W05[19&22] Spw=0 Chan=27 Freq=1.31  Corr=LL X=1.31  Y=0.197052  (1852669347/149/2103)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:48.0 BL=ea22@N04 & ea24@W05[20&22] Spw=0 Chan=37 Freq=1.33  Corr=RR X=1.33  Y=0.183905  (1852669494/149/2250)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:48.0 BL=ea24@W05 & ea27@E03[22&25] Spw=0 Chan=27 Freq=1.31  Corr=RR X=1.31  Y=0.170542  (1852669858/149/2614)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:58.0 BL=ea21@E01 & ea24@W05[19&22] Spw=0 Chan=27 Freq=1.31  Corr=LL X=1.31  Y=0.179017  (1668509051/150/2103)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:58.0 BL=ea22@N04 & ea24@W05[20&22] Spw=0 Chan=37 Freq=1.33  Corr=RR X=1.33  Y=0.177841  (1668509198/150/2250)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:21:08.0 BL=ea21@E01 & ea24@W05[19&22] Spw=0 Chan=27 Freq=1.31  Corr=LL X=1.31  Y=0.193158  (1866681459/151/2103)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:19:58.0 BL=ea19@W04 & ea24@W05[17&22] Spw=2 Chan=19 Freq=1.686 Corr=RR X=1.686 Y=0.209811  (4844192/162/1958)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:19:58.0 BL=ea19@W04 & ea24@W05[17&22] Spw=2 Chan=19 Freq=1.686 Corr=LL X=1.686 Y=0.224424  (4844193/162/1959)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:08.0 BL=ea19@W04 & ea24@W05[17&22] Spw=2 Chan=19 Freq=1.686 Corr=RR X=1.686 Y=0.231918  (1374179750/163/1958)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:08.0 BL=ea19@W04 & ea24@W05[17&22] Spw=2 Chan=19 Freq=1.686 Corr=LL X=1.686 Y=0.209708  (1374179751/163/1959)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:18.0 BL=ea08@N01 & ea24@W05[7&22]  Spw=2 Chan=19 Freq=1.686 Corr=RR X=1.686 Y=0.17744  (2898499/164/806)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:18.0 BL=ea10@N03 & ea24@W05[9&22]  Spw=2 Chan=19 Freq=1.686 Corr=LL X=1.686 Y=0.177235  (2898756/164/1063)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:18.0 BL=ea19@W04 & ea24@W05[17&22] Spw=2 Chan=19 Freq=1.686 Corr=RR X=1.686 Y=0.230346  (2899651/164/1958)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:18.0 BL=ea19@W04 & ea24@W05[17&22] Spw=2 Chan=19 Freq=1.686 Corr=LL X=1.686 Y=0.199502  (2899652/164/1959)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:28.0 BL=ea08@N01 & ea24@W05[7&22]  Spw=2 Chan=19 Freq=1.686 Corr=RR X=1.686 Y=0.177315  (6083110/165/806)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:28.0 BL=ea10@N03 & ea24@W05[9&22]  Spw=2 Chan=19 Freq=1.686 Corr=LL X=1.686 Y=0.17881  (6083367/165/1063)
Scan=190 Field=G55.7+3.4[1] Time=2010/08/23/05:20:28.0 BL=ea19@W04 & ea24@W05[17&22] Spw=2 Chan=19 Freq=1.686 Corr=RR X=1.686 Y=0.217573  (6084262/165/1958)
Found 45 points (45 unflagged) among 101376 in 0s.
</pre>


[[Image:Lband_RFI.png|400px|thumb|right|Side-by-side plots of RFI within a portion of L-Band and spectral window 0 of the observation.]]
There are three main weighting schemes (see Table 1):  


We can see that Spw 0 and 2, are the worst affected by RFI. Also, we get the corresponding frequency where the RFI is present, which looks to be 1.31 and 1.33GHz for spectral window 0, and 1.686GHz for spectral window 2. Visiting the VLA L-Band RFI [https://science.nrao.edu/facilities/vla/observing/RFI/L-Band website], we find that the 1310 and 1330 MHz 'birdies' are due to FAA ASR radars, and the one at 1686MHz is most likely due to a GOES weather satellite. We can compare the plots of the FAA ASR radar from the website, with that of spw 0, and see that the birdies (spikes) are practically identical and fall within the same frequencies.  
1) '''Natural''' weighting: This weighs data based on their rms only. Data visibility weights are gridded onto a uv-cell and summed. More visibilities in a cell will thus increase the cell's weight, which will usually emphasize the shorter baselines. Natural weighting therefore results in a larger PSF and better surface brightness sensitivity, but also a degraded resolution.  


== Priori Calibration and Flagging ==
2) '''Uniform''' weighting: The weights are first gridded as in natural weighting but then each cell is corrected such that the weights are independent of the number of visibilities inside. Compared to natural weighting, uniform weighting usually emphasizes the longer baselines. Consequently the PSF is smaller, which results in a better spatial resolution of the image. At the same time, however, the surface brightness sensitivity is reduced compared to natural weighting. The uniform weighting of the baselines is a better representation of the uv-coverage and sidelobes are more suppressed.


Normally, before we proceed with further processing, we should check the operator log for the observation to see if there were any issues noted during the run that need to be addressed. The observing log file for this observation can be found [http://www.vla.nrao.edu/operators/logs/2010/8/2010-08-23_0005_AB1345.pdf here].
3) '''Briggs''' weighting: This scheme provides a compromise between natural and uniform weighting. The transition can be controlled with the ''robust'' parameter where ''robust=-2'' is close to uniform and ''robust=2'' is close to natural weighting.  
Briggs weighting therefore offers a compromise for between spatial resolution and surface brightness sensitivity, and ''robust'' values near zero are typically used.


The log has various information, including the start/end times for the observation, frequency bands used, weather, baseline information for recently moved antennas, and any outages or issues that may have been encountered. We can see that antenna ea07 may need position corrections (calibration table SNR_G55_10s.pos will fix this), and several antennas are missing an L-Band receiver, including ea06, ea17, ea20, and ea26. We have already removed these antennas from our MS during the run of split2() and can continue with online flagging.  
Details on the weighting schemes are given in [http://www.aoc.nrao.edu/dissertations/dbriggs/ Daniel Brigg's dissertation (Chapter 3)].


=== Online Flags ===
For a visual comparison between these three weighting schemes, please see the section on "CLEAN with Weights" in this guide.


At the time of importing from the SDM-BDF raw data to a MS, we chose to process the online flags from the Flags.xml file to the FLAG_CMD sub-table within the MS. We also created a txt document which includes a list of online flags.  
[[Image:Weight_Tapering_table.png|450px|thumb|right|'''Table 1''' <br /> Table summarizing the effects of using weights and tapering.]]
'''Tapering''': In conjunction with the above weighting schemes, one can specify the ''uvtaper'' parameter within {{clean}}, which will control the radial weighting of visibilities in the uv-plane. Figure 2 graphically visualizes the uv-coverage during the observing session. The taper in {{clean}} is an elliptical Gaussian function which effectively removes long baselines and degrades the resolution. This may be desirable for extended structures when long baselines contribute a large fraction of the noise.
Tapering can, therefore, increase the surface brightness sensitivity of the data, but will decrease the point source sensitivity. Too aggressive tapering, however, may also take its toll on the surface brightness sensitivity.


The Flags.xml file holds information of flags created during the observation, such as subreflector issues, and antennas not being on source.  
We refer to the [https://casa.nrao.edu/docs/cookbook/casa_cookbook006.html CASA Cookbook Synthesis Imaging] chapter for the details of the weighting implementation in CASA's {{clean}}.
We will now want to apply these online flags to the data, but first, let's create a plot of the flags we are about to apply to get an idea of what will be flagged.
<!--
For a brief intro to the different clean algorithms, as well as other deconvolution and imaging information, please see the website kept by Urvashi R.V. [http://www.aoc.nrao.edu/~rurvashi/ImagingAlgorithmsInCasa/node2.html#SECTION00223000000000000000 here].


[[Image:flaggingreason_vs_time.png|300px|thumb|right|Online Flags]]
There are different options for <math> W(u,v) </math> (parameter ''weighting'' in {{clean}}):


We will employ the {{flagcmd}} task to apply the online flags.  
'''1. Natural''': For the CASA parameter ''weighting=’natural’'', visibilities are weighted only by the data weights, which are calculated during filling and calibration and should be equal to the inverse noise variance on that visibility. Imaging weight <math>W_i</math> of sample <math>i</math> is given by <math>W_i = \omega_i = 1/\sigma_i^2</math>


<source lang="python">
where the data weight <math>\omega_i</math> is determined from <math>\sigma_i</math>, the rms noise on visibility <math>i</math>. When data is gridded into the same uv-cell for imaging, the weights are summed, and thus a higher uv density results in higher imaging weights.
# In CASA
flagcmd(vis='SNR_G55_10s.ms', inpmode='table', reason='any', action='plot', plotfile='flaggingreason_vs_time.png')
</source>


We can see several instances of online flagging in the created image. Most notably, ea28 and ea08 had some subreflector issues througout the observation. Online flags are instances of possible missing data, including:
The weight function can be described as <math> W(u,v) = 1/ \sigma^2 </math>, where <math> \sigma^2 </math> is the noise variance. Natural weighting will maximize point source sensitivity, and provide the lowest rms noise within an image, as well as the highest signal-to-noise. It will also generally give more weight to short baselines, thus the synthesized beam will be large. This form of weighting is the default within the clean task.


* ANTENNA_NOT_ON_SOURCE
'''2. Uniform''': For ''weighting = ’uniform’'', the data weights are calculated as in ’natural’ weighting. The data is then gridded to a number of cells in the uv-plane, and after all data is gridded the uv-cells are re-weighted to have 'uniform' imaging weights. This pumps up the influence on the image of data with low weights (they are multiplied up to be the same as for the highest weighted data), which sharpens resolution and reduces the sidelobe level in the field-of-view, but increases the rms image noise.  
The JVLA antennas have slewing speeds of 20 degrees per minute in azimuth, and 40 degrees per minute in elevation. Some antennas are slower than others, and may take a few more seconds to reach the next source. The antennas can also take a few seconds to settle down due to small oscillations after having slewed.  


* SUBREFLECTOR_ERROR   
For uniform weighting, we first grid the inverse variance <math>\omega_i</math> for all selected data onto a grid with uv cell-size given by <math>2/FOV</math> where <math>FOV</math> is the specified field of view (defaults to the image field of view). This forms the gridded weights <math>W_k</math>. The weight of the <math>i</math>-th sample is then:<math>W_i=\omega_i/W_k</math>
The FRM (Focus Rotation Mount) located at the apex of the antennas, is responsible for focusing the incoming radio signal to the corresponding receiver. They can at times have issues with their focus and/or rotation axes. Being off target, so much as a few fractions of a degree can result in loss of data, depending on the frequency being observed.


Now that we've plotted the online flags, we will apply them to the MS.  
The weight function can be described as <math> W(u,v) = W(u,v) / W_k </math>, where <math> W_k </math> represents the local density of (u,v) points, otherwise known as the gridded weights. This form of weighting will increase the influence of data with lower weight, filling the (u,v) plane more uniformly, thereby reducing sidelobe levels in the field-of-view, but increasing the rms image noise. More weight is given to long baselines, therefore increasing angular resolution. Point source sensitivity is degraded due to the downweighting of some data.  


<source lang="python">
'''3. Briggs''': A flexible weighting scheme, that is a variant of uniform, and avoids giving too much weight to (u,v) points with a low natural weight. Weight function can be described as
# In CASA
<math> W(u,v) = 1/ \sqrt{1+S_N^2/S_{thresh}^2} </math>, where <math> S_N </math> is the natural weight of the cell, <math> S_{thresh} </math> is a
flagcmd(vis='SNR_G55_10s.ms', inpmode='table', reason='any', action='apply', flagbackup=False)
threshold.
</source>
A high threshold will go to a natural weight, where as a low threshold will go to a uniform weight. This form of weighting also has adjustable parameters. The ''robust'' parameter will give variation between resolution and maximum point source sensitivity. It's value can range from -2.0 (close to uniform weight) to 2.0 (close to natural weight). By default, the parameter is set to 0.0, which gives a good trade-off.
<br />
-->


The CASA logger should report the progress as the task applies these flags in chunks. Once it has finished, it will report on the percentage of data that has been flagged.  
<!--
The tapering will apodize, or filter/change the shape of the weight function (which is itself a Gaussian), which can be expressed as: <br />
<math> W(u,v) = e^{-(u^2+v^2)/t^2} </math>, where t is the adjustable tapering parameter.  
-->


=== Shadowed Antennas ===
<!--
Due to the downweight of some data, point source sensitivity can be degraded, as some antennas and baselines will suppressed and, in the extreme case, effectively removed from the data averaging. If your observation was sampled by short baselines, tapering may improve sensitivity to extended structures.
-->
== Primary and Synthesized Beam ==


Since this is the most compact JVLA configuration, there may be instances where one antenna blocks, or "shadows" another.  Therefore, we will run {{flagdata}} to remove these data (CASA Cookbook 3.4.2):
The primary beam of a single antenna defines the sensitivity across the field of view. For the VLA antennas, the main beam can be approximated by a Gaussian with a FWHM equal to <math>90*\lambda_{cm}</math> arcseconds, or <math>45/ \nu_{GHz}</math> arcminutes. Note that there are sidelobes beyond the Gaussian kernel that are sensitive to bright sources (see below). Taking our observed frequency to be the middle of L-band, 1.5 GHz, our primary beam will be around 30 arcmin.  


<source lang="python">
If your science goal is to image a source or field of view that is significantly larger than the FWHM of the VLA primary beam, then creating a mosaic from a number of pointings is usually the preferred method. For a tutorial on mosaicing, see the [https://casaguides.nrao.edu/index.php?title=VLA_Continuum_Tutorial_3C391 3C391 tutorial]. In the following, however, we will discuss methods to image large areas from single pointing data.  
# In CASA
flagdata(vis='SNR_G55_10s.ms', mode='shadow', tolerance=0.0, flagbackup=False)
</source>


In this particular observation, there does not appear to be much data affected by shadowing, as can be seen in the logger report. One reason why this may be the case, is the antennas were pointed at sources high in elevation.  
Since our data were taken in the D-configuration, we can check the [https://science.nrao.edu/facilities/vla/docs/manuals/oss/performance/resolution Observational Status Summary] section on VLA resolution to find that the synthesized beam will be around 46 arcsec. Flagging, weighting and the exact frequency, however, may result in deviations from this value. As we will see later, the synthesized beam of our data hovers around 29 arcsec and for the extreme of uniform weighting around 26"x25". Oversampling the minor axis by a factor of ~3, we will use a cell (pixel) size of 8 arcsec.


=== Zero-Amplitude Data ===
Our field contains bright point sources significantly outside the primary beam. The VLA, in particular when using multi-frequency synthesis (see below), will have significant sensitivity outside the main lobe of the primary beam. Particularly at the lower VLA frequencies, sources that are located outside the primary beam may still be bright enough to be detected with sufficient signal-to-noise in the primary beam sidelobes; which can cause artifacts to interfer with the targeted field in the main part of the primary beam. Such sources need to be cleaned to remove the dirty beam response of these interfering sources from the entire image. This can be done by either creating a very large image or by using outlier fields centered on the strongest sources (see section on outlier fields below). A large image has the added advantage of increasing the field of view for science, albeit at lower sensitivity. Other effects will start to become significant, however, like the non-coplanarity of the sky. Large image sizes will also slow down the deconvolution process. Details are provided in Sanjay Bhatnagar's presentation: [https://science.nrao.edu/science/meetings/2016/vla-data-reduction/copy_of_WF_WB.pdf "Advanced Imaging: Imaging in the Wide-band Wide-field era"] given at the [https://science.nrao.edu/science/meetings/2016/vla-data-reduction/program 2016 VLA data reduction workshop].


In addition to shadowing, there may be times during which the correlator writes out pure zero-valued data.  In order to remove this bad data, we run flagdata to remove any pure zeroes:
The calls to CLEAN within this guide will create images that are 170 arcminutes on a side ([Image Size * Cell Size]*[1arcmin / 60arsec]), or almost 6x the size of the primary beam, catching the first and second sidelobes. This is ideal for showcasing both the problems inherent in such wide-band, wide-field imaging, as well as some of the solutions currently available in CASA to deal with these issues.


<source lang="python">
Note that the execution time of {{clean}} depends on the image sizes. Large images generally take more computing time. There are some values, however, that are computationally not advisable. The logger output will then show a recommendation for the next larger but faster image size. As a rule of thumb we recommend image sizes 2<sup>n</sup> * 10 pixels (e.g., 160, 1280 pixels, etc.) for improved processing speeds.
# In CASA
flagdata(vis='SNR_G55_10s.ms', mode='clip', clipzeros=True, flagbackup=False)  
</source>


Inspecting the logger output which is generated by flagdata shows that there is a small quantity of zero-valued data (4.6%) present in this MS.


=== Quacking ===
== Clean Output Images ==
Now we can utilize the flagdata task one more time in order to run it in quaking mode.


It's common for the array to "settle down" at the start of a scan. Quacking is used to remove data at scan boundaries, and it can apply the same edit to all scans for all baselines.
As a result of the CLEAN algorithm, {{clean}} will create a number of output images. For an ''imagename='<imagename>''', these outputs are:
<source lang="python">
# In CAS
flagdata(vis='SNR_G55_10s.ms', mode='quack', quackinterval=5.0, quackmode='beg', flagbackup=False)
</source>


*quackmode='beg' : Data from the start of each scan will be flagged.
'''<imagename>.residual''': the residual after subtracting the clean model (unit: Jy/beam, where beam refers to the dirty beam).<br>
*quackinterval=5.0: Flag the first 5 seconds of every scan.
'''<imagename>.model''': the clean model, not convolved (unit: Jy/pixel).<br>
'''<imagename>.psf''': the point-spread function (dirty beam)<br>
'''<imagename>.flux''': the normalized sensitivity map. For single pointings this corresponds to the primary beam. <br>
'''<imagename>.image''': the residual + the model convolved with the clean beam. This is the final image (unit: Jy/beam, where beam refers to the clean beam).   <br>


=== Backup Data - Flagmanager ===
Additional images will be created for specific algorithms like multi-term frequency synthesis or mosaicking.


Now that we've applied online flags, clipped zero amplitude data, and removed shadowed data, we will create a backup of the MS using {{flagmanager}}.
<font color=red>'''Important:'''</font> If an image file is present in the working directory '''and''' the same name is provided in ''imagename'', {{clean}} will use that image&#151;in particular the residual and model image&#151;as a starting point for further cleaning. If you want a fresh run of {{clean}}, first remove all images of that name using 'rmtables()':


<source lang="python">
<source lang="python">
# In CASA
# In CASA
flagmanager(vis='SNR_G55_10s.ms', mode='save', versionname='after_priori_flagging')
rmtables('<imagename>.*')
</source>
</source>


From here on forward, if we make a mistake, we can always revert back to this version of the MS by setting mode='restore', and providing the version name we want to restore back to.
This method is preferable over ''rm -rf'' as it also clears the cache.


=== Hanning-Smoothing ===
<font color=red>'''Note that interrupting {{clean}} by Ctrl+C may corrupt your visibilities&#151;you may be better off choosing to let {{clean}} finish.  We are currently implementing a command that will nicely exit to prevent this from happening but, for the moment, please try to avoid Ctrl+C.'''</font>


Strong RFI sources can give rise to the Gibbs phenomenon. To remedy this ringing across the frequency channels, we can employ the Hanning smoothing algorithm via the {{hanningsmooth}} task. Note that running this task overwrites the data, and it cannot be reversed. Also, there is a loss in spectral resolution by a factor of two. In addition, the task allows for the creation of a new MS, or to directly operate on the requested column. For this tutorial, we will be directly operating on the data column. Note that hanning-smoothing will remove amplitude spikes, it is therefore not recommended for spectral analysis related science, such as HI. Also note that hanning smoothing cannot be reversed once you apply it to your data.


Let's create before and after images with the {{plotms}} task, to see the effects of hanning-smoothing the data.
== Dirty Image ==


[[Image:amp_v_freq_before_after_hanning.gif|250px|thumb|right|The effects of applying Hanning-Smoothing to our data]]
First, we will create a dirty image (Figure 3) to see the improvements as we step through several cleaning algorithms and parameters. The dirty image is the true image on the sky, convolved with the dirty beam (PSF). We will do this by running {{clean}} with parameter ''niter=0'', which will not do any cleaning.


<source lang="python">
<source lang="python">
# In CASA
# In CASA
plotms(vis='SNR_G55_10s.ms', scan='190', antenna='ea24', spw='0~2',
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.dirty',  
    xaxis='freq', yaxis='amp', coloraxis='spw', title='Before Hanning',
      imsize=1280, cell='8arcsec', interactive=False, niter=0,
    correlation='RR,LL', plotrange=[1.2,1.8,-0.01,0.25],
      stokes='I', usescratch=F)
    plotfile='amp_v_freq.beforeHanning.png')
 
hanningsmooth(vis='SNR_G55_10s.ms', datacolumn='data')
 
plotms(vis='SNR_G55_10s.ms', scan='190', antenna='ea24', spw='0~3',
    xaxis='freq', yaxis='amp', coloraxis='spw', title='After Hanning',
    correlation='RR,LL', plotrange=[1.2,1.8,-0.01,0.25],
    plotfile='amp_v_freq.afterHanning.png')
</source>
</source>


Note that the second call to plotms includes a spectral window outside of the plot range. This is to force plotms to reload the plot. By default, plotms will not redraw a plot if the inputs are unchanged. Within the GUI, checking the "reload" box next to the "Plot" button will do this as well.  
* ''imagename='SNR_G55_10s.dirty''': the root filename used for the various {{clean}} outputs.  


We can compare the generated plots and take notice that single channel RFI has spread into three channels, but it has also removed some of the worst RFI. This spreading of RFI to other channels will ultimately result in a little more data being flagged.  
* ''imsize=1280'': the image size in number of pixels. A single value will result in a square image.


== Automatic RFI excision ==
* ''cell='8arcsec''': the size of one pixel; again, entering a single value will result in a square pixel size.


Now that we're done with priori flagging, we can move on to removing some of the RFI present with auto-flagging algorithms used within {{flagdata}}. We will employ two CASA tasks, {{tfcrop}} and {{rflag}}. For further details on the two algorithm based tasks we will employ, please see the [https://science.nrao.edu/science/meetings/2016/vla-data-reduction/EVLAWorkshop_RFI_2016_UR.pdf presentation] given at the 5th VLA Data Reduction Workshop by Urvashi Rau.  
* ''niter=0'': this controls the number of iterations done in the minor cycle.


=== TFcrop ===
* ''interactive=False'': For a tutorial that covers more of an interactive clean, please see the [https://casaguides.nrao.edu/index.php?title=VLA_high_frequency_Spectral_Line_tutorial_-_IRC%2B10216 VLA high frequency Spectral Line tutorial - IRC+10216.]


{{tfcrop}} is an algorithm that detects outliers in the 2D time-frequency plane, and can operate on un-calibrated data (non bandpass-corrected). Tfcrop will iterate through chunks of time, and undergo several steps in order to find and excise different types of RFI. (CASA Cookbook 3.4.2.7)
* ''usescratch=F'': controls writing the model visibilities to the model data column. For self-calibration we currently recommend setting usescratch=T.


Step 1: Detect short-duration RFI spikes (narrow-band and broad-band). <br />
* ''stokes='I''': since we have not done any polarization calibration, we only create a total-intensity image. For using CLEAN while including various Stoke's Parameters, please see the [https://casaguides.nrao.edu/index.php?title=EVLA_Continuum_Tutorial_3C391-CASA4.5 3C391 CASA guide.]
Step 2: Search for time-persistent RFI. <br />
Step 3: Search for time-persistent, narrow-band RFI. <br />
Step 4: Search for low-level wings of very strong RFI. <br />


More details on the algorithm steps can be found on this [http://www.aoc.nrao.edu/~rurvashi/TFCrop/TFCropV1/node2.html#SECTION00023000000000000000 webpage]
<source lang="python">
 
# In CASA
We will apply the auto-flagging tfcrop algorithm for each spectral window. We will walk through the first spw and then include the rest with the same command. Let's first plot the corrected data for scan 190, with all spectral windows, so we can compare the data before and after tfcrop.
viewer('SNR_G55_10s.dirty.image')
</source><br />


<source lang="python">
<source lang="python">
# In CASA
# In CASA
plotms(vis='SNR_G55_10s.ms', scan='190', antenna='ea24', xaxis='freq', iteraxis='scan', yaxis='amp',
viewer('SNR_G55_10s.dirty.psf')
      ydatacolumn='data', plotfile='amp_v_freq_before_tfcrop.png', title='Before TFcrop',
      correlation='RR,LL', coloraxis='spw', plotrange=[1.2,2,-0.01,0.25])
</source>
</source>


Following are a set of flagdata commands which have been found to work reasonably well with these data. Please take some time to play with the parameters and the plotting capabilities. Since these runs set display='both' and action='calculate', the flags are displayed but not actually written to the MS. This allows one to try different sets of parameters before actually applying the flags to the data.
[[Image:SN_G55_10s.dirty.image.psf.png|400px|thumb|right|'''Figure 3''' <br /> A dirty image of the supernova remnant G55.7+3.4 in greyscale, with apparent sidelobes (left), and the point spread function (PSF), also known as the dirty beam (right).]]


Some representative plots are also displayed. Each column displays an individual polarization product; since we're using all four polarizations, from left to right are RR, RL, LR, and LL. The first row shows the data with current flags applied, and the second includes the flags generated by flagdata. The x-axis is channel number (the spectral window ID is displayed in the top title) and the y-axis of the first two rows is all integrations included in a time "chunk", set by the ntime parameter. These are the data considered by the tfcrop algorithm during its flagging process, and changes in ntime will have some (relatively small) affect on what data are flagged.  
The images may be easier to see in grey scale. To do this, click on Data Display Options (wrench icon, upper left corner) within the viewer and change the color map to Greyscale 1. You may also have to change the scaling power options to your liking. To change the brightness and contrast, assign a mouse button to this type of editing by clicking on the Colormap fiddling button (black/white circle icon) and click/drag the mouse over the image to change the brightness (left-right mouse movement) and contrast (up-down mouse movement).  


Each plot page displays data for a single baseline and time chunk. The buttons at the bottom allow one to step through baseline (backward as well as forward), spw, scan, and field; "Stop Display" will continue the flagging operation without the GUI, and "Quit" aborts the run.
Note that the clean beam is only defined after some clean iterations. The dirty image has therefore no beam size specified in the header and the PSF image is the representation of the response of the array to a point source. Even though it is empty with no iterations specified, {{clean}} will still produce the model file. So a continuation into actual cleaning will be possible by just restarting {{clean}} with the same image root name (and parameter ''niter>0'').


''' spw 0 '''


<source lang="python">
== Regular CLEAN & RMS Noise ==
# In CASA
flagdata(vis='SNR_G55_10s.ms', mode='tfcrop', spw='0',
        datacolumn='data', action='calculate',
        display='both', flagbackup=False)
</source>


As we iterate through the different scans and baselines, we can see the flagging we can apply, represented by the blue areas. Let's now apply these flags by changing the action parameter.
We will now create a regular clean image using mostly default values to see how deconvolution improves the image quality. The first run of {{clean}} will use a fixed number of minor cycle iterations of ''niter=1000'' (default is 500), the second will have ''niter=10000''. Note that you may have to play with the image color map and brightness/contrast to get a better view of the image details.  


<source lang="python">
<source lang="python">
# In CASA
# In CASA. Create default clean image.
flagdata(vis='SNR_G55_10s.ms', mode='tfcrop', spw='0',  
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.Reg.Clean.niter1K',  
        datacolumn='data', action='apply',  
      imsize=1280, cell='8arcsec', niter=1000, interactive=False)
        display='', flagbackup=False)
</source><br />
</source>
 
The logger will report the percentage of flagged data in the table selection. We can now apply the tfcrop algorithm to the remaining spectral windows.
 
[[Image:amp_v_freq_before_after_tfcrop_scan190.gif|250px|thumb|right|The effects of applying tfcrop to our data, for scan 190.]]


''' spw 1, 2, 3 '''
<source lang="python">
<source lang="python">
# In CASA
# In CASA.
flagdata(vis='SNR_G55_10s.ms', mode='tfcrop', spw='1~3',
viewer('SNR_G55_10s.Reg.Clean.niter1K.image')
        datacolumn='data', action='apply',
        display='both', flagbackup=False)
</source>
</source>


After tfcrop has gone through and flagged some of the worst RFI, we can inspect the log report and take note of how much has been flagged.
The logger indicates that the image was obtained in two major cycles and some improvements over the dirty image are visible. But clearly we have not cleaned deep enough yet; the image still has many sidelobes, and an inspection of the residual image shows that it still contains source flux and structure. So let's increase the ''niter'' parameter value to 10,000 and compare the images (Figure 4).


We can also use {{plotms}} to review the effects of using tfcrop and compare the image to the one created before applying tfcrop. We can see great improvements, especially for spectral window 0 and 2, which had some of the worst RFI.
<source lang="python">
# In CASA. Create default clean image with niter = 10000
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.Reg.Clean.niter10K',  
      imsize=1280, cell='8arcsec', niter=10000, interactive=False)
</source><br />


<source lang="python">
<source lang="python">
# In CASA
# In CASA.
plotms(vis='SNR_G55_10s.ms', scan='190', antenna='ea24', xaxis='freq', iteraxis='scan', yaxis='amp',
viewer('SNR_G55_10s.Reg.Clean.niter10K.image')
      ydatacolumn='data', plotfile='amp_v_freq_after_tfcrop.png', title='After TFcrop',
      correlation='RR,LL', coloraxis='spw', plotrange=[1.2,2,-0.01,0.25])
</source>
</source>


We now calculate the amount of flagged data so far in our Measurement Set by using the flagdata task with mode='summary', and apply it to a variable.  
[[Image:SNR_G55_10s.niter.1K_vs_10K.png|400px|thumb|right|'''Figure 4''' <br /> Regular run of CLEAN, with niter=1000 (left), compared to niter=10000 (right).]]


<source lang="python">
As we can see from the resulting images, increasing the niter values (minor cycles) improves our image by reducing prominent sidelobes significantly. One could now further increase the ''niter'' parameter until the residuals are down to an acceptable level. To determine the number of iterations, one needs to keep in mind that {{clean}} will fail once it starts cleaning too deeply into the noise. At that point, the cleaned flux and the peak residual flux values will start to oscillate as the number of iterations increase. This effect can be monitored on the CASA logger output. To avoid cleaning too deeply, we will set a ''threshold'' parameter that will stop minor cycle clean iterations once a peak residual value is being reached.
# In CASA
flagInfo = flagdata(vis='SNR_G55_10s.ms', mode='summary')


print("\n %2.1f%% of G55.7+3.4, %2.1f%% of 3C147, and %2.1f%% of J1925+2106 are flagged. \n" % (100.0 * flagInfo['field']['G55.7+3.4']['flagged'] / flagInfo['field']['G55.7+3.4']['total'], 100.0 * flagInfo['field']['0542+498=3C147']['flagged'] / flagInfo['field']['0542+498=3C147']['total'], 100.0 * flagInfo['field']['J1925+2106']['flagged'] / flagInfo['field']['J1925+2106']['total']))
First, we will utilize the ''SNR_G55_10s.Reg.Clean.niter1K.image'' image to give us an idea of the rms noise (your sigma value)(Figure 5).
With this image open within the viewer, click on the Rectangle Drawing button (Rectangle with R) and draw a square on the image at a position with little source or sidelobe contamination. Doing this should open up a Regions dock which holds information about the selected region, including the pixel statistics in the Statistics tab  (double clicking on the box will also bring up this information). Take notice of the rms values as you click/drag the box around empty image locations, or by drawing additional boxes at suitable positions. If the Regions dock is not displayed, click on '''View''' in the menu bar at the top of the viewer and click on the Regions check-box.


print("Spectral windows are flagged as follows:")
[[Image:SNR_G55_10s.rms.screen.png|400px|thumb|right|'''Figure 5''' <br /> Attempting to find the lowest rms value within the CLEAN'ed image using niter=1000, in order to calculate our threshold.]]


for spw in range(0,4):
The lowest rms value we found was in the order of about 4 x 10<sup>-5</sup> Jy/beam, which we will use to calculate our threshold. There really is no set standard, but fairly good threshold values can vary anywhere between 2.0&#150;4.0*sigma; using clean boxes (see the section on interactive cleaning) allows one to go to lower thresholds. For our purposes, we will choose a threshold of 2.5*sigma. Doing the math results in a value of 10 x 10<sup>-5</sup> or, equivalently, 0.10mJy/beam. Therefore, for future calls to the {{clean}} task, we will set parameter ''threshold='0.1mJy'''. The clean cycle will be stopped when the residual peak flux density is either equal to or less than the ''threshold'' value '''or''' when the maximum number of iterations ''niter'' is reached. To ensure that the stopping criterion is indeed the parameter ''threshold'', parameter ''niter'' should be set to a very high number.
    print("SPW %s: %2.1f%%" % (spw, 100.0 * flagInfo['spw'][str(spw)]['flagged'] / flagInfo['spw'][str(spw)]['total']))
</source>


The results indicate we have flagged a little over 22% of G55.7+3.4. We can now move on to the other auto-flagging algorithm, rflag.  
'''In the following, we nevertheless will use ''niter=1000'' to keep the execution times of {{clean}} on the low end as we focus on explaining different imaging methods.'''


=== RFlag ===
An alternative method to determine the approximate rms of an image is to use the [https://science.nrao.edu/facilities/vla/docs/manuals/propvla/determining VLA Exposure Calculator] and to enter the observing conditions. Make sure that the chosen bandwidth reflects the data after RFI excision.


In order to get the best possible result from the automatic RFI excision with {{rflag}}, we will first apply bandpass calibration to the MS.  Since the RFI is time-variable, using the phase calibration source to make an average bandpass over the entire observation will mitigate the amount of RFI present in the calculated bandpass.  (For the final calibration, we will use the designated bandpass source 3C147; however, since this object was only observed in the last set of scans, it doesn't sample the time variability and would not provide a good average bandpass.)


Since there are likely to be gain variations over the course of the observation, we will run {{gaincal}} to solve for an initial set of antenna-based phases over a narrow range of channels.  These will be used to create the bandpass solutions.  While amplitude variations will have little effect on the bandpass solutions, it is important to solve for these phase variations with sufficient time resolution to prevent decorrelation when vector averaging the data in computing the bandpass solutions.
== CLEAN with Weights ==


In order to choose a narrow range of channels for each spectral window which are relatively RFI-free over the course of the observation, we can look at the data with {{plotms}}. Note that it's important to only solve for phase using a narrow channel range, since an antenna-specific delay will cause the phase to vary with respect to frequency over the spectral window, perhaps by a substantial amount.
To see the effects of using different weighting schemes to the image, let's change the ''weighting'' parameter within {{clean}} and inspect the resulting images. We will be using the ''Natural, Uniform, and Briggs'' weighting algorithms. Here, we have chosen a smaller image size to mainly focus on our science target.  


<source lang="python">
<source lang="python">
# In CASA
# In CASA. Natural weighting
plotms(vis='SNR_G55_10s.ms', scan='42,65,88,11,134,157', antenna='ea24',  
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.natural', weighting='natural',
      xaxis='channel', yaxis='amp', iteraxis='spw', yselfscale=True, correlation='RR,LL')
      imsize=540, cell='8arcsec', niter=1000, interactive=False, threshold='0.1mJy')
</source>
</source>


* yselfscale=True: sets the y-scaling to be for the currently displayed spectral window, since some spectral windows have much worse RFI and will skew the scale for others.
* ''weighting'': specification of the weighting scheme. For Briggs weighting, the robust parameter will be used.


Looking at these plots, we can choose appropriate channel ranges for each SPW:
* ''threshold='0.1mJy''': threshold at which the cleaning process will halt.


<pre>
<source lang="python">
SPW 0: 20-24
# In CASA. Uniform weighting
SPW 1: 49-52
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.uniform', weighting='uniform',
SPW 2: 38-41
      imsize=540, cell='8arcsec', niter=1000, interactive=False, threshold='0.1mJy')
SPW 3: 41-44
</source><br />
</pre>


Using these channel ranges, we run {{gaincal}} to calculate phase-only solutions that will be used as input during our initial bandpass calibration. Remember - the calibration tables we are creating now are so that we can use an automatic RFI flagging algorithm. Our final calibration tables will be generated later, after automated flagging. Here are the inputs for our initial pre-bandpass phase calibration:
<source lang="python">
# In CASA. Briggs weighting, with robust set to default of 0.0
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.briggs', weighting='briggs',
      robust=0, imsize=540, cell='8arcsec', niter=1000, interactive=False, threshold='0.1mJy')
</source><br />


<source lang="python">
<source lang="python">
# In CASA
# In CASA. Open the viewer and select the created images.
gaincal(vis='SNR_G55_10s.ms', caltable='SNR_G55_10s.initPh', field='J1925+2106', solint='int',
viewer()
        spw='0:20~24,1:49~52,2:38~41,3:41~44', refant='ea24', minblperant=3,
        minsnr=3.0, calmode='p')
</source>
</source>


* caltable='SNR_G55_10s.initPh': this is the output calibration table that will be written.
[[Image:SNR_G55_11s.nat.uni.briggs.png|1000px|thumb|center|'''Figure 6''' <br /> CLEAN images created with different weighting algorithms, including natural (left), uniform (center), and briggs (right).]]
* field='J1925+2106': this is the phase calibrator we will use to calibrate the phases.
* solint='int': we request a solution for each 10-second integration.
* pw='0:20~24,1:49~52,2:38~41,3:41~44': note the syntax of this selection: a ":" is used to separate the SPW from channel selection, and "~" is used to indicate an inclusive range. 
* refant='ea24': we have chosen ea24 as the reference antenna after inspecting the antenna position diagram (see above).  It is relatively close to, but not directly in, the center of the array, which could be important in D-configuration, since you don't want the reference antenna to have a high probability of being shadowed by nearby antennas.
* minblperant=3: the minimum number of baselines which must be present to attempt a phase solution.
* minsnr=3.0: the minimum signal-to-noise a solution must have to be considered acceptable.  Note that solutions which fail this test will cause these data to be flagged downstream of this calibration step.
* calmode='p': perform phase-only solutions.
* gaintable='SNR_G55_10s.pos': use the antenna position correction for ea07 that we included in the tar file.  


Note that a number of solutions do not pass the requirements of the minimum 3 baselines (generating the terminal message "Insufficient unflagged antennas to proceed with this solve.") or minimum signal-to-noise ratio (outputting "n of x solutions rejected due to SNR being less than 3 ..."). The logger output indicates 324 solutions succeeded, out of 387 attempted.  
Figure 6 shows that the natural weighted image is most sensitive to extended emission (beam size of 46"x41"). The negative values around the extended emission, often referred to as a negative bowl, is a typical signature of unsampled visibilities near the origin of the uv-plane. That is, the flux density present at shorter baselines, or larger angular scales than measured in this observation, is not well represented. Uniform weighted data shows the highest resolution (26"x25") and Briggs parameter ''robust=0'' (default value) is a compromise with a beam of 29"x29". To image more of the extended emission, the ''robust'' parameter could be tweaked further toward more positive values.


[[Image:gain.phase_v_time.plotcal.png|250px|thumb|right| Gain Phase vs. Time for spw0, for several antennas.]]


It's always a good idea to inspect the resulting calibration table with {{plotcal}}. We will first iterate over antennas for spectral window 0.
== Multi-Scale CLEAN ==
 
Since G55.7+3.4 is an extended source with many angular scales, a more advanced form of imaging involves the use of multiple scales. MS-CLEAN is an extension of the classical CLEAN algorithm for handling extended sources. Multi-Scale CLEAN works by assuming the sky is composed of emission at different angular scales and works on them simultaneously, thereby creating a linear combination of images at different angular scales. For a more detailed description of Multi-Scale CLEAN, see T.J. Cornwell's paper [http://arxiv.org/abs/0806.2228 Multi-Scale CLEAN deconvolution of radio synthesis images].
 
We will use a set of scales&#151;expressed in units of the requested pixel or cell size&#151;which are representative of the scales that are present in the data, including a zero-scale for point sources.  


<source lang="python">
<source lang="python">
# In CASA
# In CASA
plotcal(caltable='SNR_G55_10s.initPh', xaxis='time', yaxis='phase', iteration='antenna',  
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.MultiScale',
        spw='0', plotrange=[-1,-1,-180,180])
      imsize=1280, cell='8arcsec', multiscale=[0,6,10,30,60], smallscalebias=0.9,
      interactive=False, niter=1000,  weighting='briggs', stokes='I',  
      threshold='0.1mJy', usescratch=F, imagermode='csclean')
</source>
</source>


We can see that the phase does not change much over the course of the observation. Let's now inspect the other spectral windows for a few of the antennas.  
* ''multiscale=[0,6,10,30,60]'': a set of scales on which to clean.  A good rule of thumb when using multiscale is 0, 2xbeam, 5xbeam (where beam is the synthesized beam), and larger scales up to about half the minor axis maximum scale of the mapped structure. Since these are in units of the pixel size, our chosen values will be multiplied by the requested cell size. Thus, we are requesting scales of 0 (a point source), 48, 80, 240, and 480 arcseconds (8 arcminutes).  Note that 16 arcminutes (960 arcseconds) roughly corresponds to the size of G55.7+3.4.
 
* ''smallscalebias=0.9'': This parameter is known as the '''small scale bias''', and helps with faint extended structure, by balancing the weight given to smaller structures, which tend to have higher surface brightness but lower integrated flux density. Increasing this value gives more weight to smaller scales. A value of 1.0 weights the largest scale to zero, and a value of less than 0.2 weighs all scales nearly equally. The default value is 0.6.


<source lang="python">
<source lang="python">
# In CASA
# In CASA
plotcal(caltable='SNR_G55_10s.initPh', xaxis='time', yaxis='phase', iteration='spw',
viewer('SNR_G55_10s.MultiScale.image')
        antenna='ea01,ea05,ea10,ea24', plotrange=[-1,-1,-180,180])
</source>
</source>


Now that we've determined our plots look fairly reasonable, we will create a time-averaged bandpass solutions for the phase calibration source using the {{bandpass}} task.
The logger will show how much cleaning is performed on the individual scales.


<source lang="python">
This is the fastest of the imaging techniques described here, but it's easy to see that there are artifacts in the resulting image (Figure 7). We can use the {{viewer}} to explore the point sources near the edge of the field by zooming in on them (Figure 8). Click on the Zoom button on the upper left corner and highlight an area by making a square around the portion where you would like to zoom in. Double clicking within the square will zoom in the selected area. The square can be resized by clicking/dragging the corners, or removed by pressing the Esc key. After zooming in, we can see some radio sources have prominent arcs as well as spots with a six-pointed pattern surrounding them.
# In CASA
bandpass(vis='SNR_G55_10s.ms', caltable='SNR_G55_10s.initBP', field='J1925+2106', solint='inf', combine='scan',
        refant='ea24', minblperant=3, minsnr=10.0, gaintable='SNR_G55_10s.initPh',
        interp='nearest', solnorm=False)
</source>


* solint='inf', combine='scan': the solution interval of 'inf' will automatically break by scans; this requests that the solution intervals be combined over scans, so that we will get one solution per antenna.
{|
* gaintable= 'SNR_G55_10s.initPh': we will pre-apply the initial phase solutions.
| [[Image:SN_G55.MultiScale.image.png|200px|thumb|left|'''Figure 7''' <br /> G55.7+3.4 Multi-Scale Clean]]
* interp='nearest': by default, {{gaincal}} will use linear interpolation for pre-applied calibration.  However, we want the <i>nearest</i> phase solution to be used for a given time.
| [[Image:SN_G55_MultiScale.artifacts.png|200px|thumb|center|'''Figure 8''' <br /> Artifacts around point sources]]
|}


Again, we can see that a number of solutions have been rejected by our choice of <tt>minsnr</tt>.
Next we will explore some more advanced imaging techniques to mitigate the artifacts seen towards the edges of the image.


[[Image:GainAmp_vs_Freq_bandpass.png|250px|thumb|right|Bandpasses for antennas ea01 - ea09]]
[[Image:plotcal_ea01ea05_interactive_flagging.png|250px|thumb|right|Interactive plotcal flagging of the offset points for ea01 and ea05.]]


Let us now inspect the resulting bandpass plots with plotcal.  
== Multi-Scale, Wide-Field CLEAN (w-projection) ==
[[Image:Faceting.png|250px|thumb|right|'''Figure 9''' <br /> Faceting when using wide-field gridmode, which can be used in conjunction with w-projection.]]
[[Image:SNR_G55_MS_vs_MS.wProj.png|500px|thumb|right|'''Figure 10''' <br /> Multi-Scale image of arcs around point sources far from the phase center, versus MS with w-projection (right). We can see the that combining the w-projection algorithm with the multiscale algorithm improves the resulting image by removing prominent artifacts.]]


<source lang="python">
The next {{clean}} algorithm we will employ is w-projection, which is a wide-field imaging technique that takes into account the non-coplanarity of the baselines as a function of distance from the phase center (Figure 9). For wide-field imaging, the sky curvature and non-coplanar baselines result in a non-zero w-term. The w-term, introduced by the sky and array non-coplanarity, introduces a phase term that will limit the dynamic range of the resulting image. Applying 2-D imaging to such data will result in artifacts around sources away from the phase center, as we saw in running MS-CLEAN. Note that this affects mostly the lower frequency bands where the field of view is large.
# In CASA
plotcal(caltable='SNR_G55_10s.initBP', xaxis='freq', yaxis='amp',
        iteration='antenna', subplot=331)
</source>


* subplot=331: displays 3x3 plots per screen
The w-term can be corrected by faceting (which describes the sky curvature by many smaller planes) in either the image or uv-plane. The latter is known as w-projection. A combination of the two can also be employed within {{clean}} by setting the parameter ''gridmode='widefield'''. If w-projection is employed, it will be done for each facet. Note that w-projection is an order of magnitude faster than the image-plane based faceting algorithm, but will require more memory.


We notice that antenna's ea01 and ea05 have a point that is offset from the rest. Let's plot just these two antennas, locate the point on the plot, and flag it interactively through plotcal. Please note that interactive flagging within plotcal will not create a backup, therefore it may be wise to use {{flagmanager}} before doing so.  
For more details on w-projection, as well as the algorithm itself, see [http://adsabs.harvard.edu/abs/2008ISTSP...2..647C "The Noncoplanar Baselines Effect in Radio Interferometry: The w-projection Algorithm"]. Also, the chapter on [http://www.aspbooks.org/a/volumes/article_details/?paper_id=17953 Imaging with Non-Coplanar Arrays] may be helpful.  


<source lang="python">
<source lang="python">
# In CASA
# In CASA
plotcal(caltable='SNR_G55_10s.initBP', antenna='ea01,ea05', xaxis='freq', yaxis='amp',
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.ms.wProj',
        iteration='antenna', subplot=211)
      gridmode='widefield', imsize=1280, cell='8arcsec',
      wprojplanes=-1, multiscale=[0,6,10,30,60],
      interactive=False, niter=1000,  weighting='briggs',
      stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean')
</source>
</source>


We will highlight the points by clicking on the "Mark Region" button, drawing boxes over the points, and clicking on the "Flag" button. Before doing this, one could also get information on the points by clicking on the "Locate" button, which will display details about the highlighted regions. After having flagged the offset points, your plots should update.
* ''gridmode='widefield''': Use the w-projection algorithm.


We will now apply the bandpass calibration table using {{applycal}}.
* ''wprojplanes=-1'': The number of w-projection planes to use for deconvolution. Setting to -1 forces CLEAN to utilize an acceptable number of planes for the given data.


<source lang="python">
<source lang="python">
# In CASA
# In CASA
applycal(vis='SNR_G55_10s.ms', gaintable= 'SNR_G55_10s.initBP', calwt=False)
viewer('SNR_G55_10s.ms.wProj.image')
</source>
</source>


This operation will flag data that correspond to flagged solutions, so {{applycal}} makes a backup version of the flags prior to operating on the data. Note that running {{applycal}} might take a little while.
This will take slightly longer than the previous imaging round; however, the resulting image (Figure 10) has noticeably fewer artifacts.  In particular, compare the same outlier source in the Multi-Scale w-projected image with the Multi-Scale-only image: note that the swept-back arcs have disappeared. There are still some obvious imaging artifacts remaining, though.


Now that we have bandpass-corrected data with some RFI flagged out, we will run flagdata in {{rflag}} mode (CASA Cookbook 3.4.2.8). RFlag, like TFcrop, is an autoflag algorithm which uses a sliding window statistical filter. Data is iterated through in chunks of time, where statistics are accumulated, and thresholds calculated.


Additional information on the algorithm used in rflag, as well as the statistical details that are undertaken can be found [http://www.aoc.nrao.edu/~rurvashi/FlaggerDocs/node5.html on this webpage] (sections 2.1.7).
== Multi-Scale Multi-Term Frequency Synthesis ==
[[Image:MultiFrequency_Synthesis_snapshot.png|250px|thumb|right|'''Figure 11''' <br /> Multi-Frequency Synthesis snapshot of (u,v) coverage. We can see from the image on the right, using this algorithm can greatly improve coverage, thereby improving image fidelity.]]
[[Image:SNR_G55_10s.ms.MFS.image.png|200px|thumb|right|'''Figure 12''' <br /> Multi-Scale imaged with MS-MFS and nterms=2. Artifacts around point sources close to phase center are reduced while point sources away from the phase center still show artifacts.]]
[[Image:SN_G55_MS.MFS.alpha.png|200px|thumb|right|'''Figure 13''' <br /> Spectral Index image]]


We will use flagdata with mode='rflag', and action='calculate' to first review the amount of data to be flagged. We will also change the datacolumn parameter to 'corrected', since we've applied the bandpass corrections to the MS, and in the process, created the ''corrected_data'' column in the MS table.  
A consequence of simultaneously imaging the wide fractional bandwidths available with the VLA is that the primary and synthesized beams have substantial frequency-dependent variation over the observing band. If this variation is not accounted for, it will lead to imaging artifacts and compromise the achievable image rms.


''' spw 0 '''
The coordinates of the (u,v) plane are measured in wavelengths. Observing at several frequencies results in a baseline sampling several ellipses of different sizes in the (u,v) plane. By using a method called Multi-Frequency Synthesis (MFS), which is the default cleaning mode in CLEAN, we can fill in the gaps in the single frequency (u,v) coverage to achieve much better image fidelity (Figure 11).


First, let's create a plot of our corrected data before we apply RFlag.
When observing in low-frequencies, it may prove beneficial to observe in short segments which are spread out in time. This will allow the coverage of more spatial-frequencies, permitting us to employ this algorithm more efficiently.  


<source lang="python">
The Multi-Scale Multi-Term Frequency Synthesis algorithm provides the ability to simultaneously image and fit for the intrinsic source spectrum in each pixel. The spectrum is approximated using a polynomial Taylor-term expansion in frequency, with the degree of the polynomial as a user-controlled parameter. A least-squares approach is used, along with the standard clean-type iterations.  
# In CASA
plotms(vis='SNR_G55_10s.ms', scan='190' , antenna='ea24',
xaxis='freq', yaxis='amp', ydatacolumn='corrected', coloraxis='spw', title='Before RFlag',
correlation='RR,LL', plotfile='amp_v_freq.beforeRFlag.png')
 
flagdata(vis='SNR_G55_10s.ms', mode='rflag', spw='0', datacolumn='corrected',
        action='calculate', display='both', flagbackup=False)
</source>


[[Image:rflag_default_values.png|200px|thumb|left|RFlag with Default Parameters]]
For a more detailed explanation of the MS-MFS deconvolution algorithm, please see the paper [http://arxiv.org/abs/1106.2745 A multi-scale multi-frequency deconvolution algorithm for synthesis imaging in radio interferometry]
[[Image:rflag_cutoff_1.5_sigma.png|200px|thumb|right|RFlag with a cutoff of 1.5 Sigma, shown in blue]]
 
We can see that a lot of RFI is being missed by the default timedevscale and freqdevscale parameters of 5.0. Since we know spw 0 has lots of RFI, we will change the values to 1.5 and apply them. Note that we can always decide that this still isn't good enough and click on "Quit" to change our parameter values.  


<source lang="python">
<source lang="python">
# In CASA
# In CASA
flagdata(vis='SNR_G55_10s.ms', mode='rflag', spw='0', datacolumn='corrected',  
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.ms.MFS',
        freqdevscale=1.5, timedevscale=1.5, action='apply', display='both', flagbackup=False)
      imsize=1280, cell='8arcsec', mode='mfs', nterms=2,
      multiscale=[0,6,10,30,60],  
      interactive=False, niter=1000, weighting='briggs',
      stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean')
</source>
</source>
 
Although RFlag has done a pretty good job of finding the bad data, some still remains.  One way to excise the remaining bad data is to use the mode='extend' feature in flagdata, which can extend flags along a chosen axis.  First, we will extend the flags across polarization, so if any one polarization is flagged, all data for that time / channel will be flagged:


<source lang="python">
* ''nterms=2'': the number of Taylor-terms to be used to model the frequency dependence of the sky emission.  Note that the speed of the algorithm will depend on the value used here (more terms will be slower); nterms=2 will fit a spectral index, and nterms=3 a spectral index and curvature.
# In CASA
flagdata(vis='SNR_G55_10s.ms', mode='extend', spw='0', extendpols=True,
        action='apply', display='', flagbackup=False)
</source>


Now, we will extend the flags in time and frequency, using the "growtime" and "growfreq" parameters. For the data here, the <tt>rflag</tt> algorithm seems most likely to miss RFI which should be flagged along more of the time axis, so we will try with growtime=50.0, which will flag all data for a given channel if more than 50% of that channel's time is already flagged, and growfreq=90.0, which will flag the entire spectrum for an integration if more than 90% of the channels in that integration are already flagged.
When clean is done ''<imagename>.image.tt0'' will contain a total intensity image (Figure 12), where ''tt*'' is a suffix to indicate the Taylor term: ''<imagename>.image.tt0'' is the total intensity image and ''<imagename>.image.alpha'' will contain an image of the spectral index in regions where there is sufficient signal-to-noise (Figure 13). Having this spectral index image can help convey information about the emission mechanism involved within the supernova remnant. This spectral index can also give information on the optical depth of the source.  


<source lang="python">
<source lang="python">
# In CASA
# In CASA
flagdata(vis='SNR_G55_10s.ms', mode='extend', spw='0', growtime=50.0,
viewer('SNR_G55_10s.ms.MFS.image.tt0')
        growfreq=90.0, action='apply', display='', flagbackup=False)
</source><br />
</source>
 
''' spw 1, 2, and 3 '''
 
Now, let's work on SPW's 1,2, and 3. We've chosen display='both', in order to first review the flags, and decide if too much or too litle is being flagged. We can always click on 'Quit', and decide to change the freqdevscale and timedevscale parameters, as we did with spectral window 0. Note that spw 2 has these parameters lower than spw 1/3, as it contains more RFI. To allow the flagging to continue, click on "Stop Display".


<source lang="python">
<source lang="python">
# In CASA
# In CASA
flagdata(vis='SNR_G55_10s.ms', mode='rflag',  spw='1,3', datacolumn='corrected',
viewer('SNR_G55_10s.ms.MFS.image.alpha')
        freqdevscale=4.0, timedevscale=4.0, action='apply', display='both', flagbackup=False)
</source>


flagdata(vis='SNR_G55_10s.ms', mode='extend', spw='1,3', extendpols=True,
Note: To replicate Figure 13, open the image within the viewer, click on Panel Display Options (wrench with a small P), change background color to white, and adjust your margins under the Margins button. To view the color wedge, click on the Data Display Options (wrench without the P, next to the folder icon), click on the Color Wedge button, click on Yes under display color wedge and adjust the various parameters to your liking.
        action='apply', display='', flagbackup=False)


flagdata(vis='SNR_G55_10s.ms', mode='extend', spw='1,3', growtime=50.0,
For more information on the multi-term, multi-frequency synthesis mode and its outputs, see section [https://casa.nrao.edu/docs/cookbook/casa_cookbook006.html#sec306 5.2.5.1] in the CASA cookbook.
        growfreq=90.0, action='apply', display='', flagbackup=False)


Inspect the brighter point sources in the field near the supernova remnant. You will notice that some of the artifacts, which had been symmetric around the sources themselves, are now gone. Since we did not use w-projection, however, there are still strong features related to the non-coplanar baseline effects still apparent for sources further away.


flagdata(vis='SNR_G55_10s.ms', mode='rflag',  spw='2', datacolumn='corrected',  
At this point, {{clean}} takes into account the frequency variation of the synthesized beam ''but not'' the frequency variation of the primary beam. For low frequencies and large bandwidths, this can be substantial, e.g., 1&#150;2 GHz L-band observations result in a variation of a factor of 2. One effect of such a large fractional bandwidth is that in multi-frequency synthesis, primary beam nulls will be blurred and the interferometer is sensitive everywhere in the field of view. For spectral slopes, however, a frequency-dependent primary beam causes a steepening as the higher frequencies are less sensitive at a given point away from the phase center. A correction for this effect should be made with the task {{widebandpbcor}} (see section on Primary Beam Correction below).
        freqdevscale=2.5, timedevscale=2.5, action='apply', display='both', flagbackup=False)


flagdata(vis='SNR_G55_10s.ms', mode='extend', spw='2', extendpols=True,
        action='apply', display='', flagbackup=False)


flagdata(vis='SNR_G55_10s.ms', mode='extend', spw='2', growtime=50.0,
== Multi-Scale Multi-Term Frequency, Widefield CLEAN ==
        growfreq=90.0, action='apply', display='', flagbackup=False)
</source>


Let's create an after RFlag plot to see the improvements.
We will now combine the w-projection and MS-MFS algorithms.  Be forewarned&#151;these imaging runs will take a while, and it's best to start them running and then move on to other things.


[[Image:amp_v_freq.before_after_RFlag.gif|250px|thumb|right|Before and after RFlag for scan 190.]]
Using the same parameters for the individual-algorithm images above, but combined into a single {{clean}} run, we have:


<source lang="python">
<source lang="python">
# In CASA
# In CASA
plotms(vis='SNR_G55_10s.ms', scan='190' , antenna='ea24',
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.ms.MFS.wProj',
xaxis='freq', yaxis='amp', ydatacolumn='corrected', coloraxis='spw', title='After RFlag',
      gridmode='widefield', imsize=1280, cell='8arcsec', mode='mfs',
correlation='RR,LL', plotfile='amp_v_freq.afterRFlag.png', plotrange=[1.2,2,0,2.5])
      nterms=2, wprojplanes=-1, multiscale=[0,6,10,30,60], 
</source>
      interactive=False, niter=1000, weighting='briggs',
 
      stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean')
We will now plot amplitude vs. baseline, and look for outliers which we may be able to flag.
</source><br />


<source lang="python">
<source lang="python">
# In CASA
# In CASA
plotms(vis='SNR_G55_10s.ms', scan='30,75,120,165,190,235,303', xaxis='baseline', yaxis='amp',
viewer('SNR_G55_10s.ms.MFS.wProj.image.tt0')
ydatacolumn='corrected', iteraxis='scan', correlation='RR,LL', coloraxis = 'baseline')
</source><br />
</source>
 
It appears that there are a few baselines which have consistently higher-amplitude than the others, indicating that it's probably contaminated by RFI. 
 
Use the plotms tools to identify this baseline, which turns out to be ea04&ea08 and ea04&ea21, and is mostly present for spw 1. Let's flag these baselines.


<source lang="python">
<source lang="python">
# In CASA
# In CASA
flagdata(vis='SNR_G55_10s.ms', antenna='ea04&ea08', spw='1', flagbackup=False)
viewer('SNR_G55_10s.ms.MFS.wProj.image.alpha')
flagdata(vis='SNR_G55_10s.ms', antenna='ea04&ea21', spw='1', flagbackup=False)
</source>
</source>


We will now run flagdata in summary mode to inspect how much data we have flagged thus far. We can also create a plot of the percentage of flagged data with respect to the frequency. This command will also plot the antennas, with circle size representing the percentage of flagged data per antenna.  
[[Image:SNR_G55_MS.MFS.wProj.png|1250px|thumb|center|'''Figure 14'''  Here we see the differences as the images progress through the different algorithms used: MS &rarr; MS-MFS &rarr; MS-wProjection &rarr; MS-MFS-wProjection.]]


<source lang="python">
Again, looking at the same outlier source, we can see that the major sources of error have been removed, although there are still some residual artifacts (Figure 14). One possible source of error is the time-dependent variation of the primary beam; another is the fact that we have only used ''nterms=2'', which may not be sufficient to model the spectra of some of the point sources. Some weak RFI may also show up that may need additional flagging.
# In CASA
flagInfo = flagdata(vis='SNR_G55_10s.ms', mode='summary', action='calculate', display='report', spwchan=T)


print("\n %2.1f%% of G55.7+3.4, %2.1f%% of 3C147, and %2.1f%% of J1925+2106 are flagged. \n" % (100.0 * flagInfo['field']['G55.7+3.4']['flagged'] / flagInfo['field']['G55.7+3.4']['total'], 100.0 * flagInfo['field']['0542+498=3C147']['flagged'] / flagInfo['field']['0542+498=3C147']['total'], 100.0 * flagInfo['field']['J1925+2106']['flagged'] / flagInfo['field']['J1925+2106']['total']))


print("Spectral windows are flagged as follows:")
== Imaging Outlier Fields ==


for spw in range(0,4):
When strong sources are far from the main target, but still strong enough to produce sidelobes in the main image, they should be cleaned. Sometimes, however, it is not practical to image very large fields for this purpose. An alternative is to use outlier fields. This mode will allow the user to specify a few locations that are then cleaned along with the main image. The outlier fields should be centered on strong sources that are known from sky catalogs or are identified by other means.  
    print("SPW %s: %2.1f%%" % (spw, 100.0 * flagInfo['spw'][str(spw)]['flagged'] / flagInfo['spw'][str(spw)]['total']))


</source>
In order to find outlying sources, it will help to image a large field with lower resolution and identify them within the CASA viewer. By moving the mouse cursor over the sources, we can grab their positions and use the information within an outlier file. In this example, we've called our file ''outliers.txt''. It contains the names, sizes, and positions of the outliers we want to image (see section [https://casa.nrao.edu/docs/cookbook/casa_cookbook006.html#sec356 5.3.18.1] of the CASA cookbook).


So, as a result of the flagging, we have sacrificed almost 28% of the data for G55.7+3.4, and as expected, spectral windows 0 and 2 have the most flagged data.
Open your favorite text editor and input the following:


== Interactive Flagging ==
<pre>
#content of outliers.txt
#
#outlier field1
imagename= 'Outlier1.MS.MFS'
imsize=[512,512]
phasecenter = 'J2000 19h23m27.693 22d37m37.180'
#
#outlier field2
imagename='Outlier2.MS.MFS'
imsize=[512,512]
phasecenter = 'J2000 19h25m46.888 21d22m03.365'


After plotting our data, it's obvious that there is still some amount of RFI present. We can use plotms to interactively flag points which look bad. A tutorial on flagging data interactively through plotms can be found [https://casaguides.nrao.edu/index.php?title=Data_flagging_with_plotms here]. Please note that flagging data through plotms will not create a backup, so it's important to use the flagmanager before deciding to mark your regions for flagging purposes.
</pre>
 
[[Image:corr.amp_v_freq.interactive_flagging.png|250px|thumb|right|Highlighting and interactively flagging RFI for scan 75.]]


{{clean}} will then be executed like the following:
<source lang="python">
<source lang="python">
# In CASA
# In CASA
flagmanager(vis='SNR_G55_10s.ms', mode='save',  
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR.MS.MFS-Main',
            versionname='before_interactive_flagging')
      outlierfile='outliers.txt',
 
      imsize=[512,512], cell='8arcsec', mode='mfs',
plotms(vis='SNR_G55_10s.ms', scan='30,75,120,165,190,235,303', xaxis='frequency', yaxis='amp',  
      multiscale=[0,6,10,30,60], interactive=False, niter=1000, weighting='briggs',
ydatacolumn='corrected', iteraxis='scan', correlation='RR,LL', coloraxis = 'spw')
      robust=0, stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean')
</source>
</source>


Let's iterate through scans and view scan 75. We can see that spectral window 1 has a large amplitude spike. We can use the "mark regions" button (square with green plus sign) to highlight a majority of spike, and the "locate" button (magnifying glass) to give us information on the highlighted region. We will now click on the "flag" button (red flag) and we should see the plot update, having removed the highlighted region.  
* ''outlierfile='outliers.txt''': the name of the outlier file.


== Imaging ==
We can view the images with {{viewer}} one at a time. Due to the difference in sky coordinates, they cannot be viewed on the same window display.


We will now cover topics on imaging with CASA. We will be utilizing the previous data set, which has been calibrated. A copy of the calibrated data <font color=green>(1.2GB)</font> can be downloaded from [http://casa.nrao.edu/Data/EVLA/SNRG55/SNR_G55_10s.calib.tar.gz http://casa.nrao.edu/Data/EVLA/SNRG55/SNR_G55_10s.calib.tar.gz]
<source lang="python">
# In CASA
viewer('SNR.MS.MFS-Main.image')
</source><br />


We will skip the calibration process in this guide, as examples of calibration can be found in several other guides, including the
<source lang="python">
[https://casaguides.nrao.edu/index.php/EVLA_Wide-Band_Wide-Field_Imaging:_G55.7_3.4-CASA4.4 EVLA_Wide-Band_Wide-Field_Imaging] guide.
# In CASA
viewer('Outlier1.MS.MFS.image')
</source><br />


After undergoing the calibration process, we split out the science target, G55.7+3.4, with only the RR,LL correlations, as the cross correlations were not calibrated, nor will be used for our purposes.  
<source lang="python">
# In CASA
viewer('Outlier2.MS.MFS.image')
</source>


=== The CLEAN Algorithm ===
We have changed the color map to Hot Metal 2 in order to show the different colors that can be used within the viewer. Feel free to play with other color maps that may be better suited for your screen.


[[Image:CLEAN_Cycle.png|450px|thumb|right|The CLEAN major and minor cycles, indicating the steps undertaken during gridding, projection algorithms, and creation of images.]]
[[Image:Outlier_Fields.png|1000px|thumb|center|'''Figure 15''' <br /> Images of the supernova remnant and bright sources far from the phase center (near the edge of, or beyond, the primary beam.) ]]


The CLEAN algorithm, developed by J. Högbom (1974) enabled the synthesis of complex objects, even if they have relatively poor Fourier uv-plane coverage. Poor coverage occurs with partial earth rotation synthesis, or with arrays composed of few antennas. The "dirty" image is formed by a simple Fourier inversions of the sampled visibility data, with each point on the sky being represented by a suitably scaled and centered PSF (Point Spread Function, sometimes called the dirty beam). This algorithm attempts to interpolate from the measured (u,v) points across gaps in the (u,v) coverage. It, in short, provides solutions to the convolution equation by representing radio sources by a number of point sources in an empty field.


The brightest points are found by performing a cross-correlation between the dirty image, and the PSF. The brightest parts are subtracted (minor-cycle), and the process is repeated again for the next brighter sources (major-cycle). A large part of the work in CLEAN involves shifting and scaling the dirty beam.
== Primary Beam Correction ==


The clean algorithm works well with points sources, as well as most extended objects. Where it can fall short is in speed, as convergence can be slow for extended objects, or for images containing several bright point sources. A solution to deconvolve these images would be the MEM (Maximum Entropy Method) algorithm, which has faster performance, although we can improve the CLEAN algorithm by employing other means, some of which will be mentioned below.  
In interferometry, the images formed via deconvolution are representations of the sky multiplied by the primary beam response of the antenna. The primary beam can be described as Gaussian with a size depending on the observing frequency. Images produced via {{clean}} are by default not corrected for the primary beam pattern (important for mosaics), and therefore do not have the correct flux density away from the phase center.  


''' 1. Högbom Algorithm '''<br />
Correcting for the primary beam, however, can be done during {{clean}} by using the parameter ''pbcor''. Alternatively, it can be done after imaging using the task [https://casa.nrao.edu/docs/TaskRef/impbcor-task.html impbcor] for regular data sets and {{widebandpbcorr}} for those that use Taylor-term expansion (nterms > 1). A third alternative is utilizing the task {{immath}} to manually divide the ''<imagename>.image'' by the ''<imagename>.flux'' image (''<imagename>.flux.pbcoverage'' for mosaics).  
This algorithm will initially find the strength and position of a peak in a dirty image, subtract it from the dirty image, record this position and maginitude, and repeat for further peaks. The remainder of the dirty image is known as the residuals.  


The accumulated point sources, now residing in a model, is convolved with an idealized CLEAN beam (usually a Gaussian fitted to the central lobe of the dirty beam), creating a CLEAN image. As the final step, the residuals of the dirty image are then added to the CLEAN image.  
Flux corrected images usually don't look pretty due to the noise at the edges being increased. Flux densities, however, should only be calculated from primary beam corrected images. Let's run the {{impbcor}} task to correct our multiscale image.


''' 2. Clark Algorithm '''<br />
<source lang="python">
Clark (1980), developed a FFT-based CLEAN algorithm, which more efficiently shifts and scales the dirty beam by approximating the position and strength of components using a small patch of the dirty beam. This algorithm is the default within the {{clean}} task, which involves major and minor cycles.
# In CASA
impbcor(imagename='SNR_G55_10s.MultiScale.image', pbimage='SNR_G55_10s.MultiScale.flux',
        outfile='SNR_G55_10s.MS.pbcorr.image')
</source>


The algorithm will first select a beam patch, which will include the highest exterior sidelobes. Points are then selected from the dirty image, which are up to a fraction of the image peak, and are greater than the highest exterior sidelobe of the beam. It will then conduct a list-based Högbom CLEAN, creating a model and convolution with an idealized CLEAN beam. This process is the minor cycle.
* ''imagename'': the image to be corrected


The major cycle involves transforming the point source model via a FFT (Fast-Fourier Transform), mutiplying this by the weight sampling function (more on this below), and transformed back. This is then subtracted from the dirty image, creating your CLEAN image. The process is then repeated with subsequent minor cycles. 
* ''pbimage'': the <imagename>.flux image as a representation of the primary beam (<imagename>.flux.pbcoverage for mosaics)


''' 3. Cotton-Schwab Algorithm '''<br />
<source lang="python">
This is the default imager mode (''csclean''), and is a variant of the Clark algorith in which the major cycle involves the subtraction of CLEAN components of ungridded visibility data. This allows the removal of gridding errors, as well as noise. One advantage is its ability to image and clean many seperate fields simultaneously. Fields are cleaned independently in the minor cycle, and components from all fields cleaned together in the major cycles.
# In CASA
viewer('SNR_G55_10s.MS.pbcorr.image')
</source>


This algorithm is faster than the Clark algorithm, except when dealing with a large number of visibility samples, due to the re-gridding process it undergoes. It is most useful in cleaning sensitive high-resolution images at lower frequencies where a number of confusing sources are within the primary beam.  
Let us now use the {{widebandpbcor}} task for wide-band (''nterms''>1) images. Note that for this task, we will be supplying the image name that is the prefix for the Taylor expansion images, tt0 and tt1, which must be on disk. Such files were created during the last Multi-Scale, Multi-Frequency, Widefield run of CLEAN. Task {{widebandpbcor}} will generate a set of images with a ''pbcor.image'' extension.


For more details on imaging and deconvolution, you can refer to the Astronomical Society of the Pacific Conference Series book entitled [http://www.aspbooks.org/a/volumes/table_of_contents/?book_id=292 Synthesis Imaging in Radio Astronomy II]. The chapter on [http://www.aspbooks.org/a/volumes/article_details/?paper_id=17942 Deconvolution] may prove helpful.
[[Image:SNR_G55_10s.MS.MFS.wProj.pbcor.image.png|200px|thumb|right|'''Figure 16''' <br /> Primary beam corrected image using the widebandpbcor task for the MS.MFS.wProj image created during the CLEAN process.]]
 
[[Image:SNR_G55_10s.MS.MFS.wProj.pbcor.image.alpha.png|200px|thumb|right|'''Figure 17''' <br /> Primary beam corrected spectral index image using the widebandpbcor task.]]
=== Weights and Tapering ===
 
[[Image:SNR_G55_uvcoverage.png|300px|thumb|right| u,v coverage for the 8-hour observation of the supernova remnant G055.7+3.4]]
 
When imaging data, a map is created associating the visibilities with the image. The sampling function, which is a function of the visibilities, is modified by a weight function. <math> S(u,v) \to S(u,v)W(u,v) </math>.
 
This process can be considered a convolution.  The convolution map, is the weights by which each visiblity is multiplied by before gridding is undertaken. Due to the fact that each VLA antenna performs slightly differently, different weights should be applied to each antenna. Therefore, the weight column in the data table reflects how much weight each corrected data sample should receive.
 
For a brief intro to the different clean algorithms, as well as other deconvolution and imaging information, please see the website kept by Urvashi R.V. [http://www.aoc.nrao.edu/~rurvashi/ImagingAlgorithmsInCasa/node2.html#SECTION00223000000000000000 here].
 
The following are a few of the more used forms of weighting, which can be used within the {{clean}} task. Each one has their own benefits and drawbacks. <br />
 
'''1. Natural''': The weight function can be described as <math> W(u,v) = 1/ \sigma^2 </math>, where <math> \sigma^2 </math> is the noise variance. Natural weighting will maximize point source sensitivity, and provide the lowest rms noise within an image, as well as the highest signal-to-noise. It will also generaly give more weight to short baselines, thus angular resolutions can be degraded. This form of weighting is the default within the clean task.
 
'''2. Uniform''': The weight function can be described as <math> W(u,v) = W(u,v) / W_k </math>, where <math> W_k </math> represents the local density of (u,v) points, otherwise known as the gridded weights. This form of weighting will increase the influence of data with lower weight, filling the (u,v) plane more uniformly, thereby reducing sidelobe levels in the field-of-view, but increasing the rms image noise. More weight is given to long baselines, therefore increasing angular resolution. Point source sensitivity is degraded due to the downweighting of some data.  
 
'''3. Briggs''': A flexible weighting scheme, that is a variant of uniform, and avoids giving too much weight to (u,v) points with a low natural weight. Weight function can be described as <math> W(u,v) = 1/ \sqrt{1+S_N^2/S_{thresh}^2} </math>, where <math> S_N </math> is the natural weight of the cell, <math> S_{thresh} </math> is a threshold. A high threshold will go to a natural weight, where as a low threshold will go to a uniform weight. This form of weighting also has adjustable parameters. The ''robust'' parameter will give variation between resolution and maximum point source sensitivity. It's value can range from -2.0 (close to uniform weight) to 2.0 (close to natural weight). By default, the parameter is set to 0.0, which gives a good trade-off.
<br />
 
[[Image:Weight_Tapering_table.png|450px|thumb|right| Table summarizing the effects of using weights and tapering.]]
 
'''I. Tapering''': In conjunction with weighting, we can include the ''uvtaper'' parameter within clean, which will control the radial weighting of visibilities, in the uv-plane. This in effect, reduces the visibilities, with weights decreasing as a function of uv-radius. The tapering will apodize, or filter/change the shape of the weight function (which is itself a Gaussian), which can be expressed as: <br />
<math> W(u,v) = e^{-(u^2+v^2)/t^2} </math>, where t is the adjustable tapering parameter. This process can smooth the image plane, give more weight to short baselines, but in turn degrade angular resolution. Due to the downweight of some data, point source sensitivity can be degraded. If your observation was sampled by short baselines, tapering may improve sensitivity to extended structures.
 
=== Primary and Synthesized Beam ===
 
The primary beam of the VLA antennas can be taken to be a Gaussian with FWHM equal to <math>90*\lambda_{cm}</math>  or <math>45/ \nu_{GHz}</math>. Taking our observed frequency to be the middle of the band, 1.5GHz, our primary beam will be around 30 arcmin. Note that if your science goal is to image a source, or field of view that is significantly larger than the FWHM of the VLA primary beam, then creating a mosaic from a number of pointings would be best. For a tutorial on mosaicing, see the [https://casaguides.nrao.edu/index.php?title=EVLA_Continuum_Tutorial_3C391 3C391 tutorial.]
 
Since our observation was taken in D-configuration, we can check the [https://science.nrao.edu/facilities/vla/docs/manuals/oss/performance/resolution Observational Status Summary]'s section on VLA resolution to find that the synthesized beam will be around 46 arcsec.  We want to oversample the synthesized beam by a factor of around five, so we will use a cell size of 8 arcsec. 
 
Since this field contains bright point sources significantly outside the primary beam, we will create images that are 170 arcminutes on a side ([8arcsec * Cell Size]*[1arcmin / 60arsec]), or almost 6x the size of the primary beam.  This is ideal for showcasing both the problems inherent in such wide-band, wide-field imaging, as well as some of the solutions currently available in CASA to deal with these issues.
 
First, it's worth considering why we are even interested in sources which are far outside the primary beam.  This is mainly due to the fact that the EVLA, with its wide bandwidth capabilities, is quite sensitive even far from phase center -- for example, at our observing frequencies in L-band, the primary beam gain is as much as 10% around 1 degree away.  That means that any imaging errors for these far-away sources will have a significant impact on the image rms at phase center.  The error due to a source at distance R can be parametrized as:
 
<math>
\Delta(S) = S(R) \times PB(R) \times PSF(R)
</math>
 
So, for R = 1 degree, source flux S(R) = 1 Jy, <math>\Delta(S)</math> = 1 mJy − 100 <math>{\mu}</math>Jy.  Clearly, this will be a source of significant error.
 
=== Multi-Scale Clean ===
 
Since G55.7+3.4 is an extended source with many spatial scales, the most basic (yet still reasonable) imaging procedure is to use {{clean}} with multiple scales. MS-CLEAN is an extension of the classical CLEAN algorithm for handling extended sources. It works by assuming the sky is composed of emission at different spatial scales and works on them simultaneously, thereby creating a linear combination of images at different spatials scales. For a more detailed description of Multi Scale CLEAN, see the paper by J.T. Cornwell entitled  [http://arxiv.org/abs/0806.2228 Multi-Scale CLEAN deconvolution of radio synthesis images].
 
It can also be possible to utilize [http://casa.nrao.edu/docs/TaskRef/tclean-task.html tclean], (t for test) which is a refactored version of clean, with a better interface, and provides more possible combinations of algorithms. It also allows for process computing parallelization of the imaging and deconvolution. Eventually, tclean will replace the current clean task, but for now, we will stick with the original clean, as tclean is merely experimental at the moment.
 
As is suggested, we will use a set of scales (which are expressed in units of the requested pixel, or cell, size) which are representative of the scales that are present in the data, including a zero-scale for point sources. 
 
'''Note that interrupting {{clean}} by Ctrl+C may corrupt your visibilities -- you may be better off choosing to let {{clean}} finish.  We are currently implementing a command that will nicely exit to prevent this from happening, but for the moment try to avoid Ctrl+C.'''
 
[[Image:SN_G55.MultiScale.image.png|200px|thumb|right|G55.7+3.4 Multi-Scale Clean]]
[[Image:SN_G55_MultiScale.artifacts.png|200px|thumb|right|Artifacts around point sources]]


<source lang="python">
<source lang="python">
# In CASA
# In CASA
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.MultiScale',
widebandpbcor(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.ms.MFS.wProj',  
      imsize=1280, cell='8arcsec', multiscale=[0,6,10,30,60], smallscale=0.9,
              nterms=2, action='pbcor', pbmin=0.2, spwlist=[0,1,2,3],
      interactive=False, niter=1000, pbcor=True, weighting='briggs', uvtaper='
              weightlist=[0.5,1.0,0,1.0], chanlist=[25,25,25,25], threshold='0.6mJy')
      stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean')
 
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.MultiScale',
      imsize=1280, cell='8arcsec', multiscale=[0,6,10,30,60],
      interactive=False, niter=1000,  weighting='briggs',
      stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean')
 
viewer('SN_G55_10s.MultiScale.image')
</source>
</source>


* imagename='SN_G55_10s.MultiScale': the root filename used for the various {{clean}} outputs.  These include the final image (<imagename>.image), the relative sky sensitivity over the field (<imagename>.flux), the point-spread function (also known as the dirty beam; <imagename>.psf), the clean components (<imagename>.model), and the residual image (<imagename>.residual).
Note that the spectral index image (<imagename>.alpha) will have units in Jy/beam, which is not correct. The software bug has been fixed in CASA 4.6.


* imsize=1280: the image size in number of pixels.  Note that entering a single value results in a square image with sides of this value.
* ''spwlist=[0,1,2,3]'': We want to apply this correction to all spectral windows in our calibrated measurement set.  


* cell='8arcsec': the size of one pixel; again, entering a single value will result in a square pixel size.
* ''weightlist=[0.5,1.0,0,1.0]'': Since we did not specify reference frequencies during {{clean}}, the {{widebandpbcor}} task will pick them from the provided image. Running the task, the logger reports the multiple frequencies used for the primary beam, which are 1.256, 1.429, 1.602, and 1.775 GHz. Having created an amplitude vs. frequency plot of the calibrated measurement set with colorized spectral windows using {{plotms}}, we note that the first chosen frequency lies within spectral window 0, which we know had flagged data due to RFI being present. This parameter ''weightlist'' allows us to give this chosen frequency less weight. The primary beam at 1.6 GHz lies in an area with no data, therefore we will give a weight value of zero for this frequency. The remaining frequencies 1.429 and 1.775 GHz lie within spectral windows which contained less RFI, therefore we provide a larger weight percentage.  


* multiscale=[0,6,10,30,60]: a set of scales on which to clean.  A good rule of thumb when using multiscale is [0, 2xbeam, 5xbeam] (where beam is the synthesized beam) and larger scales up to the maximum scale the interferometer can image. Since these are in units of the pixel size, our chosen values will be multiplied by the requested cell size. Thus, we are requesting scales of 0 (a point source), 48, 80, 240, and 480 arcseconds.  Note that 16 arcminutes (960 arcseconds) roughly corresponds to the size of G55.7+3.4.
* ''pbmin=0.2'': Gain level below which not to compute Taylor-coefficients or apply a primary beam correction.


*smallscale=0.9: This parameter is known as the '''small scale bias''', and helps with faint extended structure, by balancing the weight given to smaller structures which tend to be brighter, but have less flux density. Increasing this value gives more weight to smaller scales. A value of 1.0 weighs the largest scale to zero, and a value of less than 0.2 weighs all scales nearly equally. The default value is 0.6.
* ''chanlist=[25,25,25,25]'': Make primary beams at frequencies corresponding to channel 25.  


* interactive=False: we will let {{clean}} use the entire field for placing model components.  Alternatively, you could try using interactive=True, and create regions to constrain where components will be placed. However, this is a very complex field, and creating a region for every bit of diffuse emission as well as each point source can quickly become tedious. For a tutorial that covers more of an interactive clean, please see [https://casaguides.nrao.edu/index.php?title=EVLA_high_frequency_Spectral_Line_tutorial_-_IRC%2B10216 IRC+10216 tutorial.]
* ''threshold='0.6mJy''': Threshold in the intensity map, below which not to recalculate the spectral index. We adjusted the threshold to get good results: too high of a threshold resulted in errors and no spectral index image, too low of a threshold resulted in spurious values for the spectral index image.


* niter=1000: this controls the number of iterations {{clean}} will do in the minor cycle. 
<source lang="python">
 
# In CASA
* pbcor=True: We can correct for the relative sky sensitivity while running clean. This will help in creating a near circularly symmetric beam, which may be better able to image extended sources, over long observations. setting this to true, forces the raw image to be rescaled by dividing by the noise and primary beam correction image (<imagename>.flux). Note that this can also be done via the {{immath}} task if pbcor=False.
viewer('SNR_G55_10s.ms.MFS.wProj.pbcor.image.tt0')
 
</source><br />
* weighting='briggs': use Briggs weighting with a robustness parameter of 0 (halfway between uniform and natural weighting).
 
* usescratch=F: do not write the model visibilities to the model data column (only needed for self-calibration)
 
* imagermode='csclean': use the Cotton-Schwab clean algorithm
 
* stokes='I': since we have not done any polarization calibration, we only create a total-intensity image.
 
* threshold='0.1mJy': threshold at which the cleaning process will halt; i.e. no clean components with a flux less than this value will be created. This is meant to avoid cleaning what is actually noise (and creating an image with an artificially low rms).  It is advisable to set this equal to the expected rms, which can be estimated using the [http://go.nrao.edu/ect EVLA exposure calculator].  However, in our case, this is a bit difficult to do, since we have lost a hard-to-estimate amount of bandwidth due to flagging, and there is also some residual RFI present.  Therefore, we choose 0.1 mJy as a relatively conservative limit.
 
This is the fastest of the imaging techniques described here, but it's easy to see that there are artifacts in the resulting image. Note that you may have to play with the image color map and brightness/contrast to get a better view of the image details. This can be done by clicking on ''Data Display Options'' (wrench icon on top right corner), and choosing "rainbow 3" under ''basic settings''. We can use the {{viewer}} to explore the point sources near the edge of the field by zooming in on them. Some have prominent arcs, as well as spots in a six-pointed pattern surrounding them.
 
Next we will explore some more advanced imaging techniques to mitigate these artifacts.
 
=== Multi-Scale, Wide-Field Clean (w-projection) ===
 
[[Image:Faceting.png|200px|thumb|right| Faceting when using widefield gridmode, which can be used in conjunction with w-projection. ]]
[[Image:SN_G55_MS.to.MS_wProj.gif|200px|thumb|right|Multi-Scale image of arcs around point sources far from the phase center, versus MS with w-projection. We can see the that combining the w-projection algorithm with the multiscale algorithm improves the resulting image by removing prominent artifacts.]]
 
The next {{clean}} algorithm we will employ is w-projection, which is a wide-field imaging technique that takes into account the non-coplanarity of the baselines as a function of distance from the phase center. For wide-field imaging, the sky curvature and non-coplanar baselines results in a non-zero w-term. The w-term introduced by the sky and array curvature introduces a phase term that will limit the dynamic range of the resulting image. Applying 2-D imaging to such data will result in artifacts around sources away from the phase center, as we saw in running MS-CLEAN. Note that this affects mostly the lower frequency bands, especially for the more extended configurations, due to the field of view decreasing with higher frequencies.
 
The w-term can be corrected by faceting (describe the sky curvature by many smaller planes) in either the image or uv-plane, or by employing w-projection. A combination of the two can also be employed within {{clean}} by setting the parameter gridmode='widefield'. If w-projection is employed, it will be done for each facet. Note that w-projections is an order of magnitude faster than the faceting algorithm, but will require more memory.
 
For more details on w-projection, as well as the algorithm itself, see [http://adsabs.harvard.edu/abs/2008ISTSP...2..647C "The Noncoplanar Baselines Effect in Radio Interferometry: The W-Projection Algorithm"]. Also, the chapter on [http://www.aspbooks.org/a/volumes/article_details/?paper_id=17953 Imaging with Non-Coplanar Arrays] may be helpful.


<source lang="python">
<source lang="python">
# In CASA
# In CASA
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.ms.wProj',
viewer('SNR_G55_10s.ms.MFS.wProj.pbcor.image.alpha')
      gridmode='widefield', imsize=1280, cell='8arcsec',
      wprojplanes=128, multiscale=[0,6,10,30,60],
      interactive=False, niter=1000,  weighting='briggs',
      stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean')
 
viewer('SNR_G55_10s.ms.wProj.image')
</source>
</source>


* gridmode='widefield': Use the w-projection algorithm.
It's important to note that the image will cut off at about 20% of the HPBW, as we are confident of the accuracy within this percentage. Anything outside becomes less accurate, thus there is a mask associated with the creation of the corrected primary beam image (Figure 16). Spectral indices may still be unreliable below approximately 70% of the primary beam HPBW (Figure 17).  
 
* wprojplanes=128: The number of w-projection planes to use for deconvolution; 128 is the minimum recommended number.
 
This will take slightly longer than the previous imaging round; however, the resulting image has noticeably fewer artifacts. In particular, compare the same outlier source in the Multi-Scale w-projected image with the Multi-Scale-only image: note that the swept-back arcs have disappeared. There are still some obvious imaging artifacts remaining, though.


=== Multi-Scale, Multi-Frequency Synthesis ===
It would be a good exercise to use the {{viewer}} to plot both the primary beam corrected image, and the original cleaned image and compare the intensity (Jy/beam) values, which should differ slightly.
[[Image:MultiFrequency_Synthesis_snapshot.png|200px|thumb|right| Multi-Frequency Synthesis snapshot of (u,v) coverage. We can see from the image on the right, using this algorithm can greatly improve coverage, thereby improving image fidelity.]]
[[Image:SNR_G55.ms_to_mfs.gif|200px|thumb|right|Multi-Scale image artifacts versus MS-MFS artifacts near SNR, with nterms=2. We can see artifacts around point sources diminish, improving our image.]]
[[Image:SN_G55_MS.MFS.alpha.png|200px|thumb|right|Spectral Index image]]


Another consequence of simultaneously imaging the wide fractional bandwidths available with the EVLA is that the primary beam has substantial frequency-dependent variation over the observing band. If this is not accounted for, it will lead to imaging artifacts and compromise the achievable image rms.


If sources which are being imaged have intrinsically flat spectra, this will not be a problem.  However, most astronomical objects are not flat-spectrum sources, and without any estimation of the intrinsic spectral properties, the fact that the primary beam is twice as large at 2 than at 1 GHz will have substantial consequences.
== Imaging Spectral Cubes ==


Note that the dimentions of the (u,v) plane are measured in wavelengths, and therefore observing at several frequencies, a baseline can sample several ellipses in the (u,v) plane, each with different sizes. We can therefore fill in the gaps in the single frequency (u,v) coverage, hence Multi-Frequency Synthesis (MFS). Also when observing in low-frequencies, it may prove beneficial to observe in small time-chunks, which are spread out in time. This will allow the coverage of more spatial-frequencies, allowing us to employ this algorithm more efficiently.  
For spectral line imaging, CASA {{clean}} can be used in either ''mode='frequency' '' or ''mode='velocity'.'' Both modes will create a spectral axis in frequency, while the velocity mode adds an additional velocity label.  


The Multi-Scale Multi-Frequency-Synthesis (MS-MFS) algorithm provides the ability to simultaneously image and fit for the intrinsic source spectrum.  The spectrum is approximated using a polynomial in frequency, with the degree of the polynomial as a user-controlled parameter. A least-squares approach is used, along with the standard clean-type iterations. Using this method of imaging will dramatically improve our (u,v) coverage, hence improving image fidelity.
The following keywords are important for spectral modes (velocity in this example):
 
For a more detailed explanation of the MS-MFS deconvolution algorithm, please see the paper by Urvashi Rau and Tim J. Cornwell entitled [http://arxiv.org/abs/1106.2745 A multi-scale multi-frequency deconvolution algorithm for synthesis imaging in radio interferometry]


<source lang="python">
<source lang="python">
# In CASA
# In CASA
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.ms.MFS',
mode                = 'velocity'       #  Spectral gridding type (mfs, channel, velocity, frequency)
      imsize=1280, cell='8arcsec', mode='mfs', nterms=2,
    nchan          =         -1        #  Number of channels (planes) in output image; -1 = all
      multiscale=[0,6,10,30,60],
    start          =        ''       #  Velocity of first channel: e.g '0.0km/s'(''=first channel in first SpW of MS)
       interactive=False, niter=1000, weighting='briggs',
    width          =         ''       #  Channel width e.g '-1.0km/s' (''=width of first channel in first SpW of MS)
      stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean')
<snip>
    outframe       =         ''        # spectral reference frame of output image; '' =input
    veltype        =   'radio'       #  Velocity definition of output image
 
<snip>
restfreq            =         ''       #  Rest frequency to assign to image (see help)
</source>


viewer('SNR_G55_10s.ms.MFS.image.tt0')
The spectral dimension of the output cube will be defined by these parameters and {{clean}} will regrid the visibilities to the cube. Note that invoking {{cvel}} before imaging is, in most cases, not necessary even when two or more measurement sets are being provided at the same time.  
The cube is specified by a ''start'' velocity in km/s, the ''nchan'' number of channels, and a channel ''width'' (where the latter can also be negative for decreasing velocity cubes).


viewer('SNR_G55_10s.ms.MFS.image.alpha')
To correct for Doppler effects, {{clean}} also requires a ''outframe'' velocity frame where ''LSRK'' (Local Standard of Rest (kinematic)) and ''BARY'' (Sun-Earth Barycenter) are the most popular references. Parameter ''veltype'' defines whether the data will be gridded via the optical or radio velocity definition. A description of the available options and definitions can be found in the [https://science.nrao.edu/facilities/vla/docs/manuals/obsguide/modes/line VLA observing guide] and the [https://casa.nrao.edu/docs/cookbook/casa_cookbook014.html CASA cookbook]. By default, {{clean}} will produce a cube using the LRSK and radio definitions.
</source>


* nterms=2:the number of Taylor terms to be used to model the frequency dependence of the sky emission.  Note that the speed of the algorithm will depend on the value used here (more terms will be slower); of course, the image fidelity will improve with a larger number of terms (assuming the sources are sufficiently bright to be modeled more completely).
Note that {{clean}} will always work on the entire cube when searching for the highest residual fluxes in the minor cycle. To do this per channel, one can set parameter ''chaniter=T''.


This will take much longer than the two previous methods, so it would probably be a good time to have coffee or chat about EVLA data reduction with your neighbor at this point.
CASA also offers parameter ''mode='channel'.'' But we do not recommend this mode as the VLA does not observe with Doppler tracking. The fixed frequency observing of the VLA means that spectral features will drift in sky frequency over time. The Doppler correction will then be applied by {{clean}} via parameter ''mode='velocity''' or '' 'frequency' ''. Parameter ''mode='channel' '' will perform the correction only using the first timestamp in the observations. It can also be confusing when imaging some shifted, but overlapping, spectral windows.


When clean is done <imagename>.image.tt0 will contain a total intensity image, where tt0 is a suffix to indicate the Taylor term; <imagename>.image.alpha will contain an image of the spectral index in regions where there is sufficient signal-to-noise. Having this spectral index image can help convey information about the emission mechanism involved within the supernova remnant. It can also give information on the optical depth of the source. I've included a color widget on the top of the plot to give an idea of the spectral index variation.  
An example of spectral line imaging procedures is provided in the [https://casaguides.nrao.edu/index.php?title=VLA_high_frequency_Spectral_Line_tutorial_-_IRC%2B10216 VLA high frequency Spectral Line tutorial - IRC+10216].


For more information on the multi-frequency synthesis mode and its outputs, see section 5.2.5.1 in the [http://casa.nrao.edu/Doc/Cookbook/casa_cookbook.pdf CASA cookbook]. 


Inspect the brighter point sources in the field near the supernova remnant. You will notice that some of the artifacts which had been symmetric around the sources themselves are now gone; however, since we did not use W-Projection this time, there are still strong features related to the non-coplanar baseline effects still apparent for sources further away.
== Beam per Plane ==
For spectral cubes spanning relatively wide ranges in frequency, the synthesized beam can vary substantially across the channels. To account for this, CASA will calculate separate beams for each channel when the difference is more than half a pixel across the cube. All CASA image analysis tasks are capable of handling such cubes.  


=== Multi-Scale, Multi-Frequency, Widefield Clean ===
If it is desired to have a cube with a single synthesized beam, two options are available. Best is to use {{imsmooth}} with parameter ''kernel='commonbeam'.'' This task will smooth each plane to that of the lowest resolution and it will clean up all header variables such that only a single beam appears in the data. The second option is to define parameter ''resmooth=T'' directly in {{clean}}.


Finally, we will combine the W-Projection and MS-MFS algorithms to simultaneously account for both of the effects.  Be forewarned -- these imaging runs will take a while, and it's best to start them running and then move on to other things.


First, we will image the autoflagged data. Using the same parameters for the individual-algorithm images above, but combined into a single {{clean}} run, we have:
== Image Header ==


[[Image:SNR.G55.images.gif|200px|thumb|right|Here we see the differences as the images progress through the different algorithms used: MS -> MS-MFS -> MS-wProjection -> MS-MFS-wProjection.]]
The image header holds metadata associated with your CASA image. The task {{imhead}} will display this data within the casalog. We will first run {{imhead}} with parameter ''mode='summary' '':


<source lang="python">
<source lang="python">
# In CASA
# In CASA
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.ms.MFS.wProj',
imhead(imagename='SNR_G55_10s.ms.MFS.wProj.image.tt0', mode='summary')
      gridmode='widefield', imsize=1280, cell='8arcsec', mode='mfs',
      nterms=2, wprojplanes=128, multiscale=[0,6,10,30,60], 
      interactive=False, niter=1000,  weighting='briggs',
      stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean')
 
viewer('SNR_G55_10s.ms.MFS.wProj.image.tt0')
 
viewer('SNR_G55_10s.ms.MFS.wProj.image.alpha')
</source>
</source>


Again, looking at the same outlier source, we can see that the major sources of error have been removed, although there are still some residual artifacts.  One possible source of error is the time-dependent variation of the primary beam; another is the fact that we have only used nterms=2, which may not be sufficient to model the spectra of some of the point sources.
* ''mode='summary' '': gives general information about the image, including the object name, sky coordinates, image units, the telescope the data was taken with, and more.  


Ultimately, it isn't too surprising that there was still some RFI present in our auto-flagged data, since we were able to see this with {{plotms}}.  It's also possible that the auto-flagging overflagged some portions of the data, also leading to a reduction in the achievable image rms.
For further information, run parameter ''mode='list' ''. We will assign a variable'' header ''to also capture the output as a Python dictionary:
 
=== Imaging Outlier Fields ===
 
Now we will image the supernova remnant, as well as the bright sources (outliers) towards the edges of the image, creating a 512x512 pixel image for each one. For this, we will be utilizing CLEAN without widefield mode (try it to see the effects) as we will be specifying a phase center, and the images will not be too big, therefore we will not have sources very far from the phase center. We will specify a name for each image, as well as a phase center for each image. This form of imaging is useful for when you have several outlier fields you'd like to image in one go. An [https://casa.nrao.edu/Release3.3.0/docs/UserMan/UserMansu280.html outlier file] can also be created, if you have a list of sources you'd like to image.


<source lang="python">
<source lang="python">
# In CASA
# In CASA
clean(vis='SNR_G55_10s.calib.ms', imagename=['SNR.MS.MFS', 'Outlier1.MS.MFS', 'Outlier2.MS.MFS'],
header=imhead(imagename='SNR_G55_10s.ms.MFS.wProj.image.tt0', mode='list')
      imsize=[[512,512],[512,512],[512,512]], cell='8arcsec', mode='mfs',
      multiscale=[0,6,10,30,60], interactive=False, niter=1000,  weighting='briggs',
      stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean',
      phasecenter=['J2000 19h21m38.271 21d45m48.288', 'J2000 19h23m27.693 22d37m37.180', 'J2000 19h25m46.888 21d22m03.365'])
</source>
</source>


[[Image:Outlier_Fields.png|800px|thumb|center| Images of the supernova remnant and bright outlying sources. ]]
* ''mode='list' '': gives more detailed information, including beam major/minor axes, beam primary angle, and the location of the max/min intensity, and lots more. Essentially this mode displays the FITS header variables.


=== Primary Beam Correction ===


In interferomety, the images formed via deconvolution, are representations of the sky, multiplied by the primary beam response of the antenna. The primary beam can be described by a Gaussian with the size depending on the observing frequency. Correcting the primary beam can be done during clean with the '''pbcor''' parameter. It can also be done after imaging using the task [https://casa.nrao.edu/docs/TaskRef/impbcor-task.html impbcor] for regular data sets, and {{widebandpbcorr}} for those that use Taylor-term expansion (nterms > 1). The primary beam image is usually the image.flux output image given by CLEAN. The mode parameter may also be changed to divide or multiply. 
== Image Conversion ==


== Image Information ==
Finally, as an example, we convert our image from intensity to main beam brightness temperature.
 
This portions will cover topics on '''image headers''', and '''frequency reference frames'''.
We will use the standard equation
 
=== Frequency Reference Frame ===
 
The velocity within your image is calculated based on your choice of frame, velocity definition, and spectral line rest frequency. The initial frequency reference frame is initially given by the telescope, however, it can be transformed to several other frames, including:


* LSRK    - Local Standard of Rest Kinematic. Conventional LSR based on average velocity of stars in the solar neighborhood.
<math>T=1.222\times 10^{6} \frac{S}{\nu^{2} \theta_{maj} \theta_{min}} </math>
* LSRD    - Local Standard of Rest Dynamic. Velocity with respect to a frame in circular motion about the galactic center.
* BARY    - Barycentric. Referenced to JPL ephemeris DE403. Slightly different and more accurate than heliocentric.
* GEO    - Geocentric. Referenced to the Earth's center. This will just remove the observatory motion.
* TOPO    - Topocentric. Fixed observing frequency and constantly changing velocity.
* GALACTO - Galactocentric. Referenced to the dynamical center of the galaxy.
* LGROUP  - Local Group. Referenced to the mean motion of Local Group of Galaxies.
* CMB    - Cosmic Microwave Background dipole. Based on COBE measurements of dipole anisotropy.


=== Image Header ===
where the main beam brightness temperature <math>T</math> is given in K, the intensity <math>S</math> in Jy/beam, the reference frequency <math>\nu</math> in GHz, and the major and minor beam sizes <math>\theta</math> in arcseconds.


The image header holds meta data associated with your CASA image. The task {{imhead}} will display this data within the casalog. We will first run imhead with mode='summary':
For a beam of 29.30"x29.03" and a reference frequency of 1.579 GHz (as taken from the previous {{imhead}} run) we calculate the brightness temperature using {{immath}}:


<source lang="python">
<source lang="python">
# In CASA
# In CASA
imhead(imagename='SNR_G55_10s.ms.MFS.wProj.image.tt0', mode='summary')
immath(imagename='SNR_G55_10s.ms.MFS.wProj.image.tt0', mode='evalexpr',
      expr='1.222e6*IM0/1.579^2/(29.30*29.03)',
      outfile='SNR_G55_10s.ms.MFS.wProj.image.tt0-Tb')
</source>
</source>


* mode='summary': gives general information about the image, including the object name, sky coordinates, image units, the telescope the data was taken with, and more.
* ''mode='evalexpr' '': {{immath}} is used to calcuate with images
 
* ''expr'': the mathematical expression to be evaluated. The images are abbreviated as IM0, IM1, ... in the sequence of their appearance in the ''imagename'' parameter.


For further information about the image, let's now run it with mode='list':
Since {{immath}} only changes the values and not the unit of the image, we will now change the new image header'' 'bunit' ''key to 'K'. 
To do so, we will run {{imhead}} with parameter ''mode='put' '':


<source lang="python">
<source lang="python">
# In CASA
# In CASA
imhead(imagename='SNR_G55_10s.ms.MFS.wProj.image.tt0', mode='list')
imhead(imagename='SNR_G55_10s.ms.MFS.wProj.image.tt0-Tb', mode='put', hdkey='bunit', hdvalue='K')
</source>
</source>


* mode='list': gives more detailed information, including beam major/minor axes, beam primary angle, and the location of the max/min intensity, and lots more.
* ''hdkey'': the header keyword that will be changed
 
* ''hdvalue'': the new value for header keyword


We will now want to change our image header units from Jy/beam to Kelvin. To do this, we will run the task with mode='put':
Launching the viewer will now show our image in brightness temperature units:  


<source lang="python">
<source lang="python">
# In CASA
# In CASA
imhead(imagename='SNR_G55_10s.ms.MFS.wProj.image.tt0', mode='put', hdkey='bunit', hdvalue='K')
viewer('SNR_G55_10s.ms.MFS.wProj.image.tt0-Tb')
</source>
</source>


Let's also change the direction reference frame from J2000 to Galactic:
<source lang="python">
# In CASA
imhead(imagename='SNR_G55_10s.ms.MFS.wProj.image.tt0', mode='put', hdkey='equinox', hdvalue='GALACTIC')
</source>


[[Main Page | &#8629; '''CASAguides''']]
[[Main Page | &#8629; '''CASAguides''']]


-- original: Miriam Hutman <br />
<!--
--modifications: Lorant Sjouwerman (4.4.0, 2015/07/07) <br />
-- Original: Miriam Hartman <br />
--modifications: Jose Salcido (4.5.2, 2016/02/24) <br />
-- Modifications: Lorant Sjouwerman (4.4.0, 2015/07/07) <br />
{{Checked 4.5.3}}
-- Modifications: Juergen Ott (4.5.2, 2016/04/14) <br />
-- Topical Guide: Jose Salcido (4.5.2, 2016/04/18) <br />
-- Modifications: Juergen Ott (4.5.2, 2016/04/20) <br />
-- Edits: Tony Perreault (4.5.2, 2016/05/20) <br />
-->
{{Checked 4.5.2}}

Latest revision as of 15:18, 25 May 2016

Imaging

This tutorial provides guidance on imaging procedures in CASA. The tutorial covers basic continuum cleaning and the influence of image weights, as well as wide-band and wide-field imaging techniques, multi-scale clean, an outlier field setup, and primary beam correction. Spectral line imaging procedures are explained but not covered in detail. For a more thorough example of spectral line imaging, refer to the VLA high frequency Spectral Line tutorial - IRC+10216. This imaging tutorial concludes with basic image header calls and the conversion of Jy/beam to K surface brightness units through image arithmetic and header manipulation.

We will be utilizing data, taken with the Karl G. Jansky Very Large Array, of G055.7+3.4, which is a supernova remnant. The data were taken on August 23, 2010, in the first D-configuration for which the wideband capabilities of the WIDAR (Wideband Interferometric Digital ARchitecture) correlator were available. The 8 hour long observation includes all available 1 GHz of bandwidth in L-band, from 1&#150;2 GHz in frequency.

We will skip the calibration process as examples of calibration can be found in several other guides, including the VLA Continuum Tutorial 3C391 and VLA high frequency Spectral Line tutorial - IRC+10216 tutorials.

A copy of the calibrated data (1.2GB) can be downloaded from http://casa.nrao.edu/Data/EVLA/SNRG55/SNR_G55_10s.calib.tar.gz 

Your first step will be to unzip and untar the file in a terminal window (before you start CASA):

tar -xzvf  SNR_G55_10s.calib.tar.gz

Then start casa as usual via the casa command, which will bring up the ipython interface and launches the logger.


The CLEAN Algorithm

Figure 1
The CLEAN major and minor cycles, indicating the steps undertaken during gridding, projection algorithms, and creation of images.

The CLEAN algorithm, developed by J. Högbom (1974) enabled the synthesis of complex objects, even if they have relatively poor Fourier uv-plane coverage. Poor coverage occurs with partial earth rotation synthesis or with arrays composed of few antennas. The dirty image is formed by a simple Fourier inversion of the sampled visibility data, with each point on the sky being represented by a suitably scaled and centered PSF (Point Spread Function), or dirty beam, which itself is the Fourier inversion of the visibility (u,v) coverage.

The convolution with the dirty beam creates artifacts in the image and limits the dynamic range. The CLEAN algorithm attempts to remove the dirty beam pattern from the image via deconvolution. This implies that it interpolates from the measured (u,v) points across gaps in the (u,v) coverage. In short, CLEAN provides solutions to the convolution equation by representing radio sources by a number of point sources in an empty field. The brightest points are found by performing a cross-correlation between the dirty image and the PSF. The brightest parts are then subtracted, and the process is repeated again for the next brighter sources. Variants of CLEAN, such as multi-scale CLEAN, take into account extended kernels which may be better suited for extended objects.

For single pointings, CASA uses the Cotton-Schwab cleaning algorithm in the task clean (parameter imagermode='csclean'), which breaks the process into major and minor cycles (see Figure 1). To start with, the visibilities are gridded, weighted, and Fourier transformed to create a dirty image. The minor cycle then operates in the image domain to find the clean components that are added to the clean model. The model image is Fourier transformed back to the visibility domain, degridded, and subtracted from the visibilities. This creates a new residual that is then gridded, weighted, and FFT'ed again to the image domain for the next iteration. The gridding, FFT, degridding, and subtraction processes form the major cycle.

This iterative process is continued until a stopping criterion is reached, such as a maximum number of clean components or a flux threshold in the residual image.

In CASA clean, two versions of the PSF can be used: parameter psfmode = 'hogbom' uses the full-sized PSF for subtraction and is a thorough but slow method. Parameter psfmode='clark' uses a smaller beam patch, which increases the speed. The patch size and length of the minor cycle are internally chosen such that clean converges well without giving up the speed improvement. Parameter psfmode='clark' is the default option in clean.

In a final step, clean derives a Gaussian fit to the inner part of the dirty beam, which defines the clean beam. The clean model is then convolved with the clean beam and added to the last residual image to create the final image.

Note that the CASA team is currently developing a refactored clean task, called tclean. Task tclean has a better interface, provides new algorithms, and more combinations between imaging algorithms. Task tclean also includes software to parallelize the computations in a multi-processor environment. Eventually, tclean will replace the current clean task. For this guide, however, we will stick with the original clean task, as tclean is still in the development and testing phase. Nevertheless, the reader is encouraged to try tclean and send us feedback through the NRAO helpdesk.

For more details on imaging and deconvolution, we refer to the Astronomical Society of the Pacific Conference Series book entitled Synthesis Imaging in Radio Astronomy II; the chapter on Deconvolution may prove helpful. Additionally, imaging presentations are available on the Synthesis Imaging Workshop and VLA Data Reduction Workshop webpages. The CASA cookbook chapter on Synthesis Imaging provides a wealth of information on the CASA implementation of clean and related tasks.

Finally, we refer users to the VLA Observational Status Summary and the Guide to Observing with the VLA for information on the VLA capabilities and observing strategies.


Weights and Tapering

Figure 2
u,v coverage for the 8 hour observation of the supernova remnant G055.7+3.4

When imaging data, a map is created associating the visibilities with the image. The sampling function is modified by a weight function that defines the shape and size of the PSF. Weighting therefore provides some control over the spatial resolution and the surface brightness sensitivity of the map, where either direction can be emphasized.

There are three main weighting schemes (see Table 1):

1) Natural weighting: This weighs data based on their rms only. Data visibility weights are gridded onto a uv-cell and summed. More visibilities in a cell will thus increase the cell's weight, which will usually emphasize the shorter baselines. Natural weighting therefore results in a larger PSF and better surface brightness sensitivity, but also a degraded resolution.

2) Uniform weighting: The weights are first gridded as in natural weighting but then each cell is corrected such that the weights are independent of the number of visibilities inside. Compared to natural weighting, uniform weighting usually emphasizes the longer baselines. Consequently the PSF is smaller, which results in a better spatial resolution of the image. At the same time, however, the surface brightness sensitivity is reduced compared to natural weighting. The uniform weighting of the baselines is a better representation of the uv-coverage and sidelobes are more suppressed.

3) Briggs weighting: This scheme provides a compromise between natural and uniform weighting. The transition can be controlled with the robust parameter where robust=-2 is close to uniform and robust=2 is close to natural weighting. Briggs weighting therefore offers a compromise for between spatial resolution and surface brightness sensitivity, and robust values near zero are typically used.

Details on the weighting schemes are given in Daniel Brigg's dissertation (Chapter 3).

For a visual comparison between these three weighting schemes, please see the section on "CLEAN with Weights" in this guide.

Table 1
Table summarizing the effects of using weights and tapering.

Tapering: In conjunction with the above weighting schemes, one can specify the uvtaper parameter within clean, which will control the radial weighting of visibilities in the uv-plane. Figure 2 graphically visualizes the uv-coverage during the observing session. The taper in clean is an elliptical Gaussian function which effectively removes long baselines and degrades the resolution. This may be desirable for extended structures when long baselines contribute a large fraction of the noise. Tapering can, therefore, increase the surface brightness sensitivity of the data, but will decrease the point source sensitivity. Too aggressive tapering, however, may also take its toll on the surface brightness sensitivity.

We refer to the CASA Cookbook Synthesis Imaging chapter for the details of the weighting implementation in CASA's clean.


Primary and Synthesized Beam

The primary beam of a single antenna defines the sensitivity across the field of view. For the VLA antennas, the main beam can be approximated by a Gaussian with a FWHM equal to [math]\displaystyle{ 90*\lambda_{cm} }[/math] arcseconds, or [math]\displaystyle{ 45/ \nu_{GHz} }[/math] arcminutes. Note that there are sidelobes beyond the Gaussian kernel that are sensitive to bright sources (see below). Taking our observed frequency to be the middle of L-band, 1.5 GHz, our primary beam will be around 30 arcmin.

If your science goal is to image a source or field of view that is significantly larger than the FWHM of the VLA primary beam, then creating a mosaic from a number of pointings is usually the preferred method. For a tutorial on mosaicing, see the 3C391 tutorial. In the following, however, we will discuss methods to image large areas from single pointing data.

Since our data were taken in the D-configuration, we can check the Observational Status Summary section on VLA resolution to find that the synthesized beam will be around 46 arcsec. Flagging, weighting and the exact frequency, however, may result in deviations from this value. As we will see later, the synthesized beam of our data hovers around 29 arcsec and for the extreme of uniform weighting around 26"x25". Oversampling the minor axis by a factor of ~3, we will use a cell (pixel) size of 8 arcsec.

Our field contains bright point sources significantly outside the primary beam. The VLA, in particular when using multi-frequency synthesis (see below), will have significant sensitivity outside the main lobe of the primary beam. Particularly at the lower VLA frequencies, sources that are located outside the primary beam may still be bright enough to be detected with sufficient signal-to-noise in the primary beam sidelobes; which can cause artifacts to interfer with the targeted field in the main part of the primary beam. Such sources need to be cleaned to remove the dirty beam response of these interfering sources from the entire image. This can be done by either creating a very large image or by using outlier fields centered on the strongest sources (see section on outlier fields below). A large image has the added advantage of increasing the field of view for science, albeit at lower sensitivity. Other effects will start to become significant, however, like the non-coplanarity of the sky. Large image sizes will also slow down the deconvolution process. Details are provided in Sanjay Bhatnagar's presentation: "Advanced Imaging: Imaging in the Wide-band Wide-field era" given at the 2016 VLA data reduction workshop.

The calls to CLEAN within this guide will create images that are 170 arcminutes on a side ([Image Size * Cell Size]*[1arcmin / 60arsec]), or almost 6x the size of the primary beam, catching the first and second sidelobes. This is ideal for showcasing both the problems inherent in such wide-band, wide-field imaging, as well as some of the solutions currently available in CASA to deal with these issues.

Note that the execution time of clean depends on the image sizes. Large images generally take more computing time. There are some values, however, that are computationally not advisable. The logger output will then show a recommendation for the next larger but faster image size. As a rule of thumb we recommend image sizes 2n * 10 pixels (e.g., 160, 1280 pixels, etc.) for improved processing speeds.


Clean Output Images

As a result of the CLEAN algorithm, clean will create a number of output images. For an imagename='<imagename>', these outputs are:

<imagename>.residual: the residual after subtracting the clean model (unit: Jy/beam, where beam refers to the dirty beam).
<imagename>.model: the clean model, not convolved (unit: Jy/pixel).
<imagename>.psf: the point-spread function (dirty beam)
<imagename>.flux: the normalized sensitivity map. For single pointings this corresponds to the primary beam.
<imagename>.image: the residual + the model convolved with the clean beam. This is the final image (unit: Jy/beam, where beam refers to the clean beam).

Additional images will be created for specific algorithms like multi-term frequency synthesis or mosaicking.

Important: If an image file is present in the working directory and the same name is provided in imagename, clean will use that image&#151;in particular the residual and model image&#151;as a starting point for further cleaning. If you want a fresh run of clean, first remove all images of that name using 'rmtables()':

# In CASA
rmtables('<imagename>.*')

This method is preferable over rm -rf as it also clears the cache.

Note that interrupting clean by Ctrl+C may corrupt your visibilities&#151;you may be better off choosing to let clean finish. We are currently implementing a command that will nicely exit to prevent this from happening but, for the moment, please try to avoid Ctrl+C.


Dirty Image

First, we will create a dirty image (Figure 3) to see the improvements as we step through several cleaning algorithms and parameters. The dirty image is the true image on the sky, convolved with the dirty beam (PSF). We will do this by running clean with parameter niter=0, which will not do any cleaning.

# In CASA
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.dirty', 
      imsize=1280, cell='8arcsec', interactive=False, niter=0,  
      stokes='I', usescratch=F)
  • imagename='SNR_G55_10s.dirty': the root filename used for the various clean outputs.
  • imsize=1280: the image size in number of pixels. A single value will result in a square image.
  • cell='8arcsec': the size of one pixel; again, entering a single value will result in a square pixel size.
  • niter=0: this controls the number of iterations done in the minor cycle.
  • usescratch=F: controls writing the model visibilities to the model data column. For self-calibration we currently recommend setting usescratch=T.
  • stokes='I': since we have not done any polarization calibration, we only create a total-intensity image. For using CLEAN while including various Stoke's Parameters, please see the 3C391 CASA guide.
# In CASA
viewer('SNR_G55_10s.dirty.image')


# In CASA
viewer('SNR_G55_10s.dirty.psf')
Figure 3
A dirty image of the supernova remnant G55.7+3.4 in greyscale, with apparent sidelobes (left), and the point spread function (PSF), also known as the dirty beam (right).

The images may be easier to see in grey scale. To do this, click on Data Display Options (wrench icon, upper left corner) within the viewer and change the color map to Greyscale 1. You may also have to change the scaling power options to your liking. To change the brightness and contrast, assign a mouse button to this type of editing by clicking on the Colormap fiddling button (black/white circle icon) and click/drag the mouse over the image to change the brightness (left-right mouse movement) and contrast (up-down mouse movement).

Note that the clean beam is only defined after some clean iterations. The dirty image has therefore no beam size specified in the header and the PSF image is the representation of the response of the array to a point source. Even though it is empty with no iterations specified, clean will still produce the model file. So a continuation into actual cleaning will be possible by just restarting clean with the same image root name (and parameter niter>0).


Regular CLEAN & RMS Noise

We will now create a regular clean image using mostly default values to see how deconvolution improves the image quality. The first run of clean will use a fixed number of minor cycle iterations of niter=1000 (default is 500), the second will have niter=10000. Note that you may have to play with the image color map and brightness/contrast to get a better view of the image details.

# In CASA. Create default clean image.
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.Reg.Clean.niter1K', 
      imsize=1280, cell='8arcsec', niter=1000, interactive=False)


# In CASA.
viewer('SNR_G55_10s.Reg.Clean.niter1K.image')

The logger indicates that the image was obtained in two major cycles and some improvements over the dirty image are visible. But clearly we have not cleaned deep enough yet; the image still has many sidelobes, and an inspection of the residual image shows that it still contains source flux and structure. So let's increase the niter parameter value to 10,000 and compare the images (Figure 4).

# In CASA. Create default clean image with niter = 10000
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.Reg.Clean.niter10K', 
      imsize=1280, cell='8arcsec', niter=10000, interactive=False)


# In CASA.
viewer('SNR_G55_10s.Reg.Clean.niter10K.image')
Figure 4
Regular run of CLEAN, with niter=1000 (left), compared to niter=10000 (right).

As we can see from the resulting images, increasing the niter values (minor cycles) improves our image by reducing prominent sidelobes significantly. One could now further increase the niter parameter until the residuals are down to an acceptable level. To determine the number of iterations, one needs to keep in mind that clean will fail once it starts cleaning too deeply into the noise. At that point, the cleaned flux and the peak residual flux values will start to oscillate as the number of iterations increase. This effect can be monitored on the CASA logger output. To avoid cleaning too deeply, we will set a threshold parameter that will stop minor cycle clean iterations once a peak residual value is being reached.

First, we will utilize the SNR_G55_10s.Reg.Clean.niter1K.image image to give us an idea of the rms noise (your sigma value)(Figure 5). With this image open within the viewer, click on the Rectangle Drawing button (Rectangle with R) and draw a square on the image at a position with little source or sidelobe contamination. Doing this should open up a Regions dock which holds information about the selected region, including the pixel statistics in the Statistics tab (double clicking on the box will also bring up this information). Take notice of the rms values as you click/drag the box around empty image locations, or by drawing additional boxes at suitable positions. If the Regions dock is not displayed, click on View in the menu bar at the top of the viewer and click on the Regions check-box.

Figure 5
Attempting to find the lowest rms value within the CLEAN'ed image using niter=1000, in order to calculate our threshold.

The lowest rms value we found was in the order of about 4 x 10-5 Jy/beam, which we will use to calculate our threshold. There really is no set standard, but fairly good threshold values can vary anywhere between 2.0&#150;4.0*sigma; using clean boxes (see the section on interactive cleaning) allows one to go to lower thresholds. For our purposes, we will choose a threshold of 2.5*sigma. Doing the math results in a value of 10 x 10-5 or, equivalently, 0.10mJy/beam. Therefore, for future calls to the clean task, we will set parameter threshold='0.1mJy'. The clean cycle will be stopped when the residual peak flux density is either equal to or less than the threshold value or when the maximum number of iterations niter is reached. To ensure that the stopping criterion is indeed the parameter threshold, parameter niter should be set to a very high number.

In the following, we nevertheless will use niter=1000 to keep the execution times of clean on the low end as we focus on explaining different imaging methods.

An alternative method to determine the approximate rms of an image is to use the VLA Exposure Calculator and to enter the observing conditions. Make sure that the chosen bandwidth reflects the data after RFI excision.


CLEAN with Weights

To see the effects of using different weighting schemes to the image, let's change the weighting parameter within clean and inspect the resulting images. We will be using the Natural, Uniform, and Briggs weighting algorithms. Here, we have chosen a smaller image size to mainly focus on our science target.

# In CASA. Natural weighting
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.natural', weighting='natural',
      imsize=540, cell='8arcsec', niter=1000, interactive=False, threshold='0.1mJy')
  • weighting: specification of the weighting scheme. For Briggs weighting, the robust parameter will be used.
  • threshold='0.1mJy': threshold at which the cleaning process will halt.
# In CASA. Uniform weighting
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.uniform', weighting='uniform',
      imsize=540, cell='8arcsec', niter=1000, interactive=False, threshold='0.1mJy')


# In CASA. Briggs weighting, with robust set to default of 0.0
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.briggs', weighting='briggs',
      robust=0, imsize=540, cell='8arcsec', niter=1000, interactive=False, threshold='0.1mJy')


# In CASA. Open the viewer and select the created images.
viewer()
Figure 6
CLEAN images created with different weighting algorithms, including natural (left), uniform (center), and briggs (right).

Figure 6 shows that the natural weighted image is most sensitive to extended emission (beam size of 46"x41"). The negative values around the extended emission, often referred to as a negative bowl, is a typical signature of unsampled visibilities near the origin of the uv-plane. That is, the flux density present at shorter baselines, or larger angular scales than measured in this observation, is not well represented. Uniform weighted data shows the highest resolution (26"x25") and Briggs parameter robust=0 (default value) is a compromise with a beam of 29"x29". To image more of the extended emission, the robust parameter could be tweaked further toward more positive values.


Multi-Scale CLEAN

Since G55.7+3.4 is an extended source with many angular scales, a more advanced form of imaging involves the use of multiple scales. MS-CLEAN is an extension of the classical CLEAN algorithm for handling extended sources. Multi-Scale CLEAN works by assuming the sky is composed of emission at different angular scales and works on them simultaneously, thereby creating a linear combination of images at different angular scales. For a more detailed description of Multi-Scale CLEAN, see T.J. Cornwell's paper Multi-Scale CLEAN deconvolution of radio synthesis images.

We will use a set of scales&#151;expressed in units of the requested pixel or cell size&#151;which are representative of the scales that are present in the data, including a zero-scale for point sources.

# In CASA
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.MultiScale', 
      imsize=1280, cell='8arcsec', multiscale=[0,6,10,30,60], smallscalebias=0.9, 
      interactive=False, niter=1000,  weighting='briggs', stokes='I', 
      threshold='0.1mJy', usescratch=F, imagermode='csclean')
  • multiscale=[0,6,10,30,60]: a set of scales on which to clean. A good rule of thumb when using multiscale is 0, 2xbeam, 5xbeam (where beam is the synthesized beam), and larger scales up to about half the minor axis maximum scale of the mapped structure. Since these are in units of the pixel size, our chosen values will be multiplied by the requested cell size. Thus, we are requesting scales of 0 (a point source), 48, 80, 240, and 480 arcseconds (8 arcminutes). Note that 16 arcminutes (960 arcseconds) roughly corresponds to the size of G55.7+3.4.
  • smallscalebias=0.9: This parameter is known as the small scale bias, and helps with faint extended structure, by balancing the weight given to smaller structures, which tend to have higher surface brightness but lower integrated flux density. Increasing this value gives more weight to smaller scales. A value of 1.0 weights the largest scale to zero, and a value of less than 0.2 weighs all scales nearly equally. The default value is 0.6.
# In CASA
viewer('SNR_G55_10s.MultiScale.image')

The logger will show how much cleaning is performed on the individual scales.

This is the fastest of the imaging techniques described here, but it's easy to see that there are artifacts in the resulting image (Figure 7). We can use the viewer to explore the point sources near the edge of the field by zooming in on them (Figure 8). Click on the Zoom button on the upper left corner and highlight an area by making a square around the portion where you would like to zoom in. Double clicking within the square will zoom in the selected area. The square can be resized by clicking/dragging the corners, or removed by pressing the Esc key. After zooming in, we can see some radio sources have prominent arcs as well as spots with a six-pointed pattern surrounding them.

Figure 7
G55.7+3.4 Multi-Scale Clean
Figure 8
Artifacts around point sources

Next we will explore some more advanced imaging techniques to mitigate the artifacts seen towards the edges of the image.


Multi-Scale, Wide-Field CLEAN (w-projection)

Figure 9
Faceting when using wide-field gridmode, which can be used in conjunction with w-projection.
Figure 10
Multi-Scale image of arcs around point sources far from the phase center, versus MS with w-projection (right). We can see the that combining the w-projection algorithm with the multiscale algorithm improves the resulting image by removing prominent artifacts.

The next clean algorithm we will employ is w-projection, which is a wide-field imaging technique that takes into account the non-coplanarity of the baselines as a function of distance from the phase center (Figure 9). For wide-field imaging, the sky curvature and non-coplanar baselines result in a non-zero w-term. The w-term, introduced by the sky and array non-coplanarity, introduces a phase term that will limit the dynamic range of the resulting image. Applying 2-D imaging to such data will result in artifacts around sources away from the phase center, as we saw in running MS-CLEAN. Note that this affects mostly the lower frequency bands where the field of view is large.

The w-term can be corrected by faceting (which describes the sky curvature by many smaller planes) in either the image or uv-plane. The latter is known as w-projection. A combination of the two can also be employed within clean by setting the parameter gridmode='widefield'. If w-projection is employed, it will be done for each facet. Note that w-projection is an order of magnitude faster than the image-plane based faceting algorithm, but will require more memory.

For more details on w-projection, as well as the algorithm itself, see "The Noncoplanar Baselines Effect in Radio Interferometry: The w-projection Algorithm". Also, the chapter on Imaging with Non-Coplanar Arrays may be helpful.

# In CASA
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.ms.wProj',
      gridmode='widefield', imsize=1280, cell='8arcsec',
      wprojplanes=-1, multiscale=[0,6,10,30,60], 
      interactive=False, niter=1000,  weighting='briggs',
      stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean')
  • gridmode='widefield': Use the w-projection algorithm.
  • wprojplanes=-1: The number of w-projection planes to use for deconvolution. Setting to -1 forces CLEAN to utilize an acceptable number of planes for the given data.
# In CASA
viewer('SNR_G55_10s.ms.wProj.image')

This will take slightly longer than the previous imaging round; however, the resulting image (Figure 10) has noticeably fewer artifacts. In particular, compare the same outlier source in the Multi-Scale w-projected image with the Multi-Scale-only image: note that the swept-back arcs have disappeared. There are still some obvious imaging artifacts remaining, though.


Multi-Scale Multi-Term Frequency Synthesis

Figure 11
Multi-Frequency Synthesis snapshot of (u,v) coverage. We can see from the image on the right, using this algorithm can greatly improve coverage, thereby improving image fidelity.
Figure 12
Multi-Scale imaged with MS-MFS and nterms=2. Artifacts around point sources close to phase center are reduced while point sources away from the phase center still show artifacts.
Figure 13
Spectral Index image

A consequence of simultaneously imaging the wide fractional bandwidths available with the VLA is that the primary and synthesized beams have substantial frequency-dependent variation over the observing band. If this variation is not accounted for, it will lead to imaging artifacts and compromise the achievable image rms.

The coordinates of the (u,v) plane are measured in wavelengths. Observing at several frequencies results in a baseline sampling several ellipses of different sizes in the (u,v) plane. By using a method called Multi-Frequency Synthesis (MFS), which is the default cleaning mode in CLEAN, we can fill in the gaps in the single frequency (u,v) coverage to achieve much better image fidelity (Figure 11).

When observing in low-frequencies, it may prove beneficial to observe in short segments which are spread out in time. This will allow the coverage of more spatial-frequencies, permitting us to employ this algorithm more efficiently.

The Multi-Scale Multi-Term Frequency Synthesis algorithm provides the ability to simultaneously image and fit for the intrinsic source spectrum in each pixel. The spectrum is approximated using a polynomial Taylor-term expansion in frequency, with the degree of the polynomial as a user-controlled parameter. A least-squares approach is used, along with the standard clean-type iterations.

For a more detailed explanation of the MS-MFS deconvolution algorithm, please see the paper A multi-scale multi-frequency deconvolution algorithm for synthesis imaging in radio interferometry

# In CASA
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.ms.MFS',
      imsize=1280, cell='8arcsec', mode='mfs', nterms=2,
      multiscale=[0,6,10,30,60], 
      interactive=False, niter=1000,  weighting='briggs',
      stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean')
  • nterms=2: the number of Taylor-terms to be used to model the frequency dependence of the sky emission. Note that the speed of the algorithm will depend on the value used here (more terms will be slower); nterms=2 will fit a spectral index, and nterms=3 a spectral index and curvature.

When clean is done <imagename>.image.tt0 will contain a total intensity image (Figure 12), where tt* is a suffix to indicate the Taylor term: <imagename>.image.tt0 is the total intensity image and <imagename>.image.alpha will contain an image of the spectral index in regions where there is sufficient signal-to-noise (Figure 13). Having this spectral index image can help convey information about the emission mechanism involved within the supernova remnant. This spectral index can also give information on the optical depth of the source.

# In CASA
viewer('SNR_G55_10s.ms.MFS.image.tt0')


# In CASA
viewer('SNR_G55_10s.ms.MFS.image.alpha')

Note: To replicate Figure 13, open the image within the viewer, click on Panel Display Options (wrench with a small P), change background color to white, and adjust your margins under the Margins button. To view the color wedge, click on the Data Display Options (wrench without the P, next to the folder icon), click on the Color Wedge button, click on Yes under display color wedge and adjust the various parameters to your liking.

For more information on the multi-term, multi-frequency synthesis mode and its outputs, see section 5.2.5.1 in the CASA cookbook.

Inspect the brighter point sources in the field near the supernova remnant. You will notice that some of the artifacts, which had been symmetric around the sources themselves, are now gone. Since we did not use w-projection, however, there are still strong features related to the non-coplanar baseline effects still apparent for sources further away.

At this point, clean takes into account the frequency variation of the synthesized beam but not the frequency variation of the primary beam. For low frequencies and large bandwidths, this can be substantial, e.g., 1&#150;2 GHz L-band observations result in a variation of a factor of 2. One effect of such a large fractional bandwidth is that in multi-frequency synthesis, primary beam nulls will be blurred and the interferometer is sensitive everywhere in the field of view. For spectral slopes, however, a frequency-dependent primary beam causes a steepening as the higher frequencies are less sensitive at a given point away from the phase center. A correction for this effect should be made with the task widebandpbcor (see section on Primary Beam Correction below).


Multi-Scale Multi-Term Frequency, Widefield CLEAN

We will now combine the w-projection and MS-MFS algorithms. Be forewarned&#151;these imaging runs will take a while, and it's best to start them running and then move on to other things.

Using the same parameters for the individual-algorithm images above, but combined into a single clean run, we have:

# In CASA
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.ms.MFS.wProj',
      gridmode='widefield', imsize=1280, cell='8arcsec', mode='mfs',
      nterms=2, wprojplanes=-1, multiscale=[0,6,10,30,60],  
      interactive=False, niter=1000,  weighting='briggs',
      stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean')


# In CASA
viewer('SNR_G55_10s.ms.MFS.wProj.image.tt0')


# In CASA
viewer('SNR_G55_10s.ms.MFS.wProj.image.alpha')
Figure 14 Here we see the differences as the images progress through the different algorithms used: MS → MS-MFS → MS-wProjection → MS-MFS-wProjection.

Again, looking at the same outlier source, we can see that the major sources of error have been removed, although there are still some residual artifacts (Figure 14). One possible source of error is the time-dependent variation of the primary beam; another is the fact that we have only used nterms=2, which may not be sufficient to model the spectra of some of the point sources. Some weak RFI may also show up that may need additional flagging.


Imaging Outlier Fields

When strong sources are far from the main target, but still strong enough to produce sidelobes in the main image, they should be cleaned. Sometimes, however, it is not practical to image very large fields for this purpose. An alternative is to use outlier fields. This mode will allow the user to specify a few locations that are then cleaned along with the main image. The outlier fields should be centered on strong sources that are known from sky catalogs or are identified by other means.

In order to find outlying sources, it will help to image a large field with lower resolution and identify them within the CASA viewer. By moving the mouse cursor over the sources, we can grab their positions and use the information within an outlier file. In this example, we've called our file outliers.txt. It contains the names, sizes, and positions of the outliers we want to image (see section 5.3.18.1 of the CASA cookbook).

Open your favorite text editor and input the following:

#content of outliers.txt
#
#outlier field1
imagename= 'Outlier1.MS.MFS'
imsize=[512,512]
phasecenter = 'J2000 19h23m27.693 22d37m37.180'
#
#outlier field2
imagename='Outlier2.MS.MFS'
imsize=[512,512]
phasecenter = 'J2000 19h25m46.888 21d22m03.365'

clean will then be executed like the following:

# In CASA
clean(vis='SNR_G55_10s.calib.ms', imagename='SNR.MS.MFS-Main',
      outlierfile='outliers.txt',
      imsize=[512,512], cell='8arcsec', mode='mfs',
      multiscale=[0,6,10,30,60], interactive=False, niter=1000,  weighting='briggs',
      robust=0, stokes='I', threshold='0.1mJy', usescratch=F, imagermode='csclean')
  • outlierfile='outliers.txt': the name of the outlier file.

We can view the images with viewer one at a time. Due to the difference in sky coordinates, they cannot be viewed on the same window display.

# In CASA
viewer('SNR.MS.MFS-Main.image')


# In CASA
viewer('Outlier1.MS.MFS.image')


# In CASA
viewer('Outlier2.MS.MFS.image')

We have changed the color map to Hot Metal 2 in order to show the different colors that can be used within the viewer. Feel free to play with other color maps that may be better suited for your screen.

Figure 15
Images of the supernova remnant and bright sources far from the phase center (near the edge of, or beyond, the primary beam.)


Primary Beam Correction

In interferometry, the images formed via deconvolution are representations of the sky multiplied by the primary beam response of the antenna. The primary beam can be described as Gaussian with a size depending on the observing frequency. Images produced via clean are by default not corrected for the primary beam pattern (important for mosaics), and therefore do not have the correct flux density away from the phase center.

Correcting for the primary beam, however, can be done during clean by using the parameter pbcor. Alternatively, it can be done after imaging using the task impbcor for regular data sets and widebandpbcor for those that use Taylor-term expansion (nterms > 1). A third alternative is utilizing the task immath to manually divide the <imagename>.image by the <imagename>.flux image (<imagename>.flux.pbcoverage for mosaics).

Flux corrected images usually don't look pretty due to the noise at the edges being increased. Flux densities, however, should only be calculated from primary beam corrected images. Let's run the impbcor task to correct our multiscale image.

# In CASA
impbcor(imagename='SNR_G55_10s.MultiScale.image', pbimage='SNR_G55_10s.MultiScale.flux', 
        outfile='SNR_G55_10s.MS.pbcorr.image')
  • imagename: the image to be corrected
  • pbimage: the <imagename>.flux image as a representation of the primary beam (<imagename>.flux.pbcoverage for mosaics)
# In CASA
viewer('SNR_G55_10s.MS.pbcorr.image')

Let us now use the widebandpbcor task for wide-band (nterms>1) images. Note that for this task, we will be supplying the image name that is the prefix for the Taylor expansion images, tt0 and tt1, which must be on disk. Such files were created during the last Multi-Scale, Multi-Frequency, Widefield run of CLEAN. Task widebandpbcor will generate a set of images with a pbcor.image extension.

Figure 16
Primary beam corrected image using the widebandpbcor task for the MS.MFS.wProj image created during the CLEAN process.
Figure 17
Primary beam corrected spectral index image using the widebandpbcor task.
# In CASA
widebandpbcor(vis='SNR_G55_10s.calib.ms', imagename='SNR_G55_10s.ms.MFS.wProj', 
              nterms=2, action='pbcor', pbmin=0.2, spwlist=[0,1,2,3],
              weightlist=[0.5,1.0,0,1.0], chanlist=[25,25,25,25], threshold='0.6mJy')

Note that the spectral index image (<imagename>.alpha) will have units in Jy/beam, which is not correct. The software bug has been fixed in CASA 4.6.

  • spwlist=[0,1,2,3]: We want to apply this correction to all spectral windows in our calibrated measurement set.
  • weightlist=[0.5,1.0,0,1.0]: Since we did not specify reference frequencies during clean, the widebandpbcor task will pick them from the provided image. Running the task, the logger reports the multiple frequencies used for the primary beam, which are 1.256, 1.429, 1.602, and 1.775 GHz. Having created an amplitude vs. frequency plot of the calibrated measurement set with colorized spectral windows using plotms, we note that the first chosen frequency lies within spectral window 0, which we know had flagged data due to RFI being present. This parameter weightlist allows us to give this chosen frequency less weight. The primary beam at 1.6 GHz lies in an area with no data, therefore we will give a weight value of zero for this frequency. The remaining frequencies 1.429 and 1.775 GHz lie within spectral windows which contained less RFI, therefore we provide a larger weight percentage.
  • pbmin=0.2: Gain level below which not to compute Taylor-coefficients or apply a primary beam correction.
  • chanlist=[25,25,25,25]: Make primary beams at frequencies corresponding to channel 25.
  • threshold='0.6mJy': Threshold in the intensity map, below which not to recalculate the spectral index. We adjusted the threshold to get good results: too high of a threshold resulted in errors and no spectral index image, too low of a threshold resulted in spurious values for the spectral index image.
# In CASA
viewer('SNR_G55_10s.ms.MFS.wProj.pbcor.image.tt0')


# In CASA
viewer('SNR_G55_10s.ms.MFS.wProj.pbcor.image.alpha')

It's important to note that the image will cut off at about 20% of the HPBW, as we are confident of the accuracy within this percentage. Anything outside becomes less accurate, thus there is a mask associated with the creation of the corrected primary beam image (Figure 16). Spectral indices may still be unreliable below approximately 70% of the primary beam HPBW (Figure 17).

It would be a good exercise to use the viewer to plot both the primary beam corrected image, and the original cleaned image and compare the intensity (Jy/beam) values, which should differ slightly.


Imaging Spectral Cubes

For spectral line imaging, CASA clean can be used in either mode='frequency' or mode='velocity'. Both modes will create a spectral axis in frequency, while the velocity mode adds an additional velocity label.

The following keywords are important for spectral modes (velocity in this example):

# In CASA
mode                = 'velocity'        #  Spectral gridding type (mfs, channel, velocity, frequency)
     nchan          =         -1        #  Number of channels (planes) in output image; -1 = all
     start          =         ''        #  Velocity of first channel: e.g '0.0km/s'(''=first channel in first SpW of MS)
     width          =         ''        #  Channel width e.g '-1.0km/s' (''=width of first channel in first SpW of MS)
<snip>
     outframe       =         ''        #  spectral reference frame of output image; '' =input
     veltype        =    'radio'        #  Velocity definition of output image

<snip>
restfreq            =         ''        #  Rest frequency to assign to image (see help)

The spectral dimension of the output cube will be defined by these parameters and clean will regrid the visibilities to the cube. Note that invoking cvel before imaging is, in most cases, not necessary even when two or more measurement sets are being provided at the same time.

The cube is specified by a start velocity in km/s, the nchan number of channels, and a channel width (where the latter can also be negative for decreasing velocity cubes).

To correct for Doppler effects, clean also requires a outframe velocity frame where LSRK (Local Standard of Rest (kinematic)) and BARY (Sun-Earth Barycenter) are the most popular references. Parameter veltype defines whether the data will be gridded via the optical or radio velocity definition. A description of the available options and definitions can be found in the VLA observing guide and the CASA cookbook. By default, clean will produce a cube using the LRSK and radio definitions.

Note that clean will always work on the entire cube when searching for the highest residual fluxes in the minor cycle. To do this per channel, one can set parameter chaniter=T.

CASA also offers parameter mode='channel'. But we do not recommend this mode as the VLA does not observe with Doppler tracking. The fixed frequency observing of the VLA means that spectral features will drift in sky frequency over time. The Doppler correction will then be applied by clean via parameter mode='velocity' or 'frequency' . Parameter mode='channel' will perform the correction only using the first timestamp in the observations. It can also be confusing when imaging some shifted, but overlapping, spectral windows.

An example of spectral line imaging procedures is provided in the VLA high frequency Spectral Line tutorial - IRC+10216.


Beam per Plane

For spectral cubes spanning relatively wide ranges in frequency, the synthesized beam can vary substantially across the channels. To account for this, CASA will calculate separate beams for each channel when the difference is more than half a pixel across the cube. All CASA image analysis tasks are capable of handling such cubes.

If it is desired to have a cube with a single synthesized beam, two options are available. Best is to use imsmooth with parameter kernel='commonbeam'. This task will smooth each plane to that of the lowest resolution and it will clean up all header variables such that only a single beam appears in the data. The second option is to define parameter resmooth=T directly in clean.


Image Header

The image header holds metadata associated with your CASA image. The task imhead will display this data within the casalog. We will first run imhead with parameter mode='summary' :

# In CASA
imhead(imagename='SNR_G55_10s.ms.MFS.wProj.image.tt0', mode='summary')
  • mode='summary' : gives general information about the image, including the object name, sky coordinates, image units, the telescope the data was taken with, and more.

For further information, run parameter mode='list' . We will assign a variable header to also capture the output as a Python dictionary:

# In CASA
header=imhead(imagename='SNR_G55_10s.ms.MFS.wProj.image.tt0', mode='list')
  • mode='list' : gives more detailed information, including beam major/minor axes, beam primary angle, and the location of the max/min intensity, and lots more. Essentially this mode displays the FITS header variables.


Image Conversion

Finally, as an example, we convert our image from intensity to main beam brightness temperature.

We will use the standard equation

[math]\displaystyle{ T=1.222\times 10^{6} \frac{S}{\nu^{2} \theta_{maj} \theta_{min}} }[/math]

where the main beam brightness temperature [math]\displaystyle{ T }[/math] is given in K, the intensity [math]\displaystyle{ S }[/math] in Jy/beam, the reference frequency [math]\displaystyle{ \nu }[/math] in GHz, and the major and minor beam sizes [math]\displaystyle{ \theta }[/math] in arcseconds.

For a beam of 29.30"x29.03" and a reference frequency of 1.579 GHz (as taken from the previous imhead run) we calculate the brightness temperature using immath:

# In CASA
immath(imagename='SNR_G55_10s.ms.MFS.wProj.image.tt0', mode='evalexpr', 
       expr='1.222e6*IM0/1.579^2/(29.30*29.03)', 
       outfile='SNR_G55_10s.ms.MFS.wProj.image.tt0-Tb')
  • mode='evalexpr' : immath is used to calcuate with images
  • expr: the mathematical expression to be evaluated. The images are abbreviated as IM0, IM1, ... in the sequence of their appearance in the imagename parameter.

Since immath only changes the values and not the unit of the image, we will now change the new image header 'bunit' key to 'K'. To do so, we will run imhead with parameter mode='put' :

# In CASA
imhead(imagename='SNR_G55_10s.ms.MFS.wProj.image.tt0-Tb', mode='put', hdkey='bunit', hdvalue='K')
  • hdkey: the header keyword that will be changed
  • hdvalue: the new value for header keyword

Launching the viewer will now show our image in brightness temperature units:

# In CASA
viewer('SNR_G55_10s.ms.MFS.wProj.image.tt0-Tb')


CASAguides

Last checked on CASA Version 4.5.2