# This is a python script that will calibrate and image USNO 24 hour experiments using NRAOs CASA software (https://casa.nrao.edu/) # This script has been written to work with CASA 5.7. It is not guaranteed to work with any other version of CASA. import os import astropy as ap import numpy as np wd=os.getcwd() dirs=os.listdir(wd) for dr in dirs: if 'image' in dr: im_dir=wd+'/'+dr+'/' if 'x_split' in dr: x_split_dir=wd+'/'+dr+'/' # Baseline lengths for the VLBA. These are in the order of the baseline names as they appear in CASA and if you pass # msmd.baselines() for a given array it will return the array of baseline lengths for that measurement set. # you can use this to set pixel size in CASA automatically! Horray! baseline_lengths=np.array([[0,2345000,3657000,1913000,1757000,4398000,2300000,1214000,1806000,5767000], [2345000,0,3105000,744000,608000,5134000,1654000,1508000,564000,4143000], [3657000,3105000,0,3623000,3006000,7502000,1611000,3885000,3226000,2853000], [1913000,744000,3623000,0,652000,4466000,2075000,845000,417000,4839000], [1757000,608000,3006000,652000,0,4970000,1432000,1088000,236000,4458000], [4398000,5134000,7502000,4466000,4970000,0,6156000,4015000,4795000,8611000], [2300000,1654000,1611000,2075000,1432000,6156000,0,2328000,1663000,3645000], [1214000,1508000,3885000,845000,1088000,4015000,2328000,0,973000,5460000], [1806000,564000,3226000,417000,236000,4795000,1663000,973000,0,4579000], [5767000,4143000,2853000,4839000,4458000,8611000,3645000,5460000,4579000,0]]) # Task to calculate pixel size for a given measurement set. Looks at baselines in the measurement set and calculates # beam size based on maximum baseline and frequency of observations in the measurement set. Will return a string def get_unimaged(img_dir=im_dir,splt=x_split_dir): imgs=[f.split('_X')[0] for f in os.listdir(img_dir) if '_X' in f] srcs=[f.split('.ms')[0] for f in os.listdir(splt) if '.flag' not in f] imgs=np.array(imgs) srcs=np.array(srcs) unimaged=np.setdiff1d(srcs,imgs) unimaged_ms=[f+'.ms' for f in unimaged] return unimaged_ms def calc_pixel_size(msin): msmd.open(msin) spws=msmd.spwsforscan(msmd.scannumbers()[0]) mean_freq_spw_selector=int(len(spws)/2) frequency=msmd.meanfreq(spws[mean_freq_spw_selector]) wavelength_m=3e8/frequency max_baseline=max(baseline_lengths[msmd.baselines()]) hpbw_rad=wavelength_m/max_baseline hpbw_deg=hpbw_rad*180/np.pi hpbw_mas=hpbw_deg*3600*1000 pixel_size_mas=hpbw_mas/5. msmd.close() return '{pxl_size}mas'.format(pxl_size=np.round(pixel_size_mas,3)) def image_loop(ms_file, im_name, im_size=512, robust_weight=2, clean_to=8.0, start_high_mask_threshold=15.0, end_high_mask_threshold=5.0, low_mask_threshold=10.0, prllel=False): cell_size=calc_pixel_size(ms_file) old_nomask_error_line=1 current_nomask_error_line=2 with open(casalog.logfile(),'r') as fil: for cnt,line in enumerate(fil): if ('WARN' in line) & ('automasking' in line) & (cnt>current_nomask_error_line): current_nomask_error_line=cnt tmp_noisethreshold=start_high_mask_threshold while current_nomask_error_line != old_nomask_error_line: old_nomask_error_line=current_nomask_error_line os.system('rm -rf {all_im_prods}'.format(all_im_prods=im_name)) tclean(vis=ms_file, field='0', datacolumn='corrected', imagename=im_name, imsize=im_size, cell=cell_size, deconvolver='hogbom', weighting='briggs', robust=robust_weight, niter=40000, gain=0.05, nsigma=clean_to, usemask='auto-multithresh', sidelobethreshold=0.7, noisethreshold=tmp_noisethreshold, lownoisethreshold=low_mask_threshold, restart=True, savemodel='modelcolumn', parallel=prllel) with open(casalog.logfile(),'r') as fil: for cnt,line in enumerate(fil): if ('WARN' in line) & ('automasking' in line) & (cnt>current_nomask_error_line): if tmp_noisethreshold>end_high_mask_threshold: tmp_noisethreshold=tmp_noisethreshold-1.0 low_mask_threshold=tmp_noisethreshold*0.6666 current_nomask_error_line=cnt else: current_nomask_error_line=old_nomask_error_line if tmp_noisethreshold>end_high_mask_threshold: casalog.post('Image made for source {src} using noisethreshold {nth}'.format(src=im_name.split('/')[-1],nth=tmp_noisethreshold)) if tmp_noisethreshold==end_high_mask_threshold: casalog.post('Possible no image made for source {src} because masking threshold is at limit of {nth} sigma'.format(src=im_name.split('/')[-1],nth=tmp_noisethreshold)) return tmp_noisethreshold,low_mask_threshold def refant_list(msin): msmd.open(msin) antennas=msmd.antennanames(msmd.antennaids()) refant='' if 'LA' in antennas: refant+='LA,' if 'KP' in antennas: refant+='KP,' if 'PT' in antennas: refant+='PT,' if 'FD' in antennas: refant+='FD,' if 'OV' in antennas: refant+='OV,' if 'BR' in antennas: refant+='BR,' if 'NL' in antennas: refant+='NL,' if 'MK' in antennas: refant+='MK,' if 'HN' in antennas: refant+='HN,' if 'SC' in antennas: refant+='SC,' msmd.close() return refant[:-1] src_ms_list=get_unimaged(img_dir=im_dir,splt=x_split_dir) for src in src_ms_list: src_root=im_dir+src.split('.ms')[0]+'_X' os.mkdir(src_root) src_root=src_root+'/' casalog.setlogfile(src_root+src.split('.ms')[0]+'.log') src_cal_dir=src_root+'cal_dir/' src_cal_root=src_cal_dir+src.split('.ms')[0] os.mkdir(src_cal_dir) src_img_dir=src_root+'image_dir/' src_img_root=src_img_dir+src.split('.ms')[0] os.mkdir(src_img_dir) current_ms=x_split_dir+src listobs(vis=current_ms) high_noisethreshold_1,low_noisethreshold_1=image_loop(ms_file=current_ms,im_name=src_img_root+'_1',prllel=False) msmd.open(current_ms) safe_solint=[] for tme in range(20,60): safe_time=True for num in msmd.scannumbers(): time=msmd.timesforscan(num)[-1]-msmd.timesforscan(num)[0] print(time) if (time%tme<5) | (time<60): print('time={} residual={} cannot use'.format(tme,time%tme)) safe_time=False if safe_time: safe_solint.append(tme) if len(safe_solint)>1: p_solint=str(safe_solint[0])+'s' else: p_solint='inf' rf_ant=refant_list(current_ms) gaincal(vis=current_ms, caltable=src_cal_root+'.K', solint='inf', refant=rf_ant, minblperant=3, gaintype='K') gaincal(vis=current_ms, caltable=src_cal_root+'.p', solint=p_solint, refant=rf_ant, minblperant=3, gaintype='G', calmode='p',gaintable=[src_cal_root+'.K']) applycal(vis=current_ms, gaintable=[src_cal_root+'.K',src_cal_root+'.p']) high_noisethreshold_2,low_noisethreshold_2=image_loop(ms_file=current_ms, im_name=src_img_root+'_2', clean_to=5.0, start_high_mask_threshold=high_noisethreshold_1, low_mask_threshold=low_noisethreshold_1,prllel=False) if p_solint!='inf': gaincal(vis=current_ms,caltable=src_cal_root+'.ap', solint=str(safe_solint[0]*3)+'s', refant=rf_ant, minblperant=4, solnorm=True, gaintype='G', calmode='ap', gaintable=[src_cal_root+'.K',src_cal_root+'.p']) else: gaincal(vis=current_ms,caltable=src_cal_root+'.ap', solint=p_solint, refant=rf_ant, minblperant=4, solnorm=True, gaintype='G', calmode='ap', gaintable=[src_cal_root+'.K',src_cal_root+'.p']) applycal(vis=current_ms, gaintable=[src_cal_root+'.K',src_cal_root+'.p',src_cal_root+'.ap']) x,y=image_loop(ms_file=current_ms,im_name=src_img_root+'_3',clean_to=3.0,start_high_mask_threshold=5.0,end_high_mask_threshold=4.0,low_mask_threshold=3.0,prllel=False)