# 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 from astropy.io import fits import time from distutils.spawn import find_executable from recipes import tec_maps import random def calibrator_select(ms_file,cal_scan_list,num_scans,refants='LA,PT,FD,KP'): x_scan_element=0 last_line=0 sbd_pass=False num_good_antenna=0 possible_cal_scans=[] best_good_antennas={} best_num_good_antennas=0 cal_scan=0 last_line=0 with open(casalog.logfile(),'r') as fil: for cnt,line in enumerate(fil): if cnt>last_line: last_line=cnt while (not sbd_pass) & (x_scan_elementbest_num_good_antennas: print('Best num good antennas is {}'.format(best_num_good_antennas)) best_num_good_antennas=num_antennas_in_spws best_good_antennas=good_antennas possible_cal_scans=[] if num_antennas_in_spws==best_num_good_antennas: possible_cal_scans.append(cal_scan_list[x_scan_element]) if not sbd_pass: x_scan_element+=1 else: cal_scan=cal_scan_list[x_scan_element] if cal_scan: return cal_scan,best_good_antennas else: return possible_cal_scans,best_good_antennas #define working directory (should be directory that has idifits file in it) wd=os.getcwd() #making variable for the idifits file files=os.listdir(wd) for fl in files: if fl.endswith('.idifits'): idifitsfile=wd+'/'+fl #making various directories to keep files organized os.mkdir('cal_tables') os.mkdir('images') os.mkdir('ms_files') cal_dir=wd+'/cal_tables/' im_dir=wd+'/images/' ms_dir=wd+'/ms_files/' os.chmod(cal_dir,0777) os.chmod(im_dir,0777) os.chmod(ms_dir,0777) #*************************************** # Uses numpy to determine a 4th order # polynomial from the values for # gaincurve response #*************************************** def transform_poly(coeff, min_elev=0, max_elev=90): f = np.poly1d(coeff[::-1]) g = lambda x: np.sqrt(f(90 - x)) x = np.linspace(min_elev, max_elev, 64, endpoint=True) y = g(x) return np.poly1d(np.polyfit(x, y, 3)) #*************************************** # Makes dictionary for all antennas incl # in the observation #*************************************** def make_ant_dict(hdu): ant_dict={} for i in range(0,len(hdu['antenna'].data['anname'])): ant_dict[hdu['antenna'].data['antenna_no'][i]]=hdu['antenna'].data['anname'][i] return ant_dict #*************************************** # VLBA frequency ranges to find what # band the observations were taken in # needs to be expanded to include bands # other than s/x #*************************************** def find_band_letter(freq): if freq>2.2e9 and freq<2.5e9: return 'S' if freq>7.9e9 and freq<8.9e9: return 'X' if freq>3.5e9 and freq<7.9e9: return 'C' else: return 'DUH' #*************************************** # Open the idifits file in filein using # astropy (need older version of astropy) #*************************************** hdu=fits.open(idifitsfile) #*************************************** # Makes a dictionary containing all # antennas in the idifits file #*************************************** ant_dict=make_ant_dict(hdu) #*************************************** # Placeholder for timerange. Since read # directly from idifits file the GC # values should be valid for current obs #*************************************** time_range = '' if not time_range: t = time.strptime("2000y01d00h00m00s", "%Yy%jd%Hh%Mm%Ss") btime = time.mktime(t) + 40587.0 * 86400 t = time.strptime("2100y01d00h00m00s", "%Yy%jd%Hh%Mm%Ss") etime = time.mktime(t) + 40587.0 * 86400 #*************************************** # Order of columns needed to make the # CASA table # # L_Poly is the polynomial for the left # hand polarization and R_Poly is right # # Band # Low_Freq # High_Freq # Begin_Time # End_Time # Antenna # L_Poly_One # L_Poly_Two # L_Poly_Three # L_Poly_Four # R_Poly_One # R_Poly_Two # R_Poly_Three # R_Poly_Four #*************************************** row_dict={} with open(idifitsfile.split('.idifits')[0]+'.txt','w') as f: for i in range(0,len(hdu['gain_curve'].data)): for j in range(0,hdu['gain_curve'].header['no_band']): f.write(str(find_band_letter(hdu['FREQUENCY'].header['REF_FREQ']+hdu['FREQUENCY'].data['BANDFREQ'][0][j]))) f.write(' ') f.write(str(hdu['FREQUENCY'].header['REF_FREQ']+hdu['FREQUENCY'].data['BANDFREQ'][0][j])) f.write(' ') f.write(str(hdu['FREQUENCY'].header['REF_FREQ']+hdu['FREQUENCY'].data['BANDFREQ'][0][j]+hdu['frequency'].data['total_bandwidth'][0][0])) f.write(' ') f.write(str(btime)) f.write(' ') f.write(str(etime)) f.write(' ') f.write(str(ant_dict[hdu['gain_curve'].data['antenna_no'][i]])) f.write(' ') poly=transform_poly(np.split(hdu['gain_curve'].data['gain_1'][i],hdu['gain_curve'].header['no_band'])[j]) f.write(' ') f.write(str(poly[0]*np.sqrt(hdu['gain_curve'].data[i]['sens_1'][j]))) f.write(' ') f.write(str(poly[1]*np.sqrt(hdu['gain_curve'].data[i]['sens_1'][j]))) f.write(' ') f.write(str(poly[2]*np.sqrt(hdu['gain_curve'].data[i]['sens_1'][j]))) f.write(' ') f.write(str(poly[3]*np.sqrt(hdu['gain_curve'].data[i]['sens_1'][j]))) f.write(' ') f.write(str(poly[0]*np.sqrt(hdu['gain_curve'].data[i]['sens_1'][j]))) f.write(' ') f.write(str(poly[1]*np.sqrt(hdu['gain_curve'].data[i]['sens_1'][j]))) f.write(' ') f.write(str(poly[2]*np.sqrt(hdu['gain_curve'].data[i]['sens_1'][j]))) f.write(' ') f.write(str(poly[3]*np.sqrt(hdu['gain_curve'].data[i]['sens_1'][j]))) f.write('\n') # Defining column names for the CASA gaincurve table columnnames = ["BANDNAME","BFREQ","EFREQ","BTIME","ETIME","ANTENNA","GAIN"] # Defining datatypes for columns so that CASA can read the txt file datatypes = ["A","D","D","D","D","A","R4,2"] # creating CASA table to use for gaincurve to do initial amplitude calibration tb.fromascii(idifitsfile.split('.idifits')[0]+'.gc', asciifile=idifitsfile.split('.idifits')[0]+'.txt', sep=' ',columnnames=columnnames, datatypes=datatypes) #import idifits file as ms, put ms in ms_dir ms_file=ms_dir+idifitsfile.split('/')[-1].split('.')[0]+'.ms' cal_file_root=cal_dir+idifitsfile.split('/')[-1].split('.')[0] importfitsidi(fitsidifile=idifitsfile, vis=ms_file, scanreindexgap_s=15.0) # Run listobs so I can find where sources failed listobs(vis=ms_file) #for root,dirs,files in os.walk(ms_file): # for d in dirs: # os.chmod(os.path.join(root,d),0777) # for f in files: # os.chmod(os.path.join(root,f),0777) os.chmod(ms_file,0777) #create tsys and gc calibration tables gencal(vis=ms_file, caltable=cal_file_root+'.tsys', caltype='tsys', uniform=False) gencal(vis=ms_file, caltable=cal_file_root+'.gc', caltype='gc', infile=idifitsfile.split('.idifits')[0]+'.gc') #create ionosphere correction #USNO Network is down, cannot access ionosphere calibration files. This is being commented out #07/16/2020 #tec_maps.create(vis=ms_file,imname='iono') #gencal(vis=ms_file,caltable=cal_file_root+'.tec',caltype='tecim',infile='iono.IGS_TEC.im') #changing antenna table to get the diameters of the vlba antennas. #This might be something that is possible in the future just by reading the idifits file? msmd.open(ms_file) diams= len(msmd.antennaids())*[25.0] tb.open(ms_file+'/ANTENNA', nomodify=False) tb.putcol('DISH_DIAMETER', diams) tb.close() msmd.close() #DiFX correlator correction accor(vis=ms_file, caltable=cal_file_root+'.accor', solint='20s') #apply above calibrations #below is commented out because of the ionosphere calibration not working right now #applycal(vis=ms_file, #gaintable=[cal_file_root+'.tec',cal_file_root+'.accor',cal_file_root+'.gc',cal_file_root+'.tsys'], #interp=['nearest','nearest','nearest','nearest,nearest']) applycal(vis=ms_file, gaintable=[cal_file_root+'.accor',cal_file_root+'.gc',cal_file_root+'.tsys'], interp=['nearest','nearest','nearest,nearest']) flagdata(vis=ms_file,mode='quack',quackinterval=1.0,quackmode='endb') flagdata(vis=ms_file,mode='quack',quackinterval=1.0,quackmode='beg') #When doing the total fringefit I think we need to avoid scans that don't have enough #baselines. I'm going to assume we need at least 4 baselines (because I think that's #what you need for an image anyway) The code below will automatically find those scans. msmd.open(ms_file) scans=msmd.scannumbers() scn_low_bsln=[] counter=0 cal_scan_list=[] possible_cal_scans={} for scan in scans: ants=msmd.antennasforscan(scan) if len(ants)<5: scn_low_bsln.append(scan) if (scan<650) & (scan>50): if str(msmd.antennasforscan(scan))==str(msmd.antennaids()): cal_scan_list.append(scan) else: possible_cal_scans[scan]=msmd.antennasforscan(scan) num_ants_pcs=0 if not cal_scan_list: for scan in list(possible_cal_scans): if len(possible_cal_scans[scan])==num_ants_pcs: cal_scan_list.append(scan) if len(possible_cal_scans[scan])>num_ants_pcs: cal_scan_list=[] num_ants_pcs=len(possible_cal_scans[scan]) cal_scan_list.append(scan) # Below will change the array id (subarray). Maybe fringefit will read this and actually work for # Subarrays? tb.open(ms_file, nomodify=False) scan_table=tb.getcol('SCAN_NUMBER') new_arrayid=np.zeros(len(scan_table)) for scan in scn_low_bsln: new_arrayid[np.where(scan_table==scan)[0]]=1 tb.putcol('ARRAY_ID', new_arrayid) tb.close() #run aoflagger on the data! os.system('chmod -R 777 {current_directory}'.format(current_directory=wd)) ms_file_x=ms_file.split('.ms')[0]+'_x.mms' ms_file_s=ms_file.split('.ms')[0]+'_s.mms' partition(vis=ms_file,outputvis=ms_file_x,separationaxis='scan',spw='4~15',datacolumn='corrected') partition(vis=ms_file,outputvis=ms_file_s,separationaxis='scan',spw='0~3',datacolumn='corrected') os.system('chmod -R 777 {current_directory}'.format(current_directory=wd)) if find_executable('aoflagger'): print('running aoflagger!') flag_command='aoflagger {file_name}'.format(file_name=ms_file_s) print(flag_command) os.system(flag_command) flag_command='aoflagger {file_name}'.format(file_name=ms_file_x) print(flag_command) os.system(flag_command) #Flag the autocorrelations flagdata(vis=ms_file_x, autocorr=True, mode='manual') flagdata(vis=ms_file_s, autocorr=True, mode='manual') #Flag the channels at the edge of the bandpass. These channels don't calibrate well msmd.open(ms_file_x) spws='' for sp in range(0,msmd.nspw()): if sp==msmd.nspw()-1: spws=spws+'{specwin}:0~{low_flag_chan};{high_flag_chan}~{num_chans}'.format(specwin=sp,low_flag_chan=int(msmd.nchan(sp)*0.12),high_flag_chan=msmd.nchan(sp)-1-int(msmd.nchan(sp)*0.12),num_chans=msmd.nchan(sp)-1) else: spws=spws+'{specwin}:0~{low_flag_chan};{high_flag_chan}~{num_chans},'.format(specwin=sp,low_flag_chan=int(msmd.nchan(sp)*0.12),high_flag_chan=msmd.nchan(sp)-1-int(msmd.nchan(sp)*0.12),num_chans=msmd.nchan(sp)-1) flagdata(vis=ms_file_x, mode='manual', spw=spws) msmd.close() msmd.open(ms_file_s) spws='' for sp in range(0,msmd.nspw()): if sp==msmd.nspw()-1: spws=spws+'{specwin}:0~{low_flag_chan};{high_flag_chan}~{num_chans}'.format(specwin=sp,low_flag_chan=int(msmd.nchan(sp)*0.12),high_flag_chan=msmd.nchan(sp)-1-int(msmd.nchan(sp)*0.12),num_chans=msmd.nchan(sp)-1) else: spws=spws+'{specwin}:0~{low_flag_chan};{high_flag_chan}~{num_chans},'.format(specwin=sp,low_flag_chan=int(msmd.nchan(sp)*0.12),high_flag_chan=msmd.nchan(sp)-1-int(msmd.nchan(sp)*0.12),num_chans=msmd.nchan(sp)-1) flagdata(vis=ms_file_s, mode='manual', spw=spws) msmd.close() #Need to split out s and x data separately to do fringefitting and applying os.system('chmod -R 777 {current_directory}'.format(current_directory=wd)) def make_refant_list(inner_ants,outer_ants,initial_ant): refant=initial_ant inner_ants=np.setdiff1d(inner_ants,np.array(initial_ant.split(','))) for i in reversed(range(len(inner_ants))): rand_int=random.randint(0,i) refant=refant+','+inner_ants[rand_int] inner_ants=inner_ants[inner_ants!=inner_ants[rand_int]] for i in reversed(range(len(outer_ants))): rand_int=random.randint(0,i) refant=refant+','+outer_ants[rand_int] outer_ants=outer_ants[outer_ants!=outer_ants[rand_int]] return refant #Initial fringefit to correct for delay between bandpass (instrumental effect) msmd.open(ms_file_x) antennas=np.array(msmd.antennanames()) inner_antennas=np.array(['LA','PT','OV','KP','FD']) outer_antennas=np.array(['BR','HN','NL','MK','SC']) missing_inner_antennas=np.setdiff1d(inner_antennas,antennas) missing_outer_antennas=np.setdiff1d(outer_antennas,antennas) inner_antennas_available=np.setdiff1d(inner_antennas,missing_inner_antennas) outer_antennas_available=np.setdiff1d(outer_antennas,missing_outer_antennas) cal,ants=calibrator_select(ms_file=ms_file_x,refants='LA',cal_scan_list=cal_scan_list,num_scans=len(cal_scan_list)) cals={} cal_scan=0 if type(cal)==list: if 'LA' in inner_antennas_available: cal_la,ants_la=calibrator_select(ms_file=ms_file_x,refants='LA',cal_scan_list=cal_scan_list,num_scans=15) cals['LA']=cal_la if 'FD' in inner_antennas_available: cal_fd,ants_fd=calibrator_select(ms_file=ms_file_x,refants='FD',cal_scan_list=cal_scan_list,num_scans=15) cals['FD']=cal_fd if type(cal_fd)==int: cal_scan=cal_fd refants=make_refant_list(inner_ants=inner_antennas_available,outer_ants=outer_antennas_available,initial_ant='FD') if 'PT' in inner_antennas_available: cal_pt,ants_pt=calibrator_select(ms_file=ms_file_x,refants='PT',cal_scan_list=cal_scan_list,num_scans=15) cals['PT']=cal_pt if type(cal_pt)==int: cal_scan=cal_pt refants=make_refant_list(inner_ants=inner_antennas_available,outer_ants=outer_antennas_available,initial_ant='PT') if 'KP' in inner_antennas_available: cal_kp,ants_kp=calibrator_select(ms_file=ms_file_x,refants='KP',cal_scan_list=cal_scan_list,num_scans=15) cals['KP']=cal_kp if type(cal_kp)==int: cal_scan=cal_kp refants=make_refant_list(inner_ants=inner_antennas_available,outer_ants=outer_antennas_available,initial_ant='KP') if not cal_scan: ant_and_spw_for_refant={} ant_and_spw_for_refant['KP']=0 ant_and_spw_for_refant['PT']=0 ant_and_spw_for_refant['LA']=0 ant_and_spw_for_refant['FD']=0 for kee in ants_kp.keys(): ant_and_spw_for_refant['KP']+=len(ants_kp[kee]) for kee in ants_la.keys(): ant_and_spw_for_refant['LA']+=len(ants_la[kee]) for kee in ants_fd.keys(): ant_and_spw_for_refant['FD']+=len(ants_fd[kee]) for kee in ants_pt.keys(): ant_and_spw_for_refant['PT']+=len(ants_pt[kee]) refants='' for ant in sorted(ant_and_spw_for_refant, key=lambda ant: ant_and_spw_for_refant[ant], reverse=True): refants=refants+ant+',' refants=make_refant_list(inner_ants=inner_antennas_available,outer_ants=outer_antennas_available,initial_ant=refants[:-1]) cal_scan=cals[refants[0:2]][random.randint(0,len(cals[refants[0:2]])-1)] else: cal_scan=cal refants=make_refant_list(inner_ants=inner_antennas_available,outer_ants=outer_antennas_available,initial_ant='LA') fringefit(vis=ms_file_x, caltable=cal_file_root+'_x.sbd', scan=str(cal_scan), solint='inf', zerorates=True, refant=refants, minsnr=5.0, parang=True) fringefit(vis=ms_file_s, caltable=cal_file_root+'_s.sbd', scan=str(cal_scan), solint='inf', zerorates=True, refant=refants, minsnr=5.0, parang=True) #Below code to select a good solution interval is not necessary at this time #safe_solint=[] #for tme in range(15,70): # safe_time=True # for num in scans: # time=msmd.timesforscan(num)[-1]-msmd.timesforscan(num)[0] # if time%tme<5: # safe_time=False # if safe_time: # safe_solint.append(tme) fringefit(vis=ms_file_x, caltable=cal_file_root+'_x.mbd', solint='inf', refant=refants, minsnr=5.0, combine='spw', gaintable=[cal_file_root+'_x.sbd'], interp=['nearest'], parang=True) fringefit(vis=ms_file_s, caltable=cal_file_root+'_s.mbd', solint='inf', refant=refants, minsnr=5.0, combine='spw', gaintable=[cal_file_root+'_s.sbd'], interp=['nearest'], parang=True) applycal(vis=ms_file_x, gaintable=[cal_file_root+'_x.sbd',cal_file_root+'_x.mbd'], interp=['nearest','nearest','nearest','nearest,nearest','nearest','linear'], spwmap=[[],12*[0]], parang=True) applycal(vis=ms_file_s, gaintable=[cal_file_root+'_s.sbd',cal_file_root+'_s.mbd'], interp=['nearest','linear'], spwmap=[[],4*[0]], parang=True) os.mkdir('x_split_ms') x_split_root=wd+'/x_split_ms/' os.system('chmod -R 777 {current_directory}'.format(current_directory=wd)) fields=msmd.fieldnames() msmd.close() for fld in fields: split(vis=ms_file_x,field=fld,outputvis=x_split_root+fld+'.ms') os.mkdir('s_split_ms') s_split_root=wd+'/s_split_ms/' os.system('chmod -R 777 {current_directory}'.format(current_directory=wd)) for fld in fields: split(vis=ms_file_s,field=fld,outputvis=s_split_root+fld+'.ms')