{{VLA Intro}}
Category: VLACategory: CalibrationCategory: ImagingCategory: Spectral Line= 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 [http://casa.nrao.edu/Doc/Cookbook/casa_cookbook.pdf 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.
= Setting up the CASA Environment =
Start up CASA in the directory you want to use.
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.
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.
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.
'''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.
= Import the Data =
The next step is to import the multisource UVFITS data to a CASA measurement set via the importuvfits filler.
saveinputs | 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.
[[File:n5921-01.png | thumb | Logger output of listobs.]]
The output goes to the logger window; see the screenshot at right.
'''Tip:''' You can control the text size of the logger window using -A (smaller font) and -L (larger font).
A more complete listing of the listobs output follows.
INFO listobs ##### Begin Task: listobs #####
INFO listobs::::casa
================================================================================
MeasurementSet Name: /DATA/ASHLESHA_2/NGC5921/ngc5921.demo.ms MS Version 2
================================================================================
Observer: TEST Project:
Observation: VLA
Data records: 22653 Total integration time = 5280 seconds
Observed from 13-Apr-1995/09:19:00.0 to 13-Apr-1995/10:47:00.0 (TAI)
2009-12-03 16:15:38 INFO listobs::ms::summary
ObservationID = 0 ArrayID = 0
Date Timerange (TAI) Scan FldId FieldName nVis Int(s) SpwIds
13-Apr-1995/09:19:00.0 - 09:24:30.0 1 0 1331+305000* 4509 30 [0]
09:27:30.0 - 09:29:30.0 2 1 1445+099000* 1890 30 [0]
09:33:00.0 - 09:48:00.0 3 2 N5921_2 6048 30 [0]
09:50:30.0 - 09:51:00.0 4 1 1445+099000* 756 30 [0]
10:22:00.0 - 10:23:00.0 5 1 1445+099000* 1134 30 [0]
10:26:00.0 - 10:43:00.0 6 2 N5921_2 6804 30 [0]
10:45:30.0 - 10:47:00.0 7 1 1445+099000* 1512 30 [0]
(nVis = Total number of time/baseline visibilities per scan)
Fields: 3
ID Code Name RA Decl Epoch SrcId nVis
0 C 1331+305000* 13:31:08.2873 +30.30.32.9590 J2000 0 4509
1 A 1445+099000* 14:45:16.4656 +09.58.36.0730 J2000 1 5292
2 N5921_2 15:22:00.0000 +05.04.00.0000 J2000 2 12852
(nVis = Total number of time/baseline visibilities per field)
Spectral Windows: (1 unique spectral windows and 1 unique polarization setups)
SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs
0 63 LSRK 1412.66507 24.4140625 1550.19688 1413.42801 RR LL
Sources: 3
ID Name SpwId RestFreq(MHz) SysVel(km/s)
0 1331+305000* 0 1420.405752 0
1 1445+099000* 0 1420.405752 0
2 N5921_2 0 1420.405752 0
Antennas: 27:
ID Name Station Diam. Long. Lat.
0 1 VLA:N7 25.0 m -107.37.07.2 +33.54.12.9
1 2 VLA:W1 25.0 m -107.37.05.9 +33.54.00.5
2 3 VLA:W2 25.0 m -107.37.07.4 +33.54.00.9
3 4 VLA:E1 25.0 m -107.37.05.7 +33.53.59.2
4 5 VLA:E3 25.0 m -107.37.02.8 +33.54.00.5
5 6 VLA:E9 25.0 m -107.36.45.1 +33.53.53.6
6 7 VLA:E6 25.0 m -107.36.55.6 +33.53.57.7
7 8 VLA:W8 25.0 m -107.37.21.6 +33.53.53.0
8 9 VLA:N5 25.0 m -107.37.06.7 +33.54.08.0
9 10 VLA:W3 25.0 m -107.37.08.9 +33.54.00.1
10 11 VLA:N4 25.0 m -107.37.06.5 +33.54.06.1
11 12 VLA:W5 25.0 m -107.37.13.0 +33.53.57.8
12 13 VLA:N3 25.0 m -107.37.06.3 +33.54.04.8
13 14 VLA:N1 25.0 m -107.37.06.0 +33.54.01.8
14 15 VLA:N2 25.0 m -107.37.06.2 +33.54.03.5
15 16 VLA:E7 25.0 m -107.36.52.4 +33.53.56.5
16 17 VLA:E8 25.0 m -107.36.48.9 +33.53.55.1
17 18 VLA:W4 25.0 m -107.37.10.8 +33.53.59.1
18 19 VLA:E5 25.0 m -107.36.58.4 +33.53.58.8
19 20 VLA:W9 25.0 m -107.37.25.1 +33.53.51.0
20 21 VLA:W6 25.0 m -107.37.15.6 +33.53.56.4
21 22 VLA:E4 25.0 m -107.37.00.8 +33.53.59.7
23 24 VLA:E2 25.0 m -107.37.04.4 +33.54.01.1
24 25 VLA:N6 25.0 m -107.37.06.9 +33.54.10.3
25 26 VLA:N9 25.0 m -107.37.07.8 +33.54.19.0
26 27 VLA:N8 25.0 m -107.37.07.5 +33.54.15.8
27 28 VLA:W7 25.0 m -107.37.18.4 +33.53.54.8
2009-12-03 16:15:38 INFO listobs::::casa
##### End Task: listobs #####
##########################################
= Key Information from listobs | 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 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.
= Interactive Flagging =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.
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.
* 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.
* 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 Data flagging with plotms | 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 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.
* #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.''
* #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.''
* #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.
= Setting the Flux Scale =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 also looks up the radio SED for common flux calibrators and automatically assigns the total flux density.
A summary of the operation is sent to the logger window. Here's a listing of the output.
2009-12-07 15:44:22 INFO setjy::::casa
2009-12-07 15:44:22 INFO setjy ##########################################
2009-12-07 15:44:22 INFO setjy ##### Begin Task: setjy #####
2009-12-07 15:44:22 INFO setjy::::casa
2009-12-07 15:44:22 INFO setjy Adding MODEL_DATA, CORRECTED_DATA and IMAGING_WEIGHT columns
2009-12-07 15:44:23 INFO setjy Initializing MODEL_DATA (to unity) and CORRECTED_DATA (to DATA)
2009-12-07 15:44:25 INFO setjy Initialized 22653 rows.
2009-12-07 15:44:25 INFO setjy Initializing natural weights
2009-12-07 15:44:25 INFO setjy 1331+30500002_0 spwid= 0 [I=14.76, Q=0, U=0, V=0] Jy, (Perley-Taylor 99)
2009-12-07 15:44:25 INFO setjy Selecting data
2009-12-07 15:44:25 INFO setjy Selected 4509 out of 22653 visibilities.
2009-12-07 15:44:25 INFO setjy Fourier transforming: replacing MODEL_DATA column
2009-12-07 15:44:25 INFO setjy Processing after subtracting componentlist /DATA/ASHLESHA_2/NGC5921/ngc5921.demo.ms.1331+30500002_0.spw0.tempcl
2009-12-07 15:44:25 INFO setjy Performing interferometric gridding with convolution function SF
2009-12-07 15:44:25 INFO setjy::::casa
2009-12-07 15:44:25 INFO setjy ##### End Task: setjy #####
2009-12-07 15:44:25 INFO setjy ##########################################
= 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.
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''''.
= 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.]]
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.
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.
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.
If the '''scriptmode''' is set to '''False''', the plot is saved to ngc5921.demo.plotcal.png.
= Gain Calibration =
From #Inspect the Bandpass Response Curve | 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.
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 Calibrating a VLA 5 GHz continuum survey#Calibration | 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.
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.
2009-12-07 16:54:27 INFO fluxscale ##########################################
2009-12-07 16:54:27 INFO fluxscale ##### Begin Task: fluxscale #####
2009-12-07 16:54:27 INFO fluxscale::::casa
2009-12-07 16:54:27 INFO fluxscale Opening MS: ngc5921.demo.ms for calibration.
2009-12-07 16:54:27 INFO fluxscale Initializing nominal selection to the whole MS.
2009-12-07 16:54:27 INFO fluxscale Beginning fluxscale--(MSSelection version)-------
2009-12-07 16:54:27 INFO fluxscale Found reference field(s): 1331+30500002_0
2009-12-07 16:54:27 INFO fluxscale Found transfer field(s): 1445+09900002_0
2009-12-07 16:54:27 INFO fluxscale Flux density for 1445+09900002_0 in SpW=0 is: 2.48581 +/- 0.00119412 (SNR = 2081.71, nAnt= 27)
2009-12-07 16:54:27 INFO fluxscale Storing result in ngc5921.demo.fluxscale
2009-12-07 16:54:27 INFO fluxscale Writing solutions to table: ngc5921.demo.fluxscale
2009-12-07 16:54:27 INFO fluxscale::::casa
2009-12-07 16:54:27 INFO fluxscale ##### End Task: fluxscale #####
2009-12-07 16:54:27 INFO fluxscale ##########################################
= Inspect the Calibration Solutions =
[[File:ngc5921-gaincal.png|thumb|Gain calibration solutions from gaincal and fluxscale.]]
Now inspect the results of gaincal. The setup is identical to that used to plot the #Inspect the Bandpass Response Curve | 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.
= 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 Calibrating a VLA 5 GHz continuum survey|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.
In the script snippet below, the python global variables ''ftable'' and ''btable'' replace the full table names.
= Plot the Spectrum =
[[File:spectrum-plotting.png|thumb|plotms|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.
The setup 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). Following the (Tab)'''Command''' convention of the Data flagging with casaplotms|flagging tutorial, here are the settings we need.
* (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
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.
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.
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 =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 clean|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''''.
* 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 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''''.
It will take a few minutes to generate the data cube. Use imhead to look at the cube header.
The output, as follows, appears in the logger window.
INFO imhead ##########################################
INFO imhead ##### Begin Task: imhead #####
INFO imhead::::casa
INFO ImageAnalysis::summary
INFO ImageAnalysis Image name : ngc5921.demo.cleanimg.image
INFO ImageAnalysis Object name : N5921_2
INFO ImageAnalysis Image type : PagedImage
INFO ImageAnalysis Image quantity : Intensity
INFO ImageAnalysis Pixel mask(s) : None
INFO ImageAnalysis Region(s) : None
INFO ImageAnalysis Image units : Jy/beam
INFO ImageAnalysis Restoring Beam : 52.3782 arcsec, 45.7319 arcsec, -165.572 deg
INFO ImageAnalysis::summary
INFO ImageAnalysis Direction reference : J2000
INFO ImageAnalysis Spectral reference : LSRK
INFO ImageAnalysis Velocity type : RADIO
INFO ImageAnalysis Rest frequency : 1.42041e+09 Hz
INFO ImageAnalysis Pointing center : 15:22:00.000000 +05.04.00.000000
INFO ImageAnalysis Telescope : VLA
INFO ImageAnalysis Observer : TEST
INFO ImageAnalysis Date observation : 1995/04/13/00:00:00
INFO ImageAnalysis Telescope position: [-1.60119e+06m, -5.04198e+06m, 3.55488e+06m] (ITRF)
INFO ImageAnalysis::summary+
INFO ImageAnalysis Axis Coord Type Name Proj Shape Tile Coord value at pixel Coord incr Units
INFO ImageAnalysis ------------------------------------------------------------------------------------------------
INFO ImageAnalysis 0 0 Direction Right Ascension SIN 256 64 15:22:00.000 128.00 -1.500000e+01 arcsec
INFO ImageAnalysis 1 0 Direction Declination SIN 256 64 +05.04.00.000 128.00 1.500000e+01 arcsec
INFO ImageAnalysis 2 1 Stokes Stokes 1 1 I
INFO ImageAnalysis 3 2 Spectral Frequency 46 8 1.41279e+09 0.00 2.4414062e+04 Hz
INFO ImageAnalysis Velocity 1607.99 0.00 -5.152860e+00 km/s
INFO imhead::::casa
INFO imhead ##### End Task: imhead #####
INFO imhead ##########################################
= 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.
= Cube Statistics =imstat|Imstat is the tool for displaying statistics of images and cubes. The following example displays the statistics for the whole cube.
The output goes to the logger window.
INFO imstat ##########################################
INFO imstat ##### Begin Task: imstat #####
INFO imstat::::casa
INFO imstat Regions ---
INFO imstat -- bottom-left corner (pixel) [blc]: [0, 0, 0, 0]
INFO imstat -- top-right corner (pixel) [trc]: [255, 255, 0, 45]
INFO imstat -- bottom-left corner (world) [blcf]: 15:24:08.404, +04.31.59.181, I, 1.41279e+09Hz
INFO imstat -- top-right corner( world) [trcf]: 15:19:52.390, +05.35.44.246, I, 1.41389e+09Hz
INFO imstat Values ---
INFO imstat -- flux density [flux]: 4.06274 Jy
INFO imstat -- number of points [npts]: 3.01466e+06
INFO imstat -- maximum value [max]: 0.0532876
INFO imstat -- minimum value [min]: -0.010393
INFO imstat -- position of max value (pixel) [maxpos]: [134, 134, 0, 38]
INFO imstat -- position of min value (pixel) [minpos]: [117, 0, 0, 21]
INFO imstat -- position of max value (world) [maxposf]: 15:21:53.976, +05.05.29.998, I, 1.41371e+09Hz
INFO imstat -- position of min value (world) [maxposf]: 15:22:11.035, +04.31.59.966, I, 1.4133e+09Hz
INFO imstat -- Sum of pixel values [sum]: 49.0089
INFO imstat -- Sum of squared pixel values [sumsq]: 12.1519
INFO imstat::::
INFO imstat Statistics ---
INFO imstat -- Mean of the pixel values [mean]: 1.62569e-05
INFO imstat -- Variance of the pixel values : 4.03068e-06
INFO imstat -- Standard deviation of the Mean [sigma]: 0.00200765
INFO imstat -- Root mean square [rms]: 0.00200772
INFO imstat -- Median of the pixel values [median]: -1.25504e-05
INFO imstat -- Median of the deviations [medabsdevmed]: 0.00126038
INFO imstat -- Quartile [quartile]: 0.00252074
INFO imstat::::casa
INFO imstat ##### End Task: imstat #####
= The Integrated Spectrum =
[[File:spectrum-integrated.png|thumb|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 #Plot the Spectrum|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.
To generate the integrated spectrum, perform the following tasks.
* Use the player controls File:vcrNext.png 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 File:drawingSelector.png 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 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.
= Cube Moments =
[[File:spectrum-momzero.png|thumb|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|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''.
To examine the moment images, use viewer; the resulting moment zero image is displayed at right.
= 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. Note that this snippet has exportuvfits spawned to the background (async = True).
And now, the FITS cube.
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