NGC 5921: red-shifted HI emission 5.7.2: Difference between revisions

From CASA Guides
Jump to navigationJump to search
Jsalcido (talk | contribs)
351 edits
No edit summary
 
Fschinze (talk | contribs)
536 edits
 
(36 intermediate revisions by 5 users not shown)
Line 1: Line 1:
Disclaimer: Due to continuous CASA software updates, GUI images may look different on more recent versions than those shown here.
Disclaimer: Due to continuous CASA software updates, GUI images may look different on more recent versions than those shown here.


{{VLA Intro}}


[[Category: VLA]] [[Category: Calibration]] [[Category: Imaging]] [[Category: Spectral Line]]
[[Category: VLA]] [[Category: Calibration]] [[Category: Imaging]] [[Category: Spectral Line]]
{{Checked 5.7.2}}


== Overview ==
== Overview ==
Line 11: Line 12:
The following tutorial derives from an annotated script provided in the [http://casa.nrao.edu/ref_cookbook.shtml CASA Cookbook]. The script is largely reproduced and additionally annotated with figures and illustrations. It is assumed that this tutorial will be used interactively, and so interactive pauses in the original script have been removed.
The following tutorial derives from an annotated script provided in the [http://casa.nrao.edu/ref_cookbook.shtml CASA Cookbook]. The script is largely reproduced and additionally annotated with figures and illustrations. It is assumed that this tutorial will be used interactively, and so interactive pauses in the original script have been removed.


The data are included with the CASA installation.
<b>DATA</b>: The data are included with the CASA installation.
 


== Setting up the CASA Environment ==
== Setting up the CASA Environment ==
Line 21: Line 23:
mkdir NGC5921
mkdir NGC5921
cd NGC5921
cd NGC5921
casapy
casa
</source>
</source>


<!--
Now, in CASA, set up paths and global variables to facilitate the data reduction. These operations can be performed on-the-fly if you are reducing data interactively, but it's better to have them prepared from the start in the scripting environment.
Now, in CASA, set up paths and global variables to facilitate the data reduction. These operations can be performed on-the-fly if you are reducing data interactively, but it's better to have them prepared from the start in the scripting environment.


Line 38: Line 41:
imname = prefix + '.cleanimg'
imname = prefix + '.cleanimg'
</source>
</source>
-->


We'll use a python ''os'' command to get the appropriate CASA path for your installation. The use of '''os.environ.get''' is explained in [[#os.environ.get | the Appendix]].
We'll use a python ''os'' command to get the appropriate CASA path for your installation in order to import the data. The use of '''os.environ.get''' is explained in [[#os.environ.get | the Appendix]].


<source lang="python">
<source lang="python">
# In CASA
%cpaste
# Press Enter or Return, then copy/paste the following:
import os
pathname=os.environ.get('CASAPATH').split()[0]
pathname=os.environ.get('CASAPATH').split()[0]
fitsdata=pathname+'/data/demo/NGC5921.fits'
fitsdata=pathname+'/data/demo/NGC5921.fits'
--
</source>
</source>


Scripts are of course modified and repeated to the satisfaction of observer. To help clean up the bookkeeping and further avoid issues of write privileges, remove prior versions of the measurement set and calibration tables.
Scripts are of course modified and repeated to the satisfaction of observer. To help clean up the bookkeeping and further avoid issues of write privileges, remove prior versions of the measurement set and calibration tables.


This can be done with the rmtables('table_name') command.
<!--
<div style="background-color: #dddddd;">
<div style="background-color: #dddddd;">
'''Tip:''' The first command in the following code block removes the measurement set. Depending on the size of the source data, refilling a measurement set can be time-consuming, and so you may want to consider editing a separate script that preserves the measurement set but skips the '''importuvfits''' command [[#Import the Data | given below]]. This NGC5921 dataset, however, is not large and refilling goes quickly.
'''Tip:''' The first command in the following code block removes the measurement set. Depending on the size of the source data, refilling a measurement set can be time-consuming, and so you may want to consider editing a separate script that preserves the measurement set but skips the '''importuvfits''' command [[#Import the Data | given below]]. This NGC5921 dataset, however, is not large and refilling goes quickly.
Line 64: Line 77:
os.system('rm -rf '+prefix+'*')
os.system('rm -rf '+prefix+'*')
</source>
</source>
-->


== Import the Data ==
== Import the Data ==


The next step is to import the multisource UVFITS data to a CASA measurement set via the {{importuvfits}} filler.  
The next step is to import the multisource UVFITS data to a CASA measurement set via the {{importuvfits}} filler. Note that you can set each parameter for any particular task one-by-one, or you could supply the task and input parameters with one command. Here we will set each parameter value first, save them, and run the import task. Throughout the remaining tutorial, we will call upon tasks with a single command.
 
<!--
<source lang="python">
importuvfits(fitsfile=fitsdata, vis='ngc5921.demo.ms')
saveinputs('importuvfits', 'ngc5921.demo.importuvfits.saved')
</source>
-->


<source lang="python">
<source lang="python">
# Safest to start from task defaults
# Safest to start from task defaults
default('importuvfits')
default('importuvfits')
# Set up the MS filename and save as new global variable
msfile = prefix + '.ms'
# Use task importuvfits
# Use task importuvfits
fitsfile = fitsdata
fitsfile = fitsdata
vis = msfile
vis='ngc5921.demo.ms'
saveinputs('importuvfits',prefix+'.importuvfits.saved')
saveinputs('importuvfits', 'ngc5921.demo.importuvfits.saved')
importuvfits()
importuvfits()
</source>
</source>
Line 90: Line 109:
antnamescheme      =  "old"
antnamescheme      =  "old"
#importuvfits(fitsfile="/usr/lib64/casapy/30.0.9709test-001/data/demo/NGC5921.fits",vis="ngc5921.demo.ms",antnamescheme="old")
#importuvfits(fitsfile="/usr/lib64/casapy/30.0.9709test-001/data/demo/NGC5921.fits",vis="ngc5921.demo.ms",antnamescheme="old")
</pre>


</pre>


== A Summary of the Data ==
== A Summary of the Data ==
Line 112: Line 131:


<pre>
<pre>
2018-09-12 20:19:48 INFO listobs ##########################################
2018-09-12 20:19:48 INFO listobs ##### Begin Task: listobs            #####
2018-09-12 20:19:48 INFO listobs listobs(vis="ngc5921.demo.ms",selectdata=True,spw="",field="",antenna="",
2018-09-12 20:19:48 INFO listobs         uvrange="",timerange="",correlation="",scan="",intent="",
2018-09-12 20:19:48 INFO listobs         feed="",array="",observation="",verbose=True,listfile="",
2018-09-12 20:19:48 INFO listobs         listunfl=False,cachesize=50,overwrite=False)
2018-09-12 20:19:49 INFO listobs ================================================================================
2018-09-12 20:19:49 INFO listobs           MeasurementSet Name:  /lustre/aoc/sciops/akapinsk/oldVLA/ngc5921/ngc5921.demo.ms      MS Version 2
2018-09-12 20:19:49 INFO listobs ================================================================================
2018-09-12 20:19:49 INFO listobs   Observer: TEST    Project: 
2018-09-12 20:19:49 INFO listobs Observation: VLA
2018-09-12 20:19:49 INFO listobs Computing scan and subscan properties...
2018-09-12 20:19:49 INFO listobs Data records: 22653      Total elapsed time = 5310 seconds
2018-09-12 20:19:49 INFO listobs   Observed from  13-Apr-1995/09:18:45.0  to  13-Apr-1995/10:47:15.0 (TAI)
2018-09-12 20:19:49 INFO listobs  
2018-09-12 20:19:49 INFO listobs   ObservationID = 0        ArrayID = 0
2018-09-12 20:19:49 INFO listobs   Date        Timerange (TAI)          Scan  FldId FieldName            nRows    SpwIds  Average Interval(s)    ScanIntent
2018-09-12 20:19:49 INFO listobs   13-Apr-1995/09:18:45.0 - 09:24:45.0    1      0 1331+30500002_0          4509  [0]  [30]
2018-09-12 20:19:49 INFO listobs               09:27:15.0 - 09:29:45.0    2      1 1445+09900002_0          1890  [0]  [30]
2018-09-12 20:19:49 INFO listobs               09:32:45.0 - 09:48:15.0    3      2 N5921_2                  6048  [0]  [30]
2018-09-12 20:19:49 INFO listobs               09:50:15.0 - 09:51:15.0    4      1 1445+09900002_0            756  [0]  [30]
2018-09-12 20:19:49 INFO listobs               10:21:45.0 - 10:23:15.0    5      1 1445+09900002_0          1134  [0]  [30]
2018-09-12 20:19:49 INFO listobs               10:25:45.0 - 10:43:15.0    6      2 N5921_2                  6804  [0]  [30]
2018-09-12 20:19:49 INFO listobs               10:45:15.0 - 10:47:15.0    7      1 1445+09900002_0          1512  [0]  [30]
2018-09-12 20:19:49 INFO listobs           (nRows = Total number of rows per scan)
2018-09-12 20:19:49 INFO listobs Fields: 3
2018-09-12 20:19:49 INFO listobs   ID  Code Name                RA              Decl          Epoch  SrcId      nRows
2018-09-12 20:19:49 INFO listobs   0        1331+30500002_0    13:31:08.287300 +30.30.32.95900 J2000  0          4509
2018-09-12 20:19:49 INFO listobs   1        1445+09900002_0    14:45:16.465600 +09.58.36.07300 J2000  1          5292
2018-09-12 20:19:49 INFO listobs   2        N5921_2            15:22:00.000000 +05.04.00.00000 J2000  2          12852
2018-09-12 20:19:49 INFO listobs Spectral Windows:  (1 unique spectral windows and 1 unique polarization setups)
2018-09-12 20:19:49 INFO listobs   SpwID  Name  #Chans  Frame  Ch0(MHz)  ChanWid(kHz)  TotBW(kHz) CtrFreq(MHz)  Corrs 
2018-09-12 20:19:49 INFO listobs   0      none      63  LSRK    1412.665        24.414      1550.2  1413.4219  RR  LL
2018-09-12 20:19:49 INFO listobs Sources: 3
2018-09-12 20:19:49 INFO listobs   ID  Name                SpwId RestFreq(MHz)  SysVel(km/s)
2018-09-12 20:19:49 INFO listobs   0    1331+30500002_0    0    1420.405752    0           
2018-09-12 20:19:49 INFO listobs   1    1445+09900002_0    0    1420.405752    0           
2018-09-12 20:19:49 INFO listobs   2    N5921_2            0    1420.405752    0           
2018-09-12 20:19:49 INFO listobs Antennas: 27:
2018-09-12 20:19:49 INFO listobs   ID  Name  Station  Diam.    Long.        Lat.                Offset from array center (m)                ITRF Geocentric coordinates (m)       
2018-09-12 20:19:49 INFO listobs                                                                     East        North    Elevation              x              y              z
2018-09-12 20:19:49 INFO listobs   0    1    VLA:N7    25.0 m  -107.37.07.2  +33.54.12.9        -30.2623      345.7477      -0.8872 -1601155.613187 -5041783.882304  3555162.343090
2018-09-12 20:19:49 INFO listobs   1    2    VLA:W1    25.0 m  -107.37.05.9  +33.54.00.5          3.5004      -39.7725        0.9883 -1601188.991307 -5042000.530918  3554843.409670
2018-09-12 20:19:49 INFO listobs   2    3    VLA:W2    25.0 m  -107.37.07.4  +33.54.00.9        -37.1358      -25.0262        1.0383 -1601225.244615 -5041980.431775  3554855.677111
2018-09-12 20:19:49 INFO listobs   3    4    VLA:E1    25.0 m  -107.37.05.7  +33.53.59.2          6.9833      -79.6414        1.1565 -1601192.444530 -5042022.911771  3554810.411780
2018-09-12 20:19:49 INFO listobs   4    5    VLA:E3    25.0 m  -107.37.02.8  +33.54.00.5        81.5188      -37.9632        1.0246 -1601114.335629 -5042023.211477  3554844.931655
2018-09-12 20:19:49 INFO listobs   5    6    VLA:E9    25.0 m  -107.36.45.1  +33.53.53.6        536.8977    -250.3175        0.1183 -1600715.915813 -5042273.186780  3554668.167811
2018-09-12 20:19:49 INFO listobs   6    7    VLA:E6    25.0 m  -107.36.55.6  +33.53.57.7        267.7566    -124.8145        1.2815 -1600951.554888 -5042125.947753  3554772.987072
2018-09-12 20:19:49 INFO listobs   7    8    VLA:W8    25.0 m  -107.37.21.6  +33.53.53.0      -401.2640    -270.6305        2.2293 -1601614.059494 -5042001.699973  3554652.484758
2018-09-12 20:19:49 INFO listobs   8    9    VLA:N5    25.0 m  -107.37.06.7  +33.54.08.0        -16.9948      194.1215      -0.1368 -1601168.756077 -5041869.099542  3555036.914367
2018-09-12 20:19:49 INFO listobs   9    10    VLA:W3    25.0 m  -107.37.08.9  +33.54.00.1        -74.4964      -50.1921        1.1608 -1601265.132224 -5041982.597979  3554834.857504
2018-09-12 20:19:49 INFO listobs   10  11    VLA:N4    25.0 m  -107.37.06.5  +33.54.06.1        -11.7487      134.3686        0.1774 -1601173.922897 -5041902.701204  3554987.495105
2018-09-12 20:19:49 INFO listobs   11  12    VLA:W5    25.0 m  -107.37.13.0  +33.53.57.8      -179.2554    -120.8635        1.4872 -1601376.990711 -5041988.712764  3554776.381187
2018-09-12 20:19:49 INFO listobs   12  13    VLA:N3    25.0 m  -107.37.06.3  +33.54.04.8        -8.2438      94.5297        0.3947 -1601177.362708 -5041925.112425  3554954.550128
2018-09-12 20:19:49 INFO listobs   13  14    VLA:N1    25.0 m  -107.37.06.0  +33.54.01.8        -0.0030        0.0445        0.8773 -1601185.580779 -5041978.216463  3554876.396287
2018-09-12 20:19:49 INFO listobs   14  15    VLA:N2    25.0 m  -107.37.06.2  +33.54.03.5        -4.7904      54.7090        0.5774 -1601180.839839 -5041947.470902  3554921.600805
2018-09-12 20:19:49 INFO listobs   15  16    VLA:E7    25.0 m  -107.36.52.4  +33.53.56.5        348.8969    -162.6653        1.0336 -1600880.544215 -5042170.427468  3554741.431900
2018-09-12 20:19:49 INFO listobs   16  17    VLA:E8    25.0 m  -107.36.48.9  +33.53.55.1        438.6654    -204.5038        0.5027 -1600801.910482 -5042219.412805  3554706.408864
2018-09-12 20:19:49 INFO listobs   17  18    VLA:W4    25.0 m  -107.37.10.8  +33.53.59.1      -122.0163      -82.2819        1.2624 -1601315.866196 -5041985.352573  3554808.279150
2018-09-12 20:19:49 INFO listobs   18  19    VLA:E5    25.0 m  -107.36.58.4  +33.53.58.8        195.8349      -91.2758        1.2155 -1601014.427180 -5042086.300814  3554800.787928
2018-09-12 20:19:49 INFO listobs   19  20    VLA:W9    25.0 m  -107.37.25.1  +33.53.51.0      -491.1000    -331.2429        2.5539 -1601709.998072 -5042006.975455  3554602.355417
2018-09-12 20:19:49 INFO listobs   20  21    VLA:W6    25.0 m  -107.37.15.6  +33.53.56.4      -244.9704    -165.2178        1.6861 -1601447.161927 -5041992.554228  3554739.677219
2018-09-12 20:19:49 INFO listobs   21  22    VLA:E4    25.0 m  -107.37.00.8  +33.53.59.7        133.6478      -62.2829        1.0919 -1601068.773396 -5042051.970054  3554824.783566
2018-09-12 20:19:49 INFO listobs   23  24    VLA:E2    25.0 m  -107.37.04.4  +33.54.01.1        40.6649      -18.9151        0.9550 -1601150.040469 -5042000.665669  3554860.702914
2018-09-12 20:19:49 INFO listobs   24  25    VLA:N6    25.0 m  -107.37.06.9  +33.54.10.3        -23.2197      265.3902      -0.4819 -1601162.569974 -5041829.054708  3555095.873969
2018-09-12 20:19:49 INFO listobs   25  26    VLA:N9    25.0 m  -107.37.07.8  +33.54.19.0        -46.5533      532.1581      -1.8550 -1601139.422904 -5041679.082136  3555316.518142
2018-09-12 20:19:49 INFO listobs   26  27    VLA:N8    25.0 m  -107.37.07.5  +33.54.15.8        -38.0437      434.7201      -1.3387 -1601147.894127 -5041733.868915  3555235.935926
2018-09-12 20:19:49 INFO listobs   27  28    VLA:W7    25.0 m  -107.37.18.4  +33.53.54.8      -319.1171    -215.2368        1.9407 -1601526.340031 -5041996.897001  3554698.302182
2018-09-12 20:19:49 INFO listobs ##### End Task: listobs              #####
2018-09-12 20:19:49 INFO listobs ##########################################
</pre>
<!-- previous version for CASA 5.1
2011-04-25 19:52:17 INFO listobs ##########################################
2011-04-25 19:52:17 INFO listobs ##########################################
2011-04-25 19:52:17 INFO listobs ##### Begin Task: listobs            #####
2011-04-25 19:52:17 INFO listobs ##### Begin Task: listobs            #####
Line 179: Line 271:
2011-04-25 19:52:17 INFO listobs ##### End Task: listobs              #####
2011-04-25 19:52:17 INFO listobs ##### End Task: listobs              #####
2011-04-25 19:52:17 INFO listobs ##########################################
2011-04-25 19:52:17 INFO listobs ##########################################
</pre>
-->


=== Key Information from {{listobs | Listobs}} ===
=== Key Information from {{listobs | Listobs}} ===
Line 190: Line 282:
* The data were taken in a single IF (a single spectral window, SpwID = 0), divided into 63 channels.
* The data were taken in a single IF (a single spectral window, SpwID = 0), divided into 63 channels.
* Only RR and LL correlations are present; cross-pols are absent.
* Only RR and LL correlations are present; cross-pols are absent.


== Flagging ==
== Flagging ==
Line 195: Line 288:
=== Flag the autocorrelations ===
=== Flag the autocorrelations ===


We don't need the autocorrelation data, and {{flagautocorr}} gets rid of them. You shouldn't have to specify the measurement set, because the variable ''vis'' is already set, but it never hurts to be cautious.
We don't need the autocorrelation data, and we can use {{flagdata}} to get rid of them. You shouldn't have to specify the measurement set, because the variable ''vis'' is already set, but it never hurts to be cautious.


<source lang="python">
<source lang="python">
flagautocorr(vis="ngc5921.demo.ms")
flagdata(vis="ngc5921.demo.ms", autocorr=True)
</source>
</source>
CASA issues a warning when running this command, but this can be disregarded.


=== Interactive Flagging ===
=== Interactive Flagging ===


[[file:spectrum-flagging.png|thumb|Plotms settings for flagging spectral line data. Click to enlarge.]]
[[file:spectrum-flagging.png|thumb|Plotms settings for flagging spectral line data. Click to enlarge.]]


{{plotms | Plotms}} is a good tool for flagging spectral line data. Check out the [[Data flagging with casaplotms | tutorial]] that describes editing VLA continuum data. Spectral line data of course require some consideration of channels and channel averaging.  
{{plotms | Plotms}} is a good tool for flagging spectral line data. Check out the [[Data flagging with casaplotms | tutorial]] that describes editing VLA continuum data. Spectral line data of course require some consideration of channels and channel averaging.  
Line 215: Line 309:
The figure at right highlights the settings needed for effective editing of a spectral line data set. The key settings are as follows.
The figure at right highlights the settings needed for effective editing of a spectral line data set. The key settings are as follows.


* Specify the measurement set in '''File Location'''; the '''Browse''' button allows you to hunt down the measurement set.
* Specify the measurement set in '''File'''; the '''Browse''' button allows you to hunt down the measurement set.
* It's better to edit one source at a time. In the illustrated example, the flux / bandpass calibrator 1331+305* is displayed.
* It's better to edit one source at a time. In the illustrated example, the flux / bandpass calibrator 1331+305* is displayed.
* Average the channels. First, specify the central channels to remove band edge effects. Channels 6~56 in the first spectral window (IF) are appropriate (see [[#Inspect the Bandpass Response Curve]], below). In the '''Channel Averaging''' box, enter 51 channels to average over all channels in the given range.
* Average the channels. First, specify the central channels to remove band edge effects. Channels 6~56 in the first spectral window (IF) are appropriate (see [[#Inspect the Bandpass Response Curve]], below). In plotms() specify: spw=0:6~56. In the '''Channel Averaging''' box, enter 51 channels to average over all channels in the given range.
* Ideally you want the channels to have the same (''u'', ''v'') coverage (projected baseline spacings as viewed from the source); otherwise, the beam (point spread function) will be different for each channel. Therefore, if you flag data from a given channel it's usually a good idea to flag those data from all channels. Under the '''Flagging''' tab, specify '''Extend flags''' to '''Channel.'''
* Ideally you want the channels to have the same (''u'', ''v'') coverage (projected baseline spacings as viewed from the source); otherwise, the beam (point spread function) will be different for each channel. Therefore, if you flag data from a given channel it's usually a good idea to flag those data from all channels. Under the '''Flagging''' tab, specify '''Extend flags''' to '''Channel.'''


Line 227: Line 321:
Calibration of spectral line data broadly follows the approach for [[Calibrating a VLA 5 GHz continuum survey | continuum data]], except that the amplitude and phase corrections are a function of frequency and so must be corrected by '''bandpass calibration.''' The basic calibration steps follow.
Calibration of spectral line data broadly follows the approach for [[Calibrating a VLA 5 GHz continuum survey | continuum data]], except that the amplitude and phase corrections are a function of frequency and so must be corrected by '''bandpass calibration.''' The basic calibration steps follow.
* [[#Setting the Flux Scale | Set the flux scale]] of the primary calibrator, here, 1331+305 = 3C 286.
* [[#Setting the Flux Scale | Set the flux scale]] of the primary calibrator, here, 1331+305 = 3C 286.
* [[#Bandpass Calibration | Determine bandpass corrections]] based on the primary calibrator. In the script that follows, the bandpass calibration is stored in ''ngc5921.demo.bcal'', which is itself referenced by the python variable ''btable.''
* [[#Bandpass Calibration | Determine bandpass corrections]] based on the primary calibrator. In the script that follows, the bandpass calibration will be stored in ''ngc5921.demo.bcal''.
* [[#Inspect the Bandpass Response Curve | Inspect the bandpass correction]] to determine viable channels for averaging and imaging. We want to toss out end channels where the response is poor.
* [[#Inspect the Bandpass Response Curve | Inspect the bandpass correction]] to determine viable channels for averaging and imaging. We want to toss out end channels where the response is poor.
* [[#Gain Calibration | Determine the gain calibrations]] on the bandpass-corrected and channel-averaged data. In this step, we effectively turn the spectral line data into a single-channel continuum data set and calibrate accordingly. The calibration is stored in ''ngc5921.demo.gcal'', which is itself referenced by the python variable ''gtable.''
* [[#Gain Calibration | Determine the gain calibrations]] on the bandpass-corrected and channel-averaged data. In this step, we effectively turn the spectral line data into a single-channel continuum data set and calibrate accordingly. The calibration is first stored in ''ngc5921.demo.gcal''. In the second part of this step we correct the fluxscale of the .gcal table, and store the final calibration solutions with correct fluxes in ''ngc5921.demo.fluxscale'' (this is the table that needs to be applied later to the data, not the .gcal version).
* [[#Inspect the Calibration Solutions | Inspect the gain calibration solutions]] to look for any aberrant solutions that hint at bad calibrator data.
* [[#Inspect the Calibration Solutions | Inspect the gain calibration solutions]] to look for any aberrant solutions that hint at bad calibrator data.
* [[#Apply the Solutions | Apply the calibration solutions]] to the source (N5921_2). This action literally adds a new column of data to the measurement set. This new column contains the data with the gain calibration and bandpass calibration applied, but it does not overwrite the raw data in case the calibration needs revision.
* [[#Apply the Solutions | Apply the calibration solutions]] to the source (N5921_2). This action literally adds a new column of data to the measurement set. This new column contains the data with the gain calibration and bandpass calibration applied, but it does not overwrite the raw data in case the calibration needs revision.




=== Setting the Flux Scale ===
{{setjy | Setjy}} generates a source model for the primary calibrator, 1331+305 = 3C286. From CASA 5.3+ the default standard is Perley-Butler 2017, and includes resolved structure of the calibrators. This is 1.4GHz D-config and 1331+305 is sufficiently unresolved that, in principle, we don't need a model image; however, here we proceed with applying the detailed model, as a good practice.
{{setjy | Setjy}} also looks up the radio SED for common flux calibrators and automatically assigns the total flux density.
<source lang="python">
# 1331+305 = 3C286 is our primary calibrator. Use the wildcard on the end of the source name
setjy(vis='ngc5921.demo.ms', field='1331+305*', model='3C286_L.im')
</source>
A summary of the operation is sent to the logger window. Here's a listing of the output.
<pre>
2018-09-12 20:43:38 INFO setjy ##########################################
2018-09-12 20:43:38 INFO setjy ##### Begin Task: setjy              #####
2018-09-12 20:43:38 INFO setjy setjy(vis="ngc5921.demo.ms",field="1331+305*",spw="",selectdata=False,timerange="",
2018-09-12 20:43:38 INFO setjy         scan="",intent="",observation="",scalebychan=True,standard="Perley-Butler 2017",
2018-09-12 20:43:38 INFO setjy         model="3C286_L.im",modimage="",listmodels=False,fluxdensity=-1,spix=0.0,
2018-09-12 20:43:38 INFO setjy         reffreq="1GHz",polindex=[],polangle=[],rotmeas=0.0,fluxdict={},
2018-09-12 20:43:38 INFO setjy         useephemdir=False,interpolation="nearest",usescratch=False,ismms=False)
2018-09-12 20:43:38 INFO setjy {'field': '1331+305*'}
2018-09-12 20:43:38 INFO Imager Opening MeasurementSet [...]
2018-09-12 20:43:38 INFO setjy Using /home/casa/data/distro/nrao/VLA/CalModels/3C286_L.im for modimage.
2018-09-12 20:43:38 INFO setjy CASA Version 5.3.0-143 
2018-09-12 20:43:38 INFO setjy  
2018-09-12 20:43:39 INFO imager Using channel dependent flux densities
2018-09-12 20:43:39 INFO imager Selected 4509 out of 22653 rows.
2018-09-12 20:43:39 INFO imager 1331+30500002_0 (fld ind 0) spw 0  [I=15.016, Q=0, U=0, V=0] Jy @ 1.4127e+09Hz, (Perley-Butler 2017)
2018-09-12 20:43:40 INFO imager Using model image /home/casa/data/distro/nrao/VLA/CalModels/3C286_L.im
2018-09-12 20:43:40 INFO imager Scaling spw(s) [0]'s model image by channel to  I = 15.0159, 15.0118, 15.0077 Jy @(1.41265e+09, 1.41343e+09, 1.41419e+09)Hz (LSRK) for visibility prediction (a few representative values are shown).
2018-09-12 20:43:40 INFO imager The model image's reference pixel is 0.00904522 arcsec from 1331+30500002_0's phase center.
2018-09-12 20:43:40 INFO imager Will clear any existing model with matching field=1331+30500002_0 and spw=*
2018-09-12 20:43:40 INFO  Clearing model records in MS header for selected fields.
2018-09-12 20:43:40 INFO  1331+30500002_0 (id = 0) deleted.
2018-09-12 20:43:40 INFO imager Selected 4509 out of 22653 rows.
2018-09-12 20:43:40 INFO setjy ##### End Task: setjy                #####
2018-09-12 20:43:40 INFO setjy ##########################################
</pre>


=== Setting the Flux Scale ===
<!-- PREVIOUS VERSION, FOR CASA 5.1.0


{{setjy | Setjy}} generates a point source model for the primary calibrator, 1331+305 = 3C286. These data are of low angular resolution, and so the point source model is adequate for our purposes. For observations with higher angular resolution (longer baseline configurations), you can specify a model of the calibrator using the ''modimage'' parameter (see the tutorial [[Calibrating a VLA 5 GHz continuum survey#Set the Flux Scale]] for an example of how to use ''modimage'').  
{{setjy | Setjy}} generates a point source model for the primary calibrator, 1331+305 = 3C286. These data are of low angular resolution, and so the point source model is adequate for our purposes. For observations with higher angular resolution (longer baseline configurations), you can specify a model of the calibrator using the ''modimage'' parameter (see the tutorial [[Calibrating a VLA 5 GHz continuum survey#Set the Flux Scale]] for an example of how to use ''modimage'').  
Line 241: Line 375:
{{setjy | Setjy}} also looks up the radio SED for common flux calibrators and automatically assigns the total flux density.
{{setjy | Setjy}} also looks up the radio SED for common flux calibrators and automatically assigns the total flux density.


<source lang="python">
# 1331+305 = 3C286 is our primary calibrator. Use the wildcard on the end of the source name
# This is 1.4GHz D-config and 1331+305 is sufficiently unresolved that we dont need a model image. 
# For higher frequencies (particularly in A and B config) you would want to use one.
setjy(vis='ngc5921.demo.ms', field='1331+305*', modimage='')
</source>
<!--
<source lang="python">
<source lang="python">
default('setjy')
default('setjy')
Line 258: Line 400:
setjy()
setjy()
</source>
</source>
-->


A summary of the operation is sent to the logger window. Here's a listing of the output.
<!-- PREVIOUS VERSION, FOR CASA 5.1.0


<pre>
<pre>
Line 279: Line 422:
2011-04-25 19:50:55 INFO setjy ##########################################
2011-04-25 19:50:55 INFO setjy ##########################################
</pre>
</pre>
-->


=== Bandpass Calibration ===
=== Bandpass Calibration ===


The flux calibrator 1331+305 = 3C 286 now has a point-source model assigned to it. Since the point-source model doesn't change over this narrow range of frequencies, we can use the model to determine amplitude and phase (gain) corrections for each channel independently. The result is the bandpass calibration.
The flux calibrator 1331+305 = 3C 286 now has a model assigned to it. Since the bandwidth of our observations is only 1.55 MHz, the model doesn't change over this narrow range of frequencies, so we can use it to determine amplitude and phase (gain) corrections for each channel independently. The result is the bandpass calibration.


As for any antenna-based calibration scheme, we have to pick an antenna to act as the reference point for the calibration. Any antenna will do, but it's better to pick one near the center of the array. For the remainder of the calibration, we will use '''refant = '15''''.
As for any antenna-based calibration scheme, we have to pick an antenna to act as the reference point for the calibration. Any antenna will do, but it's better to pick one near the center of the array. For the remainder of the calibration, we will use '''refant = '15''''.


<source lang="python">
<source lang="python">
# We can first do the bandpass on the single 5min scan on 1331+305. At 1.4GHz phase stablility should be sufficient to do this without
# a first (rough) gain calibration. This will give us the relative antenna gain as a function of frequency.
bandpass(vis='ngc5921.demo.ms', caltable='ngc5921.demo.bcal', field='0', selectdata=False, bandtype='B', solint='inf', combine='scan', refant='15')
</source>
* field='0' : Use the flux calibrator 1331+305 = 3C286 (FIELD_ID 0) as bandpass calibrator.
* bandtype='B' : Choose bandpass solution type. Pick standard time-binned B (rather than BPOLY).
* solint='inf' and combine='scan' : Set solution interval arbitrarily long (get single bandpass).
* refant = '15' : Reference antenna Name 15 (15=VLA:N2) (Id 14)
<!--
default('bandpass')
default('bandpass')
# We can first do the bandpass on the single 5min scan on 1331+305
# We can first do the bandpass on the single 5min scan on 1331+305
Line 309: Line 464:
# antenna gain-elevation curve or atmospheric opacity
# antenna gain-elevation curve or atmospheric opacity
# (at 8GHz and above you would want these)
# (at 8GHz and above you would want these)
gaincurve = False
opacity = 0.0
opacity = 0.0
# Choose bandpass solution type
# Choose bandpass solution type
Line 321: Line 475:
saveinputs('bandpass',prefix+'.bandpass.saved')
saveinputs('bandpass',prefix+'.bandpass.saved')
bandpass()
bandpass()
</source>
-->
 


=== Inspect the Bandpass Response Curve ===
=== Inspect the Bandpass Response Curve ===


[[File:ngc5921-bandpass.png | thumb | Bandpass response curves generated by {{plotcal}}. The solutions for different antennas are indicated by differently colored plotting symbols. Plots for individual antennas can be generated by setting '''iteration = 'antenna'''' for {{plotcal}}.]]
[[File:ngc5921_bandpass_casa5.5.png | thumb | Bandpass response curves generated by {{plotms}}. The solutions for different antennas are indicated by differently colored plotting symbols. Plots for individual antennas can be generated by setting '''iteraxis = 'antenna'''' for {{plotms}}.]]
 
In the [[#Gain Calibration | gain calibration]] to follow, we will effectively convert the spectral line data into a continuum data set. Before proceeding, we need to inspect the bandpass calibration to make sure that it contains no bad values and also to inspect which channels to average to produce the continuum data. Plotms is the standard tool for plotting calibration solutions. The following commands produce the figure at right.
 
<source lang="python">
# Set up 2x1 panels - upper panel amp vs. channel
plotms(vis='ngc5921.demo.bcal', field='0', gridrows=2, gridcols=1, plotindex=0, rowindex=0, xaxis='channel', yaxis='amp', showgui=True, clearplots=False, coloraxis='antenna1')
</source>


In the [[#Gain Calibration | gain calibration]] to follow, we will effectively convert the spectral line data into a continuum data set. Before proceeding, we need to inspect the bandpass calibration to make sure that it contains no bad values and also to inspect which channels to average to produce the continuum data. {{plotcal | Plotcal}} is the standard tool for plotting calibration solutions. The following script produces the figure at right.  
<source lang="python">
# Set up 2x1 panels - lower panel phase vs. channel
plotms(vis='ngc5921.demo.bcal', field='0', gridrows=2, gridcols=1, plotindex=1, rowindex=1, xaxis='channel', yaxis='phase', showgui=True, clearplots=False, coloraxis='antenna1')
</source>


By inspection, the amplitude response curve is flat over channels 6~56; that channel range will be used to generate the continuum data for gain calibration.
By inspection, the amplitude response curve is flat over channels 6~56; that channel range will be used to generate the continuum data for gain calibration. If you want to further inspect the plots interactively and iterate over antenna, set iteraxis = 'antenna'


Notice that {{plotcal}} is run twice: once to display gain amplitudes as a function of channel (frequency), and again to plot gain phases as a function of channel.
Notice that plotms is run twice: once to display gain amplitudes as a function of channel (frequency), and again to plot gain phases as a function of channel.


<!--
If the '''scriptmode''' is set to '''False''', the plot is saved to ngc5921.demo.plotcal.png.
If the '''scriptmode''' is set to '''False''', the plot is saved to ngc5921.demo.plotcal.png.
 
-->
<source lang="python">
<!--
default('plotcal')
default('plotcal')
caltable = btable
caltable = btable
Line 356: Line 522:
# iteration = 'antenna'
# iteration = 'antenna'
plotcal()
plotcal()
</source>
-->


=== Gain Calibration ===
=== Gain Calibration ===
Line 372: Line 538:
* The default model for gain calibrations is a 1 Jy point-source. The flux scale is overridden by {{setjy}}, which has been performed for the flux calibrator. We need to transfer that flux scale to the phase calibrator using {{fluxscale}}().   
* The default model for gain calibrations is a 1 Jy point-source. The flux scale is overridden by {{setjy}}, which has been performed for the flux calibrator. We need to transfer that flux scale to the phase calibrator using {{fluxscale}}().   
* Note that fluxscale() determines the flux density of the phase calibrator and accordingly adjusts its model and calibration solutions. A report of the results are sent to the logger window.
* Note that fluxscale() determines the flux density of the phase calibrator and accordingly adjusts its model and calibration solutions. A report of the results are sent to the logger window.
* Unless you use parameter ''incremental=True'' while executing {{fluxscale}}() (the default is ''False''), the resulting .fluxscale table will replace the .gcal table at this point. This particularly important in the {{applycal}}() stage.
<source lang="python">
# Armed with the bandpass, we now solve for the time-dependent antenna gains using our previously determined bandpass.
# Note this will automatically be applied to all sources not just the one used to determine the bandpass
gaincal(vis='ngc5921.demo.ms', caltable='ngc5921.demo.gcal', gaintable=['ngc5921.demo.bcal'], interp=['nearest'], field='0,1', spw='0:6~56', gaintype='G', solint='inf', calmode='ap', refant='15')
</source>
<!-- gaincal(vis='ngc5921.demo.ms', caltable='ngc5921.demo.gcal', gaintable='ngc5921.demo.bcal', interp='nearest', field='0,1',
        spw='0:6~56', gaintype='G', solint='inf', calmode='ap', refant='15', minsnr=1.0)-->


<source lang="python">
# Now we will transfer the flux scale to the phase calibrator.
# We will be using 1331+305 (the source we did setjy on) as our flux standard reference.
# Note its extended name as in the FIELD table summary above (it has a VLA seq number appended)
fluxscale(vis='ngc5921.demo.ms', fluxtable='ngc5921.demo.fluxscale', caltable='ngc5921.demo.gcal', reference='1331*', transfer='1445*')
</source>
<!--
<source lang="python">
<source lang="python">
default('gaincal')
default('gaincal')
# Armed with the bandpass, we now solve for the
 
# time-dependent antenna gains
vis = msfile
vis = msfile
# set the name for the output gain caltable
# set the name for the output gain caltable
Line 400: Line 585:
# antenna gain-elevation curve or atmospheric opacity
# antenna gain-elevation curve or atmospheric opacity
# (at 8GHz and above you would want these)
# (at 8GHz and above you would want these)
gaincurve = False
opacity = 0.0
opacity = 0.0
# scan-based G solutions for both amplitude and phase
# scan-based G solutions for both amplitude and phase
Line 434: Line 618:
fluxscale()
fluxscale()
</source>
</source>
-->


The output from {{fluxscale}} follows. A relatively large uncertainty for the phase calibrator is a sign that something went wrong, perhaps bad solutions in {{gaincal}}. Here, the phase calibrator scaled to 2.486 &plusmn; 0.001 Jy, which looks reasonable.
The output from {{fluxscale}} follows. A relatively large uncertainty for the phase calibrator is a sign that something went wrong, perhaps bad solutions in {{gaincal}}. Here, the phase calibrator scaled to 2.486 &plusmn; 0.001 Jy, which looks reasonable.


<pre>
<pre>
2018-09-12 22:38:05 INFO fluxscale ##########################################
2018-09-12 22:38:05 INFO fluxscale ##### Begin Task: fluxscale          #####
2018-09-12 22:38:05 INFO fluxscale fluxscale(vis="ngc5921.demo.ms",caltable="ngc5921.demo.gcal",fluxtable="ngc5921.demo.fluxscale",reference="1331*",transfer="1445*",
2018-09-12 22:38:05 INFO fluxscale         listfile="",append=False,refspwmap=[-1],gainthreshold=-1.0,antenna="",
2018-09-12 22:38:05 INFO fluxscale         timerange="",scan="",incremental=False,fitorder=1,display=False)
2018-09-12 22:38:05 INFO fluxscale ****Using NEW VI2-driven calibrater tool****
2018-09-12 22:38:05 INFO fluxscale Opening MS: ngc5921.demo.ms for calibration.
2018-09-12 22:38:05 INFO fluxscale Initializing nominal selection to the whole MS.
2018-09-12 22:38:05 INFO fluxscale Beginning fluxscale--(MSSelection version)-------
2018-09-12 22:38:05 INFO fluxscale Found reference field(s): 1331+30500002_0
2018-09-12 22:38:05 INFO fluxscale Found transfer field(s):  1445+09900002_0
2018-09-12 22:38:05 INFO fluxscale Flux density for 1445+09900002_0 in SpW=0 (freq=1.41342e+09 Hz) is: 2.52609 +/- 0.00218206 (SNR = 1157.67, N = 54)
2018-09-12 22:38:05 INFO fluxscale Storing result in ngc5921.demo.fluxscale
2018-09-12 22:38:05 INFO fluxscale Writing solutions to table: ngc5921.demo.fluxscale
2018-09-12 22:38:06 INFO fluxscale CASA Version 5.3.0-143 
2018-09-12 22:38:06 INFO fluxscale  
2018-09-12 22:38:07 INFO fluxscale ##### End Task: fluxscale            #####
2018-09-12 22:38:07 INFO fluxscale ##########################################
</pre>
<!--
2011-04-25 19:54:25 INFO fluxscale ##########################################
2011-04-25 19:54:25 INFO fluxscale ##########################################
2011-04-25 19:54:25 INFO fluxscale ##### Begin Task: fluxscale          #####
2011-04-25 19:54:25 INFO fluxscale ##### Begin Task: fluxscale          #####
Line 446: Line 652:
2011-04-25 19:54:26 INFO fluxscale Found reference field(s): 1331+30500002_0
2011-04-25 19:54:26 INFO fluxscale Found reference field(s): 1331+30500002_0
2011-04-25 19:54:26 INFO fluxscale Found transfer field(s):  1445+09900002_0
2011-04-25 19:54:26 INFO fluxscale Found transfer field(s):  1445+09900002_0
2011-04-25 19:54:26 INFO fluxscale Flux density for 1445+09900002_0 in SpW=0 is: 2.49071 +/- 0.00119335 (SNR = 2087.15, nAnt= 27)
2011-04-25 19:54:26 INFO fluxscale Flux density for 1445+09900002_0 in SpW=0 (freq=1.41342e+09 Hz) is: 2.53882 +/- 0.00218946 (SNR = 1159.57, N = 54)
2011-04-25 19:54:26 INFO fluxscale Storing result in ngc5921.demo.fluxscale
2011-04-25 19:54:26 INFO fluxscale Storing result in ngc5921.demo.fluxscale
2011-04-25 19:54:26 INFO fluxscale Writing solutions to table: ngc5921.demo.fluxscale
2011-04-25 19:54:26 INFO fluxscale Writing solutions to table: ngc5921.demo.fluxscale
Line 453: Line 659:
2011-04-25 19:54:26 INFO fluxscale ##########################################
2011-04-25 19:54:26 INFO fluxscale ##########################################
</pre>
</pre>
-->


=== Inspect the Calibration Solutions ===
=== Inspect the Calibration Solutions ===
Line 461: Line 668:


<source lang="python">
<source lang="python">
# Set up 2x1 panels - upper panel amp vs. time
plotms(vis='ngc5921.demo.fluxscale', field='0,1', gridrows=2, gridcols=1, plotindex=0, rowindex=0, yaxis='amp', showgui=True, clearplots=False, coloraxis='antenna1')
</source>
<source lang="python">
# Set up 2x1 panels - lower panel phase vs. time
plotms(vis='ngc5921.demo.fluxscale', field='0,1', gridrows=2, gridcols=1, plotindex=1, rowindex=1, yaxis='phase', showgui=True, clearplots=False, coloraxis='antenna1')
</source>
The amp and phase coherence looks good. If you want to do this interactively and iterate over antenna, set iteraxis = 'antenna'.
<!--
default('plotcal')
default('plotcal')
caltable = ftable
caltable = ftable
Line 482: Line 701:
#iteration = 'antenna'
#iteration = 'antenna'
plotcal()
plotcal()
</source>
-->
 


=== Apply the Solutions ===
=== Apply the Solutions ===
Line 489: Line 709:


<source lang="python">
<source lang="python">
gaintable = ['ngc5921.demo.gcal', 'ngc5921.demo.bcal']
gaintable = ['ngc5921.demo.fluxscale', 'ngc5921.demo.bcal']
</source>
</source>


We want to correct the calibrators using themselves and transfer from 1445+099 to itself and the target N5921. Start with the fluxscale/gain and bandpass tables. We will pick the 1445+099 out of the gain table for transfer and use all of the bandpass table. Also, note that the table .fluxscale has the .gcal solutions with the correct flux scale applied, and so there is no need to invoke the .gcal again in the applycal() command below.
<source lang="python">
applycal(vis='ngc5921.demo.ms', field='1,2', gaintable=['ngc5921.demo.fluxscale','ngc5921.demo.bcal'], gainfield=['1','*'],
        interp=['linear','nearest'], spwmap=[], selectdata=False)
</source>
Now for completeness apply 1331+305 to itself.
<source lang="python">
applycal(vis='ngc5921.demo.ms', field='0', gaintable=['ngc5921.demo.fluxscale','ngc5921.demo.bcal'], gainfield=['0','*'],
        interp=['linear','nearest'], spwmap=[], selectdata=False)
</source>
<!---
In the script snippet below, the python global variables ''ftable'' and ''btable'' replace the full table names.
In the script snippet below, the python global variables ''ftable'' and ''btable'' replace the full table names.


<source lang="python">
default('applycal')
default('applycal')
vis = msfile
vis = msfile
Line 512: Line 746:
selectdata = False
selectdata = False
# as before
# as before
gaincurve = False
opacity = 0.0
opacity = 0.0
# select the fields for 1445+099 and N5921
# select the fields for 1445+099 and N5921
Line 523: Line 756:
saveinputs('applycal',prefix+'.applycal.saved')
saveinputs('applycal',prefix+'.applycal.saved')
applycal()
applycal()
</source>
 
-->


== Plot the Spectrum ==
== Plot the Spectrum ==
Line 533: Line 767:
We can do this in {{plotms}}.
We can do this in {{plotms}}.


<source lang="bash">
<source lang="python">
plotms(vis=msfile, selectdata=True, field='N5921*', spw='0:6~56', \
plotms(vis='ngc5921.demo.ms', selectdata=True, field='N5921*', spw='0:6~56', averagedata=True, avgtime='3600', avgscan=True,  
      averagedata=True, avgtime='3600', avgscan=True, avgbaseline=True, \
      avgbaseline=True, xaxis='channel', yaxis='amp', ydatacolumn='corrected')
      xaxis='channel', yaxis='amp', ydatacolumn='corrected')
</source>
</source>


Line 548: Line 781:
* (Axes)'''Y Axis''' = Amp
* (Axes)'''Y Axis''' = Amp


From inspection of this plot, it looks like channels 4~6 and 50~59 contain line-free channels, suitable to use for continuum subtraction.
Now remove 0:6~56 from the spw field to see all channels. From inspection of this plot, it looks like channels 4~6 and 50~59 contain line-free channels, suitable to use for continuum subtraction.


== Continuum Subtraction ==
== Continuum Subtraction ==
Line 555: Line 788:


<source lang="python">
<source lang="python">
split(vis='ngc5921.demo.ms', outputvis='ngc5921.demo.src.split.ms', field='N5921*', spw='', datacolumn='corrected')
</source>
<!--
default('split')
default('split')
vis = msfile
vis = msfile
Line 565: Line 802:
split()
split()
print "Created "+splitms
print "Created "+splitms
</source>
-->


This action generated a new measurement set called ''ngc5921.demo.src.split.ms'' and copied the ''calibrated'' source data (datacolumn = 'corrected') into it.
This action generated a new measurement set called ''ngc5921.demo.src.split.ms'' and copied the ''calibrated'' source data (datacolumn = 'corrected') into it.


{{uvcontsub|Uvcontsub}} subtracts the continuum from the data in the visibility (''u'', ''v'') plane.  
{{uvcontsub|Uvcontsub}} subtracts the continuum from the data in the visibility (''u'', ''v'') plane. We will be using channels 4-6 and 50-59 for continuum.


<source lang="python">
<source lang="python">
uvcontsub(vis='ngc5921.demo.src.split.ms', field='N5921*', fitspw='0:4~6;50~59', spw='0', solint=0.0, fitorder=0, want_cont=True)
</source>
<!--
default('uvcontsub')
default('uvcontsub')
vis = splitms
vis = splitms
Line 590: Line 831:
uvcontsub()
uvcontsub()
</source>
</source>
-->
Notice that {{uvcontsub}} splits the data into two new measurement sets, 'ngc5921.demo.ms.cont', which contains an average of the continuum channels, and 'ngc5921.demo.ms.contsub', which contains the continuum-subtracted spectral line data.


Notice that {{uvcontsub}} splits two new measurement sets, 'ngc5921.demo.ms.cont', which contains an average of the continuum channels, and 'ngc5921.demo.ms.contsub', which contains the continuum-subtracted spectral line data.


== Imaging ==
== Imaging ==
[[File:spectrum-uvplot.png|thumb|Plot of amplitude vs. projected baseline length (in units of the observing wavelength) produced by casaplotms. The maximum baseline is just below 5 kilo-lambda.]]
[[File:spectrum-uvplot.png|thumb|Plot of amplitude vs. projected baseline length (in units of the observing wavelength) produced by casaplotms. The maximum baseline is just below 5 kilo-lambda.]]


 
Now we can generate the primary science product, a {{tclean}} data cube (ra, dec, velocity) from the continuum-subtracted (''u'', ''v'', channel) measurement set, ''ngc5921.demo.ms.contsub''. Things to consider in using {{tclean}}:
Now we can generate the primary science product, a {{clean}} data cube (ra, dec, velocity) from the continuum-subtracted (''u'', ''v'', channel) measurement set, ''ngc5921.demo.ms.contsub''. Things to consider in using {{clean}}:
* To ensure channels aren't averaged prior to imaging, choose '''mode='channel''''.
* To ensure channels aren't averaged prior to imaging, choose '''mode='channel''''.
* Specify the channels to image using '''start = 5''', '''width = 1''' (no averaging over channels), '''nchan = 46'''; only channels 5~51 will be imaged.
* Specify the channels to image using '''start = 5''', '''width = 1''' (no averaging over channels), '''nchan = 46'''; only channels 5~51 will be imaged.
* The maximum baseline is just under 5 kilolambda (see the figure at right), and so the expected synthetic beam is roughly 1.22 &times; 206265 / 5000 = 50 arcseconds (subject to the details of ''u'', ''v'' weighting). Pixels should sample the beam better than 3 times, so 15 arcseconds is a good choice of pixel size ('''cell = ['15.0arcsec','15.0arcsec']''').
* The maximum baseline is just under 5 kilolambda (see the figure at right), and so the expected synthetic beam is roughly 1.22 &times; 206265 / 5000 = 50 arcseconds (subject to the details of ''u'', ''v'' weighting). Pixels should sample the beam better than 3 times, so 15 arcseconds is a good choice of pixel size ('''cell = ['15.0arcsec','15.0arcsec']''').
* We only want to {{clean}} down to the noise, which is easily determined by trial-and-error imaging of a single channel (choosing '''nchan=1''' and '''start''' appropriately). Here, {{clean}} stops when the maximum residual on the channel is below '''threshold='8.0mJy''''.  
* We only want to {{tclean}} down to the noise, which is easily determined by trial-and-error imaging of a single channel (choosing '''nchan=1''' and '''start''' appropriately). Here, {{tclean}} stops when the maximum residual on the channel is below '''threshold='8.0mJy''''.  


<!--
<source lang="python">
<source lang="python">
# identify the continuum subtracted measurement set
# Image the continuum subtracted measurement set
clean(vis='ngc5921.demo.src.split.ms.contsub', imagename='ngc5921.demo.clean', field='0', mode='channel', nchan=46, start=5, width=1, spw='', gain=0.1, imsize=[256,256], psfmode='clark', imagermode='', cell=['15.0arcsec','15.0arcsec'], niter=6000, threshold='8.0mJy', weighting='briggs', robust=0.5, mask = [108,108,148,148], interactive=False)
</source>
-->
 
<!-- As of CASA 5.3, the preferred and recommended deconvolution task is {{tclean}}(). The {{tclean}}() execution equivalent to the above {{clean}}() one is following: -->
 
<source lang="python">
# Image the continuum subtracted measurement set
tclean(vis='ngc5921.demo.src.split.ms.contsub', imagename='ngc5921.demo.tclean', field='0', datacolumn='data', specmode='cube', nchan=46, start=5, width=1, spw='', deconvolver='hogbom', gridder='standard', niter=6000, gain=0.1, threshold='8.0mJy', imsize=[256,256], cell=['15.0arcsec','15.0arcsec'], weighting='briggs', robust=0.5,  mask = 'box[[108pix,108pix],[148pix,148pix]]', interactive=False, pblimit=-0.2)
</source>
 
<!--
Execute only one, either {{tclean}}() or {{clean}}().
-->
 
<!--
srcsplitms = splitms + '.contsub'
srcsplitms = splitms + '.contsub'
#=====================================================================
#=====================================================================
Line 650: Line 913:
saveinputs('clean',prefix+'.clean.saved')
saveinputs('clean',prefix+'.clean.saved')
clean()
clean()
</source>
-->


Use {{imhead}} to look at the cube header:
Use {{imhead}} to look at the cube header:


<source lang="python">
<source lang="python">
imhead(imagename='ngc5921.demo.tclean.image', mode='summary')
</source>
<!--
clnimage = imname+'.image' # store the clean image in a python global for future use
clnimage = imname+'.image' # store the clean image in a python global for future use
default('imhead')
default('imhead')
Line 660: Line 927:
mode = 'summary'
mode = 'summary'
imhead()
imhead()
</source>
-->


The output, as follows, appears in the logger window.
The output, as follows, appears in the logger window.
<!--


<pre>
<pre>
2018-09-13 00:26:03 INFO imhead ##########################################
2018-09-13 00:26:03 INFO imhead ##### Begin Task: imhead            #####
2018-09-13 00:26:03 INFO imhead imhead(imagename="ngc5921.demo.clean.image",mode="summary",hdkey="",hdvalue="",verbose=False)
2018-09-13 00:26:03 INFO ImageMetaData  
2018-09-13 00:26:03 INFO ImageMetaData Image name      : ngc5921.demo.clean.image
2018-09-13 00:26:03 INFO ImageMetaData Object name      : N5921_2
2018-09-13 00:26:03 INFO ImageMetaData Image type      : PagedImage
2018-09-13 00:26:03 INFO ImageMetaData Image quantity  : Intensity
2018-09-13 00:26:03 INFO ImageMetaData Pixel mask(s)    : None
2018-09-13 00:26:03 INFO ImageMetaData Region(s)        : None
2018-09-13 00:26:03 INFO ImageMetaData Image units      : Jy/beam
2018-09-13 00:26:03 INFO ImageMetaData Restoring Beam  : 51.5763 arcsec, 47.3275 arcsec, 8.41586 deg
2018-09-13 00:26:03 INFO ImageMetaData  
2018-09-13 00:26:03 INFO ImageMetaData Direction reference : J2000
2018-09-13 00:26:03 INFO ImageMetaData Spectral  reference : LSRK
2018-09-13 00:26:03 INFO ImageMetaData Velocity  type      : RADIO
2018-09-13 00:26:03 INFO ImageMetaData Rest frequency      : 1.42041e+09 Hz
2018-09-13 00:26:03 INFO ImageMetaData Pointing center    :  15:22:00.000000  +05.04.00.000000
2018-09-13 00:26:03 INFO ImageMetaData Telescope          : VLA
2018-09-13 00:26:03 INFO ImageMetaData Observer            : TEST
2018-09-13 00:26:03 INFO ImageMetaData Date observation    : 1995/04/13/09:33:00
2018-09-13 00:26:03 INFO ImageMetaData Telescope position: [-1.60119e+06m, -5.04198e+06m, 3.55488e+06m] (ITRF)
2018-09-13 00:26:03 INFO ImageMetaData  
2018-09-13 00:26:03 INFO ImageMetaData Axis Coord Type      Name            Proj Shape Tile  Coord value at pixel    Coord incr Units
2018-09-13 00:26:03 INFO ImageMetaData ------------------------------------------------------------------------------------------------
2018-09-13 00:26:03 INFO ImageMetaData 0    0    Direction Right Ascension  SIN  256  64  15:22:00.000  128.00 -1.500000e+01 arcsec
2018-09-13 00:26:03 INFO ImageMetaData 1    0    Direction Declination      SIN  256  64 +05.04.00.000  128.00  1.500000e+01 arcsec
2018-09-13 00:26:03 INFO ImageMetaData 2    1    Stokes    Stokes                    1    1            I
2018-09-13 00:26:03 INFO ImageMetaData 3    2    Spectral  Frequency                46    8  1.41279e+09    0.00 2.4414062e+04 Hz
2018-09-13 00:26:03 INFO ImageMetaData                     Velocity                              1607.99    0.00 -5.152860e+00 km/s
2018-09-13 00:26:03 INFO imhead ##### End Task: imhead              #####
2018-09-13 00:26:03 INFO imhead ##########################################
</pre>
-->
<pre>
2019-09-03 22:25:25 INFO imhead ##### Begin Task: imhead            #####
2019-09-03 22:25:25 INFO imhead imhead(imagename="ngc5921.demo.tclean.image",mode="summary",hdkey="",hdvalue="",verbose=False)
2019-09-03 22:25:25 INFO ImageMetaData  
2019-09-03 22:25:25 INFO ImageMetaData Image name      : ngc5921.demo.tclean.image
2019-09-03 22:25:25 INFO ImageMetaData Object name      : N5921_2
2019-09-03 22:25:25 INFO ImageMetaData Image type      : PagedImage
2019-09-03 22:25:25 INFO ImageMetaData Image quantity  : Intensity
2019-09-03 22:25:25 INFO ImageMetaData Pixel mask(s)    : None
2019-09-03 22:25:25 INFO ImageMetaData Region(s)        : None
2019-09-03 22:25:25 INFO ImageMetaData Image units      : Jy/beam
2019-09-03 22:25:25 INFO ImageMetaData Restoring Beams
2019-09-03 22:25:25 INFO ImageMetaData Pol  Type Chan        Freq    Vel
2019-09-03 22:25:25 INFO ImageMetaData   I    Max    7 1.41296e+09 1571.92  51.6639 arcsec x  47.3594 arcsec pa=  8.3205 deg
2019-09-03 22:25:25 INFO ImageMetaData   I    Min  43 1.41384e+09 1386.42  51.5897 arcsec x  47.3127 arcsec pa=  8.3869 deg
2019-09-03 22:25:25 INFO ImageMetaData   I Median  13 1.41310e+09 1541.00  51.6611 arcsec x  47.3318 arcsec pa=  8.1425 deg
2019-09-03 22:25:25 INFO ImageMetaData  
2019-09-03 22:25:25 INFO ImageMetaData Direction reference : J2000
2019-09-03 22:25:25 INFO ImageMetaData Spectral  reference : LSRK
2019-09-03 22:25:25 INFO ImageMetaData Velocity  type      : RADIO
2019-09-03 22:25:25 INFO ImageMetaData Rest frequency      : 1.42041e+09 Hz
2019-09-03 22:25:25 INFO ImageMetaData Pointing center    :  15:22:00.000000  +05.04.00.000000
2019-09-03 22:25:25 INFO ImageMetaData Telescope          : VLA
2019-09-03 22:25:25 INFO ImageMetaData Observer            : TEST
2019-09-03 22:25:25 INFO ImageMetaData Date observation    : 1995/04/13/09:33:00
2019-09-03 22:25:25 INFO ImageMetaData Telescope position: [-1.60119e+06m, -5.04198e+06m, 3.55488e+06m] (ITRF)
2019-09-03 22:25:25 INFO ImageMetaData  
2019-09-03 22:25:25 INFO ImageMetaData Axis Coord Type      Name            Proj Shape Tile  Coord value at pixel    Coord incr Units
2019-09-03 22:25:25 INFO ImageMetaData ------------------------------------------------------------------------------------------------
2019-09-03 22:25:25 INFO ImageMetaData 0    0    Direction Right Ascension  SIN  256  64  15:22:00.000  128.00 -1.500000e+01 arcsec
2019-09-03 22:25:25 INFO ImageMetaData 1    0    Direction Declination      SIN  256  64 +05.04.00.000  128.00  1.500000e+01 arcsec
2019-09-03 22:25:25 INFO ImageMetaData 2    1    Stokes    Stokes                    1    1            I
2019-09-03 22:25:25 INFO ImageMetaData 3    2    Spectral  Frequency                46    8  1.41279e+09    0.00 2.4414062e+04 Hz
2019-09-03 22:25:25 INFO ImageMetaData                     Velocity                              1607.99    0.00 -5.152860e+00 km/s
2019-09-03 22:25:25 INFO imhead ##### End Task: imhead              #####
2019-09-03 22:25:25 INFO imhead ##########################################
</pre>
<!--
2011-04-25 20:10:29 INFO imhead ##########################################
2011-04-25 20:10:29 INFO imhead ##########################################
2011-04-25 20:10:29 INFO imhead ##### Begin Task: imhead            #####
2011-04-25 20:10:29 INFO imhead ##### Begin Task: imhead            #####
Line 698: Line 1,043:
2011-04-25 20:10:29 INFO imhead ##### End Task: imhead              #####
2011-04-25 20:10:29 INFO imhead ##### End Task: imhead              #####
2011-04-25 20:10:29 INFO imhead ##########################################
2011-04-25 20:10:29 INFO imhead ##########################################
</pre>
</pre>-->
 


== Additional Science Products ==
== Additional Science Products ==


If things went well, you should now have a spectral line cube (''ngc5921.demo.cleanimg.img'') as a primary science product. The demo script illustrates further how to generate cube statistics (using {{imstat}}), an integrated spectrum, and moment maps.  
If things went well, you should now have a spectral line cube (''ngc5921.demo.tclean.image'') as a primary science product. The demo script illustrates further how to generate cube statistics (using {{imstat}}), an integrated spectrum, and moment maps.  


=== Cube Statistics ===
=== Cube Statistics ===


{{imstat|Imstat}} is the tool for displaying statistics of images and cubes. The following example displays the statistics for the whole cube.
{{imstat|Imstat}} is the tool for displaying statistics of images and cubes. The following example displays the statistics for an empty region of the whole cube.


<!--
<source lang="python">
<source lang="python">
cubestat=imstat(imagename='ngc5921.demo.clean.image', box='10,10,100,100')
</source>
default('imstat')
default('imstat')
imagename = clnimage # or imagename = "ngc5921.demo.cleanimg.img"
imagename = clnimage # or imagename = "ngc5921.demo.cleanimg.img"
Line 716: Line 1,067:
#box = '108,108,148,148'
#box = '108,108,148,148'
cubestats = imstat()
cubestats = imstat()
The output goes to the logger window.
<pre>
2018-09-13 23:55:54 INFO imstat ##########################################
2018-09-13 23:55:54 INFO imstat ##### Begin Task: imstat            #####
2018-09-13 23:55:54 INFO imstat imstat(imagename="ngc5921.demo.clean.image",axes=-1,region="",box="10,10,100,100",chans="",
2018-09-13 23:55:54 INFO imstat         stokes="",listit=True,verbose=True,mask="",stretch=False,
2018-09-13 23:55:54 INFO imstat         logfile="",append=True,algorithm="classic",fence=-1,center="mean",
2018-09-13 23:55:54 INFO imstat         lside=True,zscore=-1,maxiter=-1,clmethod="auto",niter=3)
2018-09-13 23:55:54 INFO imstat Using specified box(es) 10,10,100,100
2018-09-13 23:55:54 INFO imstat Determining stats for image ngc5921.demo.clean.image
2018-09-13 23:55:54 INFO imstat Selected bounding box :
2018-09-13 23:55:54 INFO imstat     [10, 10, 0, 0] to [100, 100, 0, 45]  (15:23:58.379, +04.34.29.305, I, 1.41279e+09Hz to 15:22:28.105, +04.56.59.962, I, 1.41389e+09Hz)
2018-09-13 23:55:54 INFO imstat Statistics calculated using Classic algorithm
2018-09-13 23:55:54 INFO imstat Regions ---
2018-09-13 23:55:54 INFO imstat         -- bottom-left corner (pixel) [blc]:  [10, 10, 0, 0]
2018-09-13 23:55:54 INFO imstat         -- top-right corner (pixel) [trc]:    [100, 100, 0, 45]
2018-09-13 23:55:54 INFO imstat         -- bottom-left corner (world) [blcf]: 15:23:58.379, +04.34.29.305, I, 1.41279e+09Hz
2018-09-13 23:55:54 INFO imstat         -- top-right corner (world) [trcf]:  15:22:28.105, +04.56.59.962, I, 1.41389e+09Hz
2018-09-13 23:55:54 INFO imstat Values ---
2018-09-13 23:55:54 INFO imstat         -- flux [flux]:                            1.98066 Jy.km/s
2018-09-13 23:55:54 INFO imstat         -- number of points [npts]:                380926
2018-09-13 23:55:54 INFO imstat         -- maximum value [max]:                    0.00952374 Jy/beam
2018-09-13 23:55:54 INFO imstat         -- minimum value [min]:                    -0.0100571 Jy/beam
2018-09-13 23:55:54 INFO imstat         -- position of max value (pixel) [maxpos]: [85, 63, 0, 8]
2018-09-13 23:55:54 INFO imstat         -- position of min value (pixel) [minpos]: [30, 18, 0, 7]
2018-09-13 23:55:54 INFO imstat         -- position of max value (world) [maxposf]: 15:22:43.151, +04.47.44.907, I, 1.41298e+09Hz
2018-09-13 23:55:54 INFO imstat         -- position of min value (world) [minposf]: 15:23:38.319, +04.36.29.518, I, 1.41296e+09Hz
2018-09-13 23:55:54 INFO imstat         -- Sum of pixel values [sum]:              4.72505 Jy/beam
2018-09-13 23:55:54 INFO imstat         -- Sum of squared pixel values [sumsq]:    1.38067 Jy/beam.Jy/beam
2018-09-13 23:55:54 INFO imstat Statistics ---
2018-09-13 23:55:54 INFO imstat         -- Mean of the pixel values [mean]:        1.24041e-05 Jy/beam
2018-09-13 23:55:54 INFO imstat         -- Variance of the pixel values :          3.62436e-06 Jy/beam
2018-09-13 23:55:54 INFO imstat         -- Standard deviation of the Mean [sigma]:  0.00190378 Jy/beam
2018-09-13 23:55:54 INFO imstat         -- Root mean square [rms]:                  0.00190381 Jy/beam
2018-09-13 23:55:54 INFO imstat         -- Median of the pixel values [median]:    6.54527e-06 Jy/beam
2018-09-13 23:55:54 INFO imstat         -- Median of the deviations [medabsdevmed]: 0.00127726 Jy/beam
2018-09-13 23:55:54 INFO imstat         -- IQR [quartile]:                          0.00255404 Jy/beam
2018-09-13 23:55:54 INFO imstat         -- First quartile [q1]:                    -0.00126623 Jy/beam
2018-09-13 23:55:54 INFO imstat         -- Third quartile [q3]:                    0.0012878 Jy/beam
2018-09-13 23:55:54 INFO imstat Created Temp image  of shape [1, 1, 1, 1] with float valued pixels.
2018-09-13 23:55:54 INFO imstat Sum column unit = Jy/beam
2018-09-13 23:55:54 INFO imstat Mean column unit = Jy/beam
2018-09-13 23:55:54 INFO imstat Std_dev column unit = Jy/beam
2018-09-13 23:55:54 INFO imstat Minimum column unit = Jy/beam
2018-09-13 23:55:54 INFO imstat Maximum column unit = Jy/beam
2018-09-13 23:55:54 INFO imstat Npts          Sum          Mean          Rms          Std_dev      Minimum      Maximum   
2018-09-13 23:55:54 INFO imstat 3.809260e+05  4.725051e+00  1.240412e-05  1.903813e-03  1.903775e-03 -1.005711e-02  9.523738e-03
2018-09-13 23:55:54 INFO imstat ##### End Task: imstat              #####
2018-09-13 23:55:54 INFO imstat ##########################################
</pre>
If you executed {{tclean}}() instead of {{clean}}() in previous steps, then your image statistics may look like this:
-->
<source lang="python">
cubestat=imstat(imagename='ngc5921.demo.tclean.image', box='10,10,100,100')
</source>
</source>


The output goes to the logger window.
<!--
<pre>
2018-09-13 23:57:15 INFO imstat ##########################################
2018-09-13 23:57:15 INFO imstat ##### Begin Task: imstat            #####
2018-09-13 23:57:15 INFO imstat imstat(imagename="ngc5921.demo.tclean.image",axes=-1,region="",box="10,10,100,100",chans="",
2018-09-13 23:57:15 INFO imstat         stokes="",listit=True,verbose=True,mask="",stretch=False,
2018-09-13 23:57:15 INFO imstat         logfile="",append=True,algorithm="classic",fence=-1,center="mean",
2018-09-13 23:57:15 INFO imstat         lside=True,zscore=-1,maxiter=-1,clmethod="auto",niter=3)
2018-09-13 23:57:15 INFO imstat Using specified box(es) 10,10,100,100
2018-09-13 23:57:15 INFO imstat Determining stats for image ngc5921.demo.tclean.image
2018-09-13 23:57:15 INFO imstat Selected bounding box :
2018-09-13 23:57:15 INFO imstat     [10, 10, 0, 0] to [100, 100, 0, 45]  (15:23:58.379, +04.34.29.305, I, 1.41279e+09Hz to 15:22:28.105, +04.56.59.962, I, 1.41389e+09Hz)
2018-09-13 23:57:15 INFO imstat Statistics calculated using Classic algorithm
2018-09-13 23:57:15 INFO imstat Regions ---
2018-09-13 23:57:15 INFO imstat         -- bottom-left corner (pixel) [blc]:  [10, 10, 0, 0]
2018-09-13 23:57:15 INFO imstat         -- top-right corner (pixel) [trc]:    [100, 100, 0, 45]
2018-09-13 23:57:15 INFO imstat         -- bottom-left corner (world) [blcf]: 15:23:58.379, +04.34.29.305, I, 1.41279e+09Hz
2018-09-13 23:57:15 INFO imstat         -- top-right corner (world) [trcf]:  15:22:28.105, +04.56.59.962, I, 1.41389e+09Hz
2018-09-13 23:57:15 INFO imstat Values ---
2018-09-13 23:57:15 INFO imstat         -- flux [flux]:                            1.63176 Jy.km/s
2018-09-13 23:57:15 INFO imstat         -- number of points [npts]:                380926
2018-09-13 23:57:15 INFO imstat         -- maximum value [max]:                    0.00823404 Jy/beam
2018-09-13 23:57:15 INFO imstat         -- minimum value [min]:                    -0.00943008 Jy/beam
2018-09-13 23:57:15 INFO imstat         -- position of max value (pixel) [maxpos]: [61, 54, 0, 42]
2018-09-13 23:57:15 INFO imstat         -- position of min value (pixel) [minpos]: [30, 18, 0, 7]
2018-09-13 23:57:15 INFO imstat         -- position of max value (world) [maxposf]: 15:23:07.232, +04.45.29.778, I, 1.41381e+09Hz
2018-09-13 23:57:15 INFO imstat         -- position of min value (world) [minposf]: 15:23:38.319, +04.36.29.518, I, 1.41296e+09Hz
2018-09-13 23:57:15 INFO imstat         -- Sum of pixel values [sum]:              3.89979 Jy/beam
2018-09-13 23:57:15 INFO imstat         -- Sum of squared pixel values [sumsq]:    1.23177 Jy/beam.Jy/beam
2018-09-13 23:57:15 INFO imstat Statistics ---
2018-09-13 23:57:15 INFO imstat         -- Mean of the pixel values [mean]:        1.02377e-05 Jy/beam
2018-09-13 23:57:15 INFO imstat         -- Variance of the pixel values :          3.23353e-06 Jy/beam
2018-09-13 23:57:15 INFO imstat         -- Standard deviation of the Mean [sigma]:  0.0017982 Jy/beam
2018-09-13 23:57:15 INFO imstat         -- Root mean square [rms]:                  0.00179823 Jy/beam
2018-09-13 23:57:15 INFO imstat         -- Median of the pixel values [median]:    1.07091e-05 Jy/beam
2018-09-13 23:57:15 INFO imstat         -- Median of the deviations [medabsdevmed]: 0.00120728 Jy/beam
2018-09-13 23:57:15 INFO imstat         -- IQR [quartile]:                          0.00241453 Jy/beam
2018-09-13 23:57:15 INFO imstat         -- First quartile [q1]:                    -0.00119593 Jy/beam
2018-09-13 23:57:15 INFO imstat         -- Third quartile [q3]:                    0.00121859 Jy/beam
2018-09-13 23:57:15 INFO imstat Created Temp image  of shape [1, 1, 1, 1] with float valued pixels.
2018-09-13 23:57:15 INFO imstat Sum column unit = Jy/beam
2018-09-13 23:57:15 INFO imstat Mean column unit = Jy/beam
2018-09-13 23:57:15 INFO imstat Std_dev column unit = Jy/beam
2018-09-13 23:57:15 INFO imstat Minimum column unit = Jy/beam
2018-09-13 23:57:15 INFO imstat Maximum column unit = Jy/beam
2018-09-13 23:57:15 INFO imstat Npts          Sum          Mean          Rms          Std_dev      Minimum      Maximum   
2018-09-13 23:57:15 INFO imstat 3.809260e+05  3.899789e+00  1.023766e-05  1.798228e-03  1.798201e-03 -9.430082e-03  8.234036e-03
2018-09-13 23:57:15 INFO imstat ##### End Task: imstat              #####
2018-09-13 23:57:15 INFO imstat ##########################################
</pre>
-->


<pre>
<pre>
2019-09-03 22:40:15 INFO imstat ##########################################
2019-09-03 22:40:15 INFO imstat ##### Begin Task: imstat            #####
2019-09-03 22:40:15 INFO imstat imstat(imagename="ngc5921.demo.tclean.image",axes=-1,region="",box="10,10,100,100",chans="",
2019-09-03 22:40:15 INFO imstat         stokes="",listit=True,verbose=True,mask="",stretch=False,
2019-09-03 22:40:15 INFO imstat         logfile="",append=True,algorithm="classic",fence=-1,center="mean",
2019-09-03 22:40:15 INFO imstat         lside=True,zscore=-1,maxiter=-1,clmethod="auto",niter=3)
2019-09-03 22:40:15 INFO imstat Using specified box(es) 10,10,100,100
2019-09-03 22:40:15 INFO imstat Determining stats for image ngc5921.demo.tclean.image
2019-09-03 22:40:15 INFO imstat Selected bounding box :
2019-09-03 22:40:15 INFO imstat     [10, 10, 0, 0] to [100, 100, 0, 45]  (15:23:58.379, +04.34.29.305, I, 1.41279e+09Hz to 15:22:28.105, +04.56.59.962, I, 1.41389e+09Hz)
2019-09-03 22:40:15 INFO imstat Statistics calculated using Classic algorithm
2019-09-03 22:40:15 INFO imstat Regions ---
2019-09-03 22:40:15 INFO imstat         -- bottom-left corner (pixel) [blc]:  [10, 10, 0, 0]
2019-09-03 22:40:15 INFO imstat         -- top-right corner (pixel) [trc]:    [100, 100, 0, 45]
2019-09-03 22:40:15 INFO imstat         -- bottom-left corner (world) [blcf]: 15:23:58.379, +04.34.29.305, I, 1.41279e+09Hz
2019-09-03 22:40:15 INFO imstat         -- top-right corner (world) [trcf]:  15:22:28.105, +04.56.59.962, I, 1.41389e+09Hz
2019-09-03 22:40:15 INFO imstat Values ---
2019-09-03 22:40:15 INFO imstat         -- flux [flux]:                            1.99425 Jy.km/s
2019-09-03 22:40:15 INFO imstat         -- number of points [npts]:                380926
2019-09-03 22:40:15 INFO imstat         -- maximum value [max]:                    0.00953276 Jy/beam
2019-09-03 22:40:15 INFO imstat         -- minimum value [min]:                    -0.0100478 Jy/beam
2019-09-03 22:40:15 INFO imstat         -- position of max value (pixel) [maxpos]: [85, 63, 0, 8]
2019-09-03 22:40:15 INFO imstat         -- position of min value (pixel) [minpos]: [30, 18, 0, 7]
2019-09-03 22:40:15 INFO imstat         -- position of max value (world) [maxposf]: 15:22:43.151, +04.47.44.907, I, 1.41298e+09Hz
2019-09-03 22:40:15 INFO imstat         -- position of min value (world) [minposf]: 15:23:38.319, +04.36.29.518, I, 1.41296e+09Hz
2019-09-03 22:40:15 INFO imstat         -- Sum of pixel values [sum]:              4.76561 Jy/beam
2019-09-03 22:40:15 INFO imstat         -- Sum of squared pixel values [sumsq]:    1.38125 Jy/beam.Jy/beam
2019-09-03 22:40:15 INFO imstat Statistics ---
2019-09-03 22:40:15 INFO imstat         -- Mean of the pixel values [mean]:        1.25106e-05 Jy/beam
2019-09-03 22:40:15 INFO imstat         -- Variance of the pixel values :          3.62589e-06 Jy/beam
2019-09-03 22:40:15 INFO imstat         -- Standard deviation of the Mean [sigma]:  0.00190418 Jy/beam
2019-09-03 22:40:15 INFO imstat         -- Root mean square [rms]:                  0.00190421 Jy/beam
2019-09-03 22:40:15 INFO imstat         -- Median of the pixel values [median]:    6.37541e-06 Jy/beam
2019-09-03 22:40:15 INFO imstat         -- Median of the deviations [medabsdevmed]: 0.00127676 Jy/beam
2019-09-03 22:40:15 INFO imstat         -- IQR [quartile]:                          0.00255348 Jy/beam
2019-09-03 22:40:15 INFO imstat         -- First quartile [q1]:                    -0.00126604 Jy/beam
2019-09-03 22:40:15 INFO imstat         -- Third quartile [q3]:                    0.00128744 Jy/beam
2019-09-03 22:40:15 INFO imstat Sum column unit = Jy/beam
2019-09-03 22:40:15 INFO imstat Mean column unit = Jy/beam
2019-09-03 22:40:15 INFO imstat Std_dev column unit = Jy/beam
2019-09-03 22:40:15 INFO imstat Minimum column unit = Jy/beam
2019-09-03 22:40:15 INFO imstat Maximum column unit = Jy/beam
2019-09-03 22:40:15 INFO imstat Npts          Sum          Mean          Rms          Std_dev      Minimum      Maximum   
2019-09-03 22:40:15 INFO imstat 3.809260e+05  4.765615e+00  1.251060e-05  1.904215e-03  1.904176e-03 -1.004781e-02  9.532760e-03
2019-09-03 22:40:15 INFO imstat ##### End Task: imstat              #####
2019-09-03 22:40:15 INFO imstat ##########################################
</pre>
<!--
2011-04-25 20:12:13 INFO imstat ##########################################
2011-04-25 20:12:13 INFO imstat ##########################################
2011-04-25 20:12:13 INFO imstat ##### Begin Task: imstat            #####
2011-04-25 20:12:13 INFO imstat ##### Begin Task: imstat            #####
Line 727: Line 1,238:
2011-04-25 20:12:13 INFO imstat Using all spectral channels.
2011-04-25 20:12:13 INFO imstat Using all spectral channels.
2011-04-25 20:12:13 INFO imstat Using polarizations I
2011-04-25 20:12:13 INFO imstat Using polarizations I
2011-04-25 20:12:13 INFO imstat Determining stats for image ngc5921.demo.cleanimg.image
2011-04-25 20:12:13 INFO imstat Determining stats for image ngc5921.demo.clean.image
2011-04-25 20:12:13 INFO imstat Selected bounding box :  
2011-04-25 20:12:13 INFO imstat Selected bounding box :  
2011-04-25 20:12:13 INFO imstat     [0, 0, 0, 0] to [255, 255, 0, 45]  (15:24:08.404, +04.31.59.181, I, 1.41279e+09Hz to 15:19:52.390, +05.35.44.246, I, 1.41389e+09Hz)
2011-04-25 20:12:13 INFO imstat     [0, 0, 0, 0] to [255, 255, 0, 45]  (15:24:08.404, +04.31.59.181, I, 1.41279e+09Hz to 15:19:52.390, +05.35.44.246, I, 1.41389e+09Hz)
Line 758: Line 1,269:
2011-04-25 20:12:13 INFO imstat ##### End Task: imstat              #####
2011-04-25 20:12:13 INFO imstat ##### End Task: imstat              #####
2011-04-25 20:12:13 INFO imstat ##########################################
2011-04-25 20:12:13 INFO imstat ##########################################
</pre>
</pre> -->


=== The Integrated Spectrum ===
=== The Integrated Spectrum ===
Line 768: Line 1,279:


<source lang="python">
<source lang="python">
viewer(infile='ngc5921.demo.tclean.image')
</source>
<!--
# Store the name of the clean image into a python global
# Store the name of the clean image into a python global
clnimage = imname+'.image'
clnimage = imname+'.image'
viewer(clnimage)
viewer(clnimage)
</source>
-->


To generate the integrated spectrum, perform the following tasks.
To generate the integrated spectrum, perform the following tasks.
Line 780: Line 1,295:
* The spectrum now appears in the graphics window; see the figure at right.
* The spectrum now appears in the graphics window; see the figure at right.


You can save the integrated spectrum to a text file by clicking the [[File:save-as-text-file.png]] button on the graphics window. There are also buttons to print the figure or save the figure to disk.
You can save the integrated spectrum to a text file by clicking the [[File:save_to_text.png]] button on the graphics window. There are also buttons to print the figure or save the figure to disk.
 


=== Cube Moments ===
=== Cube Moments ===
Line 788: Line 1,304:
Cube moments are maps of weighted sums along the velocity axis. In CASA, they are generated by the task {{immoments}}. The zeroth moment (moments = 0) is a sum of intensities along the velocity axis (the integrated intensity map); the first moment (moment = 1) is the sum of velocities weighted by intensity (the ''velocity field''); the second moment (moment = 2) is a map of the velocity dispersion; see the {{immoments|helpfile}} for additional options.
Cube moments are maps of weighted sums along the velocity axis. In CASA, they are generated by the task {{immoments}}. The zeroth moment (moments = 0) is a sum of intensities along the velocity axis (the integrated intensity map); the first moment (moment = 1) is the sum of velocities weighted by intensity (the ''velocity field''); the second moment (moment = 2) is a map of the velocity dispersion; see the {{immoments|helpfile}} for additional options.


The following example produces maps of the zeroth and first moments, or the integrated intensity and velocity field. The respective measurement sets are ''ngc5921.demo.moments.integrated'' and ''ngc5921.demo.moments.weighted_coord'', stored in the python globals ''momzeroimage'' and ''momoneimage''.
The following example produces maps of the zeroth and first moments, or the integrated intensity and velocity field. The respective measurement sets are the moment zero image ''ngc5921.demo.moments.integrated'' and moment one image''ngc5921.demo.moments.weighted_coord''.


We will do the zeroth and first moments and mask out noisy pixels using hard global limits. We will also collapse along the spectral (channel) axis and include all planes.


<source lang="python">
<source lang="python">
immoments(imagename='ngc5921.demo.tclean.image', moments=[0,1], excludepix=[-100, 0.009], axis='spectral', chans='', outfile='ngc5921.demo.moments')
</source>
* moments = [0,1] : Do zeroth and first moments
* excludepix = [-100,0.009] : Mask out noisy pixels using hard global limits
* axis  = 'spectral' : Collapse along the spectral (channel) axis
* chans = '' :Include all planes
<!--
default('immoments')
default('immoments')
imagename = clnimage
imagename = clnimage
Line 810: Line 1,336:
momzeroimage = momfile + '.integrated'
momzeroimage = momfile + '.integrated'
momoneimage = momfile + '.weighted_coord'
momoneimage = momfile + '.weighted_coord'
</source>
-->


To examine the moment images, use {{viewer}}; the resulting moment zero image is displayed at right.
To examine the moment images, use {{viewer}}; the resulting moment zero image is displayed at right. Note that you may have to play with the color map (Data Display Options button in viewer) in order to replicate the image in this tutorial.


<source lang="python">
<source lang="python">
viewer(momzeroimage)
viewer(infile='ngc5921.demo.moments.integrated')
</source>
</source>


Line 822: Line 1,348:
To export the (''u'', ''v'') data and image cube as FITS files, use {{exportuvfits}} and {{exportfits}}, respectively.
To export the (''u'', ''v'') data and image cube as FITS files, use {{exportuvfits}} and {{exportfits}}, respectively.


Here's how to export the continuum-subtracted (''u'', ''v'') data. Note that this snippet has {{exportuvfits}} spawned to the background (async = True).
Here's how to export the continuum-subtracted (''u'', ''v'') data.


<source lang="python">
<source lang="python">
exportuvfits(vis='ngc5921.demo.src.split.ms.contsub', fitsfile='ngc5921.demo.contsub.uvfits', datacolumn='corrected', multisource=True)
</source>
<!--
default('exportuvfits')
default('exportuvfits')
srcuvfits = prefix + '.contsub.uvfits' # recall: prefix = 'ngc5921.demo'
srcuvfits = prefix + '.contsub.uvfits' # recall: prefix = 'ngc5921.demo'
Line 833: Line 1,363:
async = True
async = True
myhandle = exportuvfits()
myhandle = exportuvfits()
</source>
-->


And now, the FITS cube.
And now, the FITS cube.
<source lang="python">
<source lang="python">
exportfits(imagename='ngc5921.demo.tclean.image', fitsimage='ngc5921.demo.tclean.fits')
</source>
<!--
default('exportfits')
default('exportfits')
clnfits = prefix + '.cleanimg.fits'
clnfits = prefix + '.cleanimg.fits'
Line 844: Line 1,378:
saveinputs('exportfits',prefix+'.exportfits.saved')
saveinputs('exportfits',prefix+'.exportfits.saved')
myhandle2 = exportfits()
myhandle2 = exportfits()
</source>
-->


The moment maps (or any CASA images) can be similarly exported using {{exportfits}}.
The moment maps (or any CASA images) can be similarly exported using {{exportfits}}.


== Appendix: Python Notes ==
== Appendix: Python Notes ==
Line 877: Line 1,413:
</pre>
</pre>


[[VLA Tutorials | &#8629; '''VLA Tutorials''']]
[[Pre-upgrade VLA Tutorials | &#8629; '''Pre-upgrade VLA Tutorials''']]
 
{{Checked 4.7.2}}
 
--[[User:Jgallimo|Jack Gallimore]] 16:58, 7 December 2009 (UTC)


Updated: [[User:Mkrauss|Miriam Krauss]], 25 April 2011
{{Checked 5.7.2}}

Latest revision as of 17:08, 13 February 2023

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

Last checked on CASA Version 5.7.2.

Overview

The technique used to calibrate and image continuum datasets generally applies to spectral line observations, except that an additional calibration step is required. Bandpass calibration flattens the spectral response of the observations, ensuring that spectral channel images are properly calibrated in amplitude and phase.

The following tutorial derives from an annotated script provided in the CASA Cookbook. The script is largely reproduced and additionally annotated with figures and illustrations. It is assumed that this tutorial will be used interactively, and so interactive pauses in the original script have been removed.

DATA: The data are included with the CASA installation.


Setting up the CASA Environment

Start up CASA in the directory you want to use. 

# in bash
mkdir NGC5921
cd NGC5921
casa


We'll use a python os command to get the appropriate CASA path for your installation in order to import the data. The use of os.environ.get is explained in the Appendix. 

# In CASA
%cpaste

# Press Enter or Return, then copy/paste the following:
import os
pathname=os.environ.get('CASAPATH').split()[0]
fitsdata=pathname+'/data/demo/NGC5921.fits'
--

Scripts are of course modified and repeated to the satisfaction of observer. To help clean up the bookkeeping and further avoid issues of write privileges, remove prior versions of the measurement set and calibration tables.

This can be done with the rmtables('table_name') command.


Import the Data

The next step is to import the multisource UVFITS data to a CASA measurement set via the importuvfits filler. Note that you can set each parameter for any particular task one-by-one, or you could supply the task and input parameters with one command. Here we will set each parameter value first, save them, and run the import task. Throughout the remaining tutorial, we will call upon tasks with a single command. 


# Safest to start from task defaults
default('importuvfits')
# Use task importuvfits
fitsfile = fitsdata
vis='ngc5921.demo.ms'
saveinputs('importuvfits', 'ngc5921.demo.importuvfits.saved')
importuvfits()

Saveinputs saves the parameters of a given command to specified text file, handy to debug a script and see what actually was run. The parameters of importuvfits are saved to the file "ngc5921.demo.importuvfits.saved". A listing of this file follows. Notice that it is executable with execfile in CASA (remove the # commenting symbol before importuvfits to have the execfile run the command). 

CASA <71>: os.system('cat ngc5921.demo.importuvfits.saved')
taskname           = "importuvfits"
fitsfile           =  "/usr/lib64/casapy/30.0.9709test-001/data/demo/NGC5921.fits"
vis                =  "ngc5921.demo.ms"
antnamescheme      =  "old"
#importuvfits(fitsfile="/usr/lib64/casapy/30.0.9709test-001/data/demo/NGC5921.fits",vis="ngc5921.demo.ms",antnamescheme="old")


A Summary of the Data

We'll need to have a look at the observing tables to learn the calibrator and source names. The relevant command is listobs.

Logger output of listobs.
listobs(vis='ngc5921.demo.ms', verbose=True)

The output goes to the logger window; see the screenshot at right.

Tip: You can control the text size of the logger window using <ctrl>-A (smaller font) and <ctrl>-L (larger font) in Linux (<Command>-A and <Command>-L on MacOS X).

A more complete listing of the listobs output follows. 

2018-09-12 20:19:48 INFO listobs	##########################################
2018-09-12 20:19:48 INFO listobs	##### Begin Task: listobs            #####
2018-09-12 20:19:48 INFO listobs	listobs(vis="ngc5921.demo.ms",selectdata=True,spw="",field="",antenna="",
2018-09-12 20:19:48 INFO listobs	        uvrange="",timerange="",correlation="",scan="",intent="",
2018-09-12 20:19:48 INFO listobs	        feed="",array="",observation="",verbose=True,listfile="",
2018-09-12 20:19:48 INFO listobs	        listunfl=False,cachesize=50,overwrite=False)
2018-09-12 20:19:49 INFO listobs	================================================================================
2018-09-12 20:19:49 INFO listobs	           MeasurementSet Name:  /lustre/aoc/sciops/akapinsk/oldVLA/ngc5921/ngc5921.demo.ms      MS Version 2
2018-09-12 20:19:49 INFO listobs	================================================================================
2018-09-12 20:19:49 INFO listobs	   Observer: TEST     Project:   
2018-09-12 20:19:49 INFO listobs	Observation: VLA
2018-09-12 20:19:49 INFO listobs	Computing scan and subscan properties...
2018-09-12 20:19:49 INFO listobs	Data records: 22653       Total elapsed time = 5310 seconds
2018-09-12 20:19:49 INFO listobs	   Observed from   13-Apr-1995/09:18:45.0   to   13-Apr-1995/10:47:15.0 (TAI)
2018-09-12 20:19:49 INFO listobs	   
2018-09-12 20:19:49 INFO listobs	   ObservationID = 0         ArrayID = 0
2018-09-12 20:19:49 INFO listobs	  Date        Timerange (TAI)          Scan  FldId FieldName             nRows     SpwIds   Average Interval(s)    ScanIntent
2018-09-12 20:19:49 INFO listobs	  13-Apr-1995/09:18:45.0 - 09:24:45.0     1      0 1331+30500002_0           4509  [0]  [30] 
2018-09-12 20:19:49 INFO listobs	              09:27:15.0 - 09:29:45.0     2      1 1445+09900002_0           1890  [0]  [30] 
2018-09-12 20:19:49 INFO listobs	              09:32:45.0 - 09:48:15.0     3      2 N5921_2                   6048  [0]  [30] 
2018-09-12 20:19:49 INFO listobs	              09:50:15.0 - 09:51:15.0     4      1 1445+09900002_0            756  [0]  [30] 
2018-09-12 20:19:49 INFO listobs	              10:21:45.0 - 10:23:15.0     5      1 1445+09900002_0           1134  [0]  [30] 
2018-09-12 20:19:49 INFO listobs	              10:25:45.0 - 10:43:15.0     6      2 N5921_2                   6804  [0]  [30] 
2018-09-12 20:19:49 INFO listobs	              10:45:15.0 - 10:47:15.0     7      1 1445+09900002_0           1512  [0]  [30] 
2018-09-12 20:19:49 INFO listobs	           (nRows = Total number of rows per scan) 
2018-09-12 20:19:49 INFO listobs	Fields: 3
2018-09-12 20:19:49 INFO listobs	  ID   Code Name                RA               Decl           Epoch   SrcId      nRows
2018-09-12 20:19:49 INFO listobs	  0         1331+30500002_0     13:31:08.287300 +30.30.32.95900 J2000   0           4509
2018-09-12 20:19:49 INFO listobs	  1         1445+09900002_0     14:45:16.465600 +09.58.36.07300 J2000   1           5292
2018-09-12 20:19:49 INFO listobs	  2         N5921_2             15:22:00.000000 +05.04.00.00000 J2000   2          12852
2018-09-12 20:19:49 INFO listobs	Spectral Windows:  (1 unique spectral windows and 1 unique polarization setups)
2018-09-12 20:19:49 INFO listobs	  SpwID  Name   #Chans   Frame   Ch0(MHz)  ChanWid(kHz)  TotBW(kHz) CtrFreq(MHz)  Corrs  
2018-09-12 20:19:49 INFO listobs	  0      none      63   LSRK    1412.665        24.414      1550.2   1413.4219   RR  LL
2018-09-12 20:19:49 INFO listobs	Sources: 3
2018-09-12 20:19:49 INFO listobs	  ID   Name                SpwId RestFreq(MHz)  SysVel(km/s) 
2018-09-12 20:19:49 INFO listobs	  0    1331+30500002_0     0     1420.405752    0            
2018-09-12 20:19:49 INFO listobs	  1    1445+09900002_0     0     1420.405752    0            
2018-09-12 20:19:49 INFO listobs	  2    N5921_2             0     1420.405752    0            
2018-09-12 20:19:49 INFO listobs	Antennas: 27:
2018-09-12 20:19:49 INFO listobs	  ID   Name  Station   Diam.    Long.         Lat.                Offset from array center (m)                ITRF Geocentric coordinates (m)        
2018-09-12 20:19:49 INFO listobs	                                                                     East         North     Elevation               x               y               z
2018-09-12 20:19:49 INFO listobs	  0    1     VLA:N7    25.0 m   -107.37.07.2  +33.54.12.9        -30.2623      345.7477       -0.8872 -1601155.613187 -5041783.882304  3555162.343090
2018-09-12 20:19:49 INFO listobs	  1    2     VLA:W1    25.0 m   -107.37.05.9  +33.54.00.5          3.5004      -39.7725        0.9883 -1601188.991307 -5042000.530918  3554843.409670
2018-09-12 20:19:49 INFO listobs	  2    3     VLA:W2    25.0 m   -107.37.07.4  +33.54.00.9        -37.1358      -25.0262        1.0383 -1601225.244615 -5041980.431775  3554855.677111
2018-09-12 20:19:49 INFO listobs	  3    4     VLA:E1    25.0 m   -107.37.05.7  +33.53.59.2          6.9833      -79.6414        1.1565 -1601192.444530 -5042022.911771  3554810.411780
2018-09-12 20:19:49 INFO listobs	  4    5     VLA:E3    25.0 m   -107.37.02.8  +33.54.00.5         81.5188      -37.9632        1.0246 -1601114.335629 -5042023.211477  3554844.931655
2018-09-12 20:19:49 INFO listobs	  5    6     VLA:E9    25.0 m   -107.36.45.1  +33.53.53.6        536.8977     -250.3175        0.1183 -1600715.915813 -5042273.186780  3554668.167811
2018-09-12 20:19:49 INFO listobs	  6    7     VLA:E6    25.0 m   -107.36.55.6  +33.53.57.7        267.7566     -124.8145        1.2815 -1600951.554888 -5042125.947753  3554772.987072
2018-09-12 20:19:49 INFO listobs	  7    8     VLA:W8    25.0 m   -107.37.21.6  +33.53.53.0       -401.2640     -270.6305        2.2293 -1601614.059494 -5042001.699973  3554652.484758
2018-09-12 20:19:49 INFO listobs	  8    9     VLA:N5    25.0 m   -107.37.06.7  +33.54.08.0        -16.9948      194.1215       -0.1368 -1601168.756077 -5041869.099542  3555036.914367
2018-09-12 20:19:49 INFO listobs	  9    10    VLA:W3    25.0 m   -107.37.08.9  +33.54.00.1        -74.4964      -50.1921        1.1608 -1601265.132224 -5041982.597979  3554834.857504
2018-09-12 20:19:49 INFO listobs	  10   11    VLA:N4    25.0 m   -107.37.06.5  +33.54.06.1        -11.7487      134.3686        0.1774 -1601173.922897 -5041902.701204  3554987.495105
2018-09-12 20:19:49 INFO listobs	  11   12    VLA:W5    25.0 m   -107.37.13.0  +33.53.57.8       -179.2554     -120.8635        1.4872 -1601376.990711 -5041988.712764  3554776.381187
2018-09-12 20:19:49 INFO listobs	  12   13    VLA:N3    25.0 m   -107.37.06.3  +33.54.04.8         -8.2438       94.5297        0.3947 -1601177.362708 -5041925.112425  3554954.550128
2018-09-12 20:19:49 INFO listobs	  13   14    VLA:N1    25.0 m   -107.37.06.0  +33.54.01.8         -0.0030        0.0445        0.8773 -1601185.580779 -5041978.216463  3554876.396287
2018-09-12 20:19:49 INFO listobs	  14   15    VLA:N2    25.0 m   -107.37.06.2  +33.54.03.5         -4.7904       54.7090        0.5774 -1601180.839839 -5041947.470902  3554921.600805
2018-09-12 20:19:49 INFO listobs	  15   16    VLA:E7    25.0 m   -107.36.52.4  +33.53.56.5        348.8969     -162.6653        1.0336 -1600880.544215 -5042170.427468  3554741.431900
2018-09-12 20:19:49 INFO listobs	  16   17    VLA:E8    25.0 m   -107.36.48.9  +33.53.55.1        438.6654     -204.5038        0.5027 -1600801.910482 -5042219.412805  3554706.408864
2018-09-12 20:19:49 INFO listobs	  17   18    VLA:W4    25.0 m   -107.37.10.8  +33.53.59.1       -122.0163      -82.2819        1.2624 -1601315.866196 -5041985.352573  3554808.279150
2018-09-12 20:19:49 INFO listobs	  18   19    VLA:E5    25.0 m   -107.36.58.4  +33.53.58.8        195.8349      -91.2758        1.2155 -1601014.427180 -5042086.300814  3554800.787928
2018-09-12 20:19:49 INFO listobs	  19   20    VLA:W9    25.0 m   -107.37.25.1  +33.53.51.0       -491.1000     -331.2429        2.5539 -1601709.998072 -5042006.975455  3554602.355417
2018-09-12 20:19:49 INFO listobs	  20   21    VLA:W6    25.0 m   -107.37.15.6  +33.53.56.4       -244.9704     -165.2178        1.6861 -1601447.161927 -5041992.554228  3554739.677219
2018-09-12 20:19:49 INFO listobs	  21   22    VLA:E4    25.0 m   -107.37.00.8  +33.53.59.7        133.6478      -62.2829        1.0919 -1601068.773396 -5042051.970054  3554824.783566
2018-09-12 20:19:49 INFO listobs	  23   24    VLA:E2    25.0 m   -107.37.04.4  +33.54.01.1         40.6649      -18.9151        0.9550 -1601150.040469 -5042000.665669  3554860.702914
2018-09-12 20:19:49 INFO listobs	  24   25    VLA:N6    25.0 m   -107.37.06.9  +33.54.10.3        -23.2197      265.3902       -0.4819 -1601162.569974 -5041829.054708  3555095.873969
2018-09-12 20:19:49 INFO listobs	  25   26    VLA:N9    25.0 m   -107.37.07.8  +33.54.19.0        -46.5533      532.1581       -1.8550 -1601139.422904 -5041679.082136  3555316.518142
2018-09-12 20:19:49 INFO listobs	  26   27    VLA:N8    25.0 m   -107.37.07.5  +33.54.15.8        -38.0437      434.7201       -1.3387 -1601147.894127 -5041733.868915  3555235.935926
2018-09-12 20:19:49 INFO listobs	  27   28    VLA:W7    25.0 m   -107.37.18.4  +33.53.54.8       -319.1171     -215.2368        1.9407 -1601526.340031 -5041996.897001  3554698.302182
2018-09-12 20:19:49 INFO listobs	##### End Task: listobs              #####
2018-09-12 20:19:49 INFO listobs	##########################################


Key Information from listobs

Certainly the output of listobs is dense with information, but there are some particularly vital data that we'll need for the calibration.

  • The calibrators are 1331+305* (3C286, the flux and bandpass calibrator) and 1445+099* (the phase calibrator). We can use wild-cards 1331* and 1445* since they uniquely identify the sources.
  • The calibrator field indices are field='0' (1331+305) and field='1' (1445+099).
  • The name of the source in the observations list is N5921_2, or field = '2'.
  • The data were taken in a single IF (a single spectral window, SpwID = 0), divided into 63 channels.
  • Only RR and LL correlations are present; cross-pols are absent.


Flagging

Flag the autocorrelations

We don't need the autocorrelation data, and we can use flagdata to get rid of them. You shouldn't have to specify the measurement set, because the variable vis is already set, but it never hurts to be cautious. 

flagdata(vis="ngc5921.demo.ms", autocorr=True)

CASA issues a warning when running this command, but this can be disregarded.

Interactive Flagging

Plotms settings for flagging spectral line data. Click to enlarge.

plotms is a good tool for flagging spectral line data. Check out the tutorial that describes editing VLA continuum data. Spectral line data of course require some consideration of channels and channel averaging. 

plotms()


The figure at right highlights the settings needed for effective editing of a spectral line data set. The key settings are as follows.

  • Specify the measurement set in File; the Browse button allows you to hunt down the measurement set.
  • It's better to edit one source at a time. In the illustrated example, the flux / bandpass calibrator 1331+305* is displayed.
  • Average the channels. First, specify the central channels to remove band edge effects. Channels 6~56 in the first spectral window (IF) are appropriate (see #Inspect the Bandpass Response Curve, below). In plotms() specify: spw=0:6~56. In the Channel Averaging box, enter 51 channels to average over all channels in the given range.
  • Ideally you want the channels to have the same (u, v) coverage (projected baseline spacings as viewed from the source); otherwise, the beam (point spread function) will be different for each channel. Therefore, if you flag data from a given channel it's usually a good idea to flag those data from all channels. Under the Flagging tab, specify Extend flags to Channel.


With these settings, interactive flagging proceeds as for continuum data. When you're satisfied with the edits, File → Quit to return to the CASA prompt.

Calibration

Calibration of spectral line data broadly follows the approach for continuum data, except that the amplitude and phase corrections are a function of frequency and so must be corrected by bandpass calibration. The basic calibration steps follow.

  • Set the flux scale of the primary calibrator, here, 1331+305 = 3C 286.
  • Determine bandpass corrections based on the primary calibrator. In the script that follows, the bandpass calibration will be stored in ngc5921.demo.bcal.
  • Inspect the bandpass correction to determine viable channels for averaging and imaging. We want to toss out end channels where the response is poor.
  • Determine the gain calibrations on the bandpass-corrected and channel-averaged data. In this step, we effectively turn the spectral line data into a single-channel continuum data set and calibrate accordingly. The calibration is first stored in ngc5921.demo.gcal. In the second part of this step we correct the fluxscale of the .gcal table, and store the final calibration solutions with correct fluxes in ngc5921.demo.fluxscale (this is the table that needs to be applied later to the data, not the .gcal version).
  • Inspect the gain calibration solutions to look for any aberrant solutions that hint at bad calibrator data.
  • Apply the calibration solutions to the source (N5921_2). This action literally adds a new column of data to the measurement set. This new column contains the data with the gain calibration and bandpass calibration applied, but it does not overwrite the raw data in case the calibration needs revision.


Setting the Flux Scale

setjy generates a source model for the primary calibrator, 1331+305 = 3C286. From CASA 5.3+ the default standard is Perley-Butler 2017, and includes resolved structure of the calibrators. This is 1.4GHz D-config and 1331+305 is sufficiently unresolved that, in principle, we don't need a model image; however, here we proceed with applying the detailed model, as a good practice.

setjy also looks up the radio SED for common flux calibrators and automatically assigns the total flux density. 

# 1331+305 = 3C286 is our primary calibrator. Use the wildcard on the end of the source name
setjy(vis='ngc5921.demo.ms', field='1331+305*', model='3C286_L.im')


A summary of the operation is sent to the logger window. Here's a listing of the output. 

2018-09-12 20:43:38 INFO setjy	##########################################
2018-09-12 20:43:38 INFO setjy	##### Begin Task: setjy              #####
2018-09-12 20:43:38 INFO setjy	setjy(vis="ngc5921.demo.ms",field="1331+305*",spw="",selectdata=False,timerange="",
2018-09-12 20:43:38 INFO setjy	        scan="",intent="",observation="",scalebychan=True,standard="Perley-Butler 2017",
2018-09-12 20:43:38 INFO setjy	        model="3C286_L.im",modimage="",listmodels=False,fluxdensity=-1,spix=0.0,
2018-09-12 20:43:38 INFO setjy	        reffreq="1GHz",polindex=[],polangle=[],rotmeas=0.0,fluxdict={},
2018-09-12 20:43:38 INFO setjy	        useephemdir=False,interpolation="nearest",usescratch=False,ismms=False)
2018-09-12 20:43:38 INFO setjy	{'field': '1331+305*'}
2018-09-12 20:43:38 INFO Imager	Opening MeasurementSet [...]
2018-09-12 20:43:38 INFO setjy	Using /home/casa/data/distro/nrao/VLA/CalModels/3C286_L.im for modimage.
2018-09-12 20:43:38 INFO setjy	CASA Version 5.3.0-143  
2018-09-12 20:43:38 INFO setjy	   
2018-09-12 20:43:39 INFO imager	Using channel dependent flux densities
2018-09-12 20:43:39 INFO imager	Selected 4509 out of 22653 rows.
2018-09-12 20:43:39 INFO imager	1331+30500002_0 (fld ind 0) spw 0  [I=15.016, Q=0, U=0, V=0] Jy @ 1.4127e+09Hz, (Perley-Butler 2017)
2018-09-12 20:43:40 INFO imager	Using model image /home/casa/data/distro/nrao/VLA/CalModels/3C286_L.im
2018-09-12 20:43:40 INFO imager	Scaling spw(s) [0]'s model image by channel to  I = 15.0159, 15.0118, 15.0077 Jy @(1.41265e+09, 1.41343e+09, 1.41419e+09)Hz (LSRK) for visibility prediction (a few representative values are shown).
2018-09-12 20:43:40 INFO imager	The model image's reference pixel is 0.00904522 arcsec from 1331+30500002_0's phase center.
2018-09-12 20:43:40 INFO imager	Will clear any existing model with matching field=1331+30500002_0 and spw=*
2018-09-12 20:43:40 INFO  	Clearing model records in MS header for selected fields.
2018-09-12 20:43:40 INFO  	 1331+30500002_0 (id = 0) deleted.
2018-09-12 20:43:40 INFO imager	Selected 4509 out of 22653 rows.
2018-09-12 20:43:40 INFO setjy	##### End Task: setjy                #####
2018-09-12 20:43:40 INFO setjy	##########################################


Bandpass Calibration

The flux calibrator 1331+305 = 3C 286 now has a model assigned to it. Since the bandwidth of our observations is only 1.55 MHz, the model doesn't change over this narrow range of frequencies, so we can use it to determine amplitude and phase (gain) corrections for each channel independently. The result is the bandpass calibration.

As for any antenna-based calibration scheme, we have to pick an antenna to act as the reference point for the calibration. Any antenna will do, but it's better to pick one near the center of the array. For the remainder of the calibration, we will use refant = '15'. 

# We can first do the bandpass on the single 5min scan on 1331+305. At 1.4GHz phase stablility should be sufficient to do this without
# a first (rough) gain calibration. This will give us the relative antenna gain as a function of frequency.
bandpass(vis='ngc5921.demo.ms', caltable='ngc5921.demo.bcal', field='0', selectdata=False, bandtype='B', solint='inf', combine='scan', refant='15')
  • field='0' : Use the flux calibrator 1331+305 = 3C286 (FIELD_ID 0) as bandpass calibrator.
  • bandtype='B' : Choose bandpass solution type. Pick standard time-binned B (rather than BPOLY).
  • solint='inf' and combine='scan' : Set solution interval arbitrarily long (get single bandpass).
  • refant = '15' : Reference antenna Name 15 (15=VLA:N2) (Id 14)


Inspect the Bandpass Response Curve

Bandpass response curves generated by plotms. The solutions for different antennas are indicated by differently colored plotting symbols. Plots for individual antennas can be generated by setting iteraxis = 'antenna' for plotms.

In the gain calibration to follow, we will effectively convert the spectral line data into a continuum data set. Before proceeding, we need to inspect the bandpass calibration to make sure that it contains no bad values and also to inspect which channels to average to produce the continuum data. Plotms is the standard tool for plotting calibration solutions. The following commands produce the figure at right. 

# Set up 2x1 panels - upper panel amp vs. channel
plotms(vis='ngc5921.demo.bcal', field='0', gridrows=2, gridcols=1, plotindex=0, rowindex=0, xaxis='channel', yaxis='amp', showgui=True, clearplots=False, coloraxis='antenna1')

# Set up 2x1 panels - lower panel phase vs. channel
plotms(vis='ngc5921.demo.bcal', field='0', gridrows=2, gridcols=1, plotindex=1, rowindex=1, xaxis='channel', yaxis='phase', showgui=True, clearplots=False, coloraxis='antenna1')

By inspection, the amplitude response curve is flat over channels 6~56; that channel range will be used to generate the continuum data for gain calibration. If you want to further inspect the plots interactively and iterate over antenna, set iteraxis = 'antenna'

Notice that plotms is run twice: once to display gain amplitudes as a function of channel (frequency), and again to plot gain phases as a function of channel.


Gain Calibration

From inspection of the bandpass response curve, we can average channels 6~56 to produce continuum data for the calibrators. For VLA data, this averaging is specified through the spw (spectral window) parameter, which takes the form IF:Channel-range, as follows. 

spw = '0:6~56'

That is, there is only one spectral window (IF), spw = 0, and we want to average channels 6~56 within that spectral window.

Gain calibrations are otherwise determined as for continuum data.

  • gaincal() is run only on the calibrators, 1331+305 (flux calibrator) and 1445+099 (phase calibrator).
  • The default model for gain calibrations is a 1 Jy point-source. The flux scale is overridden by setjy, which has been performed for the flux calibrator. We need to transfer that flux scale to the phase calibrator using fluxscale().
  • Note that fluxscale() determines the flux density of the phase calibrator and accordingly adjusts its model and calibration solutions. A report of the results are sent to the logger window.
  • Unless you use parameter incremental=True while executing fluxscale() (the default is False), the resulting .fluxscale table will replace the .gcal table at this point. This particularly important in the applycal() stage.
# Armed with the bandpass, we now solve for the time-dependent antenna gains using our previously determined bandpass.
# Note this will automatically be applied to all sources not just the one used to determine the bandpass

gaincal(vis='ngc5921.demo.ms', caltable='ngc5921.demo.gcal', gaintable=['ngc5921.demo.bcal'], interp=['nearest'], field='0,1', spw='0:6~56', gaintype='G', solint='inf', calmode='ap', refant='15')



# Now we will transfer the flux scale to the phase calibrator. 
# We will be using 1331+305 (the source we did setjy on) as our flux standard reference.
# Note its extended name as in the FIELD table summary above (it has a VLA seq number appended)

fluxscale(vis='ngc5921.demo.ms', fluxtable='ngc5921.demo.fluxscale', caltable='ngc5921.demo.gcal', reference='1331*', transfer='1445*')


The output from fluxscale follows. A relatively large uncertainty for the phase calibrator is a sign that something went wrong, perhaps bad solutions in gaincal. Here, the phase calibrator scaled to 2.486 ± 0.001 Jy, which looks reasonable. 

2018-09-12 22:38:05 INFO fluxscale	##########################################
2018-09-12 22:38:05 INFO fluxscale	##### Begin Task: fluxscale          #####
2018-09-12 22:38:05 INFO fluxscale	fluxscale(vis="ngc5921.demo.ms",caltable="ngc5921.demo.gcal",fluxtable="ngc5921.demo.fluxscale",reference="1331*",transfer="1445*",
2018-09-12 22:38:05 INFO fluxscale	        listfile="",append=False,refspwmap=[-1],gainthreshold=-1.0,antenna="",
2018-09-12 22:38:05 INFO fluxscale	        timerange="",scan="",incremental=False,fitorder=1,display=False)
2018-09-12 22:38:05 INFO fluxscale	****Using NEW VI2-driven calibrater tool****
2018-09-12 22:38:05 INFO fluxscale	Opening MS: ngc5921.demo.ms for calibration.
2018-09-12 22:38:05 INFO fluxscale	Initializing nominal selection to the whole MS.
2018-09-12 22:38:05 INFO fluxscale	Beginning fluxscale--(MSSelection version)-------
2018-09-12 22:38:05 INFO fluxscale	 Found reference field(s): 1331+30500002_0
2018-09-12 22:38:05 INFO fluxscale	 Found transfer field(s):  1445+09900002_0
2018-09-12 22:38:05 INFO fluxscale	 Flux density for 1445+09900002_0 in SpW=0 (freq=1.41342e+09 Hz) is: 2.52609 +/- 0.00218206 (SNR = 1157.67, N = 54)
2018-09-12 22:38:05 INFO fluxscale	Storing result in ngc5921.demo.fluxscale
2018-09-12 22:38:05 INFO fluxscale	Writing solutions to table: ngc5921.demo.fluxscale
2018-09-12 22:38:06 INFO fluxscale	CASA Version 5.3.0-143  
2018-09-12 22:38:06 INFO fluxscale	   
2018-09-12 22:38:07 INFO fluxscale	##### End Task: fluxscale            #####
2018-09-12 22:38:07 INFO fluxscale	##########################################


Inspect the Calibration Solutions

Gain calibration solutions from gaincal and fluxscale.

Now inspect the results of gaincal. The setup is identical to that used to plot the bandpass response curve. The only change is that we are plotting the gaintable ngc5921.demo.gcal, and we're looking at solutions for both of the calibrator sources. The results are shown at right. 

# Set up 2x1 panels - upper panel amp vs. time
plotms(vis='ngc5921.demo.fluxscale', field='0,1', gridrows=2, gridcols=1, plotindex=0, rowindex=0, yaxis='amp', showgui=True, clearplots=False, coloraxis='antenna1')

# Set up 2x1 panels - lower panel phase vs. time
plotms(vis='ngc5921.demo.fluxscale', field='0,1', gridrows=2, gridcols=1, plotindex=1, rowindex=1, yaxis='phase', showgui=True, clearplots=False, coloraxis='antenna1')

The amp and phase coherence looks good. If you want to do this interactively and iterate over antenna, set iteraxis = 'antenna'.


Apply the Solutions

Next, apply the calibration solutions to the calibrators themselves, and finally transfer the calibration solutions by interpolation (or nearest-neighbor sampling) to the source. The relevant task is applycal, which fills out a new column (CORRECTED_DATA) of calibrated data in the measurement set without wiping out the raw data column. The application is identical to that used for continuum data, except that the bandpass table is also included in the calibration. To apply multiple calibrations at once, provide the gaintable parameter with a list of calibration tables, as follows. 

gaintable = ['ngc5921.demo.fluxscale', 'ngc5921.demo.bcal']

We want to correct the calibrators using themselves and transfer from 1445+099 to itself and the target N5921. Start with the fluxscale/gain and bandpass tables. We will pick the 1445+099 out of the gain table for transfer and use all of the bandpass table. Also, note that the table .fluxscale has the .gcal solutions with the correct flux scale applied, and so there is no need to invoke the .gcal again in the applycal() command below. 

applycal(vis='ngc5921.demo.ms', field='1,2', gaintable=['ngc5921.demo.fluxscale','ngc5921.demo.bcal'], gainfield=['1','*'], 
         interp=['linear','nearest'], spwmap=[], selectdata=False)

Now for completeness apply 1331+305 to itself. 

applycal(vis='ngc5921.demo.ms', field='0', gaintable=['ngc5921.demo.fluxscale','ngc5921.demo.bcal'], gainfield=['0','*'], 
         interp=['linear','nearest'], spwmap=[], selectdata=False)


Plot the Spectrum

plotms settings to produce the integrated spectrum from the calibrated visibilities data.

Before we attempt to image the 21 cm cube of the source, we need to subtract off the underlying continuum, which means we need to plot the integrated spectrum of the source to determine the continuum channels.

We can do this in plotms. 

plotms(vis='ngc5921.demo.ms', selectdata=True, field='N5921*', spw='0:6~56', averagedata=True, avgtime='3600', avgscan=True, 
       avgbaseline=True, xaxis='channel', yaxis='amp', ydatacolumn='corrected')

Note that we have entered all the relevent parameters via the task interface, as an alternative to entering each option into the GUI. If the symbols appear too small, the size may be increased via the Display tab: change the Unflagged Points Symbol to 'Custom' and increase the number of pixels for the plotting symbol. The resulting plot is illustrated in the figure at right. Briefly, we want to average both in time and over baselines to get the signal-to-noise necessary to reveal the 21 cm profile (see Averaging data in plotms for more details on averaging options). If you wish to enter the values directly into the GUI, you can follow the (Tab)Command convention of the flagging tutorial with the following settings :

  • (Data)field = N5921*
  • (Data)spw = 0:6~56
  • (Data)Averaging → Time = 3600 (average over some long time window)
  • (Data)Averaging → Scan = True (checkmark; average in time across scan boundaries)
  • (Data)Averaging → All Baselines = True (checkmark)
  • (Axes)X Axis = Channel
  • (Axes)Y Axis = Amp

Now remove 0:6~56 from the spw field to see all channels. From inspection of this plot, it looks like channels 4~6 and 50~59 contain line-free channels, suitable to use for continuum subtraction.

Continuum Subtraction

The next step is to split off the NGC 5921 data from the multisource measurement set and subtract the continuum. Splitting uses the split command, as follows. 

split(vis='ngc5921.demo.ms', outputvis='ngc5921.demo.src.split.ms', field='N5921*', spw='', datacolumn='corrected')


This action generated a new measurement set called ngc5921.demo.src.split.ms and copied the calibrated source data (datacolumn = 'corrected') into it.

uvcontsub subtracts the continuum from the data in the visibility (u, v) plane. We will be using channels 4-6 and 50-59 for continuum. 

uvcontsub(vis='ngc5921.demo.src.split.ms', field='N5921*', fitspw='0:4~6;50~59', spw='0', solint=0.0, fitorder=0, want_cont=True)


Notice that uvcontsub splits the data into two new measurement sets, 'ngc5921.demo.ms.cont', which contains an average of the continuum channels, and 'ngc5921.demo.ms.contsub', which contains the continuum-subtracted spectral line data.



Imaging

Plot of amplitude vs. projected baseline length (in units of the observing wavelength) produced by casaplotms. The maximum baseline is just below 5 kilo-lambda.

Now we can generate the primary science product, a tclean data cube (ra, dec, velocity) from the continuum-subtracted (u, v, channel) measurement set, ngc5921.demo.ms.contsub. Things to consider in using tclean:

  • To ensure channels aren't averaged prior to imaging, choose mode='channel'.
  • Specify the channels to image using start = 5, width = 1 (no averaging over channels), nchan = 46; only channels 5~51 will be imaged.
  • The maximum baseline is just under 5 kilolambda (see the figure at right), and so the expected synthetic beam is roughly 1.22 × 206265 / 5000 = 50 arcseconds (subject to the details of u, v weighting). Pixels should sample the beam better than 3 times, so 15 arcseconds is a good choice of pixel size (cell = ['15.0arcsec','15.0arcsec']).
  • We only want to tclean down to the noise, which is easily determined by trial-and-error imaging of a single channel (choosing nchan=1 and start appropriately). Here, tclean stops when the maximum residual on the channel is below threshold='8.0mJy'.


# Image the continuum subtracted measurement set
tclean(vis='ngc5921.demo.src.split.ms.contsub', imagename='ngc5921.demo.tclean', field='0', datacolumn='data', specmode='cube', nchan=46, start=5, width=1, spw='', deconvolver='hogbom', gridder='standard', niter=6000, gain=0.1, threshold='8.0mJy', imsize=[256,256], cell=['15.0arcsec','15.0arcsec'], weighting='briggs', robust=0.5,  mask = 'box[[108pix,108pix],[148pix,148pix]]', interactive=False, pblimit=-0.2)


Use imhead to look at the cube header: 

imhead(imagename='ngc5921.demo.tclean.image', mode='summary')


The output, as follows, appears in the logger window. 


2019-09-03 22:25:25 INFO imhead	##### Begin Task: imhead             #####
2019-09-03 22:25:25 INFO imhead	imhead(imagename="ngc5921.demo.tclean.image",mode="summary",hdkey="",hdvalue="",verbose=False)
2019-09-03 22:25:25 INFO ImageMetaData	   
2019-09-03 22:25:25 INFO ImageMetaData	Image name       : ngc5921.demo.tclean.image
2019-09-03 22:25:25 INFO ImageMetaData	Object name      : N5921_2
2019-09-03 22:25:25 INFO ImageMetaData	Image type       : PagedImage
2019-09-03 22:25:25 INFO ImageMetaData	Image quantity   : Intensity
2019-09-03 22:25:25 INFO ImageMetaData	Pixel mask(s)    : None
2019-09-03 22:25:25 INFO ImageMetaData	Region(s)        : None
2019-09-03 22:25:25 INFO ImageMetaData	Image units      : Jy/beam
2019-09-03 22:25:25 INFO ImageMetaData	Restoring Beams 
2019-09-03 22:25:25 INFO ImageMetaData	Pol   Type Chan        Freq     Vel
2019-09-03 22:25:25 INFO ImageMetaData	  I    Max    7 1.41296e+09 1571.92   51.6639 arcsec x   47.3594 arcsec pa=  8.3205 deg
2019-09-03 22:25:25 INFO ImageMetaData	  I    Min   43 1.41384e+09 1386.42   51.5897 arcsec x   47.3127 arcsec pa=  8.3869 deg
2019-09-03 22:25:25 INFO ImageMetaData	  I Median   13 1.41310e+09 1541.00   51.6611 arcsec x   47.3318 arcsec pa=  8.1425 deg
2019-09-03 22:25:25 INFO ImageMetaData	   
2019-09-03 22:25:25 INFO ImageMetaData	Direction reference : J2000
2019-09-03 22:25:25 INFO ImageMetaData	Spectral  reference : LSRK
2019-09-03 22:25:25 INFO ImageMetaData	Velocity  type      : RADIO
2019-09-03 22:25:25 INFO ImageMetaData	Rest frequency      : 1.42041e+09 Hz
2019-09-03 22:25:25 INFO ImageMetaData	Pointing center     :  15:22:00.000000  +05.04.00.000000
2019-09-03 22:25:25 INFO ImageMetaData	Telescope           : VLA
2019-09-03 22:25:25 INFO ImageMetaData	Observer            : TEST
2019-09-03 22:25:25 INFO ImageMetaData	Date observation    : 1995/04/13/09:33:00
2019-09-03 22:25:25 INFO ImageMetaData	Telescope position: [-1.60119e+06m, -5.04198e+06m, 3.55488e+06m] (ITRF)
2019-09-03 22:25:25 INFO ImageMetaData	   
2019-09-03 22:25:25 INFO ImageMetaData	Axis Coord Type      Name             Proj Shape Tile   Coord value at pixel    Coord incr Units
2019-09-03 22:25:25 INFO ImageMetaData	------------------------------------------------------------------------------------------------ 
2019-09-03 22:25:25 INFO ImageMetaData	0    0     Direction Right Ascension   SIN   256   64  15:22:00.000   128.00 -1.500000e+01 arcsec
2019-09-03 22:25:25 INFO ImageMetaData	1    0     Direction Declination       SIN   256   64 +05.04.00.000   128.00  1.500000e+01 arcsec
2019-09-03 22:25:25 INFO ImageMetaData	2    1     Stokes    Stokes                    1    1             I
2019-09-03 22:25:25 INFO ImageMetaData	3    2     Spectral  Frequency                46    8   1.41279e+09     0.00 2.4414062e+04 Hz
2019-09-03 22:25:25 INFO ImageMetaData	                     Velocity                               1607.99     0.00 -5.152860e+00 km/s
2019-09-03 22:25:25 INFO imhead	##### End Task: imhead               #####
2019-09-03 22:25:25 INFO imhead	##########################################



Additional Science Products

If things went well, you should now have a spectral line cube (ngc5921.demo.tclean.image) as a primary science product. The demo script illustrates further how to generate cube statistics (using imstat), an integrated spectrum, and moment maps.

Cube Statistics

imstat is the tool for displaying statistics of images and cubes. The following example displays the statistics for an empty region of the whole cube. 


cubestat=imstat(imagename='ngc5921.demo.tclean.image', box='10,10,100,100')



2019-09-03 22:40:15 INFO imstat	##########################################
2019-09-03 22:40:15 INFO imstat	##### Begin Task: imstat             #####
2019-09-03 22:40:15 INFO imstat	imstat(imagename="ngc5921.demo.tclean.image",axes=-1,region="",box="10,10,100,100",chans="",
2019-09-03 22:40:15 INFO imstat	        stokes="",listit=True,verbose=True,mask="",stretch=False,
2019-09-03 22:40:15 INFO imstat	        logfile="",append=True,algorithm="classic",fence=-1,center="mean",
2019-09-03 22:40:15 INFO imstat	        lside=True,zscore=-1,maxiter=-1,clmethod="auto",niter=3)
2019-09-03 22:40:15 INFO imstat	Using specified box(es) 10,10,100,100
2019-09-03 22:40:15 INFO imstat	Determining stats for image ngc5921.demo.tclean.image
2019-09-03 22:40:15 INFO imstat	Selected bounding box : 
2019-09-03 22:40:15 INFO imstat	    [10, 10, 0, 0] to [100, 100, 0, 45]  (15:23:58.379, +04.34.29.305, I, 1.41279e+09Hz to 15:22:28.105, +04.56.59.962, I, 1.41389e+09Hz)
2019-09-03 22:40:15 INFO imstat	Statistics calculated using Classic algorithm
2019-09-03 22:40:15 INFO imstat	Regions --- 
2019-09-03 22:40:15 INFO imstat	         -- bottom-left corner (pixel) [blc]:  [10, 10, 0, 0]
2019-09-03 22:40:15 INFO imstat	         -- top-right corner (pixel) [trc]:    [100, 100, 0, 45]
2019-09-03 22:40:15 INFO imstat	         -- bottom-left corner (world) [blcf]: 15:23:58.379, +04.34.29.305, I, 1.41279e+09Hz
2019-09-03 22:40:15 INFO imstat	         -- top-right corner (world) [trcf]:   15:22:28.105, +04.56.59.962, I, 1.41389e+09Hz
2019-09-03 22:40:15 INFO imstat	Values --- 
2019-09-03 22:40:15 INFO imstat	         -- flux [flux]:                            1.99425 Jy.km/s
2019-09-03 22:40:15 INFO imstat	         -- number of points [npts]:                380926
2019-09-03 22:40:15 INFO imstat	         -- maximum value [max]:                    0.00953276 Jy/beam
2019-09-03 22:40:15 INFO imstat	         -- minimum value [min]:                    -0.0100478 Jy/beam
2019-09-03 22:40:15 INFO imstat	         -- position of max value (pixel) [maxpos]: [85, 63, 0, 8]
2019-09-03 22:40:15 INFO imstat	         -- position of min value (pixel) [minpos]: [30, 18, 0, 7]
2019-09-03 22:40:15 INFO imstat	         -- position of max value (world) [maxposf]: 15:22:43.151, +04.47.44.907, I, 1.41298e+09Hz
2019-09-03 22:40:15 INFO imstat	         -- position of min value (world) [minposf]: 15:23:38.319, +04.36.29.518, I, 1.41296e+09Hz
2019-09-03 22:40:15 INFO imstat	         -- Sum of pixel values [sum]:               4.76561 Jy/beam
2019-09-03 22:40:15 INFO imstat	         -- Sum of squared pixel values [sumsq]:     1.38125 Jy/beam.Jy/beam
2019-09-03 22:40:15 INFO imstat	Statistics --- 
2019-09-03 22:40:15 INFO imstat	        -- Mean of the pixel values [mean]:         1.25106e-05 Jy/beam
2019-09-03 22:40:15 INFO imstat	        -- Variance of the pixel values :           3.62589e-06 Jy/beam
2019-09-03 22:40:15 INFO imstat	        -- Standard deviation of the Mean [sigma]:  0.00190418 Jy/beam
2019-09-03 22:40:15 INFO imstat	        -- Root mean square [rms]:                  0.00190421 Jy/beam
2019-09-03 22:40:15 INFO imstat	        -- Median of the pixel values [median]:     6.37541e-06 Jy/beam
2019-09-03 22:40:15 INFO imstat	        -- Median of the deviations [medabsdevmed]: 0.00127676 Jy/beam
2019-09-03 22:40:15 INFO imstat	        -- IQR [quartile]:                          0.00255348 Jy/beam
2019-09-03 22:40:15 INFO imstat	        -- First quartile [q1]:                     -0.00126604 Jy/beam
2019-09-03 22:40:15 INFO imstat	        -- Third quartile [q3]:                     0.00128744 Jy/beam
2019-09-03 22:40:15 INFO imstat	Sum column unit = Jy/beam
2019-09-03 22:40:15 INFO imstat	Mean column unit = Jy/beam
2019-09-03 22:40:15 INFO imstat	Std_dev column unit = Jy/beam
2019-09-03 22:40:15 INFO imstat	Minimum column unit = Jy/beam
2019-09-03 22:40:15 INFO imstat	Maximum column unit = Jy/beam
2019-09-03 22:40:15 INFO imstat	Npts          Sum           Mean          Rms           Std_dev       Minimum       Maximum     
2019-09-03 22:40:15 INFO imstat	 3.809260e+05  4.765615e+00  1.251060e-05  1.904215e-03  1.904176e-03 -1.004781e-02  9.532760e-03
2019-09-03 22:40:15 INFO imstat	##### End Task: imstat               #####
2019-09-03 22:40:15 INFO imstat	##########################################


The Integrated Spectrum

Example of the viewer rectangle selection tool on one channel of the NGC 5921 21 cm data cube. The spectral profile window is shown at right.


We saw earlier how to generate an integrated spectrum from the (u, v) measurement set. Here's how to produce the integrated spectrum from the spectral line cube. First, load the cube into viewer. 

viewer(infile='ngc5921.demo.tclean.image')


To generate the integrated spectrum, perform the following tasks.

  • Use the player controls to inspect the cube one channel at a time.
  • From the viewer Tools menu, select Spectral Profile. A new graphics window should appear.
  • By default, the rectangle selection tool is assigned to the right mouse button, and you can just right-click and drag a box over the region where you want to (spatially) integrate the spectrum. See the figure at upper right.
  • Alternatively, you can assign one of the other selection tools by right-clicking on the appropriate button.
  • The spectrum now appears in the graphics window; see the figure at right.

You can save the integrated spectrum to a text file by clicking the button on the graphics window. There are also buttons to print the figure or save the figure to disk.


Cube Moments

The moment 0 (integrated intensity) 21 cm image of NGC 5921, produced using immoments

Cube moments are maps of weighted sums along the velocity axis. In CASA, they are generated by the task immoments. The zeroth moment (moments = 0) is a sum of intensities along the velocity axis (the integrated intensity map); the first moment (moment = 1) is the sum of velocities weighted by intensity (the velocity field); the second moment (moment = 2) is a map of the velocity dispersion; see the immoments for additional options.

The following example produces maps of the zeroth and first moments, or the integrated intensity and velocity field. The respective measurement sets are the moment zero image ngc5921.demo.moments.integrated and moment one imagengc5921.demo.moments.weighted_coord.

We will do the zeroth and first moments and mask out noisy pixels using hard global limits. We will also collapse along the spectral (channel) axis and include all planes. 

immoments(imagename='ngc5921.demo.tclean.image', moments=[0,1], excludepix=[-100, 0.009], axis='spectral', chans='', outfile='ngc5921.demo.moments')
  • moments = [0,1] : Do zeroth and first moments
  • excludepix = [-100,0.009] : Mask out noisy pixels using hard global limits
  • axis = 'spectral' : Collapse along the spectral (channel) axis
  • chans = :Include all planes


To examine the moment images, use viewer; the resulting moment zero image is displayed at right. Note that you may have to play with the color map (Data Display Options button in viewer) in order to replicate the image in this tutorial. 

viewer(infile='ngc5921.demo.moments.integrated')

Export the Data

To export the (u, v) data and image cube as FITS files, use exportuvfits and exportfits, respectively.

Here's how to export the continuum-subtracted (u, v) data. 

exportuvfits(vis='ngc5921.demo.src.split.ms.contsub', fitsfile='ngc5921.demo.contsub.uvfits', datacolumn='corrected', multisource=True)


And now, the FITS cube. 

exportfits(imagename='ngc5921.demo.tclean.image', fitsimage='ngc5921.demo.tclean.fits')


The moment maps (or any CASA images) can be similarly exported using exportfits.


Appendix: Python Notes

os.system

os.system allows you to run shell commands from within python / CASA. For example: 

import os
os.system('ls -sF')

will give an OS-level listing of the current directory's contents.

os.environ.get

It's worth having a look at the output of the os.environ.get command to understand the python syntax (alternative: os.getenv). You can think of os.environ as a list of operating system environment variables, and get is a method that extracts information about the requested environment variable, here, CASAPATH. Get returns a string of whitespace separated information, and .split() turns that string into a list. The array index [0] extracts the first element of that list, which contains the path.

To illustrate, here is some example python I/O in CASA. 

CASA <12>: print os.environ.get('CASAPATH')
/usr/lib64/casapy/30.0.9709test-001 linux local el5bld64b

CASA <13>: print os.environ.get('CASAPATH').split()
['/usr/lib64/casapy/30.0.9709test-001', 'linux', 'local', 'el5bld64b']

CASA <14>: print os.environ.get('CASAPATH').split()[0]
/usr/lib64/casapy/30.0.9709test-001

↵ Pre-upgrade VLA Tutorials

Last checked on CASA Version 5.7.2.

Retrieved from "https://casaguides.nrao.edu/index.php?title=NGC_5921:_red-shifted_HI_emission_5.7.2&oldid=34800"