casa天文软件全代码记录_casa软件-优快云博客

本文链接：https://blog.youkuaiyun.com/weixin_42842694/article/details/141866260

最近的工作使用casa的频率不少，但是我经常写了代码或者弄好的脚本想要复用找起来很麻烦。我想记录一下成图所有使用过的东西，到时候我复用起来比较方便。
首先放一个doc的文档，当你知道什么代码可以实现什么功能之后，就可以自己搜着看文档了，casa的文档介绍的很清楚的
https://casadocs.readthedocs.io/en/stable/api/tt/casatasks.analysis.imregrid.html#imregrid

【importasdm 把asdm格式的casa文件格式转换成.ms 文件夹的格式】

importasdm(asdm='19A-119.sb123243.58235.79924266203',
           vis='19A-119.sb123243.58235.79924266203.ms', ocorr_mode='co',
           with_pointing_correction=True, process_flags=True)

【自动clean每个field】

msfile= '20A-346.60247-spw22.ms'
msfile_h =  msfile[:-3]+'-hanning.ms'
refantenna='ea15'
amp_cals = '3C48'
phase_cals = 'J2355+4950'
all_cals = '3C138,J0319+4130,3C48,J2355+4950'
cals_1 = '3C138'
cals_2 = 'J0319+4130'
targe_n = '*LARGE*'

cal_out_name = 'M31_spw22'
fig_dir = 'fig/'
cleansource=cal_out_name+'_hanning_cal.ms'
targetsource = cal_out_name+'_hanning_target.ms'

for i in range(0,49):
    outname = 'single-field-tclean-10cell-1500niter/M31_spw22_field-%s' %(i)
    field_name = '*LARGE_%s*' %(i)

    tclean(vis=targetsource, imagename=outname, specmode='mfs',niter=1500, gain=0.1, threshold='0mJy',gridder='standard',deconvolver='multiscale',scales=[0, 5, 15, 45], smallscalebias=0.9,interactive=False, imsize=400, cell=['10.0arcsec','10.0arcsec'], stokes='I',weighting='briggs',robust=2, pbcor=False,savemodel='modelcolumn', field=field_name, pblimit=0)

【合并多个field成为一张图】

在这里插入图片描述
例如我有好多个filed的图像，想要合并在一起。

#!/usr/bin/env python

import glob

from casatools import linearmosaic as LM
lm = LM()

"""
ss_image = 'M31_spw6_field-XX.image'
ss_pb = 'M31_spw6_field-XX.pb'
ptList = [25, 37, 38, 39, 40, 41]
images = []
weightimages = []
for pts in ptList:
    images.append(ss_image.replace('XX',str(pts)))
    weightimages.append(ss_pb.replace('XX',str(pts)))
print(images, weightimages)
"""

images = glob.glob('./*.image')
weightimages =  []

for term in images:
    weightimages.append(term.replace('image','pb'))
#需要两个文件单个filed的image文件和对应的pb文件（主波束模型）
#lm.defineoutputimage(nx=4000,ny=4000,cellx='2arcsec', imagecenter='00h41m48.00 +40d28m21.00',outputimage='test.linmos')

lm.defineoutputimage(nx=7200,ny=7200,cellx='2arcsec', imagecenter='00h42m44.00 +41d16m08.00',outputimage='m31.all.linmos')
#这个像素点数量和单个间隔，图像中心，输出文件名词都需要更改一下

lm.makemosaic(images=images,weightimages=weightimages)

【自校准】

1计算N值

在这里插入图片描述

#生成脏图(这里我只对38field做了自校准)
tclean(vis='M31_spw22_hanning_target.ms',imagename='obj.dirty', datacolumn='data', imsize=400, cell='10arcsec', pblimit=0, gridder='standard', field='*38*')

#起始模型，robust取0比较好
tclean(vis='M31_spw22_hanning_target.ms', imagename='obj.prelim_clean', datacolumn='data', imsize=400, cell='10arcsec', pblimit=0, gridder='standard', deconvolver='mtmfs', nterms=2, niter=1000, interactive=True, weighting='briggs', robust=0, savemodel='modelcolumn', field='*38*')

#打开carta看看图像的峰值和rms的比值，这个叫做动态范围。

#验证模型是否正确保存进去
plotms(vis='M31_spw22_hanning_target.ms', xaxis='UVwave', yaxis='amp', ydatacolumn='model', avgchannel='64',field='*38*')

#初始校准表
gaincal(vis='M31_spw22_hanning_target.ms', caltable='selfcal_initial.tb', solint='int', refant='ea15', calmode='p', gaintype='G', minsnr=0, field='*38*')

gaincal(vis='M31_spw22_hanning_target.ms',caltable='selfcal_combine_pol_solint_3.tb',solint='int',refant='ea15',calmode='p',gaintype='T', minsnr=0, field='*38*')
gaincal(vis='M31_spw22_hanning_target.ms',caltable='selfcal_combine_pol_solint_6.tb',solint='6s',refant='ea15',calmode='p',gaintype='T', minsnr=0, field='*38*')
gaincal(vis='M31_spw22_hanning_target.ms',caltable='selfcal_combine_pol_solint_12.tb',solint='12s',refant='ea15',calmode='p',gaintype='T', minsnr=0, field='*38*')
gaincal(vis='M31_spw22_hanning_target.ms',caltable='selfcal_combine_pol_solint_24.tb',solint='24s',refant='ea15',calmode='p',gaintype='T', minsnr=0, field='*38*')
gaincal(vis='M31_spw22_hanning_target.ms',caltable='selfcal_combine_pol_solint_48.tb',solint='48s',refant='ea15',calmode='p',gaintype='T', minsnr=0, field='*38*')
gaincal(vis='M31_spw22_hanning_target.ms',caltable='selfcal_combine_pol_solint_96.tb',solint='96s',refant='ea15',calmode='p',gaintype='T', minsnr=0, field='*38*')

#求解解的间隔选取
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
tb.open( 'selfcal_combine_pol_solint_6.tb' )
snr_6s = tb.getcol( 'SNR' ).ravel()
tb.close()
tb.open( 'selfcal_combine_pol_solint_12.tb' )
snr_12s = tb.getcol( 'SNR' ).ravel()
tb.close()
tb.open( 'selfcal_combine_pol_solint_24.tb' )
snr_24s = tb.getcol( 'SNR' ).ravel()
tb.close()
tb.open( 'selfcal_combine_pol_solint_48.tb' )
snr_48s = tb.getcol( 'SNR' ).ravel()
tb.close()
plt.hist( snr_6s, bins=50, density=True, histtype='step', label='6 seconds' )
plt.hist( snr_12s, bins=50, density=True, histtype='step', label='12 seconds' )
plt.hist( snr_24s, bins=50, density=True, histtype='step', label='24 seconds' )
plt.hist( snr_48s, bins=50, density=True, histtype='step', label='48 seconds' )
plt.legend( loc='upper right' )
plt.xlabel( 'SNR' )
print( 'P(<=6) = {0}  ({1})'.format( stats.percentileofscore( snr_6s, 6 ), '6s' ) )
print( 'P(<=6) = {0}  ({1})'.format( stats.percentileofscore( snr_12s, 6 ), '12s' ) )
print( 'P(<=6) = {0}  ({1})'.format( stats.percentileofscore( snr_24s, 6 ), '24s' ) )
print( 'P(<=6) = {0}  ({1})'.format( stats.percentileofscore( snr_48s, 6 ), '48s' ) )

#应用这个求解间隔
gaincal(vis='M31_spw22_hanning_target.ms', caltable='selfcal_combine_pol_solint_12_minsnr_6.tb', solint='12s', refant='ea15', calmode='p', gaintype='T', minsnr=6,field='*38*')

#应用
applycal(vis='M31_spw22_hanning_target.ms', gaintable='selfcal_combine_pol_solint_12_minsnr_6.tb',field='*38*')

#第一次自校准的检查
tclean(vis='M31_spw22_hanning_target.ms', imagename='obj.selfcal1_clean.3arcmin', datacolumn='corrected', imsize=400, cell='10arcsec', pblimit=0, gridder='standard', deconvolver='mtmfs', nterms=2, niter=1000, interactive=True, weighting='briggs', robust=0, savemodel='none',field='*38*')

【linearmosaic】

把下面的写出.py脚本，然后在casa里面使用execfile(‘script.py’)运行脚本就可以了。注意这个输入不能fits文件

#!/usr/bin/env python

import glob

from casatools import linearmosaic as LM
lm = LM()

images = glob.glob('./m31.spw22.f*.wsclean-MFS-image')
weightimages =  glob.glob('./M31_spw22_field-*.pb')

print(images)
print(weightimages)
#lm.defineoutputimage(nx=4000,ny=4000,cellx='2arcsec', imagecenter='00h41m48.00 +40d28m21.00',outputimage='test.linmos')

lm.defineoutputimage(nx=3000,ny=3000,cellx='5arcsec', imagecenter='00h42m44.00 +41d16m08.00',outputimage='m31.all.20A.5cell.wsclean.linmos.image')


lm.makemosaic(images=images,weightimages=weightimages)

【imregrid重新网格化】

就是把把某个图像映射到某个图像上，python也能实现。
随便贴一下模版，还可以转换坐标系，特别注意网格化带来第三和第四个轴的变化。具体参考casa doc文档

imregrid(imagename="input.image", output="output.image", template="target.image")

如果轴的顺序跟你的需要映射的图像对不上的使用imtrans改一下轴的顺序就好了

imagename = "myim.im"
outfile = "outim.im"
order = "0132"
imtrans()

获取天线位置修正因为网络带来的问题

gencal(vis=msfile_h, caltable=cal_out_name+'.antpos',caltype='antpos')  #校准天线位置

这个命令是一个联网的命令，所以所以网络不好的话大概率就会运行失败。这个尤其在你想对大量的spw做一个并行运算的时候很容易被忽略

去除延迟校准表中延迟高于5ns的值

flagdata(vis=cal_out_name+'.K0', mode='clip', clipminmax=[-5.0, 5.0], datacolumn='FPARAM', clipoutside=True, action='apply', flagbackup=False, savepars=False)