This notebook will load data, perform a tidal analyis, compare with observations, plot the results, and save the analysis in a spreadsheet. Six Tidal Constituents: M2, K1, O1, S2, P1 and N2 are considered. # imports %matplotlib inline import matplotlib.pylab as plt import numpy as np import netCDF4 as NC from scipy.optimize import curve_fit from salishsea_tools import tidetools from salishsea_tools import viz_tools from salishsea_tools import bathy_tools import collections import pandas as pd import csv import math from __future__ import division # pathname for data - all of the tide runs are stored in this directory path = '/data/nsoontie/MEOPAR/SalishSea/results/tides/' #path = '../../myResults/' #the run we want to analyze runname = 'RC5_wO1S2P1N2' #joining the two string together name = path +runname +'/' print name # grid grid = NC.Dataset('../../nemo-forcing/grid/bathy_meter_SalishSea2.nc') bathy, X, Y = tidetools.get_bathy_data(grid) filename = '/data/nsoontie/MEOPAR/analysis/compare_tides/obs_tidal_wlev_const_all.csv' filename = '../compare_tides/obs_tidal_wlev_const_all.csv' harm_obs = pd.read_csv(filename,sep=';',header=0) harm_obs = harm_obs.rename(columns={'Site': 'site', 'Lat': 'lat', 'Lon': 'lon', 'M2 amp': 'M2_amp', 'M2 phase (deg UT)': 'M2_pha', 'K1 amp': 'K1_amp', 'K1 phase (deg UT)': 'K1_pha'}) print harm_obs filename = 'other_constituents.csv' harm_other = pd.read_csv(filename,sep=',',header=0) harm_other = harm_other.rename(columns={'Site': 'site', 'Lat': 'lat', 'Lon': 'lon', 'O1 amp': 'O1_amp', 'O1 phase (deg UT)': 'O1_pha', 'P1 amp': 'P1_amp', 'P1 phase (deg UT)': 'P1_pha', 'Q1 amp': 'Q1_amp', 'Q1 phase (deg UT)': 'Q1_pha', 'S2 amp': 'S2_amp', 'S2 phase (deg UT)': 'S2_pha', 'N2 amp': 'N2_amp', 'N2 phase (deg UT)': 'N2_pha', 'K2 amp': 'K2_amp', 'K2 phase (deg UT)': 'K2_pha'}) print harm_other stations = ['PortRenfrew','SheringhamPoint','PedderBay', 'Esquimalt', 'Victoria','CloverPoint','FinnertyCove', 'FulfordHarbour', 'TumboChannel','PatosIsland','WhalerBay', 'Tsawwassen', 'Sandheads', 'PointGrey','PointAtkinson','GibsonsLanding', 'WinchelseaIs', 'HalfmoonBay','IrvinesLanding','PowellRiver', 'LittleRiver', 'Lund', 'TwinIslets','CampbellRiver','MaudeIslandE', 'NympheCove', 'SeymourNarrows','BrownBay','ChathamPoint','KelseyBay','YorkeIsland'] numsta=len(stations) #again with spaces because the text file likes that stations_obs = ['Port Renfrew','Sheringham Point','Pedder Bay', 'Esquimalt', 'Victoria','Clover Point','Finnerty Cove', 'Fulford Harbour', 'Tumbo Channel','Patos Island','Whaler Bay', 'Tsawwassen', 'Sandheads', 'Point Grey','Point Atkinson','Gibsons Landing', 'Winchelsea', 'Halfmoon Bay','Irvines Landing','Powell River', 'Little River', 'Lund', 'Twin Islets','Campbell River','Maude Island E', 'Nymphe Cove', 'Seymour Narrows','Brown Bay','Chatham Point','Kelsey Bay','Yorke Island'] fig,ax=plt.subplots(1, 1, figsize=(8, 10)) ax.pcolormesh(X,Y,bathy,cmap='winter_r') for stn in range(numsta): location = stations_obs[stn] lon=-harm_obs.lon[harm_obs.site==location] lat=harm_obs.lat[harm_obs.site==location] ax.plot(lon,lat,'.k',label=location) ax.annotate(stn, xy = (lon,lat), xytext = (5,5),ha = 'right', va = 'bottom', textcoords = 'offset points') print stn, location ax.axis([-126.1,-122,47,51]) #constants and fitting # M2 M2freq = 28.984106 # degrees per hour M2freq = M2freq*np.pi/180. # radians per hour #K1 K1freq = 15.041069*np.pi/180. #O1 O1freq = 13.943036*np.pi/180. #S2 S2freq = 30.000002*np.pi/180. #P1 P1freq = 14.958932*np.pi/180. #N2 N2freq = 28.439730*np.pi/180. # initial phase calculation # our start is currently Oct 26, 2002 # data for phase output from bdytides.F90; found in ocean.output K1ft = 1.050578 K1uvt = 296.314842 M2ft = 0.987843 M2uvt = 245.888564 O1ft = 1.081364 O1uvt = 312.950020 S2ft = 1.0 S2uvt = 0.0 P1ft = 1.0 P1uvt = 55.79460 N2ft = 0.98784 N2uvt = 353.570277 # function for fit def double(x, M2amp, M2pha, K1amp, K1pha): return (M2amp*np.cos(M2freq*x-M2pha*np.pi/180.)+ K1amp*np.cos(K1freq*x-K1pha*np.pi/180.)) # function for fitting 3 frequencies def triple(x, M2amp, M2pha, K1amp, K1pha, O1amp, O1pha): return (M2amp*np.cos(M2freq*x-M2pha*np.pi/180.)+ K1amp*np.cos(K1freq*x-K1pha*np.pi/180.)+ O1amp*np.cos(O1freq*x-O1pha*np.pi/180.)) # function for fitting 4 frequencies def quad(x, M2amp, M2pha, K1amp, K1pha, O1amp, O1pha, S2amp, S2pha): return (M2amp*np.cos(M2freq*x-M2pha*np.pi/180.)+ K1amp*np.cos(K1freq*x-K1pha*np.pi/180.)+ O1amp*np.cos(O1freq*x-O1pha*np.pi/180.)+ S2amp*np.cos(S2freq*x-S2pha*np.pi/180.)) # function for fitting 6 frequencies def sextuple(x, M2amp, M2pha, K1amp, K1pha, O1amp, O1pha, S2amp, S2pha, P1amp, P1pha, N2amp, N2pha): return (M2amp*np.cos(M2freq*x-M2pha*np.pi/180.)+ K1amp*np.cos(K1freq*x-K1pha*np.pi/180.)+ O1amp*np.cos(O1freq*x-O1pha*np.pi/180.)+ S2amp*np.cos(S2freq*x-S2pha*np.pi/180.)+ P1amp*np.cos(P1freq*x-P1pha*np.pi/180.)+ N2amp*np.cos(N2freq*x-N2pha*np.pi/180.)) #allocate space for our arrays M2_amp=[]; M2_pha=[]; K1_amp=[]; K1_pha=[] O1_amp=[]; O1_pha=[]; S2_amp=[]; S2_pha=[] P1_amp=[]; P1_pha=[]; N2_amp=[]; N2_pha=[] M2_amp_obs= np.zeros(numsta); M2_pha_obs=np.zeros(numsta) K1_amp_obs=np.zeros(numsta); K1_pha_obs=np.zeros(numsta) O1_amp_obs=np.zeros(numsta); O1_pha_obs=np.zeros(numsta) S2_amp_obs=np.zeros(numsta); S2_pha_obs=np.zeros(numsta) P1_amp_obs=np.zeros(numsta); P1_pha_obs=np.zeros(numsta) N2_amp_obs=np.zeros(numsta); N2_pha_obs=np.zeros(numsta) ts = 240 for stn in range(numsta): fT1 = NC.Dataset(name+stations[stn]+'.nc','r') time = fT1.variables["time_counter"][:]/3600. # want hours not seconds ssh = fT1.variables["sossheig"][:,0,0] fitted, cov = curve_fit(sextuple,time[ts:],ssh[ts:]) if fitted[0] < 0: fitted[0] = -fitted[0] fitted[1] = fitted[1]+180 M2_amp.append(fitted[0]*M2ft) pha = fitted[1]+M2uvt if pha > 360: pha=pha-360 M2_pha.append(pha) K1_amp.append(fitted[2]*K1ft) pha= fitted[3]+K1uvt if pha > 360: pha=pha-360 K1_pha.append(pha) if fitted[4] < 0: fitted[4] = -fitted[4] fitted[5] = fitted[5]+180 O1_amp.append(fitted[4]*O1ft) pha= fitted[5]+O1uvt if pha > 360: pha=pha-360 O1_pha.append(pha) if fitted[6] < 0: fitted[6] = -fitted[6] fitted[7] = fitted[7]+180 S2_amp.append(fitted[6]*S2ft) pha= fitted[7]+S2uvt if pha > 360: pha=pha-360 S2_pha.append(pha) if fitted[8] < 0: fitted[8] = -fitted[8] fitted[9] = fitted[9]+180 P1_amp.append(fitted[8]*P1ft) pha= fitted[9]+P1uvt if pha > 360: pha=pha-360 P1_pha.append(pha) if fitted[10] < 0: fitted[10] = -fitted[10] fitted[11] = fitted[11]+180 N2_amp.append(fitted[10]*N2ft) pha= fitted[11]+N2uvt if pha > 360: pha=pha-360 N2_pha.append(pha) #now the observations location=stations_obs[stn] M2_amp_obs[stn]=harm_obs.M2_amp[harm_obs.site==location]/100 M2_pha_obs[stn]=harm_obs.M2_pha[harm_obs.site==location] K1_amp_obs[stn]=harm_obs.K1_amp[harm_obs.site==location]/100 K1_pha_obs[stn]=harm_obs.K1_pha[harm_obs.site==location] #O1/S2/P1/N2 are in the other file if (harm_other.site==location).any(): O1_amp_obs[stn]=harm_other.O1_amp[harm_other.site==location]/100 O1_pha_obs[stn]=harm_other.O1_pha[harm_other.site==location] S2_amp_obs[stn]=harm_other.S2_amp[harm_other.site==location]/100 S2_pha_obs[stn]=harm_other.S2_pha[harm_other.site==location] P1_amp_obs[stn]=harm_other.P1_amp[harm_other.site==location]/100 P1_pha_obs[stn]=harm_other.P1_pha[harm_other.site==location] N2_amp_obs[stn]=harm_other.N2_amp[harm_other.site==location]/100 N2_pha_obs[stn]=harm_other.N2_pha[harm_other.site==location] #Mask the arrays so that we can do statistics without the 0's throwing thigns off. O1_amp_obs =np.ma.masked_values(O1_amp_obs, 0) O1_pha_obs =np.ma.masked_values(O1_pha_obs, 0) S2_amp_obs =np.ma.masked_values(S2_amp_obs, 0) S2_pha_obs =np.ma.masked_values(S2_pha_obs, 0) P1_amp_obs =np.ma.masked_values(P1_amp_obs, 0) P1_pha_obs =np.ma.masked_values(P1_pha_obs, 0) N2_amp_obs =np.ma.masked_values(N2_amp_obs, 0) N2_pha_obs =np.ma.masked_values(N2_pha_obs, 0) #Plotting M2 labels=['JdF/Islands','SoG','North'] split1=8; split2=20 fig=tidetools.plot_scatter_pha_amp(M2_amp,M2_amp_obs,M2_pha,M2_pha_obs,'M2',figsize=(14,6), split1=split1,split2=split2, labels=labels) ax_amp,ax_pha = fig.axes min_value, max_value = ax_amp.set_xlim(0, 1.2) ax_amp.plot([min_value, max_value], [min_value, max_value], color='red',lw=2) min_value, max_value = ax_pha.set_xlim(0, 360) ax_pha.plot([min_value, max_value], [min_value, max_value], color='red',lw=2) #Plotting - K1 fig=tidetools.plot_scatter_pha_amp(K1_amp,K1_amp_obs,K1_pha,K1_pha_obs,'K1',figsize=(14,6), split1=split1, split2=split2, labels=labels) ax_amp,ax_pha = fig.axes min_value, max_value = ax_amp.set_xlim(0, 1.2) ax_amp.plot([min_value, max_value], [min_value, max_value], color='red',lw=2) min_value, max_value = ax_pha.set_xlim(0, 360) ax_pha.plot([min_value, max_value], [min_value, max_value], color='red',lw=2) #Plotting - O1 fig=tidetools.plot_scatter_pha_amp(O1_amp,O1_amp_obs,O1_pha,O1_pha_obs,'O1',figsize=(14,6), split1=split1, split2=split2, labels=labels) ax_amp,ax_pha = fig.axes min_value, max_value = ax_amp.set_xlim(0, 1.2) ax_amp.plot([min_value, max_value], [min_value, max_value], color='red',lw=2) min_value, max_value = ax_pha.set_xlim(0, 360) ax_pha.plot([min_value, max_value], [min_value, max_value], color='red',lw=2) #Plotting - S2 fig=tidetools.plot_scatter_pha_amp(S2_amp,S2_amp_obs,S2_pha,S2_pha_obs,'S2',figsize=(14,6), split1=split1, split2=split2, labels=labels) ax_amp,ax_pha = fig.axes min_value, max_value = ax_amp.set_xlim(0, 1.2) ax_amp.plot([min_value, max_value], [min_value, max_value], color='red',lw=2) min_value, max_value = ax_pha.set_xlim(0, 360) ax_pha.plot([min_value, max_value], [min_value, max_value], color='red',lw=2) #Plotting - P1 fig=tidetools.plot_scatter_pha_amp(P1_amp,P1_amp_obs,P1_pha,P1_pha_obs,'P1',figsize=(14,6), split1=split1, split2=split2, labels=labels) ax_amp,ax_pha = fig.axes min_value, max_value = ax_amp.set_xlim(0, 1.2) ax_amp.plot([min_value, max_value], [min_value, max_value], color='red',lw=2) min_value, max_value = ax_pha.set_xlim(0, 360) ax_pha.plot([min_value, max_value], [min_value, max_value], color='red',lw=2) #Plotting - N2 fig=tidetools.plot_scatter_pha_amp(N2_amp,N2_amp_obs,N2_pha,N2_pha_obs,'N2',figsize=(14,6), split1=split1, split2=split2, labels=labels) ax_amp,ax_pha = fig.axes min_value, max_value = ax_amp.set_xlim(0, 1.2) ax_amp.plot([min_value, max_value], [min_value, max_value], color='red',lw=2) min_value, max_value = ax_pha.set_xlim(0, 360) ax_pha.plot([min_value, max_value], [min_value, max_value], color='red',lw=2) def mean(diff): return np.mean(abs(diff)) def rms(diff): return np.sqrt(np.mean(diff**2)) def complex_diff(Ao,go,Am,gm): #calculates complex differences between observations and model #Ao, go - amplitude and phase from observations #Am, gm - amplitude and phase from model D = np.sqrt((Ao*np.cos(np.pi*go/180)-Am*np.cos(np.pi*gm/180))**2 + (Ao*np.sin(np.pi*go/180)-Am*np.sin(np.pi*gm/180))**2) return D #R R_M2 = M2_amp/M2_amp_obs R_K1 = K1_amp/K1_amp_obs #delta phi (adjust so between -180, 180) Dphi_M2=M2_pha-M2_pha_obs; Dphi_M2 = Dphi_M2 -360*(Dphi_M2>180) + 360*(Dphi_M2<-180) Dphi_K1=K1_pha-K1_pha_obs Dphi_K1 = Dphi_K1 -360*(Dphi_K1>180) + 360*(Dphi_K1<-180) #Complex differences D_M2= complex_diff(np.array(M2_amp_obs),np.array(M2_pha_obs), np.array(M2_amp),np.array(M2_pha)) D_K1= complex_diff(np.array(K1_amp_obs),np.array(K1_pha_obs), np.array(K1_amp),np.array(K1_pha)) D_O1= complex_diff(np.ma.array(O1_amp_obs),np.ma.array(O1_pha_obs), np.ma.array(O1_amp),np.ma.array(O1_pha)) D_S2= complex_diff(np.ma.array(S2_amp_obs),np.ma.array(S2_pha_obs), np.ma.array(S2_amp),np.ma.array(S2_pha)) D_P1= complex_diff(np.ma.array(P1_amp_obs),np.ma.array(P1_pha_obs), np.ma.array(P1_amp),np.ma.array(P1_pha)) D_N2= complex_diff(np.ma.array(N2_amp_obs),np.ma.array(N2_pha_obs), np.ma.array(N2_amp),np.ma.array(N2_pha)) outfile = runname+'.csv' with open(outfile, 'wb') as csvfile: writer = csv.writer(csvfile, delimiter=',') writer.writerow([ 'Station Name', 'R (M2)', 'Delta phi (M2)', 'D (M2)', 'R (K1)', 'Delta phi (K1)', 'D (K1)' ]) for stn in range(numsta): location = stations_obs[stn] writer.writerow([stations_obs[stn], R_M2[stn], Dphi_M2[stn], D_M2[stn], R_K1[stn], Dphi_K1[stn], D_K1[stn]]) #write averages and rms writer.writerow(['Mean Difference', mean(M2_amp-M2_amp_obs),mean(Dphi_M2),mean(D_M2), mean(K1_amp-K1_amp_obs),mean(Dphi_K1),mean(D_K1)]) writer.writerow(['RMS Difference', rms(M2_amp-M2_amp_obs),rms(Dphi_M2),rms(D_M2), rms(K1_amp-K1_amp_obs),rms(Dphi_K1),rms(D_K1)]) #without the north writer.writerow(['Mean Difference no North no PR', mean(M2_amp[1:split2]-M2_amp_obs[1:split2]),mean(Dphi_M2[1:split2]),mean(D_M2[1:split2]), mean(K1_amp[1:split2]-K1_amp_obs[1:split2]),mean(Dphi_K1[1:split2]),mean(D_K1[1:split2])]) writer.writerow(['RMS Difference no North no PR', rms(M2_amp[1:split2]-M2_amp_obs[1:split2]),rms(Dphi_M2[1:split2]),rms(D_M2[1:split2]), rms(K1_amp[1:split2]-K1_amp_obs[1:split2]),rms(Dphi_K1[1:split2]),rms(D_K1[1:split2])]) plt.figure(figsize=(20,12)) plt.subplot(3,2,1) plt.plot(M2_amp, '-bo', label = 'model') plt.plot(M2_amp_obs, 'r-o', label = 'observation') plt.title('M2 Amplitude') plt.legend( loc='upper left' ) plt.subplot(3,2,2) plt.plot(K1_amp, '-bo', label = 'model') plt.plot(K1_amp_obs, 'r-o', label = 'observation') plt.title('K1 Amplitude') plt.subplot(3,2,3) # use the un-wrap function to plot the M2 phase more smoothly pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha)*np.pi/180.) plt.plot(pha_uw, '-bo', label = 'model') pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha_obs)*np.pi/180.) plt.plot(pha_uw, 'r-o', label = 'observation') plt.title('M2 Phase') plt.subplot(3,2,4) pha_uw = 180./np.pi * np.unwrap(np.array(K1_pha)*np.pi/180.) plt.plot(pha_uw, '-bo', label = 'model') plt.plot(K1_pha_obs, 'r-o', label = 'observation') plt.title('K1 Phase') plt.subplot(3,2,5) plt.plot(D_M2, '-mo', label = 'M2') plt.plot(D_K1, '-go', label = 'K1') plt.plot((0,30),(0.05,0.05),'k') plt.plot((0,30),(0.10,0.10),'r') plt.title('D error') plt.legend( loc='upper left' ) cmap = plt.get_cmap('PuBu') cmap.set_bad('burlywood') fig,axs=plt.subplots(3, 2, figsize=(8,15)) constituent = ('M2', 'K1', 'O1', 'S2', 'P1', 'N2') error_D = (D_M2, D_K1, D_O1, D_S2, D_P1, D_N2) for row in range(3): for ax, error_D1, const in zip(axs[row], error_D[row*2:row*2+2], constituent[row*2:row*2+2]): ax.pcolormesh(X,Y,bathy,cmap='PuBu') for stn in range(numsta): location = stations_obs[stn] lon=-harm_obs.lon[harm_obs.site==location] lat=harm_obs.lat[harm_obs.site==location] if error_D1 [stn] <= 0.05: ax.plot(lon,lat,'og',label=location,markersize=10,markeredgecolor='g') if error_D1 [stn] > 0.1: ax.plot(lon,lat,'or',label=location,markersize=10,markeredgecolor='r') if 0.1 >= error_D1[stn] > 0.05: ax.plot(lon,lat,'oy',label=location,markersize=10,markeredgecolor='y') ax.annotate(stn, xy = (lon,lat), xytext = (5,5),ha = 'right', va = 'bottom', textcoords = 'offset points') ax.set_title(const) ax.axis([-126.1,-122,47,51]) fig, axs = plt.subplots(4,2,figsize=(10,10)) axs[0,0].plot(np.array(O1_amp)/np.array(K1_amp), '-bo', label = 'model') axs[0,0].plot((0,28),(0.553,0.553), 'r-', label = 'observation') axs[0,0].set_title('O1/K1 Amplitude') axs[0,1].plot(np.array(O1_pha)-np.array(K1_pha), '-bo', label = 'model') axs[0,1].plot((0,28),(-23.1,-23.1), 'r-', label = 'observation') axs[0,1].set_title('O1-K1 Phase') axs[1,0].plot(np.array(S2_amp)/np.array(M2_amp), '-bo', label = 'model') axs[1,0].plot((0,28),(0.250,0.250), 'r-', label = 'observation') axs[1,0].set_title('S2/M2 Amplitude') pha_uw = 180./np.pi * np.unwrap((np.array(S2_pha)-np.array(M2_pha))*np.pi/180.) axs[1,1].plot(pha_uw, '-bo', label = 'model') axs[1,1].plot((0,28),( 28.7, 28.7), 'r-', label = 'observation') axs[1,1].set_title('S2-M2 Phase') axs[2,0].plot(np.array(P1_amp)/np.array(K1_amp), '-bo', label = 'model') axs[2,0].plot((0,28),(0.31,0.31), 'r-', label = 'observation') axs[2,0].set_title('P1/K1 Amplitude') pha_uw = 180./np.pi * np.unwrap((np.array(P1_pha)-np.array(K1_pha))*np.pi/180.) axs[2,1].plot(pha_uw, '-bo', label = 'model') axs[2,1].plot((0,28),(-2.5,-2.5), 'r-', label = 'observation') axs[2,1].set_title('P1-K1 Phase') axs[3,0].plot(np.array(N2_amp)/np.array(M2_amp), '-bo', label = 'model') axs[3,0].plot((0,28),(0.213,0.213), 'r-', label = 'observation') axs[3,0].set_title('N2/M2 Amplitude') pha_uw = 180./np.pi * np.unwrap((np.array(N2_pha)-np.array(M2_pha))*np.pi/180.) axs[3,1].plot(pha_uw, '-bo', label = 'model') axs[3,1].plot((0,28),(-27.2, -27.2), 'r-', label = 'observation') axs[3,1].set_title('N2-M2 Phase') #allocate space for our arrays M2_amp=np.zeros((numsta,3)); M2_pha=np.zeros((numsta,3)) K1_amp=np.zeros((numsta,3)); K1_pha=np.zeros((numsta,3)) O1_amp=np.zeros((numsta,3)); O1_pha=np.zeros((numsta,3)) S2_amp=np.zeros((numsta,3)); S2_pha=np.zeros((numsta,3)) P1_amp=np.zeros((numsta,3)); P1_pha=np.zeros((numsta,3)) N2_amp=np.zeros((numsta,3)); N2_pha=np.zeros((numsta,3)) for it,ts in zip(range(3),(240,300,360)): for stn in range(numsta): fT1 = NC.Dataset(name+stations[stn]+'.nc','r') time = fT1.variables["time_counter"][:]/3600. # want hours not seconds ssh = fT1.variables["sossheig"][:,0,0] fitted, cov = curve_fit(sextuple,time[ts:],ssh[ts:]) if fitted[0] < 0: fitted[0] = -fitted[0] fitted[1] = fitted[1]+180 M2_amp[stn,it] = fitted[0]*M2ft pha = fitted[1]+M2uvt if pha > 360: pha=pha-360 M2_pha[stn,it] = pha K1_amp[stn,it] =fitted[2]*K1ft pha= fitted[3]+K1uvt if pha > 360: pha=pha-360 K1_pha[stn,it]= pha if fitted[4] < 0: fitted[4] = -fitted[4] fitted[5] = fitted[5]+180 O1_amp[stn,it] =fitted[4]*O1ft pha= fitted[5]+O1uvt if pha > 360: pha=pha-360 O1_pha[stn,it]= pha if fitted[6] < 0: fitted[6] = -fitted[6] fitted[7] = fitted[7]+180 S2_amp[stn,it] =fitted[6]*S2ft pha= fitted[7]+S2uvt if pha > 360: pha=pha-360 S2_pha[stn,it]= pha if fitted[8] < 0: fitted[8] = -fitted[8] fitted[9] = fitted[9]+180 P1_amp[stn,it] = fitted[8]*P1ft pha= fitted[9]+P1uvt if pha > 360: pha=pha-360 P1_pha[stn,it] =pha if fitted[10] < 0: fitted[10] = -fitted[10] fitted[11] = fitted[11]+180 N2_amp[stn,it] = fitted[10]*N2ft pha= fitted[11]+N2uvt if pha > 360: pha=pha-360 N2_pha[stn,it] = pha print 'M2' print np.mean(M2_amp[29:31]),2*np.std(M2_amp[29:31]) print (M2_amp_obs[29]+M2_amp_obs[29])/2., np.mean(M2_amp[29:31])-(M2_amp_obs[29]+M2_amp_obs[29])/2. print 'S2' print ' South' print np.mean(S2_amp[14:18]/M2_amp[14:18]), 2*np.std(S2_amp[14:18]/M2_amp[14:18]) print np.mean(S2_pha[14:18]-M2_pha[14:18]) print ' North' print np.mean(S2_amp[29:31]/M2_amp[29:31]), 2*np.std(S2_amp[29:31]/M2_amp[29:31]) print np.mean(S2_pha[29:31]-M2_pha[29:31])+360., 2*np.std(S2_pha[29:31]-M2_pha[29:31]) print 'N2' print ' South' print np.mean(N2_amp[14:18]/M2_amp[14:18]), 2*np.std(N2_amp[14:18]/M2_amp[14:18]) print M2_pha[14:18] for i in range(14,18): for j in range (3): if N2_pha[i,j] >300: N2_pha[i,j] = N2_pha[i,j]-360 print N2_pha[14:18] print np.mean(N2_pha[14:18]-M2_pha[14:18]),2*np.std(N2_pha[14:18]-M2_pha[14:18]) print ' North' print np.mean(N2_amp[29:31]/M2_amp[29:31]), 2*np.std(N2_amp[29:31]/M2_amp[29:31]) print np.mean(N2_pha[29:31]-M2_pha[29:31]), 2*np.std(N2_pha[29:31]-M2_pha[29:31]) fig,axs = plt.subplots(6,2,figsize=(15,20)) pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha[:,0])*np.pi/180.) axs[0,1].plot(pha_uw, '-bo', label = 'model') pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha[:,1])*np.pi/180.) axs[0,1].plot(pha_uw, '-go', label = 'model') pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha[:,2])*np.pi/180.) axs[0,1].plot(pha_uw, '-ko', label = 'model') pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha_obs)*np.pi/180.) axs[0,1].plot(pha_uw, 'r-*', label = 'observation') axs[0,1].set_title('M2 Phase') axs[0,0].plot(M2_amp[:,0], '-bo', label = 'model') axs[0,0].plot(M2_amp[:,1], '-go', label = 'model') axs[0,0].plot(M2_amp[:,2], '-ko', label = 'model') axs[0,0].plot(M2_amp_obs, 'r-*', label = 'observation') axs[0,0].set_title('M2 Amp') pha_uw = 180./np.pi * np.unwrap(np.array(K1_pha[:,0])*np.pi/180.) axs[1,1].plot(pha_uw, '-bo', label = 'model') pha_uw = 180./np.pi * np.unwrap(np.array(K1_pha[:,1])*np.pi/180.) axs[1,1].plot(pha_uw, '-go', label = 'model') pha_uw = 180./np.pi * np.unwrap(np.array(K1_pha[:,2])*np.pi/180.) axs[1,1].plot(pha_uw, '-ko', label = 'model') axs[1,1].plot(K1_pha_obs, 'r-*', label = 'observation') axs[1,1].set_title('K1 Phase') axs[1,0].plot(K1_amp[:,0], '-bo', label = 'model') axs[1,0].plot(K1_amp[:,1], '-go', label = 'model') axs[1,0].plot(K1_amp[:,2], '-ko', label = 'model') axs[1,0].plot(K1_amp_obs, 'r-*', label = 'observation') axs[1,0].set_title('K1 Amp') axs[2,1].plot(O1_pha[:,0], '-bo', label = 'model') axs[2,1].plot(O1_pha[:,1], '-go', label = 'model') axs[2,1].plot(O1_pha[:,2], '-ko', label = 'model') axs[2,1].plot(O1_pha_obs, 'r-*', label = 'observation') axs[2,1].set_title('O1 Phase') axs[2,0].plot(O1_amp[:,0], '-bo', label = 'model') axs[2,0].plot(O1_amp[:,1], '-go', label = 'model') axs[2,0].plot(O1_amp[:,2], '-ko', label = 'model') axs[2,0].set_title('O1 Amp') pha_uw = 180./np.pi * np.unwrap(np.array(S2_pha[:,0])*np.pi/180.) axs[3,1].plot(pha_uw, '-bo', label = 'model') pha_uw = 180./np.pi * np.unwrap(np.array(S2_pha[:,1])*np.pi/180.) axs[3,1].plot(pha_uw, '-go', label = 'model') pha_uw = 180./np.pi * np.unwrap(np.array(S2_pha[:,2])*np.pi/180.) axs[3,1].plot(pha_uw, '-ko', label = 'model') axs[3,1].set_title('S2 Phase') axs[3,0].plot(S2_amp[:,0], '-bo', label = 'model') axs[3,0].plot(S2_amp[:,1], '-go', label = 'model') axs[3,0].plot(S2_amp[:,2], '-ko', label = 'model') axs[3,0].set_title('S2 Amp') axs[4,0].plot(P1_amp[:,0], '-bo', label = 'model') axs[4,0].plot(P1_amp[:,1], '-go', label = 'model') axs[4,0].plot(P1_amp[:,2], '-ko', label = 'model') axs[4,0].set_title('P1 Amp') pha_uw = 180./np.pi * np.unwrap(np.array(P1_pha[:,0])*np.pi/180.) axs[4,1].plot(pha_uw, '-bo', label = 'model') pha_uw = 180./np.pi * np.unwrap(np.array(P1_pha[:,1])*np.pi/180.) axs[4,1].plot(pha_uw, '-go', label = 'model') pha_uw = 180./np.pi * np.unwrap(np.array(P1_pha[:,2])*np.pi/180.) axs[4,1].plot(pha_uw, '-ko', label = 'model') axs[4,1].set_title('P1 Phase') axs[5,0].plot(N2_amp[:,0], '-bo', label = 'model') axs[5,0].plot(N2_amp[:,1], '-go', label = 'model') axs[5,0].plot(N2_amp[:,2], '-ko', label = 'model') axs[5,0].set_title('N2 Amp') pha_uw = 180./np.pi * np.unwrap(np.array(N2_pha[:,0])*np.pi/180.) axs[5,1].plot(pha_uw, '-bo', label = 'model') pha_uw = 180./np.pi * np.unwrap(np.array(N2_pha[:,1])*np.pi/180.) axs[5,1].plot(pha_uw, '-go', label = 'model') pha_uw = 180./np.pi * np.unwrap(np.array(N2_pha[:,2])*np.pi/180.) axs[5,1].plot(pha_uw, '-ko', label = 'model') axs[5,1].set_title('N2 Phase') fig,axs = plt.subplots(1,2,figsize=(15,10)) K1mean=np.average(K1_pha, axis=1) K1max = np.max(K1_pha,axis=1) K1min = np.min(K1_pha,axis=1) asymmetric_error = [ K1mean-K1min, K1max-K1mean] axs[0].errorbar(range(31),K1mean, yerr = asymmetric_error) axs[0].plot(K1_pha_obs, 'r-*', label = 'observation') K1mean=np.average(K1_amp, axis=1) K1max = np.max(K1_amp,axis=1) K1min = np.min(K1_amp,axis=1) asymmetric_error = [ K1mean-K1min, K1max-K1mean] axs[1].errorbar(range(31),K1mean, yerr = asymmetric_error) axs[1].plot(K1_amp_obs, 'r-*', label = 'observation')