First Look at Line Imaging: Difference between revisions

From CASA Guides
Jump to navigationJump to search
Tag: New redirect
 
(31 intermediate revisions by 5 users not shown)
Line 1: Line 1:
This lesson steps you through imaging of spectral lines.
#REDIRECT [[First_Look_at_Line_Imaging_CASA_6.5.4]]
 
You should have downloaded the data package as part of the first imaging tutorial.  If you haven't done that yet, check the [http://casaguides.nrao.edu/index.php?title=First_Look_at_Imaging First Look At Imaging] Guide for instructions.
 
We start here
with the self-calibrated data from the previous tutorial. We could
just as easily begin with the non self-calibrated data, but we would
not then benefit from the improvements that the self calibration
provided.
 
Note that it often makes good sense to self-calibrate on the continuum source, even when the main science goal is to image spectral line emission.  This is especially true when the continuum is bright and not overly complex.  Sometimes, though, it can make better sense to self-calibrate on a spectral line, and apply those solutions to the entire data set.  For now, we will apply the continuum self-cal solutions to our spectral line data to optimize our spectral line images.
 
First, copy the self-calibrated measurement set to the
current directory.  This measurement set contains both spectral line and continuum data for our science target.
 
<source lang="python">
# In CASA
os.system("rm -rf sis14_twhya_selfcal.ms")
os.system("cp -r ../working_data/sis14_twhya_selfcal.ms .")
</source>
 
== Continuum subtraction ==
 
Our first step is to remove the continuum emission from the UV data. Here, continuum emission means broadband emission that spans all spectral channels.  The continuum may have some mild frequency dependence, and that's OK. 
The cleanest way to do the continuum subtraction is to operate on the UV data, fitting the amplitudes using only those spectral channels that are free from spectral line emission. 
We will do so using the CASA task {{uvcontsub}}.
 
The main subtlety in the UV continuum subtraction is that we want
to be sure to fit the continuum using only channels that contain
mostly continuum (and not line) emission. To this end, the first
thing that we will do is plot the integrated spectrum of the science target
using the {{plotms}} task.  We specify the science source (field 5) and
plot the amplitude as a function channel number, averaging the data over
time.  (The parameter avgscan=True tells CASA to allow time averaging across scans.)
 
<figure id="Imaging-tutorial-line-plotms-1.png">
[[File:Imaging-tutorial-line-plotms-1.png|thumb|<caption>In this plot of amplitude vs. channel number, we can see several high amplitudes that represent detected spectral line emission.  Note that the amplitudes away from the spectral line are well above zero, showing the continuum emission.</caption>]]
</figure>
 
<source lang="python">
# In CASA
plotms(vis='sis14_twhya_selfcal.ms',
xaxis='channel',
yaxis='amp',
field='5',
avgspw=False,
avgtime='1e9',
avgscan=True,
avgbaseline=True,
showgui = True)
</source>
 
Inside the {{plotms}} tabs, you can toggle the axis to frequency and
then back to channel - the N2H+ line is at 372.67249 GHz and you can see that there
are several high amplitudes around channel 250 representing a weak but still visible spectral line.
The rest of the spectrum seems reasonably clear for measuring the continuum. To be sure that we exclude all spectral line emission,
let's avoid channels 240-280 in our continuum fit.
 
In our call to uvcontsub, we specify the line-free regions by setting
fitspw='0:240~280' and then set excludechans=True to tell CASA that
these are the channels NOT to fit.  In this case we'll fit only a constant amplitude (flat)
continuum by setting the fitorder to 0.  We target only field '5',
the science source. Note that we want to attempt this subtraction
for every individual UV data point and so we request a solution
interval equal to the integration time.
 
On a larger ALMA dataset, this step may take some time depending on your computing resources!
 
<source lang="python">
# In CASA
os.system('rm -rf sis14_twhya_selfcal.ms.contsub')
uvcontsub(vis = 'sis14_twhya_selfcal.ms',
field = '5',
fitspw = '0:240~280',
excludechans = True,
fitorder = 0,
solint='int')
</source>
 
In CASA 5.1, you will see a warning that uvcontsub is using the original Visibility Iterator, which is the underlying iteration and retrieval mechanism used when processing visibility data.
A new version of the Visibility Iterator has been developed and deployed in several other calibration tasks in CASA 5.0 and later (see [https://science.nrao.edu/enews/casa_005/ CASA eNews] for more information).
We can ignore this warning.
 
<pre style="background-color: #fffacd;">
WARN calibrater Forcing use of OLD VisibilityIterator.
</pre>
 
The output is a continuum-subtracted UV data set that has the
file extension ".contsub". We can now plot that using {{plotms}},
similar to what we did before, and look at the spectrum, which
should now contain only the spectral line. Note that because we only
ran uvcontsub on field '5' above, and created a new MS, the science data is now labeled as field '0' (the first
and only field) in the new measurement set.  (Run {{listobs}} on the continuum-subtracted data set to see how the field is now labeled, if you like.)
 
<figure id="Imaging-tutorial-line-plotms-2.png">
[[File:Imaging-tutorial-line-plotms-2.png|thumb|<caption>A plot of amplitude vs. channel number, after continuum subtraction.  Note that because the amplitude is positive definite, you don't expect the line-free channels to go entirely down to zero.  They should just
average to a very small value that reflects the noise.</caption>]]
</figure>
 
<source lang="python">
# In CASA
plotms(vis='sis14_twhya_selfcal.ms.contsub',
xaxis='channel',
yaxis='amp',
field='0',
avgspw=False,
avgtime='1e9',
avgscan=True,
avgbaseline=True,
showgui = True)
</source>
 
The continuum subtraction looks reasonably good.  We'll see better when we image the cube.
 
== Imaging the spectral line ==
 
Line imaging is like continuum imaging with (as you would expect)
the additional dimension of spectral information. In the call to {{tclean}}, we now need to set specmode="cube".
Further, we need to specify the third axis onto which the data will be gridded by
setting appropriate values (and units) for the start, width, and nchan parameters.
This third axis can be specified in units of channel, velocity, or frequency. Any of these choices would be valid for these data.
Here, we image 15 channels in the output data cube, each 0.5 km/s in width, starting at 0 km/s in the LSRK reference frame by setting start="0.0km/s", nchan=15, and width="0.5km/s".
The velocity is defined relative to the N2H+ rest frequency, which can be looked up in [http://www.cv.nrao.edu/php/splat/ splatalogue].
 
Set the N2H+ rest frequency in a python variable as follows:
 
<source lang="python">
# In CASA
restfreq = '372.67249GHz'
</source>
 
Next we will use {{tclean}} to make the cube, targeting the continuum-subtracted data.  Set field to '0'
and spw to '0' (the only options in this data set), and use specmode="cube"
with 15 channels of 0.5 km/s starting at 0.0 km/s defined in the
LSRK frame (a common frame for Galactic and extragalactic work).
We also set restoringbeam='common' to ensure a constant beam size for each plane in the cube.
Otherwise, this {{tclean}} will resemble our previous calls when we imaged the continuum.
 
<figure id="Imaging-tutorial-line-clean.png">
[[File:Imaging-tutorial-line-clean.png|thumb|<caption>Here we show the clean viewer, after using the tape deck control to step up to channel 5, which shows some line emission.</caption>]]
</figure>
 
<source lang="python">
# In CASA
os.system('rm -rf twhya_n2hp.*')
tclean(vis = 'sis14_twhya_selfcal.ms.contsub',
imagename = 'twhya_n2hp',
field = '0',
spw = '0',
specmode = 'cube',
nchan = 15,
start = '0.0km/s',
width = '0.5km/s',
outframe = 'LSRK',
restfreq = restfreq,
deconvolver= 'hogbom',
gridder = 'standard',
imsize = [250, 250],
cell = '0.08arcsec',
phasecenter = 0,
weighting = 'briggs',
robust = 0.5,
restoringbeam='common',
interactive = True,
niter=5000)
</source>
 
We specified interactive mode, so the clean GUI will come up. When it does, step through the
channels in the cube using the tape deck controls. We see emission over a few
channels near the center of the cube.  Notice the location of the emission
shifts from channel to channel -- Cool!  This is showing velocity structure from
the rotating disk!
 
Since the line emission has fairly low signal-to-noise, we don't expect cleaning to make a dramatic improvement here.  But, let's give it a shot to try to improve the image.  For spectral cubes, it's possible to specify a clean mask individually for each velocity channel, or you can specify a single mask that gets applied to all channels.  There is a button at the top of the {{tclean}} GUI that toggles between these options.  Both are viable approaches.  For a high fidelity telescope like ALMA it's often (but not always) fine to use a single mask across all channels.
 
The output of the {{tclean}} task is a data cube.  You can use the viewer to
inspect the cube, plot spectra, estimate the noise, and overlay this
result with the earlier continuum image.
 
<figure id="Imaging-tutorial-line-profile.png">
[[File:Imaging-tutorial-line-profile.png|thumb|<caption>Here we show a spectral line profile plotted using the viewer.</caption>]]
</figure>
 
<source lang="python">
# In CASA
imview("twhya_n2hp.image")
</source>
 
At this point you should play with viewer.  Change the velocity axis
definitions, measure statistics from the image planes, and so on. 
 
The viewer gives you the capability to make spectral profile plots.  On the viewer GUI, select "Tools -> Spectral Profile ...".  Then draw a box or ellipse around a region of interest, or set a single point using the "point marking" button to get a spectrum from a cut through a single pixel on the image.  The figure here shows a spectrum from a single-pixel cut through the cube, with the pixel set on a lobe of bright spectral line emission.
 
== Primary beam correction ==
 
As with continuum imaging, the line data cubes produced by CASA do not have the primary beam correction by default. Remember that to form an astronomically correct image of the sky, the output of tclean needs to be divided by this primary beam (or, in the case of mosaics, the combination of primary beam patterns used to make the mosaic). For final calculations you will therefore want to produce a primary beam-corrected version of your cube.
You can create the primary beam-corrected cube in {{tclean}} by setting the parameter pbcor = True. This will produce an additional file with the extension .image.pbcor, which is the cleaned cube corrected for the primary beam.
 
If you've already cleaned your data and don't want to rerun {{tclean}}, CASA stores the primary beam information needed to make this correction in the file with the .pb extension. The CASA task impbcor can be used to combine the .pb image with the output image from tclean to produce a primary-beam corrected image. First remove the old primary beam corrected image if it exists:
 
<source lang="python">
# In CASA
os.system('rm -rf twhya_n2hp.pbcor.image')
</source>
 
Now correct the image:
 
<source lang="python">
# In CASA
impbcor(imagename='twhya_n2hp.image',
pbimage='twhya_n2hp.pb',
outfile='twhya_n2hp.pbcor.image')
</source>

Latest revision as of 17:16, 12 October 2023