""" Plot the total power spectrum for the files in a given folder Given a folder, finds all the spectra on it, parse information from their names and header and plot them. wmax, Nov 6, 2014 """ import scipy # Remember to use pyplot, because this is the matplotlib recommendation import matplotlib.pyplot as pyplot import os import sys import imp import time # load script configuration parameters from PWD # this will ensure that only the configredrfi.py file on current directory is # loaded PWD = os.environ['PWD'] try: configredrfi = imp.load_source('configredrfi', os.path.join(PWD, 'configredrfi.py')) from configredrfi import * except IOError: error_message = """No configredrfi.py file was found in %s use 'redrfi --help' for help\n""" % (PWD) sys.exit(error_message) #------------------------------------------------------------------------------- def load_spectrum(filename): """ Load the spectrum file created with POSSM in AIPS input ----- filename the name of the file output ------ data a dictionary with the data notes ----- The file has a header. The first line of data comes after a line with column names as the example below shows Channel IF Polar Frequency Velocity Real(Jy) Imag(Jy) 1 1 RR 2090.7500 0.00000 40.8612 0.00000 2 1 RR 2090.7656 -2.24047 0.990043 0.00000 3 1 RR 2090.7812 -4.48093 0.531951 0.00000 wmax, Nov 6, 2014 """ fileObj = open(filename) # Arrays to save data points CHAN = [] IF = [] POL = [] FREQ = [] VEL = [] REAL = [] IMAG = [] # Reading state. Goes to True after column description reading_data = False for line in fileObj: # Matches the line with column names, just before the data start if line[0:7] == 'Channel': reading_data = True continue elif reading_data == True: # Parse the line and save to arrays fields = line.strip().split() CHAN.append(float(fields[0])) IF.append(int(fields[1])) POL.append(fields[2]) FREQ.append(float(fields[3])) VEL.append(float(fields[4])) REAL.append(float(fields[5])) IMAG.append(float(fields[6])) # Save in dictionary as scipy.arrays for easier manipulation data = {'chan': scipy.array(CHAN), 'if': scipy.array(IF), 'pol': scipy.array(POL), 'freq': scipy.array(FREQ), 'vel': scipy.array(VEL), 'real': scipy.array(REAL), 'imag': scipy.array(IMAG) } return data #------------------------------------------------------------------------------- def plot_total_power_single_scan(data, ylim=None, filename='total_power.pdf'): """ Plot total power output from AIPS input ----- data Dictionary with data, output of load_data() ylim Limit for y-axis. If none uses max and min from data filename Output plot file name output ------ A file with 'filename' with the plot with requested normalization notes ----- wmax, Nov 8, 2014 """ # Verify the file is single polarizarition. Abort if not assert len(scipy.unique(data['pol'])) == 1,\ 'More than one polarization in file' # Determine the number of IFs in the file n_ifs = len(scipy.unique(data['if'])) # Configure plot parameters pyplot.rcParams['xtick.labelsize'] = 6 pyplot.rcParams['ytick.labelsize'] = 6 pyplot.rcParams['figure.subplot.wspace'] = 0.0 # Find a common scale for the plots if not given if ylim == None: ymax = max(data['real']) ymin = min(data['real']) else: ymin, ymax = ylim fig = pyplot.figure(figsize=(4.0 * n_ifs, 2)) # Loop over IFs for i in range(1, n_ifs + 1): if_index = scipy.where(data['if'] == i)[0] pyplot.subplot(1, n_ifs, i) pyplot.plot(data['freq'][if_index], data['real'][if_index], '-') pyplot.ylim(ymin, ymax) pyplot.ticklabel_format(axis='x', style='sci', scilimits=(0, 0)) pyplot.tick_params(axis='x', which='major') # remove yticks for all plots except the first one if i > 1 and i < n_ifs: pyplot.tick_params(axis='y', which='both', labelleft='off', labelright='off') # Add right yticks for last plot elif i == n_ifs: pyplot.tick_params(axis='y', which='both', labelleft='off', labelright='on') pyplot.savefig(filename, bbox_inches='tight') # Close figure to avoid filling up memory pyplot.close() # Restore rcParams pyplot.rcdefaults() return #------------------------------------------------------------------------------- def compute_spectra_envelope_if(data_list, normalize=False, method='minrms'): """ Takes a list of scans with as given antenna and returns the spectral total power envelope. This envelope is the minimum, maximum total power observed for all the scans. The envelope is obtained for each frequency channel and is meant to give a sense of the range of observed total powers input ----- data_list A list of data for scans. Each element like the output of load_spectrum() normalize Normalize the spectrum by the median to get values around 1 method Chooses a normlaization method. Two are currently implemented, but more can be implemented later output ------ data_env Same structure as data dictionaty, but with max and min for real and imag notes ----- wmax, Nov 8, 2014 """ # Check the chan, if, pol and freq are the same. Take first as reference # If not quit chan_ref = data_list[0]['chan'] if_ref = data_list[0]['if'] pol_ref = data_list[0]['pol'] freq_ref = data_list[0]['freq'] n_spectra = len(data_list) # Normalize if requested # In the normalization for S-band, only consider band in 2.2 - 2.4 GHz band # Use mean for bands with no overlap with that band if normalize == True: # Find common frequencies with 2200 - 2400 MHz band, not in strong RFI # band in 2315 - 2345 MHz i_freq_in_band =\ scipy.where((2200.0 <= freq_ref) & (freq_ref <= 2400.0))[0] i_freq_rfi = scipy.where((2300.0 <= freq_ref) & (freq_ref <= 2360.0))[0] i_freq_in_band = scipy.array([i for i in i_freq_in_band if\ i not in i_freq_rfi]) # If empty, it means there is no overlap, normalize with all band if len(i_freq_in_band) == 0: if method == 'mean': norm_factor = scipy.mean(data_list[0]['real']) elif method == 'median': norm_factor = scipy.median(data_list[0]['real']) elif method == 'minrms': norm_factor = scipy.sum(data_list[0]['real']**2) /\ scipy.sum(data_list[0]['real']) # Or normalize using the channels on the passband else: if method == 'mean': norm_factor = scipy.mean(data_list[0]['real'][i_freq_in_band]) elif method == 'median': norm_factor = scipy.median(data_list[0]['real'][i_freq_in_band]) elif method == 'minrms': norm_factor =\ scipy.sum(data_list[0]['real'][i_freq_in_band]**2) /\ scipy.sum(data_list[0]['real'][i_freq_in_band]) else: norm_factor = 1.0 # Initialize arrays and apply normalization real_min = data_list[0]['real'] / norm_factor real_max = data_list[0]['real'] / norm_factor real_ave = 0 # Sums for RMS estimate real_rms_s2 = 0 real_rms_s = 0 for data in data_list: # Check the chan, if, pol and freq are the same. Take first as reference # If not quit if scipy.array_equal(chan_ref, data['chan']) and\ scipy.array_equal(if_ref, data['if']) and\ (pol_ref == data['pol']).all and\ scipy.array_equal(freq_ref, data['freq']): # Normalize if requested if normalize == True: # If empty, it means there is no overlap, normalize with all band if len(i_freq_in_band) == 0: if method == 'mean': norm_factor = scipy.mean(data['real']) elif method == 'median': norm_factor = scipy.median(data['real']) elif method == 'minrms': norm_factor = scipy.sum(data['real']**2) /\ scipy.sum(data['real']) # Or normalize using the channels on the passband else: if method == 'mean': norm_factor = scipy.mean(data['real'][i_freq_in_band]) elif method == 'median': norm_factor =\ scipy.median(data['real'][i_freq_in_band]) elif method == 'minrms': norm_factor =\ scipy.sum(data['real'][i_freq_in_band]**2) /\ scipy.sum(data['real'][i_freq_in_band]) ################ # FOR DEVELOPMENT #pyplot.clf() #pyplot.plot(data['freq'], data['real'], label='raw') #pyplot.axhline(1, color='k') #pyplot.axhline(norm_factor, color='r', label='normalization') #pyplot.plot(data['freq'], # data['real'] / norm_factor, # color='r', # label='normalized') #pyplot.axhline(1, color='k') #pyplot.legend(loc=0, frameon=False) #pyplot.ylim(0, 5) #pyplot.show() ################# else: norm_factor = 1.0 # Normalize real_norm_t = data['real'] / norm_factor real_min_t = real_norm_t real_max_t = real_norm_t real_ave += real_norm_t real_rms_s2 += real_norm_t**2 real_rms_s += real_norm_t # Update arrays real_min = scipy.minimum(real_min, real_min_t) real_max = scipy.maximum(real_max, real_max_t) else: print 'Not all scans have same observing parameters' return {} # Estimate the mean and RMS real_ave = real_ave / n_spectra real_rms =\ scipy.sqrt((real_rms_s2 - 2.0 * real_ave * real_rms_s +\ n_spectra * real_ave**2) / n_spectra) # Initialize data_env data_env = {'chan': chan_ref, 'if': if_ref, 'pol': pol_ref, 'freq': freq_ref, 'real_min': real_min, 'real_max': real_max, 'real_ave': real_ave, 'real_rms': real_rms } return data_env #------------------------------------------------------------------------------- def compute_spectra_envelope(data_list, normalize=False, method='mean'): """ Takes a list of scans with as given antenna and returns the spectral total power envelope. This envelope is the minimum, maximum total power observed for all the scans. The envelope is obtained for each frequency channel and is meant to give a sense of the range of observed total powers input ----- data_list A list of data for scans. Each element like the output of load_spectrum() normalize Normalize the spectrum by the median to get values around 1 method Chooses a normlaization method. Two are currently implemented, but more can be implemented later output ------ data_env Same structure as data dictionaty, but with max and min for real and imag notes ----- wmax, Nov 8, 2014 """ # Determine the number of IFs and its numbers in the file ifs_unique = scipy.unique(data_list[0]['if']) # Empty dictionary initialization data_env = None for if_i in ifs_unique: # Create a data_list with data for a single IF data_list_if = [] for i in xrange(len(data_list)): # Get elements for a given IF if_index = scipy.where(data_list[i]['if'] == if_i)[0] # Fill a temporal element with a single IF data_list_if_t = {} for key in data_list[i].keys(): data_list_if_t[key] = data_list[i][key][if_index] # Append to list with single IF data_list_if.append(data_list_if_t) # Compute envepole for this IF data_env_if = compute_spectra_envelope_if(data_list_if, normalize=normalize, method=method) # Append to other IFs or create if empty if data_env == None: data_env = data_env_if else: for key in data_env.keys(): data_env[key] = scipy.concatenate((data_env[key], data_env_if[key])) return data_env #------------------------------------------------------------------------------- def plot_total_power_envelope(data, ylim=None, filename='total_power_envelope.pdf', xyratio_panels=2.0, plot_type='envelope'): """ Plot total power output from AIPS input ----- data_env Dictionary with data envelope, output of compute_spectra_envelope(data_list) ylim Limit for y-axis. If none uses max and min from data filename Output plot file name xyratio_panels plot_type output ------ A file with 'filename' with the plot and another one with the data notes ----- wmax, Nov 8, 2014 """ # Determine the number of IFs in the file ifs_unique = scipy.unique(data['if']) n_ifs = len(ifs_unique) # Configure plot parameters pyplot.rcParams['xtick.labelsize'] = 12 pyplot.rcParams['ytick.labelsize'] = 12 pyplot.rcParams['figure.subplot.wspace'] = 0.0 # Find a common scale for the plots of not given if ylim == None: ymax = max(data['real_max']) ymin = min(data['real_min']) else: ymin, ymax = ylim fig = pyplot.figure(figsize=(xyratio_panels * 5.0 * n_ifs, 5.0)) # Sort the IFs by the first frequency min_freq = [] for if_i in ifs_unique: if_index = scipy.where(data['if'] == if_i)[0] min_freq.append(min(data['freq'][if_index])) ifs_unique_sorted = ifs_unique[scipy.argsort(min_freq)] # Loop over IFs for i, if_i in zip(range(1, n_ifs + 1), ifs_unique_sorted): # Find the data for a given IF if_index = scipy.where(data['if'] == if_i)[0] # Arrays to plot freq = data['freq'][if_index] real_min = data['real_min'][if_index] real_max = data['real_max'][if_index] real_ave = data['real_ave'][if_index] real_rms = data['real_rms'][if_index] # Plot pyplot.subplot(1, n_ifs, i) if plot_type == 'envelope': pyplot.fill_between(freq, real_min, real_max, facecolor='0.75', edgecolor='0.75') pyplot.plot(freq, real_ave, color='0.0') elif plot_type == 'rms': pyplot.fill_between(freq, real_ave - real_rms, real_ave + real_rms, facecolor='0.75', edgecolor='0.75') pyplot.plot(freq, real_ave, color='0.0') pyplot.ylim(ymin, ymax) pyplot.ticklabel_format(axis='x', style='plain', useOffset=True) pyplot.tick_params(axis='x', which='major') pyplot.xticks(scipy.linspace(min(freq), max(freq), 5)) # remove yticks for all plots except the first one if i > 1 and i < n_ifs: pyplot.tick_params(axis='y', which='both', labelleft='off', labelright='off') # Add right yticks for last plot elif i == n_ifs: pyplot.tick_params(axis='y', which='both', labelleft='off', labelright='on') pyplot.savefig('%s.%s' % (filename, 'pdf'), bbox_inches='tight') # Close figure to avoid filling up memory pyplot.close() # Save data to text file for later plotting on the web scipy.savetxt('%s.%s' % (filename, 'csv'), scipy.transpose((data['freq'], data['if'], data['real_min'], data['real_max'], data['real_ave'], data['real_rms'])), fmt=['%f', '%i', '%f', '%f', '%f', '%f'], delimiter=',', header='freq,if,real_min,real_max,real_ave,real_rms') # Restore rcParams pyplot.rcdefaults() return #------------------------------------------------------------------------------- def load_envelope_data(filename): """ Load the total power envelope data saved by plot_total_power_envelope() wmax, nov 14, 2014 """ freq, ifn, real_min, real_max, real_ave, real_rms =\ scipy.loadtxt(filename, delimiter=',', comments='#', unpack=True) data_env = {'freq': freq, 'if': ifn, 'real_min': real_min, 'real_max': real_max, 'real_ave': real_ave, 'real_rms': real_rms } return data_env #------------------------------------------------------------------------------- # Run #------------------------------------------------------------------------------- for name_prefix in name_prefixes: aips_filename = name_prefix print '--------------------------------------------------' print 'Working on %s' % (aips_filename) for antenna in antennas: print '' print 'Plotting total power for %s' % (antenna) for pol in stokes: # Find a clever way to get files for a given site given polarization listdir = os.listdir(aips_output) files = [filename for filename in listdir if\ filename.find('%s-%s-%s' % (aips_filename, antenna, pol)) >= 0 and\ os.path.splitext(filename)[-1] == '.txt'] files.sort() # Check that there is data for this antenna, if not go to next one if len(files) == 0: print 'No data found for %s in %s polarization' % (antenna, pol) continue # List of arrays with total power data for each scan data_list = [] for filename in files: data = load_spectrum(os.path.join(aips_output, filename)) data_list.append(data) if plot_single_scans == True: plot_total_power_single_scan(data, ylim=None, filename=os.path.join(plots_output, filename.replace('txt', 'png'))) data_env = compute_spectra_envelope(data_list, normalize=normalize, method=norm_method) plot_total_power_envelope(data_env, ylim=(0.0, 4.0), filename=os.path.join(plots_output, 'tpower-env-%s-%s-%s' %\ (aips_filename, antenna, pol)), plot_type=plot_type) # Plot the envelope for all the setups combined on the same site print '--------------------------------------------------' print 'Total power envelope for all setups combined\n' for antenna in antennas: print 'Getting data for %s\n' % (antenna) for pol in stokes: # Find a clever way to get files for a given site given polarization listdir = os.listdir(plots_output) files = [filename for filename in listdir if\ filename.find('%s-%s' % (antenna, pol)) >= 0 and\ filename.find('tpower-env-') >= 0 and\ os.path.splitext(filename)[-1] == '.csv'] files.sort() # Check that there is data for this antenna, if not go to next one if len(files) == 0: print 'No data found for %s in %s polarization' % (antenna, pol) continue # List of arrays with total power envelope data for each scan data_env_all_ifs = {} data_env_all_ifs['freq'] = scipy.array([]) data_env_all_ifs['if'] = scipy.array([]) data_env_all_ifs['real_min'] = scipy.array([]) data_env_all_ifs['real_max'] = scipy.array([]) data_env_all_ifs['real_ave'] = scipy.array([]) data_env_all_ifs['real_rms'] = scipy.array([]) # IF index maximum n_if_last = 1 for filename in files: data_env = load_envelope_data(os.path.join(plots_output, filename)) # Fix the IF number ifs = scipy.unique(data_env['if']) index_ifs = [] # Get index arrays for all IFs before modifying # this is to avoid confusion for if_n in ifs: index_ifs.append(scipy.where(data_env['if'] == if_n)[0]) # Change the IF numbers starting from n_if_last for index_if in index_ifs: data_env['if'][index_if] = n_if_last n_if_last += 1 # Add data to IF envelope all ifs data_env_all_ifs['freq'] =\ scipy.concatenate((data_env_all_ifs['freq'], data_env['freq'])) data_env_all_ifs['if'] =\ scipy.concatenate((data_env_all_ifs['if'], data_env['if'])) data_env_all_ifs['real_min'] =\ scipy.concatenate((data_env_all_ifs['real_min'], data_env['real_min'])) data_env_all_ifs['real_max'] =\ scipy.concatenate((data_env_all_ifs['real_max'], data_env['real_max'])) data_env_all_ifs['real_ave'] =\ scipy.concatenate((data_env_all_ifs['real_ave'], data_env['real_ave'])) data_env_all_ifs['real_rms'] =\ scipy.concatenate((data_env_all_ifs['real_rms'], data_env['real_rms'])) # Plot envelope for an antenna and polarization plot_total_power_envelope(data_env_all_ifs, ylim=(0.0, 3.0), filename=os.path.join(plots_output, 'tpower-all-env-ifs-%s-%s' %\ (antenna, pol)), xyratio_panels=1.0, plot_type=plot_type)