This is a python program to inspect GBT data under the Obit fits file format.
This is the “manual” of the software:
#!/usr/bin/python # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see <http://www.gnu.org/licenses/>. # # inspectGBT Version 1.0 # Author: Francesco de Gasperin (astro@voo.it) # Program to interactively plot GBT data organized in Obit fits format import sys, os import pyfits import numpy as np import matplotlib import matplotlib.pyplot as plt import matplotlib.colors as clr from matplotlib.widgets import Cursor, Button, RadioButtons, CheckButtons, Slider from mpl_toolkits.axes_grid1 import make_axes_locatable # for decent colorbar plt.rcParams['font.size']=12 class viewer_2d(object): def __init__(self, valON, valOFF, scans, RAs, DECs, time=None, freq=None): """ Shows a given array in a 2d-viewer. Input: z, an 2d array. x,y coordinters are optional. """ self.valON = valON self.valOFF = valOFF self.scans = scans self.RAs = RAs self.DECs = DECs self.z = np.zeros(shape = valON.shape[1:]) # to inizialize the overview if time == None: self.alltime = np.arange(self.valON.shape[1]) else: self.alltime = time if freq == None: self.allfreq = np.arange(self.valON.shape[2]) else: self.allfreq = freq self.fig = plt.figure() # Set parameters for initial overview self.pol = 0 # XX is the default self.cal = 'diff' # ON-OFF is the default self.fft_time = False self.fft_freq = False self.log_time = False self.log_freq = False self.time = self.alltime self.freq = self.allfreq # scan slide (select the first plotted scan and how many scans) slide_ax1 = plt.subplot2grid((20,20),(19,13),rowspan=1,colspan=3) slide_ax2 = plt.subplot2grid((20,20),(19,17),rowspan=1,colspan=3) self.slide_beg_scan = Slider(slide_ax1, 'From', min(self.scans), max(self.scans), valinit=min(self.scans), valfmt='%1i') self.slide_end_scan = Slider(slide_ax2, 'To', min(self.scans), max(self.scans), valinit=max(self.scans), valfmt='%1i') self.slide_beg_scan.slidermax = self.slide_end_scan self.slide_end_scan.slidermin = self.slide_beg_scan # Doing some layout with subplots: self.fig.subplots_adjust(0.05,0.05,0.98,0.98,0.1) self.overview = plt.subplot2grid((20,20),(0,0),rowspan=20,colspan=10) self.overview.set_xlabel('Channel') self.overview.set_ylabel('Sample') divider = make_axes_locatable(self.overview) self.cax = divider.append_axes("bottom", "3%", pad="5%") # creating overview image self.overview_im = self.overview.imshow(self.z, interpolation='nearest', aspect='auto') self.draw_overview() # time/freq subplots self.time_subplot = plt.subplot2grid((20,20),(0,11),rowspan=9,colspan=8) self.freq_subplot = plt.subplot2grid((20,20),(10,11),rowspan=9,colspan=8) #Adding widgets, to not be garbage-collected they are put in a list: # cursor for selecting time/freq cursor = Cursor(self.overview, useblit=True, color='black', linewidth=2 ) # reset button but_ax = plt.subplot2grid((20,20),(15,19),rowspan=1,colspan=1) reset_button = Button(but_ax,'Reset') # legend button but_ax2 = plt.subplot2grid((20,20),(16,19),rowspan=1,colspan=1) legend_button = Button(but_ax2,'Legend') # point button but_ax3 = plt.subplot2grid((20,20),(17,19),rowspan=1,colspan=1) point_button = Button(but_ax3,'Pointings') # quit button but_ax4 = plt.subplot2grid((20,20),(18,19),rowspan=1,colspan=1) quit_button = Button(but_ax4,'Quit') # polarizaitions radio buttons butpol_ax = plt.subplot2grid((20,20),(0,19),rowspan=2,colspan=1) pol_selector = RadioButtons(butpol_ax, labels=['XX','YY','XY','YX'], active=0, activecolor='black') # ON/OFF/diff button butcal_ax = plt.subplot2grid((20,20),(2,19),rowspan=2,colspan=1) cal_selector = RadioButtons(butcal_ax, labels=['diff','ON','OFF'], active=0, activecolor='black') # fft buttons butfft_ax1 = plt.subplot2grid((20,20),(9,11),rowspan=1,colspan=1) butfft_ax2 = plt.subplot2grid((20,20),(19,11),rowspan=1,colspan=1) fft_time_switch = CheckButtons(butfft_ax1, ['Time FFT'], [False]) fft_freq_switch = CheckButtons(butfft_ax2, ['Freq FFT'], [False]) # log buttons butlog_ax1 = plt.subplot2grid((20,20),(9,12),rowspan=1,colspan=1) butlog_ax2 = plt.subplot2grid((20,20),(19,12),rowspan=1,colspan=1) log_time_switch = CheckButtons(butlog_ax1, ['Time log'], [False]) log_freq_switch = CheckButtons(butlog_ax2, ['Freq log'], [False]) self._widgets = [cursor,reset_button,legend_button,point_button,quit_button,pol_selector,fft_time_switch, fft_freq_switch,cal_selector,log_time_switch,log_freq_switch] #connect events reset_button.on_clicked(self.reset_plots) legend_button.on_clicked(self.show_legend) point_button.on_clicked(self.show_point) quit_button.on_clicked(self.quit) pol_selector.on_clicked(self.select_pol) cal_selector.on_clicked(self.select_cal) fft_time_switch.on_clicked(self.fft_data) fft_freq_switch.on_clicked(self.fft_data) log_time_switch.on_clicked(self.log_data) log_freq_switch.on_clicked(self.log_data) self.fig.canvas.mpl_connect('button_press_event',self.click) def quit(self, event): """Kill all windows open""" print "Done." plt.close('all') def get_time_idx(self): """Get time index based on slide selection""" beg_scan = self.slide_beg_scan.val end_scan = self.slide_end_scan.val idx1 = np.where(self.scans >= beg_scan)[0] idx2 = np.where(self.scans <= end_scan)[0] return list(set(idx1) & set(idx2)) def draw_overview(self): """Draw the main plot""" idx = self.get_time_idx() self.time = self.alltime[idx] if self.cal == 'ON': self.z = self.valON[self.pol][idx] if self.cal == 'OFF': self.z = self.valOFF[self.pol][idx] if self.cal == 'diff': self.z = self.valON[self.pol][idx] - self.valOFF[self.pol][idx] self.overview.set_title('cal: '+str(self.cal)+' - pol: '+str(self.pol)) self.overview_im.set_data(self.z) if np.amin(self.z) > 0: self.overview_im.set_norm(clr.LogNorm(vmin=np.amin(self.z), vmax=np.amax(self.z))) else: self.overview_im.set_norm(clr.Normalize(vmin=-np.std(self.z), vmax=np.std(self.z))) self.overview_im.set_extent([np.amin(self.freq),np.amax(self.freq),np.amin(self.time),np.amax(self.time)]) plt.colorbar(self.overview_im, cax=self.cax, orientation='horizontal') def draw_subplot_time(self, x, y, label): """Draw a subplot content""" if self.fft_time: y = abs(np.fft.rfft(y)) x = range(len(y)) if self.log_time: self.time_subplot.set_yscale('log') y = abs(y) else: self.time_subplot.set_yscale('linear') c = self.time_subplot.plot(x, y, label = label) self.time_subplot.relim() self.time_subplot.autoscale_view() return c def draw_subplot_freq(self, x, y, label): """Draw a subplot content""" if self.fft_freq: y = abs(np.fft.rfft(y)) x = range(len(y)) if self.log_freq: self.freq_subplot.set_yscale('log') y = abs(y) else: self.freq_subplot.set_yscale('linear') c = self.freq_subplot.plot(x, y, label = label) self.freq_subplot.relim() self.freq_subplot.autoscale_view() return c def show_legend(self, event): """Shows legend for the plots""" for pl in [self.time_subplot,self.freq_subplot]: if len(pl.lines)>0: pl.legend() plt.draw() def show_point(self, event): """Shows pointings in another plot""" timeidx = self.get_time_idx() self.fig2 = plt.figure() ax = self.fig2.add_subplot(111) ax.set_title('Pointings') ax.set_xlabel('RA [deg]') ax.set_ylabel('DEC [deg]') scansidx = [[i for i, x in enumerate(self.scans[timeidx]) if x == y] for y in np.unique(self.scans[timeidx])] for scanidx in scansidx: ax.plot(self.RAs[timeidx][scanidx], self.DECs[timeidx][scanidx], marker='o') self.fig2.show() def reset_plots(self, event=None, plot=None): """Clears the subplots.""" if plot == None: plots = [self.time_subplot,self.freq_subplot,self.overview] if plot == 'time': plots = [self.time_subplot] if plot == 'freq': plots = [self.freq_subplot] for j in plots: j.lines=[] j.legend_ = None plt.draw() self.draw_overview() def fft_data(self, axis): """Set the FFT in one of the two plots""" if axis == 'Time FFT': self.fft_time = not self.fft_time self.reset_plots('change fft', 'time') if axis == 'Freq FFT': self.fft_freq = not self.fft_freq self.reset_plots('change fft', 'freq') def log_data(self, axis): """Set the Log-scale in one of the two plots""" if axis == 'Time log': self.log_time = not self.log_time self.reset_plots(plot='time') if axis == 'Freq log': self.log_freq = not self.log_freq self.reset_plots(plot='freq') def select_pol(self, pol): """Change the displayed polarization.""" self.pol = {'XX':0,'YY':1,'XY':2,'YX':3}[pol] self.reset_plots() self.draw_overview() def select_cal(self, cal): """Change the displayed polarization.""" self.cal = cal self.reset_plots() self.draw_overview() def click(self,event): """ What to do, if a click on the figure happens: 1. Check which axis 2. Get data coord's. 3. Plot resulting data. 4. Update Figure """ if event.inaxes==self.overview: # Get nearest data freqpos = np.argmin(np.abs(event.xdata-self.freq)) timepos = np.argmin(np.abs(event.ydata-self.time)) # Check which mouse button (1: left, 2: central, 3: right): if event.button == 1 or event.button == 2: c, = self.draw_subplot_freq(self.freq, self.z[timepos,:], label=str(self.time[timepos])) self.overview.axhline(self.time[timepos],color=c.get_color(), lw=2) if event.button == 3 or event.button == 2: c, = self.draw_subplot_time(self.time, self.z[:,freqpos], label=str(self.freq[freqpos])) self.overview.axvline(self.freq[freqpos],color=c.get_color(), lw=2) # Click on onlt subplot will activate the other if event.inaxes == self.freq_subplot: freqpos = np.argmin(np.abs(event.xdata-self.freq)) c, = self.draw_subplot_time(self.time, self.z[:,freqpos],label=str(self.freq[freqpos])) # self.overview.axhline(self.freq[freqpos],color=c.get_color(), lw=2) if event.inaxes == self.time_subplot: timepos = np.argmin(np.abs(event.xdata-self.time)) c, = self.draw_subplot_freq(self.freq, self.z[timepos,:], label=str(self.time[timepos])) # self.overview.axvline(self.time[timepos],color=c.get_color(), lw=2) #Show it plt.draw() if __name__=='__main__': if len(sys.argv) != 2: print "Usage: inspectGBT.py filename.fits" sys.exit(0) try: fitsfile = pyfits.open(sys.argv[1]) except: print "ERROR opening the fits file." sys.exit(1) datatab = fitsfile[4] maskcal = datatab.data.field('CAL') == 1 dataON = datatab.data[maskcal].field('DATA') dataOFF = datatab.data[~maskcal].field('DATA') scans = datatab.data.field('SCAN')[maskcal] RAs = datatab.data.field('RA')[maskcal] DECs = datatab.data.field('DEC')[maskcal] dataON = np.rollaxis(dataON, 1)[0] # move+remove feed dataON = np.rollaxis(dataON, 3)[0] # move+remove weigth dataON = np.rollaxis(dataON, 2) # move pol dataOFF = np.rollaxis(dataOFF, 1)[0] # move+remove feed dataOFF = np.rollaxis(dataOFF, 3)[0] # move+remove weigth dataOFF = np.rollaxis(dataOFF, 2) # move pol #Put it in the viewer inspector = viewer_2d(dataON, dataOFF, scans, RAs, DECs) #Show it plt.show()
Categories: Astronomy, Coding, Python, Science, Telescopes