MRK 6: red-shifted HI absorption 5.5.0: Difference between revisions

From CASA Guides
Jump to navigationJump to search
No edit summary
(Updating to CASA 5.5.0)
Line 6: Line 6:
== Overview ==
== Overview ==


This example demonstrates a trickier spectral line data set. In this case, several radio bright active galaxies were observed to look for redshifted 21 cm absorption. The data were obtained using 4IF mode, which means that the full range of velocities were split into two spectral windows (in this example, CASA identifies them as spectral windows 14 and 15). Producing the data cube will require eventually stitching these windows back together. Happily, {{clean}} includes flexible regridding and interpolation in the spectral domain. {{clean|Clean}} further allows you to adjust the velocity reference frame.
This example demonstrates a trickier spectral line data set. In this case, several radio bright active galaxies were observed to look for redshifted 21 cm absorption. The data were obtained using 4IF mode, which means that the full range of velocities were split into two spectral windows (in this example, CASA identifies them as spectral windows 10 and 11). Producing the data cube will require eventually stitching these windows back together. Happily, {{tclean}} includes flexible regridding and interpolation in the spectral domain. {{tclean}} further allows you to adjust the velocity reference frame.


== Retrieve the data from the NRAO archive ==
== Retrieve the data from the NRAO archive ==
Line 67: Line 67:


<source lang="python">
<source lang="python">
# In CASA
# generate a listobs file with information about the observations
# generate a listobs file with information about the observations
listobs(vis='ab658.ms', listfile='ab658.listobs')
listobs(vis='ab658.ms', listfile='ab658.listobs')
Line 74: Line 75:


<source lang="python">
<source lang="python">
# In CASA
# plot antenna positions
# plot antenna positions
plotants(vis="ab658.ms", figfile="ab658_plotants.png")
plotants(vis="ab658.ms", figfile="ab658_plotants.png")
Line 81: Line 83:
For the purposes of this tutorial, we are interested only in the source [http://en.wikipedia.org/wiki/B._E._Markarian Markarian] [http://nedwww.ipac.caltech.edu/cgi-bin/nph-objsearch?objname=mrk+6&extend=no&hconst=73&omegam=0.27&omegav=0.73&corr_z=1&out_csys=Equatorial&out_equinox=J2000.0&obj_sort=RA+or+Longitude&of=pre_text&zv_breaker=30000.0&list_limit=5&img_stamp=YES 6], but the reduction techniques could be applied to any of the sources in the measurement set.  
For the purposes of this tutorial, we are interested only in the source [http://en.wikipedia.org/wiki/B._E._Markarian Markarian] [http://nedwww.ipac.caltech.edu/cgi-bin/nph-objsearch?objname=mrk+6&extend=no&hconst=73&omegam=0.27&omegav=0.73&corr_z=1&out_csys=Equatorial&out_equinox=J2000.0&obj_sort=RA+or+Longitude&of=pre_text&zv_breaker=30000.0&list_limit=5&img_stamp=YES 6], but the reduction techniques could be applied to any of the sources in the measurement set.  


[[File:MRK6-shot3.png|thumb|{{listobs|Listobs}} output for AB658. The data for MRK 6 and its calibrators are highlighted.]]
[[File:MKN6_Listobs.png|thumb|{{listobs|Listobs}} output for AB658. The data for MKN 6 and its calibrators are highlighted.]]


Using {{listobs}} and {{plotants}} we learn the source names, spectral windows, and viable reference antennas. A screenshot of the {{listobs}} output is given at right.
Using {{listobs}} and {{plotants}} we learn the source names, spectral windows, and viable reference antennas. A screenshot of the {{listobs}} output is given at right.
Line 88: Line 90:
| '''Source''' || MKN6
| '''Source''' || MKN6
|-
|-
| '''Amp Cal''' || 0134+329 = 3C48
| '''Flux Density Scale Cal''' || 0134+329 = 3C48
|-
|-
| '''Phase Cal''' || 1003+830
| '''Complex Gain Cal''' || 1003+830
|-  
|-  
| '''Central Antennas''' || VA27, VA09, VA23
| '''Central Antennas''' || VA27, VA09, VA23
Line 102: Line 104:
<source lang="python">
<source lang="python">
sourceName = "MKN6"
sourceName = "MKN6"
phaseCal = "1003+830"
CG_Cal = "1003+830"
ampCal = "0134+329"
FDS_Cal = "0134+329"
refAnt = "VA27" # or perhaps "VA09" or "VA23"
refAnt = "VA27" # or perhaps "VA09" or "VA23"
spwStr = '10, 11'
spwStr = '10, 11'
Line 116: Line 118:
== Editing the Data ==
== Editing the Data ==


[[File:MRK6-shot1.png|thumbnail|Plotms display of (uncalibrated) visibility amplitudes vs. time for the source MKN6.]]
[[File:MKN6_Amp_vs_Time.png|thumbnail|Plotms display of (uncalibrated) visibility amplitudes vs. time for the source MKN 6.]]


First inspect the data using {{plotms}}.  
First inspect the data using {{plotms}}.  


<source lang="python">
<source lang="python">
# In CASA
plotms(vis='ab658.ms', field='MKN6', xaxis='time', yaxis='amp')
plotms(vis='ab658.ms', field='MKN6', xaxis='time', yaxis='amp')
</source>
</source>
Line 139: Line 142:
Calibration of HI absorption spectra follows the basic techniques outlined in the [[Calibrating a VLA 5 GHz continuum survey|continuum tutorial]] and [[NGC 5921: red-shifted HI emission|HI emission tutorial]].
Calibration of HI absorption spectra follows the basic techniques outlined in the [[Calibrating a VLA 5 GHz continuum survey|continuum tutorial]] and [[NGC 5921: red-shifted HI emission|HI emission tutorial]].


First, set the flux of the amplitude calibrator. Notice that we are being careful to keep track of the spectral windows and not using the default ''spw''.
First, set the flux denisty of the flux density scale calibrator. Notice that we are being careful to keep track of the spectral windows and not using the default ''spw''.


<source lang="python">
<source lang="python">
# set the flux standard
# In CASA
# set the flux denisty scale standard
setjy(vis='ab658.ms', field='0134+329', spw='10,11', modimage='')
setjy(vis='ab658.ms', field='0134+329', spw='10,11', modimage='')
</source>
</source>
Line 149: Line 153:


<source lang="python">
<source lang="python">
# In CASA
# bandpass calibration
# bandpass calibration
bandpass(vis='ab658.ms', spw='10,11', caltable='ab658.bcal', gaintable='', gainfield='', interp='', field='0134+329', selectdata=False, bandtype='B', solint='inf', combine='scan', refant='VA27')
bandpass(vis='ab658.ms', spw='10,11', caltable='ab658.bcal', gaintable='', gainfield='', interp='', field='0134+329', selectdata=False, bandtype='B', solint='inf', combine='scan', refant='VA27')
Line 155: Line 160:
<!--
<!--
<source lang="python">
<source lang="python">
# set the flux standard
# set the flux density scale standard
default('setjy')
default('setjy')
vis = measurementSet
vis = measurementSet
field = ampCal
field = FDS_Cal
spw = spwStr # this is key to getting the calibration right!
spw = spwStr # this is key to getting the calibration right!
modimage = ''
modimage = ''
Line 170: Line 175:
btable = 'ab658.bcal' # bandpass
btable = 'ab658.bcal' # bandpass
gtable = 'ab658.gcal' # gain
gtable = 'ab658.gcal' # gain
ftable = sourceName + '.fluxscale' # flux
ftable = sourceName + '.fluxscale' # flux density scale
#
#
if os.path.exists(btable): rmtables(btable)
if os.path.exists(btable): rmtables(btable)
Line 186: Line 191:
gainfield = ''
gainfield = ''
interp = ''
interp = ''
field = ampCal
field = FDS_Cal
selectdata = False
selectdata = False
opacity = 0.0
opacity = 0.0
Line 197: Line 202:
-->
-->


[[File:MRK6-shot2.png|thumb|Bandpass calibration curves. Note that the response curve is flat between channels 2 and 27.]]
[[File:MKN6_PlotCal_Bandpass_Cal_Curves.png|thumb|Bandpass calibration curves. Note that the response curve is flat between channels 2 and 27.]]


In order to perform the amplitude and phase calibrations, we'll need to inspect the bandpass response curve to select appropriate channels to average. We want to choose those channels where the response is relatively flat. The result of plotcal is shown at right: the bandpass is flat over channels 2:27.  
In order to perform the flux density scale and complex gain calibrations, we'll need to inspect the bandpass response curve to select appropriate channels to average. We want to choose those channels where the response is relatively flat. The result of plotms is shown at right: the bandpass is flat over channels 2:27.  


<source lang="python">
<source lang="python">
Line 207: Line 212:
# Press Enter or Return, then copy/paste all of the following:
# Press Enter or Return, then copy/paste all of the following:
# plot the bandpass calibration, amp vs channel
# plot the bandpass calibration, amp vs channel
plotcal(caltable='ab658.bcal', spw='10,11', field='0134+329', subplot=211, yaxis='amp', showgui=True)
plotms(vis='ab658.bcal', spw='10~11', field='0134+329', gridrows=2, gridcols=1, rowindex=0, plotindex=0, yaxis='gainamp', xaxis='channel', coloraxis='Antenna1')


# plot the bandpass calibration, gain phase vs channel
# plot the bandpass calibration, gain phase vs channel
plotcal(caltable='ab658.bcal', spw='10,11', field='0134+329', subplot=212, yaxis='phase', showgui=True)
plotms(vis='ab658.bcal', spw='10~11', field='0134+329', gridrows=2, gridcols=1, rowindex=1, plotindex=1, yaxis='gainphase', xaxis='channel', coloraxis='Antenna1', clearplots=False)
--
--
</source>
</source>
Line 220: Line 225:
spw = spwStr
spw = spwStr
caltable = btable
caltable = btable
field = ampCal
field = FDS_Cal
subplot = 211
subplot = 211
yaxis = 'amp'
yaxis = 'amp'
Line 231: Line 236:
-->
-->


The following steps generate amplitude and phase solutions for the calibrators, and the fluxscale is bootstrapped from the amplitude calibrator onto the phase calibrator, effectively turning the phase calibrator into a local flux calibrator for the science target.
The following steps generate amplitude and complex gain solutions for the calibrators, and the fluxscale is bootstrapped from the flux density scale calibrator onto the complex gain calibrator, effectively turning the complex gain calibrator into a local flux density scale calibrator for the science target.


<source lang="python">
<source lang="python">
# generate amplitude and gain solutions for phase and amplitude calibrators
# In CASA
# generate amplitude and gain solutions for complex gain and flux density scale calibrators
gaincal(vis='ab658.ms', caltable='ab658.gcal', gaintable='ab658.bcal', gainfield='', interp='nearest', field='0134+329, 1003+830',  
gaincal(vis='ab658.ms', caltable='ab658.gcal', gaintable='ab658.bcal', gainfield='', interp='nearest', field='0134+329, 1003+830',  
         spw='10:2~27, 11:2~27', selectdata=False, gaintype='G', solint='inf', combine='', calmode='ap', minsnr=1.0, refant='VA27')
         spw='10:2~27, 11:2~27', selectdata=False, gaintype='G', solint='inf', combine='', calmode='ap', minsnr=1.0, refant='VA27')
Line 240: Line 246:


<source lang="python">
<source lang="python">
# Bootstrap flux scale
# In CASA
# Bootstrap flux density scale
fluxscale(vis='ab658.ms', fluxtable='MKN6.fluxscale', caltable='ab658.gcal', reference='0134+329', transfer='1003+830')
fluxscale(vis='ab658.ms', fluxtable='MKN6.fluxscale', caltable='ab658.gcal', reference='0134+329', transfer='1003+830')
</source>
</source>
Line 252: Line 259:
gainfield = ''
gainfield = ''
interp = 'nearest'
interp = 'nearest'
field = ampCal + ',' + phaseCal
field = FDS_Cal + ',' + CG_Cal
# specify data only from the flat region of the bandpass response
# specify data only from the flat region of the bandpass response
spw = '10:2~27,11:2~27'
spw = '10:2~27,11:2~27'
Line 266: Line 273:
#=====================================================================
#=====================================================================
#
#
# Bootstrap flux scale
# Bootstrap flux density scale
#
#
vis = measurementSet
vis = measurementSet
Line 272: Line 279:
fluxtable = ftable
fluxtable = ftable
caltable = gtable
caltable = gtable
reference = ampCal
reference = FDS_Cal
transfer = phaseCal
transfer = CG_Cal
fluxscale()
fluxscale()
</source>
</source>
-->
-->


The task {{fluxscale}} will give the bootstrapped flux density of the phase calibrator. If all goes well, something like the following should appear in the log window.
The task {{fluxscale}} will give the bootstrapped flux density of the complex gain calibrator. If all goes well, something like the following should appear in the log window.


<pre>
<pre>
Line 284: Line 291:
  Found reference field(s): 0134+329
  Found reference field(s): 0134+329
  Found transfer field(s):  1003+830
  Found transfer field(s):  1003+830
  Flux density for 1003+830 in SpW=10 (freq=1.39215e+09 Hz) is: 0.492631 +/- 0.00341773 (SNR = 144.14, N = 52)
...
  Flux density for 1003+830 in SpW=11 (freq=1.3943e+09 Hz) is: 0.489476 +/- 0.00286404 (SNR = 170.904, N = 52)
  Flux density for 1003+830 in SpW=10 (freq=1.3921e+09 Hz) is: 0.492452 +/- 0.00331148 (SNR = 148.711, N = 52)
  Flux density for 1003+830 in SpW=11 (freq=1.39425e+09 Hz) is: 0.489206 +/- 0.00287348 (SNR = 170.249, N = 52)
...
  Storing result in ab658.fluxscale
  Storing result in ab658.fluxscale
</pre>
</pre>
Line 292: Line 301:


<source lang="python">
<source lang="python">
# In CASA
# apply calibrations to science target
# apply calibrations to science target
applycal(vis='ab658.ms', gaintable=['MKN6.fluxscale','ab658.bcal'], gainfield=['1003+830', '*'], interp=['linear','nearest'], field='MKN6', spw='10,11', spwmap=[], selectdata=False)
applycal(vis='ab658.ms', gaintable=['MKN6.fluxscale','ab658.bcal'], gainfield=['1003+830', '*'], interp=['linear','nearest'], field='MKN6', spw='10,11', spwmap=[], selectdata=False)
Line 297: Line 307:


<source lang="python">
<source lang="python">
#In CASA
# split out the corrected data column (which is placed in the data column in the new MS) for our science target.
# split out the corrected data column (which is placed in the data column in the new MS) for our science target.
split(vis='ab658.ms', outputvis='ab658-MKN6.split.ms', field='MKN6', spw='10:2~27, 11:2~27', datacolumn='corrected')
split(vis='ab658.ms', outputvis='ab658-MKN6.split.ms', field='MKN6', spw='10:2~27, 11:2~27', datacolumn='corrected')
Line 308: Line 319:
vis = measurementSet
vis = measurementSet
gaintable = [ftable,btable]
gaintable = [ftable,btable]
gainfield = [phaseCal,'*']
gainfield = [CG_Cal,'*']
interp = ['linear','nearest']
interp = ['linear','nearest']
field = sourceName
field = sourceName
Line 335: Line 346:
== Continuum Subtraction ==
== Continuum Subtraction ==


Ideally, we would use {{plotms}} to identify the channels where HI absorption is present (see [[NGC 5921: red-shifted HI emission#Plot the Spectrum|the HI emission tutorial.]]) The absorption line is however weak and isolated against one of a number of radio continuum sources in the data; as such, the absorption is not apparent in averaged visibilities. Instead, we'll have to produce a dirty cube using {{clean}} and inspect the spectrum in the image domain.
Ideally, we would use {{plotms}} to identify the channels where HI absorption is present (see [[NGC 5921: red-shifted HI emission#Plot the Spectrum|the HI emission tutorial.]]) The absorption line is however weak and isolated against one of a number of radio continuum sources in the data; as such, the absorption is not apparent in averaged visibilities. Instead, we'll have to produce a dirty cube using {{tclean}} and inspect the spectrum in the image domain.


=== Search for the Absorption Line ===
=== Search for the Absorption Line ===
Line 342: Line 353:


<source lang="python">
<source lang="python">
# In CASA
# image cube for spectral window 0
# image cube for spectral window 0
clean(vis='ab658-MKN6.split.ms', imagename='temp0', spw='0:0~25', mode='channel', nchan=-1, start=0, width=1, interpolation='nearest', niter=2000,  
tclean(vis='ab658-MKN6.split.ms', imagename='temp0', spw='0:0~25', specmode='cube', interpolation='nearest', niter=2000,  
       cell=['0.4arcsec', '0.4arcsec'], phasecenter='J2000 06h52m12.2 +74d25m37')
       cell=['0.4arcsec', '0.4arcsec'], phasecenter='J2000 06h52m12.2 +74d25m37')
</source>
</source>


<source lang="python">
<source lang="python">
# In CASA
# make a new image cube with spectral window 1
# make a new image cube with spectral window 1
clean(vis='ab658-MKN6.split.ms', imagename='temp1', spw='1:0~25', mode='channel', nchan=-1, start=0, width=1, interpolation='nearest', niter=2000,  
tclean(vis='ab658-MKN6.split.ms', imagename='temp1', spw='1:0~25', specmode='cube', interpolation='nearest', niter=2000,  
       cell=['0.4arcsec', '0.4arcsec'], phasecenter='J2000 06h52m12.2 +74d25m37')
       cell=['0.4arcsec', '0.4arcsec'], phasecenter='J2000 06h52m12.2 +74d25m37')
</source>
</source>
Line 375: Line 388:
-->
-->


It's not ordinarily necessary to give the phase center, but this is a good opportunity to use the regridding capabilities of {{clean}}. The data were originally taken in B1950 coordinates, and the observations followed the good practice of offsetting the pointing center slightly away from the source coordinates. Here, we specify the J2000 coordinates (looked up on [http://nedwww.ipac.caltech.edu NED]) and allow {{clean}} to handle the precession and centering that we require.
It's not ordinarily necessary to give the phase center, but this is a good opportunity to use the regridding capabilities of {{tclean}}. The data were originally taken in B1950 coordinates, and the observations followed the good practice of offsetting the pointing center slightly away from the source coordinates. Here, we specify the J2000 coordinates (looked up on [http://nedwww.ipac.caltech.edu NED]) and allow {{tclean}} to handle the precession and centering that we require.


The output cubes will go to '''temp0.image''' and '''temp1.image'''.
The output cubes will go to '''temp0.image''' and '''temp1.image'''.
Line 382: Line 395:


<source lang="python">
<source lang="python">
# In CASA
viewer(infile='temp0.image', displaytype='raster')
viewer(infile='temp0.image', displaytype='raster')
</source>
</source>
Line 393: Line 407:
-->
-->


[[File:MRK6-shot4-annotated.png|thumb|Spectral profile for MRK 6 in the 2nd spectral window, '''spw = 1'''. The location of the weak absorption line is marked by green arrow.]]
[[File:MKN6_SpectralProfile.png|thumb|Spectral profile for MKN 6 in the 2nd spectral window, '''spw = 1'''. The weak absorption line can be seen around channel 15.]]


Now, use {{viewer}} to generate a spectrum. The source spectrum in '''spw = 1''' is shown in the figure at right.
Now, use {{viewer}} to generate a spectrum. The source spectrum in '''spw = 1''' is shown in the figure at right.
Line 401: Line 415:
* Right-click and drag a box around the source; the spectrum will appear in the '''Image Profile''' window.
* Right-click and drag a box around the source; the spectrum will appear in the '''Image Profile''' window.
* On the '''Image Profile''' window, change the velocity convention to '''Channel''' using the pull-down menu.
* On the '''Image Profile''' window, change the velocity convention to '''Channel''' using the pull-down menu.
* ''Optional:'' Save the numerical data to a file. Click the save-as-text button [[File:Save-as-text-file.png]] and write the data to a file; for this example, we'll use '''temp0.txt'''.
* ''Optional:'' Save the numerical data to a file. Click the Save button [[File:Save.png]] and write the data to a file; for this example, we'll use '''temp0.txt'''.


Repeat for '''spw = 1'''.
Repeat for '''spw = 1'''.
Line 422: Line 436:


<source lang="python">
<source lang="python">
# In CASA
uvcontsub(vis='ab658-MKN6.split.ms', fitspw='0:0~25, 1:0~9; 21~25', solint='int', fitorder=1)
uvcontsub(vis='ab658-MKN6.split.ms', fitspw='0:0~25, 1:0~9; 21~25', solint='int', fitorder=1)
</source>
</source>
Line 445: Line 460:
== Imaging the Continuum Subtracted Cube ==
== Imaging the Continuum Subtracted Cube ==


[[File:MRK6-shot6.png|thumb|Viewer windows displaying the HI (21cm) absorption line cube of MRK6. The viewer window shows the channel containing the strongest absorption (at 5456 km/s), and the spectral profile window shows the average spectrum around the position of the strongest absorption.]]
[[File:MKN6_Viewer_and_HI_Line.png|thumb|Viewer windows displaying the HI (21cm) absorption line cube of MKN 6. The viewer window shows the channel containing the strongest absorption (at 5479 km/s), and the spectral profile window shows the average spectrum around the position of the strongest absorption.]]


Inspecting the preliminary (continuum-containing) cubes, temp0.image and temp1.image, that we generated [[#Search for the Absorption Line|above]], we can determine good parameters for the velocity axis. Here are the velocity options I came up with.
Inspecting the preliminary (continuum-containing) cubes, temp0.image and temp1.image, that we generated [[#Search for the Absorption Line|above]], we can determine good parameters for the velocity axis.


<source lang="python">
<source lang="python">
clean(vis='ab658-MKN6.split.ms.contsub', spw='', imagename='MKN6-HIABS', mode='velocity', nchan=-1, start='-252.349km/s', interpolation='linear', width='',           
# In CASA
       outframe='bary', veltype='optical', restfreq='1420.405752MHz', niter=2000, threshold='1.0mJy', weighting='natural', cell=['0.4arcsec', '0.4arcsec'],  
tclean(vis='ab658-MKN6.split.ms.contsub', spw='', imagename='MKN6-HIABS', specmode='cube', start='-252.349km/s', interpolation='linear',           
       outframe='BARY', veltype='optical', restfreq='1420.405752MHz', niter=2000, threshold='1.0mJy', weighting='natural', cell=['0.4arcsec', '0.4arcsec'],  
       phasecenter = 'J2000 06h52m12.2 +74d25m37')
       phasecenter = 'J2000 06h52m12.2 +74d25m37')
</source>
</source>


* width='' : Set the channel width equal to the first channel in the first spectral window.
* threshold='1.0mJy' : avoid over-cleaning since it appears that the channel rms is about 0.5 mJy/beam
* threshold='1.0mJy' : avoid over-cleaning since it appears that the channel rms is about 0.5 mJy/beam
* cell = ['0.4arcsec', '0.4arcsec'] : VLA-A configuration expects 1.5--2 arcsec beams. 0.4arcsec pixels ensure sampling the beam by better than 3x
* cell = ['0.4arcsec', '0.4arcsec'] : VLA-A configuration expects 1.5--2 arcsec beams. 0.4arcsec pixels ensure sampling the beam by better than 3x
Line 495: Line 510:
-->
-->


The ({{viewer}}) results are illustrated in the figure at upper right.
The ({{viewer}}) results are illustrated in the figure at right.


=== Imaging the Continuum Dataset ===
=== Imaging the Continuum Dataset ===


[[File:MRK6-shot7.png|thumb|21cm continuum image of MRK 6. The image has been scaled, using (Data)(Adjust)(Basic Settings)'''Scaling Power Cycles''' = -1.9, to bring out artifacts apparent in the background. The 6-pointed star pattern around the source is symptomatic of residual phase-errors, which can be healed using [[Calibrating a VLA 5 GHz continuum survey#Self-calibration|self-calibration]], and the arcuate stripes across the field are sidelobes of neighboring sources that need to be cleaned using [[Imaging Flanking Fields|flanking fields]].]]
[[File:MKN6_Continuum.png|thumb|21cm continuum image of MKN 6. The image has been scaled, using (Data)(Adjust)(Basic Settings) '''Scaling Power Cycles''' = -1.9, to bring out artifacts apparent in the background. The 6-pointed star pattern around the source is symptomatic of residual phase-errors, which can be healed using [[Calibrating a VLA 5 GHz continuum survey#Self-calibration|self-calibration]], and the arcuate stripes across the field are sidelobes of neighboring sources that need to be cleaned using [[Imaging Flanking Fields|flanking fields]].]]


We can similarly image the continuum data with just a few changes of the {{clean}} parameters. Here, we call clean by selecting the spectral windows and channels pertaining to only the continuum.
We can similarly image the continuum data with just a few changes of the {{tclean}} parameters. Here, we call tclean by selecting the spectral windows and channels pertaining to only the continuum.


<source lang="python">
<source lang="python">
clean(vis='ab658-MKN6.split.ms', mode='mfs', spw='0:0~25,1:0~9;21~25', imagename='MKN6-CONT', restfreq='1420.405752MHz', niter=2000,  
# In CASA
tclean(vis='ab658-MKN6.split.ms', specmode='mfs', spw='0:0~25,1:0~9;21~25', imagename='MKN6-CONT', restfreq='1420.405752MHz', niter=2000,  
       threshold='1.0mJy', weighting='natural', cell=['0.4arcsec', '0.4arcsec'], phasecenter = 'J2000 06h52m12.2 +74d25m37')
       threshold='1.0mJy', weighting='natural', cell=['0.4arcsec', '0.4arcsec'], phasecenter = 'J2000 06h52m12.2 +74d25m37')
</source>
</source>
Line 523: Line 539:
== Extract the Absorption Spectrum ==
== Extract the Absorption Spectrum ==


{{viewer|Viewer}} provides a handy first-look at the spectrum, but the extracted spectrum is an average over a spatial box. What we really want is the ''integrated'' spectrum, not the average. To do this, we'll use CASA's {{imstat}}, {{imval}}, and {{imhead}} tools to extract cube information to python variables and then let python build up the spectrum.
The {{viewer|Viewer}} provides a handy first-look at the spectrum, but the extracted spectrum is an average over a spatial box. What we really want is the ''integrated'' spectrum, not the average. To do this, we'll use CASA's {{imstat}}, {{imval}}, and {{imhead}} tools to extract cube information to python variables and then let python build up the spectrum.


=== Extract the Spectrum: load the cube data into python arrays ===
=== Extract the Spectrum: load the cube data into python arrays ===
Line 530: Line 546:


<source lang="python">
<source lang="python">
# In CASA
cubeStat=imstat("MKN6-HIABS.image")
cubeStat=imstat("MKN6-HIABS.image")
</source>
</source>
Line 536: Line 553:


<source lang="python">
<source lang="python">
# In CASA
cubeStat
cubeStat
</source>
</source>
Line 543: Line 561:
<pre>
<pre>
{'blc': array([0, 0, 0, 0], dtype=int32),
{'blc': array([0, 0, 0, 0], dtype=int32),
  'blcf': '06:52:24.903, +74.24.45.777, I, 1.39522e+09Hz',
  'blcf': '06:52:17.165, +74.25.16.997, I, 1.39542e+09Hz',
'flux': array([-0.05243921]),
  'max': array([ 0.00286199]),
  'max': array([ 0.00279196]),
  'maxpos': array([27, 68,  0,  6], dtype=int32),
  'maxpos': array([17382,  0,  33], dtype=int32),
  'maxposf': '06:52:14.485, +74.25.44.199, I, 1.39484e+09Hz',
  'maxposf': '06:52:07.732, +74.25.18.597, I, 1.39199e+09Hz',
  'mean': array([ -5.30860147e-06]),
  'mean': array([ -3.96712512e-07]),
  'medabsdevmed': array([ 0.00031676]),
  'medabsdevmed': array([ 0.00029433]),
  'median': array([ -3.58635123e-06]),
  'median': array([ 0.]),
  'min': array([-0.01634453]),
  'min': array([-0.01497807]),
  'minpos': array([49, 53, 0, 10], dtype=int32),
  'minpos': array([127, 131,   0,   7], dtype=int32),
  'minposf': '06:52:12.299, +74.25.38.200, I, 1.39444e+09Hz',
  'minposf': '06:52:12.299, +74.25.38.200, I, 1.39453e+09Hz',
  'npts': array([ 470000.]),
  'npts': array([ 3145728.]),
'q1': array([-0.00032042]),
  'quartile': array([ 0.00058867]),
'q3': array([ 0.00031312]),
  'rms': array([ 0.00047121]),
  'quartile': array([ 0.00063354]),
  'sigma': array([ 0.00047121]),
  'rms': array([ 0.00050688]),
  'sum': array([-1.24794966]),
  'sigma': array([ 0.00050685]),
  'sumsq': array([ 0.69847239]),
  'sum': array([-2.49504269]),
  'trc': array([255, 255,   0, 47], dtype=int32),
  'sumsq': array([ 0.12075419]),
  'trcf': '06:51:59.574, +74.26.27.778, I, 1.39063e+09Hz'}
  'trc': array([99, 99, 0, 46], dtype=int32),
  'trcf': '06:52:07.331, +74.25.56.597, I, 1.39093e+09Hz'}
</pre>
</pre>


Line 567: Line 586:


<source lang=python>
<source lang=python>
# In CASA
%cpaste
%cpaste


Line 578: Line 598:


<source lang=python>
<source lang=python>
# In CASA
%cpaste
%cpaste


Line 603: Line 624:


<div style="background-color: #dddddd;">
<div style="background-color: #dddddd;">
'''Comment:''' The absorption line source is unresolved by this observation, or, to put it another way, is entirely contained within a synthesized beam. The extraction at a single spatial pixel gives the integrated spectrum and therefore changes units from surface brightness (mJy/beam) to flux density (mJy). However, the data have been forced into a rectangular grid, and there's no guarantee that the absorption minimum centers within a pixel. Better to use {{imfit}} to determine more accurately the position of the minimum, and perhaps re-{{clean}} with a new phasecenter. To avoid complicating the present tutorial, we'll just stick with the minimum in the current {{clean}} grid.
'''Comment:''' The absorption line source is unresolved by this observation, or, to put it another way, is entirely contained within a synthesized beam. The extraction at a single spatial pixel gives the integrated spectrum and therefore changes units from surface brightness (mJy/beam) to flux density (mJy). However, the data have been forced into a rectangular grid, and there's no guarantee that the absorption minimum centers within a pixel. Better to use {{imfit}} to determine more accurately the position of the minimum, and perhaps re-{{tclean}} with a new phasecenter. To avoid complicating the present tutorial, we'll just stick with the minimum in the current {{tclean}} grid.
</div>
</div>


Line 639: Line 660:
=== Plot the spectrum using matplotlib ===
=== Plot the spectrum using matplotlib ===


[[File:MRK6-shot8.png|thumb|Absorption spectrum of MRK 6 in frequency space.]]
[[File:MKN6_Absorption_Spectrum.png|thumb|Absorption spectrum of MKN 6 in frequency space.]]


Now we have the spectrum stored in <tt>freqSpec</tt>, <tt>cubeSpec['data']</tt>. Use [http://matplotlib.sourceforge.net/api/pyplot_api.html matplotlib.pyplot] to plot the spectrum (see figure at right).
Now we have the spectrum stored in <tt>freqSpec</tt>, <tt>cubeSpec['data']</tt>. Use [http://matplotlib.sourceforge.net/api/pyplot_api.html matplotlib.pyplot] to plot the spectrum (see figure at right).
Line 693: Line 714:
-->
-->


[[File:MRK6-shot9.png|thumb|Publication quality spectrum of MRK 6.]]
[[File:MKN6_Absorption_Spectrum_Publication_Ready.png|thumb|Publication quality spectrum of MKN 6.]]


Plot the spectrum again, but this time try for publication-ready quality (see figure at right). This example is somewhat more clever to highlight the use of TeX typography in <tt>matplotlib</tt>.  Note that you will have to have a compatible version of LaTeX installed for this to be successful.
Plot the spectrum again, but this time try for publication-ready quality (see figure at right). This example is somewhat more clever to highlight the use of TeX typography in <tt>matplotlib</tt>.  Note that you will have to have a compatible version of LaTeX installed for this to be successful.


<source lang="python">
<source lang="python">
# In CASA
# set up fonts
# set up fonts
%cpaste
%cpaste
Line 741: Line 763:


<source lang="python">
<source lang="python">
# In CASA
# First, combine the spectrum into a single array
# First, combine the spectrum into a single array
%cpaste
%cpaste
Line 746: Line 769:
# Press Enter or Return, then copy/paste the following:
# Press Enter or Return, then copy/paste the following:
spectrum = numpy.vstack((velocity, cubeSpec['data']))
spectrum = numpy.vstack((velocity, cubeSpec['data']))
numpy.savetxt("mrk6.txt", numpy.transpose(spectrum))
numpy.savetxt("mkn6.txt", numpy.transpose(spectrum))
--
--
</source>
</source>
Line 754: Line 777:
# First, combine the spectrum into a single array
# First, combine the spectrum into a single array
spectrum = numpy.vstack((velocity, cubeSpec['data']))
spectrum = numpy.vstack((velocity, cubeSpec['data']))
numpy.savetxt("mrk6.txt", numpy.transpose(spectrum))
numpy.savetxt("mkn6.txt", numpy.transpose(spectrum))
</source>
</source>
-->
-->


Here's a listing of the first few lines of ''mrk6.txt''.  
Here's a listing of the first few lines of ''mkn6.txt''.  


<pre>
<pre>
5.379999999489542461e+03 -5.940478877164423466e-04
5.367264765747175261e+03 -2.500247210264205933e-04
5.399939886933506386e+03 -3.919676819350570440e-04
5.388621623358961187e+03 -2.288779069203883410e-04
5.419882380281690530e+03 -3.093948180321604013e-04
5.409981470534922664e+03 4.886248498223721981e-04
5.439827480044931690e+03 -1.092202146537601948e-03
5.431344307902921173e+03 3.546753214322961867e-05
5.459775186734200361e+03 -4.310846794396638870e-03
5.452710136090820924e+03 -1.000610413029789925e-03
5.479725500860732609e+03 -3.458642400801181793e-03
5.474078955726684399e+03 -2.162145450711250305e-03
5.499678422935699928e+03 -7.933050510473549366e-04
5.495450767438905132e+03 -1.413836842402815819e-03
5.519633953470603956e+03 -2.004600130021572113e-03
5.516825571855946691e+03 -9.065763442777097225e-04
5.539592092976883578e+03 -8.163747377693653107e-03
5.538203369606402703e+03 -3.753502853214740753e-03
5.559552841966238702e+03 -1.676503755152225494e-02
5.559584161319202394e+03 -1.479859836399555206e-02
</pre>
</pre>


Line 776: Line 799:


<source lang="python">
<source lang="python">
# In CASA
%cpaste
%cpaste


# Press Enter or Return, then copy/paste the following:
# Press Enter or Return, then copy/paste the following:
# import numpy # if you haven't already
# import numpy # if you haven't already
spectrum = numpy.loadtxt("mrk6.txt",dtype="float")
spectrum = numpy.loadtxt("mkn6.txt",dtype="float")
x = spectrum[:,0]
x = spectrum[:,0]
y = spectrum[:,1]
y = spectrum[:,1]
Line 786: Line 810:
</source>
</source>


{{Checked 5.1.1}}
{{Checked 5.5.0}}


<!--
<!--

Revision as of 13:47, 26 August 2019

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

Overview

This example demonstrates a trickier spectral line data set. In this case, several radio bright active galaxies were observed to look for redshifted 21 cm absorption. The data were obtained using 4IF mode, which means that the full range of velocities were split into two spectral windows (in this example, CASA identifies them as spectral windows 10 and 11). Producing the data cube will require eventually stitching these windows back together. Happily, tclean includes flexible regridding and interpolation in the spectral domain. tclean further allows you to adjust the velocity reference frame.

Retrieve the data from the NRAO archive

Download your data from the VLA Archive; in this example we'll use the publicly available survey AB658. These data were published in Gallimore et al. (1999).

Download all files associated with AB658. With the present archive defaults, you should have obtained the following VLA archive files:

'AB658_A921122.xp1'
'AB658_A921122.xp2'
'AB658_A921122.xp3'
'AB658_A921122.xp4'

This tutorial broadly follows the techniques described in the continuum survey tutorial and the 21 cm emission tutorial, and the basic calibration steps are presented only in outline.

Loading the Data

Recall that the python file globber glob is your friend here!

Note: IPython's (5.1.0) latest updates prevent copy/pasting several lines of code. To remedy this, we will use the %cpaste command to allow multiple lines to be placed in the terminal.

# In CASA
%cpaste

# Press Enter or Return, then copy/paste the following:
# sort the list of data files and place their names in a list variable
from glob import glob

fileList = sorted(glob("AB658*.xp?"))

# only import if ab658.ms doesn't yet exist, otherwise importvla will throw an error
importvla(archivefiles=fileList, vis='ab658.ms') 
--


We will now generate a listing of the details for the observations using the listobs task.

# In CASA
# generate a listobs file with information about the observations
listobs(vis='ab658.ms', listfile='ab658.listobs')

Next we will generate a plot of antenna positions. This will help identify a good reference antenna.

# In CASA
# plot antenna positions
plotants(vis="ab658.ms", figfile="ab658_plotants.png")

For the purposes of this tutorial, we are interested only in the source Markarian 6, but the reduction techniques could be applied to any of the sources in the measurement set.

listobs output for AB658. The data for MKN 6 and its calibrators are highlighted.

Using listobs and plotants we learn the source names, spectral windows, and viable reference antennas. A screenshot of the listobs output is given at right.

Source MKN6
Flux Density Scale Cal 0134+329 = 3C48
Complex Gain Cal 1003+830
Central Antennas VA27, VA09, VA23
Spectral Window IDs 10, 11


The selection of the spectral windows is key to getting the calibration right. Each science target in the multisource measurement set has a tuning respective of the source redshift, and the calibration observations are matched to those spectral windows. Were the spectral windows parameter (usually spw) left to its default value (all available spectral windows), many calibration routines would necessarily fail.


Editing the Data

Plotms display of (uncalibrated) visibility amplitudes vs. time for the source MKN 6.

First inspect the data using plotms.

# In CASA
plotms(vis='ab658.ms', field='MKN6', xaxis='time', yaxis='amp')


There is some obvious junk isolated in time; see the screenshot at right. Follow the tutorial to flag discrepant data for the source and calibrators.

Calibration

Calibration of HI absorption spectra follows the basic techniques outlined in the continuum tutorial and HI emission tutorial.

First, set the flux denisty of the flux density scale calibrator. Notice that we are being careful to keep track of the spectral windows and not using the default spw.

# In CASA
# set the flux denisty scale standard
setjy(vis='ab658.ms', field='0134+329', spw='10,11', modimage='')

Bandpass calibration comes next.

# In CASA
# bandpass calibration
bandpass(vis='ab658.ms', spw='10,11', caltable='ab658.bcal', gaintable='', gainfield='', interp='', field='0134+329', selectdata=False, bandtype='B', solint='inf', combine='scan', refant='VA27')


Bandpass calibration curves. Note that the response curve is flat between channels 2 and 27.

In order to perform the flux density scale and complex gain calibrations, we'll need to inspect the bandpass response curve to select appropriate channels to average. We want to choose those channels where the response is relatively flat. The result of plotms is shown at right: the bandpass is flat over channels 2:27.

# In CASA
%cpaste

# Press Enter or Return, then copy/paste all of the following:
# plot the bandpass calibration, amp vs channel
plotms(vis='ab658.bcal', spw='10~11', field='0134+329', gridrows=2, gridcols=1, rowindex=0, plotindex=0, yaxis='gainamp', xaxis='channel', coloraxis='Antenna1')

# plot the bandpass calibration, gain phase vs channel
plotms(vis='ab658.bcal', spw='10~11', field='0134+329', gridrows=2, gridcols=1, rowindex=1, plotindex=1, yaxis='gainphase', xaxis='channel', coloraxis='Antenna1', clearplots=False)
--


The following steps generate amplitude and complex gain solutions for the calibrators, and the fluxscale is bootstrapped from the flux density scale calibrator onto the complex gain calibrator, effectively turning the complex gain calibrator into a local flux density scale calibrator for the science target.

# In CASA
# generate amplitude and gain solutions for complex gain and flux density scale calibrators
gaincal(vis='ab658.ms', caltable='ab658.gcal', gaintable='ab658.bcal', gainfield='', interp='nearest', field='0134+329, 1003+830', 
        spw='10:2~27, 11:2~27', selectdata=False, gaintype='G', solint='inf', combine='', calmode='ap', minsnr=1.0, refant='VA27')
# In CASA
# Bootstrap flux density scale
fluxscale(vis='ab658.ms', fluxtable='MKN6.fluxscale', caltable='ab658.gcal', reference='0134+329', transfer='1003+830')


The task fluxscale will give the bootstrapped flux density of the complex gain calibrator. If all goes well, something like the following should appear in the log window.

Beginning fluxscale--(MSSelection version)-------
 Found reference field(s): 0134+329
 Found transfer field(s):  1003+830
 ...
 Flux density for 1003+830 in SpW=10 (freq=1.3921e+09 Hz) is: 0.492452 +/- 0.00331148 (SNR = 148.711, N = 52)
 Flux density for 1003+830 in SpW=11 (freq=1.39425e+09 Hz) is: 0.489206 +/- 0.00287348 (SNR = 170.249, N = 52)
 ...
 Storing result in ab658.fluxscale

We should inspect the calibration solutions to make sure they are reasonable, but we'll be cavalier and tread forward. Apply the calibrations to the science target and split that target into its own measurement set.

# In CASA
# apply calibrations to science target
applycal(vis='ab658.ms', gaintable=['MKN6.fluxscale','ab658.bcal'], gainfield=['1003+830', '*'], interp=['linear','nearest'], field='MKN6', spw='10,11', spwmap=[], selectdata=False)
#In CASA
# split out the corrected data column (which is placed in the data column in the new MS) for our science target.
split(vis='ab658.ms', outputvis='ab658-MKN6.split.ms', field='MKN6', spw='10:2~27, 11:2~27', datacolumn='corrected')
  • spw='10:2~27, 11:2~27 : Note that we have selected only the inner channels, 2 through 27 for spectral windows 10, 11.


Continuum Subtraction

Ideally, we would use plotms to identify the channels where HI absorption is present (see the HI emission tutorial.) The absorption line is however weak and isolated against one of a number of radio continuum sources in the data; as such, the absorption is not apparent in averaged visibilities. Instead, we'll have to produce a dirty cube using tclean and inspect the spectrum in the image domain.

Search for the Absorption Line

uvcontsub currently identifies continuum ranges by channel rather than frequency or velocity. We'll make two cubes, one for each spectral window. Note that the act of using split re-indexes the spectral windows and channels. Therefore spw='2~27' have been reindexed and are now spw='0~25' for both channels 10 & 11 (now reindexed as channels 0 & 1).

# In CASA
# image cube for spectral window 0
tclean(vis='ab658-MKN6.split.ms', imagename='temp0', spw='0:0~25', specmode='cube', interpolation='nearest', niter=2000, 
      cell=['0.4arcsec', '0.4arcsec'], phasecenter='J2000 06h52m12.2 +74d25m37')
# In CASA
# make a new image cube with spectral window 1
tclean(vis='ab658-MKN6.split.ms', imagename='temp1', spw='1:0~25', specmode='cube', interpolation='nearest', niter=2000, 
      cell=['0.4arcsec', '0.4arcsec'], phasecenter='J2000 06h52m12.2 +74d25m37')


It's not ordinarily necessary to give the phase center, but this is a good opportunity to use the regridding capabilities of tclean. The data were originally taken in B1950 coordinates, and the observations followed the good practice of offsetting the pointing center slightly away from the source coordinates. Here, we specify the J2000 coordinates (looked up on NED) and allow tclean to handle the precession and centering that we require.

The output cubes will go to temp0.image and temp1.image.

Now we can inspect the cube in viewer.

# In CASA
viewer(infile='temp0.image', displaytype='raster')


Spectral profile for MKN 6 in the 2nd spectral window, spw = 1. The weak absorption line can be seen around channel 15.

Now, use viewer to generate a spectrum. The source spectrum in spw = 1 is shown in the figure at right.

  • If necessary. use the player buttons VcrNext.png to scan through the channels until you see the source.
  • (Tools)Spectral Profile
  • Right-click and drag a box around the source; the spectrum will appear in the Image Profile window.
  • On the Image Profile window, change the velocity convention to Channel using the pull-down menu.
  • Optional: Save the numerical data to a file. Click the Save button Save.png and write the data to a file; for this example, we'll use temp0.txt.

Repeat for spw = 1.

viewer(infile='temp1.image', displaytype='raster')


Perform the Continuum Subtraction

Looking over the plots and text files, only channels 0~25 in either spectral window contain usable data. The absorption line is located in 1:15 and surrounding channels. Just to play it safe, we'll use the channels outside of 1:10~20 (spectral window 1, channels 10-20) to define the continuum for subtraction.

# In CASA
uvcontsub(vis='ab658-MKN6.split.ms', fitspw='0:0~25, 1:0~9; 21~25', solint='int', fitorder=1)


uvcontsub produces the uv subtracted measurement set:

  • ab658-MKN6.split.ms.contsub (continuum-subtracted visibilities).

Imaging the Continuum Subtracted Cube

Viewer windows displaying the HI (21cm) absorption line cube of MKN 6. The viewer window shows the channel containing the strongest absorption (at 5479 km/s), and the spectral profile window shows the average spectrum around the position of the strongest absorption.

Inspecting the preliminary (continuum-containing) cubes, temp0.image and temp1.image, that we generated above, we can determine good parameters for the velocity axis.

# In CASA
tclean(vis='ab658-MKN6.split.ms.contsub', spw='', imagename='MKN6-HIABS', specmode='cube', start='-252.349km/s', interpolation='linear',          
      outframe='BARY', veltype='optical', restfreq='1420.405752MHz', niter=2000, threshold='1.0mJy', weighting='natural', cell=['0.4arcsec', '0.4arcsec'], 
      phasecenter = 'J2000 06h52m12.2 +74d25m37')
  • threshold='1.0mJy' : avoid over-cleaning since it appears that the channel rms is about 0.5 mJy/beam
  • cell = ['0.4arcsec', '0.4arcsec'] : VLA-A configuration expects 1.5--2 arcsec beams. 0.4arcsec pixels ensure sampling the beam by better than 3x


The (viewer) results are illustrated in the figure at right.

Imaging the Continuum Dataset

21cm continuum image of MKN 6. The image has been scaled, using (Data)(Adjust)(Basic Settings) Scaling Power Cycles = -1.9, to bring out artifacts apparent in the background. The 6-pointed star pattern around the source is symptomatic of residual phase-errors, which can be healed using self-calibration, and the arcuate stripes across the field are sidelobes of neighboring sources that need to be cleaned using flanking fields.

We can similarly image the continuum data with just a few changes of the tclean parameters. Here, we call tclean by selecting the spectral windows and channels pertaining to only the continuum.

# In CASA
tclean(vis='ab658-MKN6.split.ms', specmode='mfs', spw='0:0~25,1:0~9;21~25', imagename='MKN6-CONT', restfreq='1420.405752MHz', niter=2000, 
      threshold='1.0mJy', weighting='natural', cell=['0.4arcsec', '0.4arcsec'], phasecenter = 'J2000 06h52m12.2 +74d25m37')


The result is shown at right. The continuum image needs some more help. Phase errors produce the 6-pointed star pattern surrounding the source and can be healed using self-calibration. The striping across the field is actually the combination of sidelobes from neighboring sources; these sidelobes can be suppressed by cleaning flanking fields.

Extract the Absorption Spectrum

The viewer provides a handy first-look at the spectrum, but the extracted spectrum is an average over a spatial box. What we really want is the integrated spectrum, not the average. To do this, we'll use CASA's imstat, imval, and imhead tools to extract cube information to python variables and then let python build up the spectrum.

Extract the Spectrum: load the cube data into python arrays

First, we need to determine where in the cube to extract the spectrum; equivalently, we need to find the location of the minimum value in the cube. imstat is the tool for this.

# In CASA
cubeStat=imstat("MKN6-HIABS.image")

The variable cubeStat is now a python dictionary. To see what it contains, just enter the variable name:

# In CASA
cubeStat

and the following summary of cubeStat appears.

{'blc': array([0, 0, 0, 0], dtype=int32),
 'blcf': '06:52:17.165, +74.25.16.997, I, 1.39542e+09Hz',
 'max': array([ 0.00286199]),
 'maxpos': array([27, 68,  0,  6], dtype=int32),
 'maxposf': '06:52:14.485, +74.25.44.199, I, 1.39484e+09Hz',
 'mean': array([ -5.30860147e-06]),
 'medabsdevmed': array([ 0.00031676]),
 'median': array([ -3.58635123e-06]),
 'min': array([-0.01634453]),
 'minpos': array([49, 53,  0, 10], dtype=int32),
 'minposf': '06:52:12.299, +74.25.38.200, I, 1.39444e+09Hz',
 'npts': array([ 470000.]),
 'q1': array([-0.00032042]),
 'q3': array([ 0.00031312]),
 'quartile': array([ 0.00063354]),
 'rms': array([ 0.00050688]),
 'sigma': array([ 0.00050685]),
 'sum': array([-2.49504269]),
 'sumsq': array([ 0.12075419]),
 'trc': array([99, 99,  0, 46], dtype=int32),
 'trcf': '06:52:07.331, +74.25.56.597, I, 1.39093e+09Hz'}

These entries can be accessed individually like an array variable. For example, the following commands copy the x and y positions of the absorption line to scalar variables.

# In CASA
%cpaste

# Press Enter or Return, then copy/paste the following:
x2extract = cubeStat['minpos'][0]
y2extract = cubeStat['minpos'][1]
--

Now, extract the spectrum using imval.

# In CASA
%cpaste

# Press Enter or Return, then copy/paste the following:
extractBox = "%d,%d" % (x2extract, y2extract) # copy the position to a string
cubeSpec = imval("MKN6-HIABS.image", box=extractBox, stokes='I')
--


cubeSpec is another python dictionary, and the spectrum is stored in cubeSpec['data'].

Comment: The absorption line source is unresolved by this observation, or, to put it another way, is entirely contained within a synthesized beam. The extraction at a single spatial pixel gives the integrated spectrum and therefore changes units from surface brightness (mJy/beam) to flux density (mJy). However, the data have been forced into a rectangular grid, and there's no guarantee that the absorption minimum centers within a pixel. Better to use imfit to determine more accurately the position of the minimum, and perhaps re-tclean with a new phasecenter. To avoid complicating the present tutorial, we'll just stick with the minimum in the current tclean grid.

Generate the frequency axis

We can derive the frequency axis from the header.

# In CASA
%cpaste

# Press Enter or Return, then copy/paste the following:
import numpy # make sure vector and array arithmetic options are loaded
cubeHead = imhead("MKN6-HIABS.image", mode="list")
nSpec = cubeStat['trc'][3] + 1 # get the number of frequency channels. 
f0 = float(cubeHead['crval4']) # reference freq in Hz
df = float(cubeHead['cdelt4']) # channel width in Hz
i0 = cubeHead['crpix4'] # reference pixel
freqSpec = (numpy.arange(nSpec) - i0)*df + f0
--


Plot the spectrum using matplotlib

Absorption spectrum of MKN 6 in frequency space.

Now we have the spectrum stored in freqSpec, cubeSpec['data']. Use matplotlib.pyplot to plot the spectrum (see figure at right).

# In CASA
%cpaste


# Press Enter or Return, then copy/paste the following:
import matplotlib.pyplot as plt
plt.clf() # clear the plot (figure)
plt.plot(freqSpec/1.e9, cubeSpec['data'], 'k-')
plt.xlabel("Freq (GHz)")
plt.ylabel("Flux Density (mJy)")
plt.show()
--


A publication-quality figure: plot the spectrum against velocity

Nice, now let's make a velocity axis instead. We'll use the optical convention to convert frequency to velocity; the radio convention is offered as an alternative but is commented out.

# In CASA
%cpaste

# Press Enter or Return, then copy/paste the following:
import matplotlib.pyplot as plt
fRest = cubeHead['restfreq'][0]
velocity = 299792.458 * (fRest/freqSpec - 1.0) # optical convention, km/s
# velocity = 299792.458 * (1.0 - freqSpec/fRest) # radio convention, km/s
--


Publication quality spectrum of MKN 6.

Plot the spectrum again, but this time try for publication-ready quality (see figure at right). This example is somewhat more clever to highlight the use of TeX typography in matplotlib. Note that you will have to have a compatible version of LaTeX installed for this to be successful.

# In CASA
# set up fonts
%cpaste


# Press Enter or Return, then copy/paste the following:
from matplotlib import rc
rc('font',**{'family':'serif','serif':['Times']})
rc('text', usetex=True)
plt.clf() # clear the plot (figure)
plt.plot(velocity, cubeSpec['data'], 'k-')
# the r prefix to the following format strings prevents \ from triggering an escape sequence, 
# like \n for line feed
plt.xlabel(r"$cz$ (km s$^{-1}$)",fontsize=16) 
plt.ylabel(r"$S_{\nu}$ (mJy)",fontsize=16)
plt.ylim(-0.023,0.005)
plt.savefig("mkn6-abs.eps") # save the figure as encapsulated postscript, suitable for journals
plt.savefig("mkn6-abs.png") # save the figure as a bitmap, better suited for astro-ph
--


Save the spectrum to a text file

It's handy to have a simple, columnar text file to inspect and analyze later. Or perhaps you might want to analyze the spectrum outside of CASA. Numpy makes it easy to save and restore text files.

# In CASA
# First, combine the spectrum into a single array
%cpaste

# Press Enter or Return, then copy/paste the following:
spectrum = numpy.vstack((velocity, cubeSpec['data']))
numpy.savetxt("mkn6.txt", numpy.transpose(spectrum))
--


Here's a listing of the first few lines of mkn6.txt.

5.367264765747175261e+03 -2.500247210264205933e-04
5.388621623358961187e+03 -2.288779069203883410e-04
5.409981470534922664e+03  4.886248498223721981e-04
5.431344307902921173e+03  3.546753214322961867e-05
5.452710136090820924e+03 -1.000610413029789925e-03
5.474078955726684399e+03 -2.162145450711250305e-03
5.495450767438905132e+03 -1.413836842402815819e-03
5.516825571855946691e+03 -9.065763442777097225e-04
5.538203369606402703e+03 -3.753502853214740753e-03
5.559584161319202394e+03 -1.479859836399555206e-02

Here's how to restore that text file back into CASA (or python, for that matter).

# In CASA
%cpaste

# Press Enter or Return, then copy/paste the following:
# import numpy # if you haven't already
spectrum = numpy.loadtxt("mkn6.txt",dtype="float")
x = spectrum[:,0]
y = spectrum[:,1]
--

Last checked on CASA Version 5.5.0.