%matplotlib inline import pandas as pd import scipy.signal as sig import scipy.stats as stats import matplotlib.pyplot as plt import os import numpy as np from datetime import datetime def iterate_months(path, title=''): """ Changes to the specificed path assumed to be a buoy year directory and loads the displacement DataFrame for each month plotting the spectra for each method using Welch's method """ print "path", path os.chdir(path) months_data = os.listdir('.') max_psd = {} for month_dir in months: month_path = path + month_dir try: welch_df_path = month_path + '/month_welch_df' if os.path.isfile(welch_df_path): welch_df = pd.read_pickle(welch_df_path) else: dataframe = read_dataframe(month_path) os.chdir(month_path) welch_df = create_welch_df(dataframe) welch_df.to_pickle(welch_df_path) plt, data, max_psd = calc_welch(welch_df, psd_dic = max_psd, save = True, line_name = month_dir ) except: print "No Data", month_path os.chdir(path) plt.title(title) return max_psd def create_displacements_df(path): displacements_df = read_dataframe(path) try: os.chdir(path) except: print "No directory", path return displacements_df def create_welch_df(dataframe): """Takes a single column displacement DataFrame in metres""" x,y = sig.welch(dataframe, 1.28, noverlap=None) dic = {'x':x,'y':y} welch_df = pd.DataFrame(dic) return welch_df def calc_welch(welch_df, psd_dic = {}, save = False, line_name = ''): """ Given a dataframe will run welch's methods on the values, printing out the maxmimum power spectral density and the frequency it occurs at, plots the output and returns the data """ psd_dic[line_name] = [welch_df.x.max(), welch_df.y.ix[welch_df.y==welch_df.y.max()]] plt.plot(welch_df.x,welch_df.y, label=line_name) plt.ylabel('S(f) m^2/Hz') plt.xlabel('Frequency (Hz)') plt.legend() return plt, welch_df, psd_dic def read_dataframe(path, start = None, end = None, filter_signal_error=True): """ Given a path to a folder containing a DataFrame file called raw_plus_std created by hebtools returns the a subset of the heave column converted from centimetres to metres """ try: displacement_df = pd.read_pickle(path + '/' + 'raw_plus_std') if filter_signal_error: displacement_df = displacement_df[displacement_df.signal_error == 0] return displacement_df.heave[start:end]/100 except: return pd.Series() def read_raw(raw_file_path, seg_overlap = None, filter = False): """ Given a path to a raw file reads it using Pandas and returns a DataFrame object optionally filtered. """ raw_df = pd.read_csv(raw_file_path, names=['signal_quality', 'heave', 'north', 'west']) if filter: raw_df = raw_df[raw_df.signal_quality<1] return raw_df def plot_raw_file(raw_file_path, seg_overlap = None, filter = False, column = 'heave'): """ Takes path to a raw file as input with option to plot filter/unfiltered or any of the directions returns output of calc_welch the plt object and data in chart. """ raw_df = read_raw(raw_file_path, seg_overlap, filter) print raw_df.describe() print "Skewness", stats.skew(raw_df[column]) print "Kurtosis", stats.kurtosis(raw_df[column]) welch_df = create_welch_df(raw_df[column]/100) return calc_welch(welch_df, line_name = raw_file_path + ' Filter ' + str(filter) ) if os.name == 'posix': path_prefix = '/media/Data/Datawell/' c_prefix = '/media/System/' else: path_prefix = 'D:/Datawell/' c_prefix = 'C:/' months = ['January','February','March','April','May','June','July','August','September','October','November','December'] buoys = ['Siadar_HebMarine1', 'Bragar_HebMarine2', 'Roag_Wavegen'] years = [2011, 2012, 2013] def iterate_buoys_months(path_prefix, buoys, year): title_suffix = ' Buoy monthly Power Spectral Density ' for buoy in buoys: fig = plt.figure(figsize=(20,10), dpi=100) year_path = path_prefix + buoy + '/' + year + '/' title = buoy + title_suffix + year max_psd = iterate_months(year_path, title) iterate_buoys_months(path_prefix, buoys, '2012') def iterate_buoys_one_month(buoys, year, month): for buoy in buoys: try: path = path_prefix + buoy + '/' + year + '/' + month os.chdir(path) #dataframe = read_dataframe(path) welch_df = pd.read_pickle('month_welch_df') #welch_df = create_welch_df(dataframe) plt, january_data, max_psd = calc_welch(welch_df, line_name = path) except: print "No data" for month in months: fig = plt.figure(figsize=(20,10), dpi=100) iterate_buoys_one_month(buoys, '2012', month) fig = plt.figure(figsize=(20,10), dpi=100) plt, data, psd_dic = plot_raw_file(path_prefix + 'Bragar_HebMarine2/2012/February/Bragar_HebMarine2}2012-02-01T01h00Z.raw') plt.title('First half hour of February 2012 Bragar buoy') from IPython.display import Image Image(filename='D:/image.png') fig = plt.figure(figsize=(20,10), dpi=100) plt, heave_spectra_unfiltered_2013_02_04_13_00, psd_dic = plot_raw_file(path_prefix + 'Bragar_HebMarine2/2013/February/Bragar_HebMarine2}2013-02-04T13h00Z.raw', filter = False) plt, data_filtered_02_04_13_00, psd_dic = plot_raw_file(path_prefix + 'Bragar_HebMarine2/2013/February/Bragar_HebMarine2}2013-02-04T13h00Z.raw', filter = True) plt.title('Comparing filtered and unfiltered spectra at the peak of the 4th Feb 2013 Storm') plt.legend() # Find the raw file with the worst Signal Quality in February 2013 to see the effect of filtering os.chdir(path_prefix + 'Bragar_HebMarine2/2013/February/') raw_displacements = pd.read_pickle('raw_plus_std') grouped_raw_displacements = raw_displacements.groupby('file_name') sig_qual_files_sum = grouped_raw_displacements.sig_qual.sum() sig_qual_files_sum[sig_qual_files_sum.values == sig_qual_files_sum.max()] fig = plt.figure(figsize=(20,10), dpi=100) plt, filtered_2013_02_04_00_30, psd_filtered = plot_raw_file(path_prefix + 'Bragar_HebMarine2/2013/February/Bragar_HebMarine2}2013-02-04T00h30Z.raw', filter = True) plt, unfiltered_2013_02_04_00_30, psd_unfiltered = plot_raw_file(path_prefix + 'Bragar_HebMarine2/2013/February/Bragar_HebMarine2}2013-02-04T00h30Z.raw', filter = False) plt.title('File with worst Signal Quality February 2013 Bragar buoy') # There is a buoy computed spectra dated 2013/02/04 17:33 this computes a similar spectra for the same # period from the DataFrame produced by hebtools path = 'Bragar_HebMarine2/2013/February' os.chdir(path_prefix + path) feb_raw = pd.read_pickle('raw_plus_std') subset = feb_raw[datetime(2013,02,03,23,47,33):datetime(2013,02,04,00,17,33)] fig = plt.figure(figsize=(20,10), dpi=100) welch_df = create_welch_df(subset.heave/100) plt, data, psd = calc_welch(welch_df, line_name = path) plt.title("Dataframe computed spectra 3rd/4th Feb 23:57-00:17 to match Buoy spectra") # Example W@ves21 Spectra plot from IPython.display import Image Image(filename='D:\Profiles\le12jm\My Documents\My Pictures\plot_02_04_00_30.png') def dat_to_x_y_list(data_list_strings): """ Function which takes the strings from an open dat file and outputs a nested list of x and y values for plotting dat files can be output by right clicking on spectra plot in W@ves21 and selecting 'Save Chart Data' and changing the default file from OC to dat """ data = ' '.join(data_list_strings).split(' ') data_values = [] values = [] for i in data: if i == 'Y': data_values.append(values) values = [] if '.' in i: values.append(i) data_values.append(values) return data_values # Comparing a Datawell buoy and pc spectra taken over a similar period ( PC 00:00-00:30, Buoy 23:47-00:17 ) # Also comparing with a spectra produced from a hebtools dataframe for the period of the buoy spectra os.chdir(path_prefix + "Bragar_HebMarine2/2013/February/") pc_spectra = dat_to_x_y_list(open('20140904110659.dat')) buoy_spectra = dat_to_x_y_list(open('2013_02_04_00_47_33_buoy_spectra.dat')) fig = plt.figure(figsize=(20,10), dpi=100) plt.plot(pc_spectra[0],pc_spectra[1], label='PC Spectra 00:30:00') plt.ylabel('S(f) m^2/Hz') plt.xlabel('Frequency (Hz)') plt.title('Plot of Spectra 4th February 2013 - 00:30') plt.plot(buoy_spectra[0],buoy_spectra[1], label='Buoy Spectra 00:17:33') welch_df = create_welch_df(subset.heave/100) plt, data, psd = calc_welch(welch_df, line_name = 'Dataframe Spectra 23:47-00:17') plt.legend() fig = plt.figure(figsize=(20,10), dpi=100) plt, west_filtered_2013_02_04_00_30, psd = plot_raw_file(path_prefix + 'Bragar_HebMarine2/2013/February/Bragar_HebMarine2}2013-02-04T00h30Z.raw', filter = True, column='west') plt.title('West Displacement - Peak of Storm February 2013 Bragar buoy') fig = plt.figure(figsize=(20,10), dpi=100) plt.plot(filtered_2013_02_04_00_30.x, filtered_2013_02_04_00_30.y, label = "Dataframe Filtered") plt.plot(unfiltered_2013_02_04_00_30.x, unfiltered_2013_02_04_00_30.y, label = "Dataframe Unfiltered") plt.plot(pc_spectra[0], pc_spectra[1], label = "Datawell PC Spectra") plt.plot(west_filtered_2013_02_04_00_30.x, west_filtered_2013_02_04_00_30.y, label = "Dataframe West Displacement Filtered") plt.title('Spectra Comparison Bragar February 4th 2014 00:00-00:30') plt.legend() from IPython.display import Image Image(filename='D:\Profiles\le12jm\My Documents\My Pictures\wave_stats_02_04_00_30.png') fig = plt.figure(figsize=(20,10), dpi=100) plt, data, psd = plot_raw_file(path_prefix + 'Bragar_HebMarine2/2013/February/Bragar_HebMarine2}2013-02-04T00h30Z.raw', filter = True, column='north') plt.title('North Displacement - Peak of Storm February 2013 Bragar buoy') # Percentage of records returned with an upper filter bound of 1. 2120/2303.0 from IPython.display import Image Image(filename='D:\Profiles\le12jm\My Documents\My Pictures\Spectra_Summary_00_30_W@ves21.png')