Stacking Multiple Spectral Lines at Same Position: Difference between revisions

From CASA Guides
Jump to navigationJump to search
(Created page with "{{Using CASA Hints Tips Tricks}} Category: CASA Analysis This guide presents a method in CASA for averaging together multiple transitions of the same line. This process incr...")
 
No edit summary
 
(6 intermediate revisions by the same user not shown)
Line 2: Line 2:
[[Category: CASA Analysis]]
[[Category: CASA Analysis]]


This guide presents a method in CASA for averaging together multiple transitions of the same line. This process increases the signal to noise of the line observations. It may not be suitable for all science cases though.
This guide presents a method in CASA for averaging together multiple transitions of the same line. This process increases the signal to noise of the line observations. It may not be suitable for all science cases.


The observations in this guide were of 8 hydrogen radio recombination lines in the center of galaxy. The lines are all emitted from the same region. Therefore, the strategy is to use {{clean}} to put each spectral line cube on the same velocity scale. Then the cubes are smoothed to the same resolution as the cube with the coarsest resolution using {{imsmooth}}. Finally, the spectral line cubes are averaged together using {{immath}}.
The observations in this guide were of 5 different hydrogen radio recombination lines. The correlator was configured to observe one spectral line per spectral window. The data were calibrated and split off so that each spectral line had its own visibility set.
 
The lines are all emitted from the same region in the galaxy. Therefore, a strategy to combine them is to use {{clean}} to put each spectral line cube on the same velocity scale. Then the cubes are smoothed to the same resolution as the cube with the coarsest resolution using {{imsmooth}}. Finally, the spectral line cubes are averaged together using {{immath}}. It is not necessary to regrid the visibilities using {{cvel}} prior to imaging because {{clean}} will regrid them regardless.


This script assumes that you have already continuum subtracted the data using {{uvcontsub}}. If you need to continuum subtract in the image domain using {{imcontsub}}, it should be done after the clean, not after the spectral line cubes are averaged together. Be careful with your continuum subtraction -- small errors in continuum subtraction can lead to spurious lines.
This script assumes that you have already continuum subtracted the data using {{uvcontsub}}. If you need to continuum subtract in the image domain using {{imcontsub}}, it should be done after the clean, not after the spectral line cubes are averaged together. Be careful with your continuum subtraction -- small errors in continuum subtraction can lead to spurious lines.


Instead of smoothing after cleaning, another option would be to have all the images use the same restoring beam size. You can set this parameter in {{clean}}.
Below we do a detailed walkthrough of the script.
First, we define which spectral line visibilities we are going to use and their rest frequencies.
<source lang='python'>
linems = ['ic342_cband_line_if2_cal_spw0.ms.contsub',
          'ic342_cband_line_if2_cal_spw1.ms.contsub',
          'ic342_cband_line_if2_cal_spw2.ms.contsub',
          'ic342_cband_line_if2_cal_spw3.ms.contsub',
          'ic342_cband_line_if2_cal_spw4.ms.contsub']
restfreq = ['6.2891GHz',
            '6.4788GHz',
            '6.6761GHz',
            '6.8815GHz',
            '7.0954GHz']
</source>
Now we define the parameters that we are going to use to image the spectral line data. Depending on your data set, you may want to change some of these parameters.
<source lang='python'>
velstart = '-250km/s'    # velocity to start the cube at
chanwidth = '10km/s'      # width of channels in the cube
nchan = 50                # number of channels in the cube
outframe = 'bary'        # use a barycentric reference frame
veltype = 'optical'      # use the optical velocity definition
imsize = [1536,1536]
cell = ['0.75arcsec']
</source>
Now set the parameters for a shallow clean on the data.
<source lang='python'>
threshold = '20mJy'
niter = 300
</source>
Finally, we set the name of the output image
<source lang='python'>
output_stacked_image = 'RRL_if2_stack.image'
</source>
Now we run the script.
<source lang='python'>
# creating the output image names
output_line_image = [linems[i].replace('.ms.contsub','.contsub') for i in range(len(linems))]
# creating the images for each spectral line data set
for i in range(len(linems)):
   
    clean(vis=linems[i],
          imagename=output_line_image[i],
          mode='velocity',
          nchan=nchan,
          start=velstart,
          width=chanwidth,
          outframe=outframe,
          veltype=veltype,
          niter=niter,
          threshold=threshold,
          interactive=False,
          imsize=imsize,
          cell=cell,
        restfreq=restfreq[i])
# Getting the sizes of the beam.
beammajor = []
beamminor = []
beampa = []
for image in output_line_image :
    myimagename = image + '.image'
    beammajor.append(imhead(imagename=myimagename,mode='get',hdkey='beammajor'))
    beamminor.append(imhead(imagename=myimagename,mode='get',hdkey='beamminor'))
    beampa.append(imhead(imagename=myimagename,mode='get',hdkey='beampa'))
# determining what size the spectral line images should be smoothed to.
targetbeammajor = max(beammajor)
targetidx = beammajor.index(targetbeammajor)
targetbeamminor = beamminor[targetidx]
targetbeampa = beampa[targetidx]
# smoothing the images
for image in output_line_image:
    inputimage = image + '.image'
    outputimage = image + '.image' + '.smooth'
    if image != output_line_image[targetidx]:
        imsmooth(imagename=inputimage,outfile=outputimage,kernel='gauss',
                major=str(targetbeammajor['value'])+targetbeammajor['unit'],
                minor=str(targetbeamminor['value'])+targetbeamminor['unit'],
                pa=str(targetbeampa['value'])+targetbeampa['unit'],
                targetres=True)
    else:
        # need to skip the one with the largest resolution or imsmooth fails.
        os.system('cp -ir '+ inputimage + ' ' + outputimage)
# creating the expression for immath. It doesn't appear to know about mean.
output_line_image = [output_line_image[i] + '.image' for i in range(len(output_line_image))]
myim = ['IM' + str(i) for i in range(len(output_line_image))]
myexp =  '(' + '+'.join(myim) + ')/' + str(float(len(output_line_image)))
# combining the images
immath(imagename=output_line_image,outfile=output_stacked_image,mode='evalexpr',expr=myexp)
</source>
--[[User:Akepley|Amanda Kepley]] 19:03, 23 April 2012 (PDT)


[[Main Page | &#8629; '''CASAguides''']] <br>
[[Main Page | &#8629; '''CASAguides''']] <br>

Latest revision as of 22:04, 23 April 2012

Hints, Tips, and Tricks

This guide presents a method in CASA for averaging together multiple transitions of the same line. This process increases the signal to noise of the line observations. It may not be suitable for all science cases.

The observations in this guide were of 5 different hydrogen radio recombination lines. The correlator was configured to observe one spectral line per spectral window. The data were calibrated and split off so that each spectral line had its own visibility set.

The lines are all emitted from the same region in the galaxy. Therefore, a strategy to combine them is to use clean to put each spectral line cube on the same velocity scale. Then the cubes are smoothed to the same resolution as the cube with the coarsest resolution using imsmooth. Finally, the spectral line cubes are averaged together using immath. It is not necessary to regrid the visibilities using cvel prior to imaging because clean will regrid them regardless.

This script assumes that you have already continuum subtracted the data using uvcontsub. If you need to continuum subtract in the image domain using imcontsub, it should be done after the clean, not after the spectral line cubes are averaged together. Be careful with your continuum subtraction -- small errors in continuum subtraction can lead to spurious lines.

Instead of smoothing after cleaning, another option would be to have all the images use the same restoring beam size. You can set this parameter in clean.

Below we do a detailed walkthrough of the script.

First, we define which spectral line visibilities we are going to use and their rest frequencies.

linems = ['ic342_cband_line_if2_cal_spw0.ms.contsub',
          'ic342_cband_line_if2_cal_spw1.ms.contsub',
          'ic342_cband_line_if2_cal_spw2.ms.contsub',
          'ic342_cband_line_if2_cal_spw3.ms.contsub',
          'ic342_cband_line_if2_cal_spw4.ms.contsub']
restfreq = ['6.2891GHz',
            '6.4788GHz',
            '6.6761GHz',
            '6.8815GHz',
            '7.0954GHz']

Now we define the parameters that we are going to use to image the spectral line data. Depending on your data set, you may want to change some of these parameters.

velstart = '-250km/s'     # velocity to start the cube at
chanwidth = '10km/s'      # width of channels in the cube
nchan = 50                # number of channels in the cube
outframe = 'bary'         # use a barycentric reference frame
veltype = 'optical'       # use the optical velocity definition
imsize = [1536,1536]
cell = ['0.75arcsec']

Now set the parameters for a shallow clean on the data.

threshold = '20mJy'
niter = 300

Finally, we set the name of the output image

output_stacked_image = 'RRL_if2_stack.image'

Now we run the script.

# creating the output image names
output_line_image = [linems[i].replace('.ms.contsub','.contsub') for i in range(len(linems))]

# creating the images for each spectral line data set
for i in range(len(linems)):
    
    clean(vis=linems[i],
          imagename=output_line_image[i],
          mode='velocity',
          nchan=nchan,
          start=velstart,
          width=chanwidth,
          outframe=outframe,
          veltype=veltype,
          niter=niter,
          threshold=threshold,
          interactive=False,
          imsize=imsize,
          cell=cell,
         restfreq=restfreq[i])

# Getting the sizes of the beam.
beammajor = []
beamminor = []
beampa = []

for image in output_line_image :
    myimagename = image + '.image'
    beammajor.append(imhead(imagename=myimagename,mode='get',hdkey='beammajor'))
    beamminor.append(imhead(imagename=myimagename,mode='get',hdkey='beamminor'))
    beampa.append(imhead(imagename=myimagename,mode='get',hdkey='beampa'))

# determining what size the spectral line images should be smoothed to.
targetbeammajor = max(beammajor)
targetidx = beammajor.index(targetbeammajor)
targetbeamminor = beamminor[targetidx]
targetbeampa = beampa[targetidx]

# smoothing the images
for image in output_line_image:
    inputimage = image + '.image'
    outputimage = image + '.image' + '.smooth'
    if image != output_line_image[targetidx]:
        imsmooth(imagename=inputimage,outfile=outputimage,kernel='gauss',
                 major=str(targetbeammajor['value'])+targetbeammajor['unit'],
                 minor=str(targetbeamminor['value'])+targetbeamminor['unit'],
                 pa=str(targetbeampa['value'])+targetbeampa['unit'],
                 targetres=True)
    else:
        # need to skip the one with the largest resolution or imsmooth fails.
        os.system('cp -ir '+ inputimage + ' ' + outputimage)


# creating the expression for immath. It doesn't appear to know about mean.
output_line_image = [output_line_image[i] + '.image' for i in range(len(output_line_image))]
myim = ['IM' + str(i) for i in range(len(output_line_image))]
myexp =  '(' + '+'.join(myim) + ')/' + str(float(len(output_line_image)))

# combining the images
immath(imagename=output_line_image,outfile=output_stacked_image,mode='evalexpr',expr=myexp)


--Amanda Kepley 19:03, 23 April 2012 (PDT)

CASAguides