|
|
Line 1: |
Line 1: |
| = Overview = | | = Overview = |
|
| |
|
| |
| = Combine and Image the 7m+12m =
| |
|
| |
|
| |
| == 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')
| |
| </source>
| |
|
| |
| 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 =
| |