#!/usr/bin/env ParselTongue # reads an Aips data file and plots a dynamic spectrum for all baselines # on a single diagram # # changelog # 07jul2010: reads in multiple IFs # reads in only every nth data point to speed up reading the data # Enno Middelberg 2009 import sys if len(sys.argv)==1: print """\n dynspec-flagger.py written by Enno Middelberg 2009 Usage: dynspec-flagger.py user.name.class.disk.seq [-minmax n,m -h -noint -noaips] where user = your Aips user number name = Aips file name class = Aips file class disk = Aips disk number seq = Aips file sequence optional arguments: -minmax n,m min and max of colour range in plots -h make a PNG hardcopy of dynamic spectra, go into interactive mode afterwards unless -noint is specified -noint don't go into interactive mode (useful for scripting purposes) -noaips read data from pickled file created in a previous run. Aips file needs to exist nevertheless to determine some data properties such as number of antennas etc. -sigmacut n flag baseline when its amplitude standard deviation is larger than n -meancut n flag baseline when its amplitude mean is larger than n -read n read only every nth data point (speeds up reading in the data) \n""" sys.exit(0) from AIPS import AIPS from Wizardry.AIPSData import AIPSUVData import numpy as n import pylab as p import matplotlib.cm as cm import math as m import pprint import cPickle import time as tm import pdb import os from matplotlib.ticker import FuncFormatter # Function to convert seconds to hh:mm:ss.ss format, returns a string def time2hms(seconds): h=int(seconds/3600) m=int(seconds % 3600)/60 s=seconds-(h*3600)-(m*60) h=`h` m=`m` s="%4.2f" % s hms=h.zfill(2)+":"+m.zfill(2)+":"+s.zfill(4) return hms # function to explain usage of interactive display def printhelp(): print """\n You have the following options for the interactive display: left click - open up a window with an enlarged view of the selected subplot shift + left click - flag the selected baseline right click - flag all baselines to the antenna labelled \"Ant XX\" shift + right click - flag all baselines to the antenna labelled \"Bsl XX\" \n""" # this makes meaningful x-axis tic labels # see http://matplotlib.sourceforge.net/examples/pylab_examples/custom_ticker1.html def xticlabels(x, pos): 'The two args are the value and tick position' return '%1.3f' % ((reffreq+(x+refpix)*refinc)/1e9) # this function evaluates the clicks on the plot def onclick(event): # each plot is identified via a dictionary # the key is the subplot object selplt=plots[event.__dict__['inaxes']][0] selant=plots[event.__dict__['inaxes']][1] selbsl=plots[event.__dict__['inaxes']][2] selpol=plots[event.__dict__['inaxes']][3] # left click makes a new plot if event.__dict__['button']==1 and event.__dict__['key']==None: zoom=p.figure() bslname="%i-%i" % (selant, selbsl) data=alldata[selpol][bslname] im = p.imshow(data[:,1:], origin='lower', vmin=vmin, vmax=vmax, aspect='auto', interpolation='nearest') p.xlabel('Channel Frequency / GHz') p.ylabel('Visibility number') ax=p.gca() formatter = FuncFormatter(xticlabels) ax.xaxis.set_major_formatter(formatter) p.title('Data on baseline %s' % bslname) p.show() # shift left click flags a baseline if event.__dict__['button']==1 and event.__dict__['key']=='shift': outstring="antennas=%i baseline=%i stokes=\'%s\' reason=\'dynspec-flagger.py\' /" % (selant, selbsl, polnames[selpol]) print outstring outfile=open('%s.flag' % aipsfile, 'a') outfile.write('%s\n' % outstring) outfile.close() #pprint.pprint(event.__dict__) #temp = p.subplot(nx, ny, selplt) #im = p.imshow(n.array([[-1],[-1]]), origin='lower', vmin=vmin, vmax=vmax, aspect='auto') #fig.add_subplot(temp) #event.__dict__['inaxes'].cla() #event.__dict__['inaxes'].set_ticklabels([], visible=False) #pprint.pprint(fig.__dict__) #fig.__dict__['axes'].pop(fig.__dict__['axes'].index(event.__dict__['inaxes'])) #pprint.pprint(fig.__dict__) #p.draw() # right click flags an antenna if event.__dict__['button']==3 and event.__dict__['key']==None: outstring="antennas=%i stokes=\'%s\' reason=\'dynspec-flagger.py\' /" % (selant, polnames[selpol]) print outstring outfile=open('%s.flag' % aipsfile, 'a') outfile.write('%s\n' % outstring) outfile.close() # shift right click flags a baseline if event.__dict__['button']==3 and event.__dict__['key']=='shift': outstring="antennas=%i stokes=\'%s\' reason=\'dynspec-flagger.py\' /" % (selbsl, polnames[selpol]) print outstring outfile=open('%s.flag' % aipsfile, 'a') outfile.write('%s\n' % outstring) outfile.close() # evaluate keyboard events def on_key_event(event): #pprint.pprint(event.__dict__) # quit if q is pressed if event.__dict__['key']=='q': print " q pressed - exiting.\n" sys.exit(0) # print help if needed if event.__dict__['key']=='h': printhelp() # refresh the plot if event.__dict__['key']=='r': p.draw() ###################################### # # evaluate the command line # aipsfile=sys.argv[1] AIPS.userno=int(aipsfile.split('.')[0]) infile =aipsfile.split('.')[1] inclass=aipsfile.split('.')[2] indisk =int(aipsfile.split('.')[3]) inseq =int(aipsfile.split('.')[4]) print "\n Reading data from file/class/disk/seq %s.%s.%i.%i, user number %i" % \ (infile, inclass, indisk, inseq, int(aipsfile.split('.')[0])) # do we make a hardcopy? if '-h' in sys.argv: print " Storing a hardcopy in %s.png" % aipsfile # do we go into interactive mode? if '-noint' in sys.argv: print " Not going into interactive mode" # generate the UVFLG file to which flags are appended later #pdb.set_trace() now=tm.asctime(tm.localtime()) if os.path.exists('%s.flag' % aipsfile): print " File %s.flag exists - appending new flags to it." % aipsfile outfile=open('%s.flag' % aipsfile, 'a') outfile.write('!\n! Flags appended by dynspec-flagger.py on %s\n!\n' % now) outfile.close() else: print " Creating %s.flag to store the flags." % aipsfile outfile=open('%s.flag' % aipsfile, 'w') outfile.write('! Flags generated by dynspec-flagger.py on %s\n' % now) outfile.write('opcode = \'FLAG\'\n') outfile.close() # flag based on sigma? if "-sigmacut" in sys.argv: index=sys.argv.index('-sigmacut') sigmacut=float(sys.argv[index+1]) print " Flagging data if sigma exceeds %4.3f" % sigmacut outfile=open('%s.flag' % aipsfile, 'a') outfile.write('! Flagging data if sigma exceeds %4.3f\n' % sigmacut) outfile.close() # flag based on mean? if "-meancut" in sys.argv: index=sys.argv.index('-meancut') meancut=float(sys.argv[index+1]) print " Flagging data if mean exceeds %4.3f" % meancut outfile=open('%s.flag' % aipsfile, 'a') outfile.write('! Flagging data if sigma exceeds %4.3f\n' % meancut) outfile.close() # read only every nth data point? read=1 if "-read" in sys.argv and not ('-noaips' in sys.argv): index=sys.argv.index('-read') read=int(sys.argv[index+1]) print " Reading every %i data point(s)" % read if "-read" in sys.argv and '-noaips' in sys.argv: print " ### WARNING -read does have no effect when data is restored from previous run using the -noaips switch" # # ######################################## # open up the Aips data file uvdata=AIPSUVData(infile, inclass, indisk, inseq) # figure out how many plots to draw nant=len(uvdata.antennas) nx=nant ny=nant print " Number of stations is %i, going to draw %i x %i plots" % (nant, nx, ny) # figure out number of channels nchan=uvdata.header['naxis'][uvdata.header['ctype'].index('FREQ')] print " Number of channels: %i" % nchan # store frequency of reference channel and channel increment reffreq=uvdata.header['crval'][uvdata.header['ctype'].index('FREQ')] refinc=uvdata.header['cdelt'][uvdata.header['ctype'].index('FREQ')] refpix=uvdata.header['crpix'][uvdata.header['ctype'].index('FREQ')] # this could be needed to determine the frequencies very accurately #refsign=cmp(refinc, 0.0) #print "refsign=", refsign # figure out number of IFs if 'IF' in uvdata.header['ctype']: nif=uvdata.header['naxis'][uvdata.header['ctype'].index('IF')] print " Number of IFs: %i" % nif else: print " Keyword IF not found in header, assuming number of IFs is 1" nif=1 # try to guess polarisation names AN=uvdata.table('AN', 1) for x in AN: polnames={0: 2*x['poltya'][0], 1: 2*x['poltyb'][0]} print " Polarisation labels of last antenna in AN table are %s, %s, using this as plot labels" % (polnames[0], polnames[1]) # generate a dictionary of baselines data_RR={} data_LL={} for x in range(1, len(uvdata.antennas)): for y in range(x+1, len(uvdata.antennas)+1): code="%i-%i" % (x,y) data_RR[code]=n.zeros([0,nchan*nif+1]) data_LL[code]=n.zeros([0,nchan*nif+1]) # read in the data count=int(0) if not '-noaips' in sys.argv: for row in uvdata: if (count % read)==0: # ignore autocorrelations #if row.baseline[0]==row.baseline[1]: # continue #amp_RR=n.sqrt(row.visibility[0][:,0][:,0]**2 + row.visibility[0][:,0][:,1]**2) #amp_LL=n.sqrt(row.visibility[0][:,1][:,0]**2 + row.visibility[0][:,1][:,1]**2) # # When Aips data is flagged in UVCOP, weights are negated, but data are kept. Here, data is only copied # when weights are positive, otherwise amplitudes are set to -1 # this fills an array upon a condition # n.where(condition, value_if_true, value_if_false) amp_RR=n.zeros((nchan*nif)) amp_LL=n.zeros((nchan*nif)) for j in range(nif): amp_RR[j*nchan:(j*nchan+nchan)]=n.where(row.visibility[j][:,0][:,2]>0.0, (n.sqrt(row.visibility[j][:,0][:,0]**2 + row.visibility[j][:,0][:,1]**2)), -1.0) amp_LL[j*nchan:(j*nchan+nchan)]=n.where(row.visibility[j][:,1][:,2]>0.0, (n.sqrt(row.visibility[j][:,1][:,0]**2 + row.visibility[j][:,1][:,1]**2)), -1.0) #print row #print "baseline: ", row.baseline #print "visibility:", row.visibility #print "amp_RR: ", amp_RR amp_I=(amp_RR+amp_LL)/2.0 t=n.array(row.time) amp_I=n.hstack((t, amp_I)) amp_RR=n.hstack((t, amp_RR)) amp_LL=n.hstack((t, amp_LL)) index="%i-%i" % (row.baseline[0], row.baseline[1]) data_RR[index]=n.vstack((data_RR[index], amp_RR)) data_LL[index]=n.vstack((data_LL[index], amp_LL)) print ' Time, baseline, nvis: %s %5s %i \r' % (time2hms(86400*t), index, data_RR[index].shape[0]), count=count+1 print # put the data into a list so we can loop over it alldata=[data_RR, data_LL] print " Saving data to %s.data..." % aipsfile file=open('%s.data' % aipsfile, 'w') cPickle.dump(alldata, file, 1) file.close() print " ...done." else: print " Reading data from %s.data..." % aipsfile file=open('%s.data' % aipsfile, 'r') alldata=cPickle.load(file) file.close() print " ...done." # initialize a directory for the subplot properties plots={} # initialise the figure fig=p.figure(figsize=(12,12)) fig.subplots_adjust(wspace=0.02, hspace=0.02) fig.suptitle('Data from Aips file %s' % aipsfile) # get min and max for plots if '-minmax' in sys.argv: index=sys.argv.index('-minmax') vmin=float(sys.argv[index+1].split(',')[0]) vmax=float(sys.argv[index+1].split(',')[1]) print " Found min, max = %2.1f, %2.1f for the colour transfer function" % (vmin, vmax) else: maxamps=[] for pol in range(2): for ant in range(1, (len(uvdata.antennas)+1)): for bsl in range((ant+1), (len(uvdata.antennas)+1)): # get the data and only do something if data exists for the actual baseline bslname="%i-%i" % (ant, bsl) data=alldata[pol][bslname] try: maxamps.append(data.max()) except: pass print "Maximum amplitudes per baseline and polarisation: " print n.sort(n.array(maxamps)) vmax=n.array(maxamps).max() vmin=0.0 # vmax=1.0 median=[] # make a colourful plot for pol in range(2): for ant in range(1, (len(uvdata.antennas)+1)): for bsl in range((ant+1), (len(uvdata.antennas)+1)): # get the data and only do something if data exists for the actual baseline bslname="%i-%i" % (ant, bsl) data=alldata[pol][bslname] # the following lines pickle the data from a single array to a file so that it can be analysed off-line #if bslname=='10-27' and pol==0: # print " Saving data on baseline %s to %s.data..." % (bslname, bslname) # file=open('%s.data' % bslname, 'w') # cPickle.dump(data[:,1:], file, 1) # file.close() # print " ...done." # the following lines can be useful for debugging - fill data array with known distribution and analyse it #data=n.random.normal(1,0.1, data.shape) #print data.mean(), data.std() #data[:,0:20]=-1.0 # create a masked array (ignoring flagged=negative data) to calculate mean and stddev masked_data=n.ma.masked_array(data[:,1:], data[:,1:]<0.0) sigma=n.std(masked_data) mean=n.mean(masked_data) median.append([sigma, mean]) output='' print " Now analysing baseline %s, polarisation %s: sigma = %4.3f, mean = %4.3f " % (bslname, polnames[pol], sigma, mean), # do nothing when there is no data if data.shape[0]==0: output=output+"... no data found on this baseline" # when there is data, flag it if requested else: is_flagged=0 reason='' if "-sigmacut" in sys.argv and sigma>sigmacut: data[:,1:]=data[:,1:]*0.0-1.0 output=output+"... flagged due high sigma" reason='sigma=%4.3f' % sigma is_flagged=1 if "-meancut" in sys.argv and mean>meancut: data[:,1:]=data[:,1:]*0.0-1.0 output=output+"... flagged due high mean" is_flagged=1 reason='%s, mean=%4.3f' % (reason, mean) if is_flagged==1: outstring="antennas=%i baseline=%i stokes=\'%s\' reason=\'%s\' /" % (ant, bsl, polnames[pol], reason) outfile=open('%s.flag' % aipsfile, 'a') outfile.write('%s\n' % outstring) outfile.close() # plotting can only be skipped if interactive mode is switched off and no hardcopy is desired: if not ('-noint' in sys.argv and (not '-h' in sys.argv)): # determine the location of the subplot # RR goes into the lower left triangle if pol==0: plotnum=ant+(bsl-1)*nant output=output+"... subplot number %i" % plotnum temp=p.subplot(nx, ny, plotnum) plots[temp]=[plotnum, ant, bsl, pol] # only plot labels at the edges of the diagram if ant==1: p.ylabel('%s Bsl %i' % (polnames[0], bsl), rotation=0) if bsl==nant: p.xlabel('%s Ant %i' % (polnames[0], ant), rotation=90) # LL goes into the upper right triangle else: plotnum=bsl+(ant-1)*nant output=output+"... subplot number %i" % plotnum temp = p.subplot(nx, ny, plotnum) plots[temp]=[plotnum, ant, bsl, pol] temp.yaxis.set_label_position("right") temp.xaxis.set_label_position("top") if ant==1: p.xlabel('%s Bsl %i' % (polnames[1], bsl), rotation=90) if bsl==nant: p.ylabel('%s Ant %i' % (polnames[1], ant), rotation=0) # until I figure out how to draw tick labels only # on the outer axes, no tick labels are drawn whatsoever temp.xaxis.set_ticklabels([], visible=False) temp.yaxis.set_ticklabels([], visible=False) im = p.imshow(data[:,1:], origin='lower', vmin=vmin, vmax=vmax, aspect='auto', interpolation='nearest') #im = p.imshow(data[:,1:], origin='lower', vmin=vmin, vmax=vmax, aspect='auto') #im = p.imshow(data[:,1:], origin='lower', aspect='auto') fig.add_subplot(temp) print output if '-h' in sys.argv: print "\n Writing the hardcopy" fig.set_size_inches((nant*2,nant*2)) p.savefig('%s.png' % aipsfile) median=n.array(median)# print "\n The median of the standard deviations is %4.3f and the median of the means is %4.3f" % (n.median(median[:,0]), n.median(median[:,1])) # plot the distribution of standard deviation and mean #maxbin=median.max() #bins=n.arange(0,maxbin, (maxbin/10)) #hist_std=n.histogram(median[:,0], bins) #hist_mean=n.histogram(median[:,1], bins) #histograms=p.figure() #n, bins_std, patches = p.hist(median[:,0], 10, normed=1, facecolor='green', alpha=0.75) #temp=p.subplot(2,1,1) #p.plot(bins_std) #histograms.add_subplot(temp) # #m, bins_mean, patches = p.hist(median[:,1], 10, normed=1, facecolor='blue', alpha=0.75) #temp=p.subplot(2,1,2) #p.plot(bins_mean) #histograms.add_subplot(temp) #l=p.plot(bins_std, bins_mean) if not '-noint' in sys.argv: fig.set_size_inches((12,12)) printhelp() cid = fig.canvas.mpl_connect('button_press_event', onclick) cld = fig.canvas.mpl_connect('key_press_event', on_key_event) stats=p.figure() p.scatter(median[:,0], median[:,1], marker='+') p.xlabel('Standard deviation') p.ylabel('Mean') p.title('Distribution of the standard deviations and means') p.show()