Difference between revisions of "M100 Band3 Combine 4.2.2"

From CASA Guides
Jump to navigationJump to search
(Created page with "M100 Combination of ACA+12m and Feathering")
 
Line 1: Line 1:
M100 Combination of ACA+12m and Feathering
+
== Overview ==
 +
 
 +
 
 +
== Split off CO spws ==
 +
 
 +
<source lang="python">
 +
# In CASA
 +
os.system('m100_*m.ms.listobs')
 +
listobs('M100_Band3_12m_CalibratedData.ms',listfile='m100_12m.ms.listobs')
 +
listobs('M100_Band3_7m_CalibratedData.ms',listfile='m100_7m.ms.listobs')
 +
 
 +
Investigation shows that the CO is spw=0 for the 12m data and in
 +
spws 1,3 for the 7m data. There are two for the 7m data because some
 +
of the data was taken with a slightly different correlator setup.
 +
 
 +
Also note that the integratin time per visibility is different:
 +
10.1s for the 7m and 6.05s for the 12m data as this will be
 +
important to correctly weighting the data.
 +
 
 +
<source lang="python">
 +
# In CASA
 +
os.system('rm -rf m100_12m_CO.ms')
 +
split(vis='M100_Band3_12m_CalibratedData.ms',
 +
      outputvis='m100_12m_CO.ms',spw='0',field='M100',
 +
      datacolumn='data',keepflags=False)
 +
os.system('rm -rf m100_7m_CO.ms')
 +
split(vis='M100_Band3_7m_CalibratedData.ms',
 +
      outputvis='m100_7m_CO.ms',spw='1,3',field='M100',
 +
      datacolumn='data',keepflags=False)
 +
</source>
 +
 
 +
Because the continuum is too weak to contribute to a 5km/s
 +
channel we will forgo the the continuum subtraction step.
 +
 
 +
== Concat with 1/sigma**2 scaling of Weights ==
 +
 
 +
When combining data with disparate properties it is very important
 +
that the relative weights of each visibility be in the correct
 +
proportion according to the radiometer equation. Formally the
 +
visibility weights should be proportinal to 1/sigma**2 where sigma
 +
is the variance or rms noise of a given visibility.
 +
 
 +
CASA currently scales the weights by 1/[(Tsys(i) * Tsys(j)] if
 +
calwt=True for # the Tsys table applycal.
 +
 
 +
# Plot the weights
 +
<source lang="python">
 +
# In CASA
 +
os.system('rm -rf 7m_WT.png 12m_WT.png')
 +
plotms(vis='m100_12m_CO.ms',yaxis='wt',xaxis='uvdist',spw='0~2:200',
 +
      coloraxis='spw',plotfile='12m_WT.png')
 +
 
 +
plotms(vis='m100_7m_CO.ms',yaxis='wt',xaxis='uvdist',spw='0~2:200',
 +
      coloraxis='spw',plotfile='7m_WT.png')
 +
</source>
 +
 
 +
Thus 7m WT = (7./12.)**4*(10.1/6.05) = 0.193 to
 +
account for the difference in telescope size and integration time
 +
per visibility
 +
 
 +
<source lang="python">
 +
# In CASA
 +
os.system('rm -rf M100_Intcombo_0.193.ms')
 +
concat(vis=['m100_12m_CO.ms','m100_7m_CO.ms'],
 +
      concatvis='M100_Intcombo_0.193.ms',
 +
      visweightscale=[1,0.193])
 +
</source>
 +
 
 +
<source lang="python">
 +
# In CASA
 +
os.system('rm -rf Intcombo_0.193_WT.png')
 +
plotms(vis='M100_Intcombo_0.193.ms',yaxis='wt',xaxis='uvdist',spw='0~2:200',
 +
      coloraxis='spw',plotfile='Intcombo_0.193_WT.png')
 +
</source>
 +
 
 +
 
 +
<source lang="python">
 +
# In CASA
 +
os.system('rm -rfM100_Intcombo_0.193 _uvdist.png')
 +
plotms(vis='M100_Intcombo_0.193.ms',yaxis='amp',xaxis='uvdist',spw='',
 +
      avgchannel='5000',
 +
      coloraxis='spw',plotfile='M100_Intcombo_0.193_uvdist.png')
 +
 
 +
os.system('rm -rfM100_Intcombo_0.193 _vel.png')
 +
plotms(vis='M100_Intcombo_0.193.ms',yaxis='amp',xaxis='velocity',spw='',
 +
      avgtime='1e8',coloraxis='spw',
 +
      transform=True,freqframe='LSRK',restfreq='115.271201800GHz',
 +
      plotfile='M100_Intcombo_0.193_vel.png')
 +
 
 +
os.system('rm -rfM100_Intcombo_0.193 _chan.png')
 +
plotms(vis='M100_Intcombo_0.193.ms',yaxis='amp',xaxis='channel',spw='',
 +
      avgtime='1e8',coloraxis='spw',plotfile='M100_Intcombo_0.193_chan.png')
 +
</source>
 +
 
 +
 +
== Image Using Automasking Technique ==
 +
 
 +
This script attempts to do an iterative automasking proceedure down
 +
to a user specified stop*rms where stop is typcially 2-3.  It takes
 +
some care to ensure that the generated masks (i) only have values of
 +
0 or 1; (ii) are themselves masked at the minpb level. It also
 +
removes very small masked regions that are consistent with noise
 +
bumps using a function in scipy. A typical seeting is to remove mask
 +
regions that are 1/2 the beam area in pixels.
 +
 
 +
=== Define Parameters ===
 +
<source lang="python">
 +
# In CASA
 +
 
 +
# Initialize
 +
import scipy.ndimage
 +
 
 +
# Define clean parameters
 +
vis='M100_Intcombo_0.193.ms'
 +
prename='M100_Intcombo_0.193_cube'
 +
myimage=prename+'.image'
 +
myflux=prename+'.flux'
 +
mymask=prename+'.mask'
 +
myresidual=prename+'.residual'
 +
imsize=800
 +
cell='0.5arcsec'
 +
minpb=0.2
 +
restfreq='115.271201800GHz'
 +
outframe='LSRK'
 +
spw='0~2'
 +
width='5km/s'
 +
start='1400km/s'
 +
nchan=70
 +
robust=0.5
 +
phasecenter='J2000 12h22m54.9 +15d49m10'
 +
scales=[0]
 +
smallscalebias=0.6
 +
 
 +
# Setup stopping criteria
 +
stop=3. # multiplier for rms
 +
 
 +
# Minimum size multiplier for beam area for removing very
 +
# small mask regions.
 +
pixelmin=0.5  # reasonable default is 1/2 the beam area
 +
</source>
 +
 
 +
=== Make Initial Dirty Image and find rms noise ===
 +
 
 +
<source lang="python">
 +
# In CASA
 +
# Make initial dirty image
 +
os.system('rm -rf '+prename+'.* ' +prename+'_*')
 +
clean(vis=vis,imagename=prename,
 +
      imagermode='mosaic',ftmachine='mosaic',minpb=minpb,
 +
      imsize=imsize,cell=cell,spw=spw,
 +
      weighting='briggs',robust=robust,phasecenter=phasecenter,
 +
      mode='velocity',width=width,start=start,nchan=nchan,     
 +
      restfreq=restfreq,outframe=outframe,veltype='radio',
 +
      mask='',
 +
      niter=0,interactive=F)
 +
</source>
 +
 
 +
<source lang="python">
 +
# In CASA
 +
# Find properties of the dirty image
 +
myimage=prename+'.image'
 +
bigstat=imstat(imagename=myimage)
 +
peak= bigstat['max'][0]
 +
print 'peak in cube = '+str(peak)
 +
# sets threshold of first loop, try 2-4. Subsequent
 +
# loops are set thresh/2
 +
thresh = peak /4.
 +
 
 +
 
 +
if True:
 +
    # If True: find the rms in two line-free channels   
 +
    chanstat=imstat(imagename=myimage,chans='4')
 +
    rms1= chanstat['rms'][0]
 +
    chanstat=imstat(imagename=myimage,chans='66')
 +
    rms2= chanstat['rms'][0]
 +
    rms=0.5*(rms1+rms2)       
 +
else:
 +
    # Set rms by hand
 +
    rms=0.013
 +
 
 +
 
 +
print 'rms in a channel = '+str(rms)
 +
</source>
 +
 
 +
=== Define minimum size of masked region  ===
 +
 
 +
<source lang="python">
 +
# In CASA
 +
# Deterimine the beam area in pixels for later removal of very small
 +
# mask regions
 +
major=imhead(imagename=myimage,mode='get',hdkey='beammajor')['value']
 +
minor=imhead(imagename=myimage,mode='get',hdkey='beamminor')['value']
 +
pixelsize=float(cell.split('arcsec')[0])
 +
beamarea=(major*minor*pi/(4*log(2)))/(pixelsize**2)
 +
print 'beamarea in pixels =', beamarea
 +
</source>
 +
 
 +
=== Automasking Loop ===
 +
 
 +
<source lang="python">
 +
# In CASA
 +
n=-1
 +
while (thresh >= stop*rms): 
 +
    n=n+1
 +
    print 'clean threshold this loop is', thresh
 +
    threshmask = prename+'_threshmask' +str(n)
 +
    maskim = prename+'_fullmask' +str(n)
 +
    immath(imagename = [myresidual],
 +
          outfile = threshmask,
 +
          expr = 'iif(IM0 > '+str(thresh) +',1.0,0.0)',
 +
          mask=myflux+'>'+str(minpb))
 +
    if (n==0):
 +
        os.system('cp -r '+threshmask+' '+maskim+'.pb')
 +
        print 'This is the first loop'
 +
    else:
 +
        makemask(mode='copy',inpimage=myimage,
 +
                inpmask=[threshmask,mymask],
 +
                output=maskim)
 +
        imsubimage(imagename=maskim, mask=myflux+'>'+str(minpb),
 +
                  outfile=maskim+'.pb')   
 +
    print 'Combined mask ' +maskim+' generated.'
 +
 
 +
    # Remove small masks
 +
    os.system('cp -r '+maskim+'.pb ' +maskim+'.pb.min')
 +
    maskfile=maskim+'.pb.min'
 +
    ia.open(maskfile)
 +
    mask=ia.getchunk()         
 +
    labeled,j=scipy.ndimage.label(mask)                   
 +
    myhistogram = scipy.ndimage.measurements.histogram(labeled,0,j+1,j+1)
 +
    object_slices = scipy.ndimage.find_objects(labeled)
 +
    threshold=beamarea*pixelmin
 +
    for i in range(j):
 +
        if myhistogram[i+1]<threshold:
 +
            mask[object_slices[i]] = 0
 +
 
 +
 
 +
    ia.putchunk(mask)
 +
    ia.done()
 +
    print 'Small masks removed and ' +maskim +'.pb.min generated.'
 +
 
 +
    os.system('rm -rf '+mymask+'')
 +
    clean(vis=vis,imagename=prename,
 +
          imagermode='mosaic',ftmachine='mosaic',minpb=minpb,
 +
          imsize=imsize,cell=cell,spw=spw,
 +
          weighting='briggs',robust=robust,phasecenter=phasecenter,
 +
          mode='velocity',width=width,start=start,nchan=nchan,     
 +
          restfreq=restfreq,outframe=outframe,veltype='radio',
 +
          mask = maskim+'.pb.min',
 +
          multiscale=scales,smallscalebias=smallscalebias,
 +
          interactive = F,
 +
          niter = 10000,
 +
          threshold = str(thresh) +'Jy/beam')
 +
 
 +
 
 +
    if thresh==stop*rms: break
 +
    thresh = thresh/2.
 +
    # Run a final time with stop*rms if more than a little above
 +
    # stop*rms. Also make a back-up of next to last image
 +
    if (thresh < stop*rms and thresh*2.>1.05*stop*rms):
 +
        thresh=stop*rms 
 +
        os.system('cp -r '+myimage+' '+myimage+str(n))
 +
</source>
 +
 
 +
== Moment Maps for 7m+12m ==
 +
 
 +
<source lang="python">
 +
# In CASA
 +
os.system('rm -rf M100_Intcombo_0.193_cube.image.mom0')
 +
immoments(imagename = 'M100_Intcombo_0.193_cube.image',
 +
        moments = [0],
 +
        axis = 'spectral',chans = '10~61',
 +
        mask='M100_Intcombo_0.193_cube.flux>0.3',
 +
        outfile = 'M100_Intcombo_0.193_cube.image.mom0')
 +
</source>
 +
 
 +
<source lang="python">
 +
# In CASA
 +
myimage='M100_Intcombo_0.193_cube.image'
 +
chanstat=imstat(imagename=myimage,chans='4')
 +
rms1= chanstat['rms'][0]
 +
chanstat=imstat(imagename=myimage,chans='66')
 +
rms2= chanstat['rms'][0]
 +
rms=0.5*(rms1+rms2)
 +
print 'rms in a channel = '+str(rms)   
 +
os.system('rm -rf M100_Intcombo_0.193_cube.image.mom1')
 +
immoments(imagename = 'M100_Intcombo_0.193_cube.image',
 +
        moments = [1],
 +
        axis = 'spectral',chans = '10~61',
 +
        mask='M100_Intcombo_0.193_cube.flux>0.3',
 +
        includepix = [rms*6,100.],
 +
        outfile = 'M100_Intcombo_0.193_cube.image.mom1')
 +
</source>
 +
 
 +
 
 +
<source lang="python">
 +
# In CASA
 +
os.system('rm -rf M100_Intcombo_0.193_cube.flux.1ch')
 +
imsubimage(imagename='M100_Intcombo_0.193_cube.flux',
 +
          outfile='M100_Intcombo_0.193_cube.flux.1ch',
 +
          chans='35')
 +
</source>
 +
 
 +
<source lang="python">
 +
# In CASA
 +
os.system('rm -rf M100_Intcombo_0.193_cube.image.mom*.png')
 +
imview (raster=[{'file': 'M100_Intcombo_0.193_cube.image.mom0',
 +
                'range': [-1,25.],'scaling': -0.5,'colorwedge': T}],
 +
        zoom={'blc': [150,150],'trc': [650,650]},
 +
        out='M100_Intcombo_0.193_cube.image.mom0.png')
 +
imview (raster=[{'file': 'M100_Intcombo_0.193_cube.image.mom1',
 +
                'range': [1455,1695],'colorwedge': T}],
 +
        zoom={'blc': [150,150],'trc': [650,650]},
 +
        out='M100_Intcombo_0.193_cube.image.mom1.png')
 +
</source>
 +
 
 +
<source lang="python">
 +
# In CASA
 +
os.system('rm -rf M100_Intcombo_0.193_cube.image.mom0.pbcor')
 +
immath(imagename=['M100_Intcombo_0.193_cube.image.mom0', \
 +
                      'M100_Intcombo_0.193_cube.flux.1ch'],
 +
        expr='IM0/IM1',
 +
        outfile='M100_Intcombo_0.193_cube.image.mom0.pbcor')
 +
</source>
 +
 
 +
 
 +
<source lang="python">
 +
# In CASA
 +
imview (raster=[{'file': 'M100_Intcombo_0.193_cube.image.mom0',
 +
                'range': [-1,25.],'scaling': -0.5},
 +
                {'file': 'M100_Intcombo_0.193_cube.image.mom0.pbcor',
 +
                'range': [-1,25.],'scaling': -0.5}],
 +
        zoom={'blc': [150,150],'trc': [650,650]})
 +
</source>
 +
 
 +
 
 +
== Feathering TP Image with Interferrometric Image==

Revision as of 14:04, 17 June 2013

Overview

Split off CO spws

# In CASA
os.system('m100_*m.ms.listobs')
listobs('M100_Band3_12m_CalibratedData.ms',listfile='m100_12m.ms.listobs')
listobs('M100_Band3_7m_CalibratedData.ms',listfile='m100_7m.ms.listobs')

Investigation shows that the CO is spw=0 for the 12m data and in
spws 1,3 for the 7m data. There are two for the 7m data because some
of the data was taken with a slightly different correlator setup.

Also note that the integratin time per visibility is different:
10.1s for the 7m and 6.05s for the 12m data as this will be
important to correctly weighting the data.

<source lang="python">
# In CASA
os.system('rm -rf m100_12m_CO.ms')
split(vis='M100_Band3_12m_CalibratedData.ms',
      outputvis='m100_12m_CO.ms',spw='0',field='M100',
      datacolumn='data',keepflags=False)
os.system('rm -rf m100_7m_CO.ms')
split(vis='M100_Band3_7m_CalibratedData.ms',
      outputvis='m100_7m_CO.ms',spw='1,3',field='M100',
      datacolumn='data',keepflags=False)

Because the continuum is too weak to contribute to a 5km/s channel we will forgo the the continuum subtraction step.

Concat with 1/sigma**2 scaling of Weights

When combining data with disparate properties it is very important that the relative weights of each visibility be in the correct proportion according to the radiometer equation. Formally the visibility weights should be proportinal to 1/sigma**2 where sigma is the variance or rms noise of a given visibility.

CASA currently scales the weights by 1/[(Tsys(i) * Tsys(j)] if calwt=True for # the Tsys table applycal.

  1. Plot the weights
# In CASA
os.system('rm -rf 7m_WT.png 12m_WT.png')
plotms(vis='m100_12m_CO.ms',yaxis='wt',xaxis='uvdist',spw='0~2:200',
       coloraxis='spw',plotfile='12m_WT.png')

plotms(vis='m100_7m_CO.ms',yaxis='wt',xaxis='uvdist',spw='0~2:200',
       coloraxis='spw',plotfile='7m_WT.png')

Thus 7m WT = (7./12.)**4*(10.1/6.05) = 0.193 to account for the difference in telescope size and integration time per visibility

# In CASA
os.system('rm -rf M100_Intcombo_0.193.ms')
concat(vis=['m100_12m_CO.ms','m100_7m_CO.ms'],
       concatvis='M100_Intcombo_0.193.ms',
       visweightscale=[1,0.193])
# In CASA
os.system('rm -rf Intcombo_0.193_WT.png')
plotms(vis='M100_Intcombo_0.193.ms',yaxis='wt',xaxis='uvdist',spw='0~2:200',
       coloraxis='spw',plotfile='Intcombo_0.193_WT.png')


# In CASA
os.system('rm -rfM100_Intcombo_0.193 _uvdist.png')
plotms(vis='M100_Intcombo_0.193.ms',yaxis='amp',xaxis='uvdist',spw='',
       avgchannel='5000',
       coloraxis='spw',plotfile='M100_Intcombo_0.193_uvdist.png') 

os.system('rm -rfM100_Intcombo_0.193 _vel.png')
plotms(vis='M100_Intcombo_0.193.ms',yaxis='amp',xaxis='velocity',spw='',
       avgtime='1e8',coloraxis='spw',
       transform=True,freqframe='LSRK',restfreq='115.271201800GHz',
       plotfile='M100_Intcombo_0.193_vel.png')

os.system('rm -rfM100_Intcombo_0.193 _chan.png')
plotms(vis='M100_Intcombo_0.193.ms',yaxis='amp',xaxis='channel',spw='',
       avgtime='1e8',coloraxis='spw',plotfile='M100_Intcombo_0.193_chan.png')


Image Using Automasking Technique

This script attempts to do an iterative automasking proceedure down to a user specified stop*rms where stop is typcially 2-3. It takes some care to ensure that the generated masks (i) only have values of 0 or 1; (ii) are themselves masked at the minpb level. It also removes very small masked regions that are consistent with noise bumps using a function in scipy. A typical seeting is to remove mask regions that are 1/2 the beam area in pixels.

Define Parameters

# In CASA

# Initialize 
import scipy.ndimage 

# Define clean parameters
vis='M100_Intcombo_0.193.ms'
prename='M100_Intcombo_0.193_cube'
myimage=prename+'.image'
myflux=prename+'.flux'
mymask=prename+'.mask'
myresidual=prename+'.residual'
imsize=800
cell='0.5arcsec'
minpb=0.2
restfreq='115.271201800GHz'
outframe='LSRK'
spw='0~2'
width='5km/s'
start='1400km/s'
nchan=70
robust=0.5
phasecenter='J2000 12h22m54.9 +15d49m10'
scales=[0]
smallscalebias=0.6

# Setup stopping criteria 
stop=3. # multiplier for rms 

# Minimum size multiplier for beam area for removing very 
# small mask regions. 
pixelmin=0.5  # reasonable default is 1/2 the beam area

Make Initial Dirty Image and find rms noise

# In CASA
# Make initial dirty image
os.system('rm -rf '+prename+'.* ' +prename+'_*')
clean(vis=vis,imagename=prename,
      imagermode='mosaic',ftmachine='mosaic',minpb=minpb,
      imsize=imsize,cell=cell,spw=spw,
      weighting='briggs',robust=robust,phasecenter=phasecenter,
      mode='velocity',width=width,start=start,nchan=nchan,      
      restfreq=restfreq,outframe=outframe,veltype='radio',
      mask='',
      niter=0,interactive=F)
# In CASA
# Find properties of the dirty image
myimage=prename+'.image'
bigstat=imstat(imagename=myimage)
peak= bigstat['max'][0]
print 'peak in cube = '+str(peak)
# sets threshold of first loop, try 2-4. Subsequent
# loops are set thresh/2
thresh = peak /4.


if True:
    # If True: find the rms in two line-free channels    
    chanstat=imstat(imagename=myimage,chans='4')
    rms1= chanstat['rms'][0]
    chanstat=imstat(imagename=myimage,chans='66')
    rms2= chanstat['rms'][0]
    rms=0.5*(rms1+rms2)        
else:
    # Set rms by hand
    rms=0.013


print 'rms in a channel = '+str(rms)

Define minimum size of masked region

# In CASA
# Deterimine the beam area in pixels for later removal of very small
# mask regions
major=imhead(imagename=myimage,mode='get',hdkey='beammajor')['value']
minor=imhead(imagename=myimage,mode='get',hdkey='beamminor')['value']
pixelsize=float(cell.split('arcsec')[0])
beamarea=(major*minor*pi/(4*log(2)))/(pixelsize**2)
print 'beamarea in pixels =', beamarea

Automasking Loop

# In CASA
n=-1
while (thresh >= stop*rms):   
    n=n+1
    print 'clean threshold this loop is', thresh
    threshmask = prename+'_threshmask' +str(n)
    maskim = prename+'_fullmask' +str(n)
    immath(imagename = [myresidual],
           outfile = threshmask,
           expr = 'iif(IM0 > '+str(thresh) +',1.0,0.0)',
           mask=myflux+'>'+str(minpb))
    if (n==0):
        os.system('cp -r '+threshmask+' '+maskim+'.pb')
        print 'This is the first loop'
    else:
        makemask(mode='copy',inpimage=myimage,
                 inpmask=[threshmask,mymask],
                 output=maskim)
        imsubimage(imagename=maskim, mask=myflux+'>'+str(minpb),
                   outfile=maskim+'.pb')     
    print 'Combined mask ' +maskim+' generated.'

    # Remove small masks
    os.system('cp -r '+maskim+'.pb ' +maskim+'.pb.min')
    maskfile=maskim+'.pb.min'
    ia.open(maskfile)
    mask=ia.getchunk()           
    labeled,j=scipy.ndimage.label(mask)                     
    myhistogram = scipy.ndimage.measurements.histogram(labeled,0,j+1,j+1)
    object_slices = scipy.ndimage.find_objects(labeled)
    threshold=beamarea*pixelmin
    for i in range(j):
        if myhistogram[i+1]<threshold:
            mask[object_slices[i]] = 0


    ia.putchunk(mask)
    ia.done()
    print 'Small masks removed and ' +maskim +'.pb.min generated.'

    os.system('rm -rf '+mymask+'')
    clean(vis=vis,imagename=prename,
          imagermode='mosaic',ftmachine='mosaic',minpb=minpb,
          imsize=imsize,cell=cell,spw=spw,
          weighting='briggs',robust=robust,phasecenter=phasecenter,
          mode='velocity',width=width,start=start,nchan=nchan,      
          restfreq=restfreq,outframe=outframe,veltype='radio',
          mask = maskim+'.pb.min',
          multiscale=scales,smallscalebias=smallscalebias,
          interactive = F,
          niter = 10000,
          threshold = str(thresh) +'Jy/beam')


    if thresh==stop*rms: break
    thresh = thresh/2.
    # Run a final time with stop*rms if more than a little above
    # stop*rms. Also make a back-up of next to last image
    if (thresh < stop*rms and thresh*2.>1.05*stop*rms):
        thresh=stop*rms  
        os.system('cp -r '+myimage+' '+myimage+str(n))

Moment Maps for 7m+12m

# In CASA
os.system('rm -rf M100_Intcombo_0.193_cube.image.mom0')
immoments(imagename = 'M100_Intcombo_0.193_cube.image',
         moments = [0],
         axis = 'spectral',chans = '10~61',
         mask='M100_Intcombo_0.193_cube.flux>0.3',
         outfile = 'M100_Intcombo_0.193_cube.image.mom0')
# In CASA
myimage='M100_Intcombo_0.193_cube.image'
chanstat=imstat(imagename=myimage,chans='4')
rms1= chanstat['rms'][0]
chanstat=imstat(imagename=myimage,chans='66')
rms2= chanstat['rms'][0]
rms=0.5*(rms1+rms2)
print 'rms in a channel = '+str(rms)    
os.system('rm -rf M100_Intcombo_0.193_cube.image.mom1')
immoments(imagename = 'M100_Intcombo_0.193_cube.image',
         moments = [1],
         axis = 'spectral',chans = '10~61',
         mask='M100_Intcombo_0.193_cube.flux>0.3',
         includepix = [rms*6,100.],
         outfile = 'M100_Intcombo_0.193_cube.image.mom1')


# In CASA
os.system('rm -rf M100_Intcombo_0.193_cube.flux.1ch')
imsubimage(imagename='M100_Intcombo_0.193_cube.flux',
           outfile='M100_Intcombo_0.193_cube.flux.1ch',
           chans='35')
# In CASA
os.system('rm -rf M100_Intcombo_0.193_cube.image.mom*.png')
imview (raster=[{'file': 'M100_Intcombo_0.193_cube.image.mom0',
                 'range': [-1,25.],'scaling': -0.5,'colorwedge': T}],
         zoom={'blc': [150,150],'trc': [650,650]},
         out='M100_Intcombo_0.193_cube.image.mom0.png')
imview (raster=[{'file': 'M100_Intcombo_0.193_cube.image.mom1',
                 'range': [1455,1695],'colorwedge': T}],
         zoom={'blc': [150,150],'trc': [650,650]},
         out='M100_Intcombo_0.193_cube.image.mom1.png')
# In CASA
os.system('rm -rf M100_Intcombo_0.193_cube.image.mom0.pbcor')
immath(imagename=['M100_Intcombo_0.193_cube.image.mom0', \
                       'M100_Intcombo_0.193_cube.flux.1ch'],
        expr='IM0/IM1',
        outfile='M100_Intcombo_0.193_cube.image.mom0.pbcor')


# In CASA
imview (raster=[{'file': 'M100_Intcombo_0.193_cube.image.mom0',
                 'range': [-1,25.],'scaling': -0.5},
                {'file': 'M100_Intcombo_0.193_cube.image.mom0.pbcor',
                 'range': [-1,25.],'scaling': -0.5}],
         zoom={'blc': [150,150],'trc': [650,650]})


Feathering TP Image with Interferrometric Image