In [12]:
# 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
Run Details¶
First, let's define the run that we will be analyzing. We can analyze a different run by changing runname in the cell below. A spreadsheet called tide_runs.ods contains a list of runs that we can look at.
In [13]:
# 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 = 'RC4_wO1S2'
#joining the two string together
name = path +runname +'/'
print name
/data/nsoontie/MEOPAR/SalishSea/results/tides/RC4_wO1S2/
We'll also load the bathymetry data in case we want to look at that. The package tidetools has a function get_SS_bathy_data() that returns bathymetry and grid data.
In [14]:
# grid
grid = NC.Dataset('../../nemo-forcing/grid/bathy_meter_SalishSea2.nc')
bathy, X, Y = tidetools.get_bathy_data(grid)
Observations¶
Next, we can load some observations from a text file: /data/nsoontie/MEOPAR/analysis/compare_tides/obs_tidal_wlev_const_all.csv Note: This file contains a mix of M2/K1 measurements from Foreman et al (1995), US tidal harmonics, Foreman et al (2004) and Foreman et al (2012) (for Northern tides).
In [15]:
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
site lat lon M2_amp M2_pha K1_amp K1_pha
0 Sooke 48.36700 123.7330 43.8 282.7 56.9 266.4
1 Port Angeles 48.12500 123.4400 51.8 307.4 66.9 261.4
2 Pedder Bay 48.33100 123.5490 34.2 308.0 62.7 269.0
3 Esquimalt 48.43300 123.4330 36.7 317.1 64.3 268.1
4 Clover Point 48.40500 123.3470 40.3 320.3 64.2 269.8
5 Victoria 48.41700 123.3670 37.3 316.1 62.7 269.2
6 Finnerty Cove 48.47300 123.2950 44.7 357.7 70.8 277.5
7 Port Townsend 48.14500 122.7550 65.2 350.0 75.0 270.8
8 Sidney 48.65000 123.4000 55.4 5.9 76.7 277.6
9 Patricia Bay 48.65000 123.4500 60.3 14.4 76.0 281.3
10 Maple Bay 48.81700 123.6170 68.5 17.0 79.3 281.2
11 Fulford Harbour 48.76700 123.4500 58.2 12.7 75.3 280.0
12 Ladysmith 48.98300 123.8000 70.8 16.3 79.8 281.8
13 Patos Island 48.78300 122.9670 68.0 25.0 79.0 285.6
14 Tumbo Channel 48.79200 123.1080 72.6 31.0 81.1 286.9
15 Whaler Bay 48.88500 123.3250 83.4 32.9 84.7 287.5
16 Silva Bay 49.15300 123.7000 92.2 32.0 86.5 286.7
17 Ferndale 48.83300 122.7170 72.3 23.8 80.1 283.6
18 Blaine 48.99000 122.7600 77.4 25.1 82.3 284.3
19 Tsawwassen 48.99000 123.1330 81.1 27.8 83.4 284.8
20 Sandheads 49.10000 123.3000 86.9 30.9 83.7 286.5
21 Point Grey 49.25000 123.2670 94.5 33.9 90.6 287.0
22 Point Atkinson 49.33300 123.2500 91.8 31.2 86.2 286.1
23 Squamish 49.70000 123.1500 94.2 31.2 87.4 286.8
24 Gibsons Landing 49.40000 123.5000 94.7 30.1 87.2 285.2
25 Halfmoon Bay 49.51700 123.9170 96.4 31.5 88.0 285.8
26 Irvines Landing 49.63300 124.0500 98.8 31.9 88.0 286.7
27 Winchelsea 49.30000 124.0830 95.2 32.6 87.5 286.7
28 Northwest Bay 49.30000 124.2000 95.6 32.7 87.2 286.7
29 Cherry Point 48.86300 122.7570 73.2 21.8 81.5 281.9
30 Friday Harbour 48.54500 123.0100 56.4 11.0 76.3 279.4
31 Rosario 48.64700 122.8700 58.1 4.3 75.3 277.3
32 North Beach 48.71200 122.9080 68.4 18.1 79.2 281.8
33 Bellingham 48.74500 122.4950 65.3 15.2 77.0 281.1
34 Anacortes 48.52300 122.6130 62.2 5.2 77.3 277.2
35 Reservation Bay 48.41500 122.6520 58.1 344.7 73.8 269.5
36 Hanbury Point 48.58000 123.1720 52.8 357.3 75.3 275.1
37 Sheringham Point 48.37500 123.9180 48.7 263.3 54.8 261.5
38 Neah Bay 48.36670 124.6117 78.9 246.2 49.6 248.5
39 Sekiu Clallam Bay 48.26330 124.2967 66.9 258.0 54.2 252.2
40 Port Angeles 48.12500 123.4400 51.4 307.1 66.7 261.8
41 Port Townsend 48.11170 122.7567 68.5 349.7 76.6 270.7
42 Budd Inlet 47.09830 122.8950 145.9 30.3 87.7 289.6
43 Seattle 47.60170 122.3383 107.2 10.6 83.4 277.0
44 Bangor 47.74830 122.7267 102.4 4.1 83.5 274.8
45 Foulweather Bluff 47.92670 122.6167 88.3 3.4 76.5 275.9
46 Everett 47.98000 122.2217 100.9 11.6 80.9 276.9
47 Sneeoosh Point 48.40000 122.5467 102.6 18.3 78.4 282.0
48 Turner Bay 48.44500 122.5550 94.4 16.7 75.4 281.4
49 Armitage Island 48.53500 122.7967 57.3 0.5 75.6 276.4
50 Friday Harbour 48.54670 123.0100 56.5 9.7 75.8 278.8
51 Richardson 48.44670 122.9000 52.2 340.1 71.3 270.9
52 Cherry Point 48.86330 122.7567 73.4 22.8 81.7 282.8
53 Blaine 48.99167 122.7650 76.3 24.8 78.4 286.3
54 Port Renfrew 48.55000 124.4300 70.8 241.1 45.3 254.1
55 Little River 49.74000 124.9200 99.4 32.9 90.2 287.0
56 Twin Islets 50.03000 124.9300 101.3 35.4 90.4 287.5
57 Campbell River 50.04000 125.2400 82.5 18.4 84.6 284.0
58 Seymour Narrows 50.13000 125.3400 94.6 320.1 69.2 272.1
59 Owen Bay 50.31000 125.2200 85.0 319.9 67.8 272.7
... ... ... ... ... ... ...
[77 rows x 7 columns]
This is a list of observations that we can compare with our model output. Now we have a struc object called harm_obs that contains the data printed above.
Model¶
We don't have model output at all of the above locations. The model outputs are listed below. There is a location.nc file in the run directory for each of the stations listed below.
In [16]:
stations = ['PortRenfrew','SheringhamPoint','PedderBay', 'Esquimalt',
'Victoria','CloverPoint','FinnertyCove', 'FulfordHarbour',
'TumboChannel','PatosIsland','WhalerBay', 'Tsawwassen',
'Sandheads', 'PointGrey','PointAtkinson','GibsonsLanding',
'HalfmoonBay','IrvinesLanding','PowellRiver', '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',
'Halfmoon Bay','Irvines Landing','Powell River', 'Lund',
'Twin Islets','Campbell River','Maude Island E', 'Nymphe Cove',
'Seymour Narrows','Brown Bay','Chatham Point','Kelsey Bay','Yorke Island']
Next, we can plot these locations on a map of our domain.
In [17]:
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])
0 Port Renfrew 1 Sheringham Point 2 Pedder Bay 3 Esquimalt 4 Victoria 5 Clover Point 6 Finnerty Cove 7 Fulford Harbour 8 Tumbo Channel 9 Patos Island 10 Whaler Bay 11 Tsawwassen 12 Sandheads 13 Point Grey 14 Point Atkinson 15 Gibsons Landing 16 Halfmoon Bay 17 Irvines Landing 18 Powell River 19 Lund 20 Twin Islets 21 Campbell River 22 Maude Island E 23 Nymphe Cove 24 Seymour Narrows 25 Brown Bay 26 Chatham Point 27 Kelsey Bay 28 Yorke Island
Out[17]:
[-126.1, -122, 47, 51]
Note: Some day it would be worthwhile to place the numbers more carefully so that they don't overlap.
Tidal Harmonics¶
We need a way of determing the amplitude and phase of M2/K1/O1/S2 from our model output. We will do this by fitting our model water levels to cosine curves with the known frequency of M2/K1/O1/S2.
In [18]:
#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.
# 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
In [19]:
# 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.))
In [20]:
# 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.))
In [21]:
# 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.))
Now we can apply this fit to our model output.
In [22]:
#allocate space for our arrays
M2_amp=[]
M2_pha=[]
K1_amp=[]
K1_pha=[]
O1_amp=[]
O1_pha=[]
S2_amp=[]
S2_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)
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(quad,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)
#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]
The model data is saved in lists M2_amp, M2_pha, K1_amp, K1_pha. We have also saved the observations in M2_amp_obs, etc.
We can compare model and observations by plotting.
In [23]:
#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)
Out[23]:
[<matplotlib.lines.Line2D at 0x52ceb10>]
In [24]:
#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)
Out[24]:
[<matplotlib.lines.Line2D at 0x595c810>]
The model performs well when the dots are close to the red line.
Statistics¶
We would like to save some statistics so that we can determine which runs give us the best match with observations. So, we will define some functions that will help us calculate statistics.
Mean Error (absolute value)¶
In [25]:
def mean(diff):
return np.mean(abs(diff))
RMS Error¶
In [26]:
def rms(diff):
return np.sqrt(np.mean(diff**2))
Complex differences¶
This is a way of measuring distances in the complex plane. We can think of our tidal amplitude and phase as a point on the complex plane. So we would like to measure the distance between a point given by the model and a point given by the observations. The function below does this.
In [27]:
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
Some other things we will look at are
$R = \frac{A_m}{A_o}$, the ratio of modelled to observed amplitude and
$\Delta \phi = \phi_m - \phi_o$, the difference betwen modelled and observed phase.
In [28]:
#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))
Saving the results¶
We will now save these statistics in a spreadsheet
In [29]:
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])])
Now there is a csv file in this directory with data about this run. It should be called runname.csv (where runname is the string we defined at the beginning of the notebook).
Things to try:
- Add the complex differences information printed in the .csv file to tide_runs.ods. Also, if you notice any discrepancies, you can correct them. (Check M2 amplitude at Yorke Island)
- Work through the notebook with a different run. There is a list of runs in tide_runs.ods
- Commit and push any changes you've made to this notebook and tide_runs.odt.
Try this: * hg status (see what changes have been made) * hg in * hg commit mynotebook.ipynb (write a commit message and then save and exit) * hg commit tide_runs.odt * hg pull --rebase * pg push
- Add any new csv files to the repository.
Try this: * hg add filename.csv * hg commit filename.csv * hg pull --rebase * hg push
- Repeat with all the runs listed in tide_runs.odt
Plots comparing measured and model amplitudes and phases¶
In [30]:
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
M2_pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha)*np.pi/180.)
plt.plot(M2_pha_uw, '-bo', label = 'model')
M2_pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha_obs)*np.pi/180.)
plt.plot(M2_pha_uw, 'r-o', label = 'observation')
plt.title('M2 Phase')
plt.subplot(3,2,4)
plt.plot(K1_pha, '-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, '-bo', 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' )
Out[30]:
<matplotlib.legend.Legend at 0x6652490>
In [31]:
plt.figure(figsize=(6,6))
plt.plot(D_M2, '-co', label = 'M2')
plt.plot(D_K1, '-mo', 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(1, 2, figsize=(20,12))
error_D = (D_M2, D_K1)
for ax,error_D1 in zip(axs, error_D):
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.axis([-126.1,-122,47,51])
First Plot: M2 D error. Second Plot: K1 D error. Markers increase in size as D error increases.¶
Green: D error <= 0.05, Yellow: 0.05 < D error <= 0.1, Red: D error > 0.1
In [32]:
fig, axs = plt.subplots(2,2,figsize=(10,5))
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')
Out[32]:
<matplotlib.text.Text at 0x8b32f10>
In [37]:
print np.mean(O1_pha[14:18]-K1_pha[14:18])
-24.2412463451
In [33]:
#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))
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(quad,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
In [34]:
fig,axs = plt.subplots(4,2,figsize=(15,15))
M2_pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha[:,0])*np.pi/180.)
axs[0,1].plot(M2_pha_uw, '-bo', label = 'model')
M2_pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha[:,1])*np.pi/180.)
axs[0,1].plot(M2_pha_uw, '-go', label = 'model')
M2_pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha[:,2])*np.pi/180.)
axs[0,1].plot(M2_pha_uw, '-ko', label = 'model')
M2_pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha_obs)*np.pi/180.)
axs[0,1].plot(M2_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')
axs[1,1].plot(K1_pha[:,0], '-bo', label = 'model')
axs[1,1].plot(K1_pha[:,1], '-go', label = 'model')
axs[1,1].plot(K1_pha[:,2], '-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].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')
S2_pha_uw = 180./np.pi * np.unwrap(np.array(S2_pha[:,0])*np.pi/180.)
axs[3,1].plot(S2_pha_uw, '-bo', label = 'model')
S2_pha_uw = 180./np.pi * np.unwrap(np.array(S2_pha[:,1])*np.pi/180.)
axs[3,1].plot(S2_pha_uw, '-go', label = 'model')
S2_pha_uw = 180./np.pi * np.unwrap(np.array(S2_pha[:,2])*np.pi/180.)
axs[3,1].plot(S2_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')
Out[34]:
<matplotlib.text.Text at 0xe111bd0>
In [35]:
fig,axs = plt.subplots(2,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,0].errorbar(range(29),K1mean, yerr = asymmetric_error)
axs[0,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[0,1].errorbar(range(29),K1mean, yerr = asymmetric_error)
axs[0,1].plot(K1_amp_obs, 'r-*', label = 'observation')
Out[35]:
[<matplotlib.lines.Line2D at 0x7dcd310>]
In [35]: