This notebook will load data, perform a tidal analyis, compare with observations, plot the results, and save the analysis in a spreadsheet. Eight Tidal Constituents: M2, K1, O1, S2, P1, N2, Q1 and K2 are considered.
- Links consistuents that cannot be separated over 40 days
- Corrects nodal calculations
In [57]:
# imports
%matplotlib inline
import matplotlib.pyplot 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
from salishsea_tools import nc_tools
#from nowcast import figures
import scipy.io as sio
import collections
import pandas as pd
import csv
import math
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 [58]:
# pathname for data - all of the tide runs are stored in this directory
#path = '/data/nsoontie/MEOPAR/SalishSea/results/tides/'
path = '../../myResults/NEMO36_Tides/'
#the run we want to analyze
runname = 'GmO_N36_D'
nemo = '3.6'
start = 'apr2003' # 'oct2002' or 'apr2003' 'dec2015'
#joining the two string together
name = path +runname +'/'
print (name)
if nemo == '3.4':
timeinc = 0
elif nemo == '3.6':
if start == 'dec2015':
timeinc = 3659558400
elif start == 'apr2003':
timeinc = 3259872000.
TSout = {'West': {'M2': {'Amp': 0,
'Pha': 0},
'K1': {'Amp': 0,
'Pha': 0},
'O1': {'Amp': 0,
'Pha': 0},
'S2': {'Amp': 0,
'Pha': 0},
'N2': {'Amp': 0,
'Pha': 0},
'Q1': {'Amp': 0,
'Pha': 0}
},
'North': {'M2': {'Amp': 0,
'Pha': 0},
'K1': {'Amp': 0,
'Pha': 0},
'O1': {'Amp': 0,
'Pha': 0},
'S2': {'Amp': 0,
'Pha': 0},
'N2': {'Amp': 0,
'Pha': 0},
'Q1': {'Amp': 0,
'Pha': 0}
}
}
../../myResults/NEMO36_Tides/GmO_N36_D/
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 [59]:
# grid
grid = NC.Dataset('../../nemo-forcing/grid/bathy_downonegrid.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 [60]:
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 .. ... ... ... ... ... ... ... 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 60 Big Bay 50.36000 125.1300 75.5 14.9 83.3 283.5 61 Chatham Point 50.33000 125.4400 90.3 305.1 65.4 270.5 62 Yorke Island 50.44000 125.9700 117.1 271.8 55.8 260.0 63 Powell River 49.86000 124.5500 100.7 34.3 90.4 286.6 64 Lund 49.98000 124.7600 102.2 35.4 88.9 287.9 65 Nymphe Cove 50.13000 125.3600 61.5 350.4 77.0 279.9 66 Brown Bay 50.16000 125.3700 93.5 315.9 67.9 270.1 67 Maude Island E 50.13000 125.3300 55.6 7.4 81.1 283.9 68 Welsford Island 50.22000 125.1300 99.4 35.1 91.1 286.9 69 Redonda Bay 50.26000 124.9900 97.5 36.7 87.1 287.4 70 Channel Islands 50.31000 124.7500 102.6 35.9 89.9 288.0 71 Turnback Point 50.42000 125.1200 102.0 37.0 91.7 287.6 72 Orford Bay 50.59000 124.8600 101.5 37.2 90.3 288.1 73 Waddington Harbour 50.87000 124.8700 103.4 38.0 89.2 288.2 74 Shoal Bay 50.46000 125.3600 89.9 307.5 66.6 269.6 75 Kelsey Bay 50.39000 125.9600 117.0 276.3 57.7 261.4 76 Tacoma 47.26670 122.4133 113.9 11.8 83.8 277.9 [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.
In [61]:
filename = '../Idalia/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)
site lat lon O1_amp O1_pha P1_amp P1_pha \
0 Neah Bay 48.385 -124.616 30.90 231.50 15.50 244.60
1 Port Renfrew 48.537 -124.476 28.30 234.80 14.07 250.60
2 Port Angeles 48.129 -123.400 39.10 241.60 20.70 259.40
3 Victoria 48.413 -123.399 37.00 247.80 19.70 264.60
4 Port Townsend 48.112 -122.758 45.00 249.90 23.90 268.40
5 Bangor 47.748 -122.727 46.60 251.90 26.00 273.90
6 Seattle 47.605 -122.338 45.80 255.40 25.20 274.50
7 Tacoma 47.267 -122.413 45.90 255.10 25.50 277.20
8 Cherry Point 48.863 -122.758 45.60 260.00 25.60 281.40
9 Friday Harbor 48.540 -123.010 42.30 256.40 23.60 274.90
10 Hanbury Point 48.580 -123.172 43.60 253.60 23.40 271.40
11 Sidney 48.658 -123.383 44.40 255.80 24.20 275.20
12 Fulford Harbour 48.765 -123.453 43.00 257.80 23.40 277.80
13 Patos Island 48.783 -122.967 45.50 262.10 24.50 284.60
14 Tsawwassen 48.991 -123.137 47.20 261.80 25.90 282.60
15 Point Atkinson 49.334 -123.250 48.30 263.20 26.80 283.10
16 Winchelsea Islands 49.300 -124.083 47.70 263.50 27.40 286.20
17 Little River 49.744 -124.918 49.26 263.94 28.62 285.67
18 Twin Islets 50.029 -124.936 49.29 264.24 28.62 286.97
19 Campbell River 50.042 -125.247 48.46 263.74 24.60 280.57
20 Seymour Narrows 50.135 -125.347 41.27 254.54 21.28 271.47
21 Owen Bay 50.311 -125.223 38.19 251.34 20.97 267.47
22 Big Bay 50.394 -125.136 46.63 262.44 25.33 282.07
23 Chatham Point 50.332 -125.441 37.46 249.04 20.39 265.97
24 Yorke Island 50.444 -125.975 32.16 241.04 17.10 257.67
25 Alert Bay 50.588 -126.937 30.60 239.84 16.00 251.77
26 Port Hardy 50.720 -127.476 29.70 233.50 15.40 245.50
27 Montagu Point 50.639 -126.213 31.10 237.60 16.60 251.30
28 Siwash Bay 50.680 -125.763 31.30 239.40 17.10 253.20
29 Winter Harbour 50.490 -128.044 27.26 231.20 13.39 242.90
30 Bella Bella 52.177 -128.111 27.80 236.20 14.20 247.20
31 Tofino 49.144 -125.937 24.50 227.20 12.30 237.90
Q1_amp Q1_pha S2_amp S2_pha N2_amp N2_pha K2_amp K2_pha
0 5.50 222.10 22.80 272.6 16.60 222.80 6.00 266.40
1 5.04 225.90 21.04 268.7 15.15 217.30 4.92 263.10
2 6.60 232.80 14.60 326.4 11.60 280.10 2.70 332.70
3 6.10 236.00 10.20 332.8 9.10 292.00 2.00 341.90
4 7.40 243.60 16.80 13.0 14.20 321.80 5.00 18.30
5 8.00 247.20 25.70 29.5 20.80 333.50 7.30 28.50
6 7.50 250.60 25.80 37.9 21.20 340.20 7.20 36.50
7 7.60 250.60 28.20 37.8 22.50 341.20 8.20 39.60
8 7.60 253.20 17.90 50.3 15.40 354.50 5.00 50.50
9 6.80 244.00 13.30 34.9 12.20 341.30 3.50 40.60
10 7.50 247.00 12.70 18.0 11.30 324.90 3.80 37.90
11 7.50 247.00 13.20 26.8 12.00 334.60 3.80 37.90
12 7.00 251.60 13.90 37.2 11.90 342.60 3.90 40.00
13 7.80 253.20 16.70 54.8 14.30 354.20 4.90 58.50
14 6.90 258.50 20.00 55.0 17.20 0.20 5.60 59.40
15 7.70 258.80 22.90 59.9 18.40 2.90 6.20 59.90
16 8.00 257.40 23.60 62.0 20.60 5.60 6.40 64.60
17 8.38 257.20 25.02 61.6 21.64 5.42 6.80 62.56
18 7.89 258.59 25.82 64.8 21.82 9.12 6.92 63.66
19 8.08 252.39 20.27 43.6 19.20 2.82 5.42 49.76
20 7.25 244.99 30.27 339.6 20.48 290.52 8.29 333.06
21 6.37 244.89 27.52 339.6 17.89 290.92 6.89 335.26
22 8.20 224.79 19.29 35.3 15.94 346.02 4.72 35.56
23 5.82 243.69 29.44 326.8 19.57 276.22 8.05 322.36
24 5.33 234.89 38.56 301.2 25.73 248.12 10.70 293.76
25 5.18 231.09 40.63 290.0 26.97 237.72 11.19 279.96
26 5.00 224.30 42.00 281.4 27.30 227.80 10.90 276.20
27 5.20 230.40 49.60 292.7 31.60 238.70 12.50 285.50
28 5.20 232.50 50.60 296.7 32.70 242.50 14.00 290.00
29 4.89 224.50 29.55 273.1 20.74 219.00 7.87 265.80
30 4.90 225.60 40.10 280.0 27.10 227.50 10.90 271.10
31 4.40 219.60 27.90 269.5 20.30 215.60 7.60 261.60
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 [62]:
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']
Next, we can plot these locations on a map of our domain.
In [63]:
fig,ax=plt.subplots(1, 1, figsize=(16, 20))
cmap = plt.get_cmap('winter_r')
cmap.set_bad('teal')
PNW_coastline = sio.loadmat('/ocean/rich/more/mmapbase/bcgeo/PNW.mat')
#figures.plot_map(ax, PNW_coastline, lat_range=(46.5, 52), lon_range=(-128.5, -121.5))
ax.pcolormesh(X,Y,bathy,cmap=cmap)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(14)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(14)
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,'ow',label=location)
ax.annotate(stn, xy = (lon,lat), xytext = (5,5),ha = 'right', va = 'bottom',
textcoords = 'offset points', fontsize=14, color='#F87902')
ax.text(
-127, 47.7,
'Pacific Ocean',
fontsize=30, color='blue')
ax.text(
-124., 51,
' British Columbia',
fontsize=30, color='black')
ax.text(
-123.8, 46.6, 'Washington State',
fontsize=30, color='black')
ax.text(
-122.3, 47.65, ' Puget\nSound',
fontsize=20, color='yellow')
ax.text(
-124.5, 48.1,
' Juan de Fuca Strait',
fontsize=20, color='yellow', rotation=-9)
ax.text(
-124.2, 49.4,
'Strait of Georgia',
fontsize=20, color='yellow', rotation=-30)
ax.text(
-127.5, 50.4,
'Vancouver Island',
fontsize=25, color='black', rotation=-25)
Out[63]:
<matplotlib.text.Text at 0x7f2694462fd0>
In [64]:
# set the split that defines Juan de Fuca and the North
split1=8; split2=21
print (split1, split2)
8 21
Tidal Harmonics¶
In [65]:
sday=2-1; eday = 23-1
fig, ax = plt.subplots(1,1,figsize=(12,5))
for stn in range(12, 16):
print (name)
print (stations[stn])
fT1 = NC.Dataset(name+stations[stn]+'.nc','r')
time = (fT1.variables["time_counter"][:]-timeinc)/3600. # want hours not seconds
ssh = fT1.variables["sossheig"][:,0,0]
print (ssh.shape)
print (ssh[0])
ax.plot(time[(sday-1)*24*2:eday*24*2],ssh[(sday-1)*24*2:eday*24*2])
print (time[0], fT1.variables['time_counter'][0], timeinc)
../../myResults/NEMO36_Tides/GmO_N36_D/ Sandheads (1920,) -1.59517 ../../myResults/NEMO36_Tides/GmO_N36_D/ PointGrey (1920,) -1.56449 ../../myResults/NEMO36_Tides/GmO_N36_D/ PointAtkinson (1920,) -1.57639 ../../myResults/NEMO36_Tides/GmO_N36_D/ GibsonsLanding (1920,) -1.59548 0.25 3259872900.0 3259872000.0
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 [66]:
#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.
#Q1
Q1freq = 13.398661*np.pi/180.
#K2
K2freq = 30.082138*np.pi/180.
# initial phase calculation
# our start is currently Oct 26, 2002
# data for phase output from bdytides.F90; found in ocean.output
# NEMO 3.4 writes out in degrees, NEMO 3.6 writes out in radians
# NEMO 3.6 ut phase of tide due to nodes, vt phase of tide relative to Greenwich
# first one printed out is stable, second one is varying fast
# first on is ut second one is vt
if start == 'oct2002':
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
Q1ft = 1.081364
Q1uvt = 60.631733
K2ft = 1.114095
K2uvt = 52.129248
# for start of Apr 21, 2003
elif start == 'apr2003':
K1ft = 1.065505
K1uvt = 111.481741
M2ft = 0.982328
M2uvt = 250.506179
O1ft = 1.105495
O1uvt = 142.040782
S2ft = 1.000000
S2uvt = 0.000000
P1ft = 1.000000
P1uvt = 241.335269
N2ft = 0.982328
N2uvt = 205.684028
Q1ft = 1.105495
Q1uvt = 97.218631
K2ft = 1.159036
K2uvt = 42.361669
elif start == 'dec2015':
M2ft = 1.03774639446235
M2uvt = (-12.5687951531251+15.0747977989126)*180./np.pi
S2ft = 1.00000000000000
S2uvt = (0.000000000000000+6.28609418926624)*180./np.pi
N2ft = 1.03774639446235
N2uvt = (-12.5687951531251+15.0140965178508)*180./np.pi
K1ft = 0.882044509440140
K1uvt = (-1.198932264051458E-002+6.25441043983728)*180./np.pi
O1ft = 0.806243380430995
O1uvt = (-12.5500821639088 + 8.82038735907528)*180./np.pi
Q1ft = 0.806243380430995
Q1uvt = (-12.5500821639088 + 8.75968607801354)*180./np.pi
K2ft = 0.746702716634735
K2uvt = (-2.156203094515553E-002 + 15.6504135332644)*180./np.pi
P1ft = 1.00000000000000
P1uvt = (0.000000000000000E+000 + 3.168374942895902E-002)*180./np.pi
However, with our typical 40-day output, we cannot separate K1 from P1 and K2 from S2. So we want to work ratios with these.
In [67]:
K1_amp_obs=np.zeros(numsta); K1_pha_obs=np.zeros(numsta)
K2_amp_obs=np.zeros(numsta); K2_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)
P1K1_amp = np.zeros(numsta); P1K1_pha = np.zeros(numsta)
K2S2_amp = np.zeros(numsta); K2S2_pha = np.zeros(numsta)
for stn in range(numsta):
location=stations_obs[stn]
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]
if (harm_other.site==location).any():
K2_amp_obs[stn]=harm_other.K2_amp[harm_other.site==location]/100
K2_pha_obs[stn]=harm_other.K2_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]
P1K1_amp[stn] = P1_amp_obs[stn] / K1_amp_obs[stn]
P1K1_pha[stn] = P1_pha_obs[stn] - K1_pha_obs[stn]
K2S2_amp[stn] = K2_amp_obs[stn] / S2_amp_obs[stn]
K2S2_pha[stn] = K2_pha_obs[stn] - S2_pha_obs[stn]
In [68]:
def interpolate_ratios(amp_ratio, pha_diff):
for stn in range(numsta):
if (amp_ratio[stn] == 0) or (np.isnan(amp_ratio[stn])):
if (amp_ratio[stn-1] != 0) and not(np.isnan(amp_ratio[stn-1])):
if amp_ratio[stn+1] != 0 and not(np.isnan(amp_ratio[stn+1])):
amp_ratio[stn] = 0.5*(amp_ratio[stn-1] + amp_ratio[stn+1])
pha_diff[stn] = 0.5*(pha_diff[stn-1] + pha_diff[stn+1])
else:
amp_ratio[stn] = amp_ratio[stn-1]
pha_diff[stn] = pha_diff[stn-1]
elif (amp_ratio[stn+1] != 0) and not(np.isnan(amp_ratio[stn+1])):
amp_ratio[stn] = amp_ratio[stn+1]
pha_diff[stn] = pha_diff[stn+1]
print (stations[stn], amp_ratio[stn], pha_diff[stn])
In [69]:
interpolate_ratios(P1K1_amp, P1K1_pha)
SheringhamPoint 0.31059602649 -3.5 PedderBay 0.31059602649 -3.5 Esquimalt 0.312395301921 -4.05 CloverPoint 0.314194577352 -4.6 FinnertyCove 0.312475774732 -3.4 TumboChannel 0.310441777195 -1.6 WhalerBay 0.310339070516 -1.6 Sandheads 0.310551558753 -2.2 PointGrey 0.310728215571 -2.6 GibsonsLanding 0.31090487239 -3.0 HalfmoonBay 0.31090487239 -3.0 IrvinesLanding 0.31090487239 -3.0 PowellRiver 0.31090487239 -3.0 Lund 0.313748896372 -1.765 MaudeIslandE 0.290780141844 -3.43 NympheCove 0.299147296356 -2.03 BrownBay 0.309644075586 -2.58 KelseyBay 0.309112656605 -3.43
In [70]:
interpolate_ratios(K2S2_amp, K2S2_pha)
SheringhamPoint 0.233840304183 -5.6 PedderBay 0.233840304183 -5.6 Esquimalt 0.214959367778 1.75 CloverPoint 0.196078431373 9.1 FinnertyCove 0.23832698547 5.95 TumboChannel 0.286994356611 3.25 WhalerBay 0.286706586826 4.05 Sandheads 0.28 4.4 PointGrey 0.275371179039 2.2 GibsonsLanding 0.270742358079 0.0 HalfmoonBay 0.270742358079 0.0 IrvinesLanding 0.270742358079 0.0 PowellRiver 0.270742358079 0.0 Lund 0.269375826599 -0.57 MaudeIslandE 0.26739023187 6.16 NympheCove 0.270629374276 -0.19 BrownBay 0.273653008342 -5.49 KelseyBay 0.275463563278 -5.94
In [71]:
fig,ax = plt.subplots(2, 2, figsize=(10,10))
ax[0,0].plot(P1K1_amp, '*')
ax[0,1].plot(P1K1_pha, 'o')
ax[1,0].plot(K2S2_amp, '*')
ax[1,1].plot(K2S2_pha, 'o')
Out[71]:
[<matplotlib.lines.Line2D at 0x7f2696b96cf8>]
In [72]:
# function for fitting 6 frequencies and inferring 2 with constant
def sixplustwo(xplus, M2amp, M2pha, K1amp, K1pha, O1amp, O1pha, S2amp, S2pha,
N2amp, N2pha, Q1amp, Q1pha, Constant):
x = xplus[0]
P1K1amp = xplus[1]
P1K1pha = xplus[2]
K2S2amp = xplus[3]
K2S2pha = xplus[4]
return (M2amp * M2ft * np.cos(M2freq*x - (M2pha - M2uvt) *np.pi/180.) +
K1amp * K1ft * np.cos(K1freq*x - (K1pha - K1uvt) *np.pi/180.) +
O1amp * O1ft * np.cos(O1freq*x - (O1pha - O1uvt) *np.pi/180.) +
S2amp * S2ft * np.cos(S2freq*x - (S2pha - S2uvt) *np.pi/180.) +
P1K1amp * K1amp *P1ft * np.cos(P1freq*x - (K1pha + P1K1pha - P1uvt) *np.pi/180.) +
N2amp * N2ft * np.cos(N2freq*x - (N2pha - N2uvt) *np.pi/180.) +
Q1amp * Q1ft * np.cos(Q1freq*x - (Q1pha - Q1uvt) *np.pi/180.) +
K2S2amp * S2amp *K2ft * np.cos(K2freq*x - (S2pha + K2S2pha - K2uvt) *np.pi/180.) + Constant )
In [73]:
# Testing to Find the Time Shift
fPA = NC.Dataset(path+runname+'/PointAtkinson.nc')
print (path+runname)
#nc_tools.show_dataset_attrs(fPA)
nc_tools.show_variables(fPA)
nc_tools.show_variable_attrs(fPA,'time_counter')
print (fPA.variables['time_counter'][0]-900.)
print (fPA.variables['time_counter'][1]-fPA.variables['time_counter'][0])
print (fPA.variables['time_centered'][0])
fPO = NC.Dataset('/ocean/sallen/allen/research/MEOPAR/myResults/oldtopog/PointAtkinson.nc')
nc_tools.show_variable_attrs(fPO,'time_counter')
print (fPO.variables['time_counter'][0])
../../myResults/NEMO36_Tides/GmO_N36_D
odict_keys(['nav_lat', 'nav_lon', 'sossheig', 'time_centered', 'time_centered_bounds', 'time_counter', 'time_counter_bounds'])
<class 'netCDF4._netCDF4.Variable'>
float64 time_counter(time_counter)
axis: T
standard_name: time
long_name: Time axis
calendar: gregorian
units: seconds since 1900-01-01 00:00:00
time_origin: 1900-01-01 00:00:00
bounds: time_counter_bounds
unlimited dimensions: time_counter
current shape = (1920,)
filling off
3259872000.0
1800.0
3259872900.0
<class 'netCDF4._netCDF4.Variable'>
float64 time_counter(time_counter)
axis: T
standard_name: time
units: seconds since 2003-04-21 00:00:00
calendar: gregorian
title: Time
long_name: Time axis
time_origin: 2003-APR-21 00:00:00
bounds: time_counter_bnds
unlimited dimensions: time_counter
current shape = (3840,)
filling on, default _FillValue of 9.969209968386869e+36 used
450.0
Now we can apply this fit to our model output.
Use Full Time Series for Best Estimate of Amplitudes/Phases¶
In [74]:
#allocate space for our arrays
M2_amp=[]; M2_pha=[]; K1_amp=[]; K1_pha=[]
O1_amp=[]; O1_pha=[]; S2_amp=[]; S2_pha=[]
N2_amp=[]; N2_pha=[]
Q1_amp=[]; Q1_pha=[]
M2_amp_obs=np.zeros(numsta); M2_pha_obs=np.zeros(numsta)
N2_amp_obs=np.zeros(numsta); N2_pha_obs=np.zeros(numsta)
Q1_amp_obs=np.zeros(numsta); Q1_pha_obs=np.zeros(numsta)
O1_amp_obs=np.zeros(numsta); O1_pha_obs=np.zeros(numsta)
ts = 240
te = ssh.shape[0]
timeplus = np.zeros([5,te-ts])
for stn in range(numsta):
print (name+stations[stn])
fT1 = NC.Dataset(name+stations[stn]+'.nc','r')
time = (fT1.variables["time_counter"][ts:te]-timeinc)/3600. # want hours not seconds
ssh = fT1.variables["sossheig"][:,0,0]
timeplus[0] = time
timeplus[1] = P1K1_amp[stn]*np.ones_like(time); timeplus[2] = P1K1_pha[stn]*np.ones_like(time)
timeplus[3] = K2S2_amp[stn]*np.ones_like(time); timeplus[4] = K2S2_pha[stn]*np.ones_like(time)
fitted, cov = curve_fit(sixplustwo,timeplus,ssh[ts:te])
if fitted[0] < 0:
fitted[0] = -fitted[0]
fitted[1] = fitted[1]+180
M2_amp.append(fitted[0])
pha = fitted[1]
if pha > 360:
pha=pha-360
elif pha < 0:
pha = pha+360
if stn == 6:
print (pha)
M2_pha.append(pha)
if fitted[2] < 0:
fitted[2] = - fitted[2]
fitted[3] = fitted[3] + 180
K1_amp.append(fitted[2])
pha = fitted[3]
if pha > 360:
pha = pha-360
elif pha < 0:
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])
pha= fitted[5]
if pha > 360:
pha=pha-360
elif pha < 0:
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])
pha= fitted[7]
if pha > 360:
pha=pha-360
elif pha < 0:
pha = pha + 360
S2_pha.append(pha)
if fitted[8] < 0:
fitted[8] = -fitted[8]
fitted[9] = fitted[9]+180
N2_amp.append(fitted[8])
pha= fitted[9]
if pha > 360:
pha=pha-360
elif pha < 0:
pha = pha + 360
N2_pha.append(pha)
if fitted[10] < 0:
fitted[10] = -fitted[10]
fitted[11] = fitted[11]+180
Q1_amp.append(fitted[10])
pha= fitted[11]
if pha > 360:
pha=pha-360
Q1_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/Q1/K2 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]
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]
Q1_amp_obs[stn]=harm_other.Q1_amp[harm_other.site==location]/100
Q1_pha_obs[stn]=harm_other.Q1_pha[harm_other.site==location]
#Mask the arrays so that we can do statistics without the 0's throwing things 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)
N2_amp_obs =np.ma.masked_values(N2_amp_obs, 0)
N2_pha_obs =np.ma.masked_values(N2_pha_obs, 0)
Q1_amp_obs =np.ma.masked_values(Q1_amp_obs, 0)
Q1_pha_obs =np.ma.masked_values(Q1_pha_obs, 0)
../../myResults/NEMO36_Tides/GmO_N36_D/PortRenfrew ../../myResults/NEMO36_Tides/GmO_N36_D/SheringhamPoint ../../myResults/NEMO36_Tides/GmO_N36_D/PedderBay ../../myResults/NEMO36_Tides/GmO_N36_D/Esquimalt ../../myResults/NEMO36_Tides/GmO_N36_D/Victoria ../../myResults/NEMO36_Tides/GmO_N36_D/CloverPoint ../../myResults/NEMO36_Tides/GmO_N36_D/FinnertyCove 354.123857399 ../../myResults/NEMO36_Tides/GmO_N36_D/FulfordHarbour ../../myResults/NEMO36_Tides/GmO_N36_D/TumboChannel ../../myResults/NEMO36_Tides/GmO_N36_D/PatosIsland ../../myResults/NEMO36_Tides/GmO_N36_D/WhalerBay ../../myResults/NEMO36_Tides/GmO_N36_D/Tsawwassen ../../myResults/NEMO36_Tides/GmO_N36_D/Sandheads ../../myResults/NEMO36_Tides/GmO_N36_D/PointGrey ../../myResults/NEMO36_Tides/GmO_N36_D/PointAtkinson ../../myResults/NEMO36_Tides/GmO_N36_D/GibsonsLanding ../../myResults/NEMO36_Tides/GmO_N36_D/HalfmoonBay ../../myResults/NEMO36_Tides/GmO_N36_D/IrvinesLanding ../../myResults/NEMO36_Tides/GmO_N36_D/PowellRiver ../../myResults/NEMO36_Tides/GmO_N36_D/Lund ../../myResults/NEMO36_Tides/GmO_N36_D/TwinIslets ../../myResults/NEMO36_Tides/GmO_N36_D/CampbellRiver ../../myResults/NEMO36_Tides/GmO_N36_D/MaudeIslandE ../../myResults/NEMO36_Tides/GmO_N36_D/NympheCove ../../myResults/NEMO36_Tides/GmO_N36_D/SeymourNarrows ../../myResults/NEMO36_Tides/GmO_N36_D/BrownBay ../../myResults/NEMO36_Tides/GmO_N36_D/ChathamPoint ../../myResults/NEMO36_Tides/GmO_N36_D/KelseyBay ../../myResults/NEMO36_Tides/GmO_N36_D/YorkeIsland
In [75]:
stn = 20
timeplus[0] = time
timeplus[1] = P1K1_amp[stn]*np.ones_like(time); timeplus[2] = P1K1_pha[stn]*np.ones_like(time)
timeplus[3] = K2S2_amp[stn]*np.ones_like(time); timeplus[4] = K2S2_pha[stn]*np.ones_like(time)
thefit = sixplustwo(timeplus, M2_amp[stn], M2_pha[stn], K1_amp[stn], K1_pha[stn], O1_amp[stn], O1_pha[stn],
S2_amp[stn], S2_pha[stn], N2_amp[stn], N2_pha[stn], Q1_amp[stn], Q1_pha[stn], 0.)
theobs = sixplustwo(timeplus, M2_amp_obs[stn], M2_pha_obs[stn], K1_amp_obs[stn], K1_pha_obs[stn],
O1_amp_obs[stn], O1_pha_obs[stn],
S2_amp_obs[stn], S2_pha_obs[stn], N2_amp_obs[stn], N2_pha_obs[stn],
Q1_amp_obs[stn], Q1_pha_obs[stn], 0.)
fT1 = NC.Dataset(name+stations[stn]+'.nc','r')
ssh = fT1.variables["sossheig"][:,0,0]
fig, ax = plt.subplots(3,1,figsize=(15,15))
ax[0].plot(time, thefit, linewidth=3)
ax[0].plot(time, theobs, linewidth=2)
ax[0].plot(time, ssh[ts:te]-0.1)
ax[0].set_title('Model Output')
ax[1].plot(time, ssh[ts:te]-thefit-0.1)
ax[1].set_title('Model Ouput - Fit')
ax[2].plot(time, theobs-thefit)
ax[2].set_title("Observations Predictions - Fit")
Out[75]:
<matplotlib.text.Text at 0x7f2696aaaef0>
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 [76]:
print (M2_pha[4], M2_pha[14], M2_pha[14]-M2_pha[4]+360)
315.54791804 31.3253692208 75.7774511806
In [77]:
print (M2_pha_obs[4], M2_pha_obs[14], M2_pha_obs[14]+360-M2_pha_obs[4])
316.1 31.2 75.1
In [78]:
print (M2_amp_obs)
[ 0.708 0.487 0.342 0.367 0.373 0.403 0.447 0.582 0.726 0.68 0.834 0.811 0.869 0.945 0.918 0.947 0.964 0.988 1.007 1.022 1.013 0.825 0.556 0.615 0.946 0.935 0.903 1.17 1.171]
In [79]:
#Plotting M2
labels=['JdF/Islands','SoG','North']
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)
SoGratio = np.mean(M2_amp_obs[split1:split2]/M2_amp[split1:split2])
print (np.mean(M2_amp_obs[split1:split2]), 'mean amp obs')
print ("SoG", M2_pha_obs[split1:split2]-M2_pha[split1:split2])
TSout['West']['M2']['Amp'] = SoGratio
TSout['West']['M2']['Pha'] = np.mean(M2_pha_obs[split1:split2]-M2_pha[split1:split2])
TSout['North']['M2']['Amp'] = np.mean(M2_amp_obs[-2:]/M2_amp[-2:])
TSout['North']['M2']['Pha'] = np.mean(M2_pha_obs[-2:]-M2_pha[-2:])
print (SoGratio, 'Amp')
print (TSout['West']['M2']['Pha'], 'Pha')
print ("North", M2_pha_obs[-2:]-M2_pha[-2:])
print (TSout['North']['M2']['Amp'])
print (TSout['North']['M2']['Pha'])
0.901846153846 mean amp obs SoG [-2.19542223 -1.22479499 0.07586193 -1.74975972 0.16507876 3.12732116 -0.12536922 -1.6678781 -0.52894291 -0.25039887 0.99400532 2.07805186 1.36689132] 0.999946570494 Amp 0.00497263834354 Pha North [ 1.76690119 -1.70555429] 0.999294951295 0.0306734513582
In [80]:
#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)
TSout['West']['K1']['Amp'] = np.mean(K1_amp_obs[split1:split2]/K1_amp[split1:split2])
TSout['West']['K1']['Pha'] = np.mean(K1_pha_obs[split1:split2]-K1_pha[split1:split2])
TSout['North']['K1']['Amp'] = np.mean(K1_amp_obs[-2:]/K1_amp[-2:])
TSout['North']['K1']['Pha'] = np.mean(K1_pha_obs[-2:]-K1_pha[-2:])
print ("SoG", K1_pha_obs[split1:split2]-K1_pha[split1:split2])
print (TSout['West']['K1']['Amp'])
print (TSout['West']['K1']['Pha'])
print ("North", K1_pha_obs[-2:]-K1_pha[-2:])
print (TSout['North']['K1']['Amp'] )
print (TSout['North']['K1']['Pha'] )
SoG [-0.87284528 1.43708679 0.34042665 -0.65854925 0.1025621 0.38410467 0.36231416 -1.40094529 -0.59264421 -0.35142789 -0.37853586 0.73817389 0.03054662] 0.990025093955 -0.0661332994153 North [ 7.39314778 5.91381236] 0.909192405015 6.65348006746
In [81]:
#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)
TSout['West']['O1']['Amp'] = np.mean(O1_amp_obs[split1:split2]/O1_amp[split1:split2])
TSout['West']['O1']['Pha'] = np.mean(O1_pha_obs[split1:split2]-O1_pha[split1:split2])
TSout['North']['O1']['Amp'] = np.mean(O1_amp_obs[-2:]/O1_amp[-2:])
TSout['North']['O1']['Pha'] = np.mean(O1_pha_obs[-2:]-O1_pha[-2:])
print ("SoG", O1_pha_obs[split1:split2]-O1_pha[split1:split2])
print (TSout['West']['O1']['Amp'])
print (TSout['West']['O1']['Pha'])
print ("North", O1_pha_obs[-2:]-O1_pha[-2:])
print (TSout['North']['O1']['Amp'] )
print (TSout['North']['O1']['Pha'] )
SoG [-- 0.5162589643491629 -- -1.0728417350749169 -- -- 0.25128443051977456 -- -- -- -- -- -0.12375121139046996] 0.993814504425 -0.107262387899 North [-- 7.098203504888517] 0.931566748215 7.09820350489
In [82]:
#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)
TSout['West']['S2']['Amp'] = np.mean(S2_amp_obs[split1:split2]/S2_amp[split1:split2])
TSout['West']['S2']['Pha'] = np.mean(S2_pha_obs[split1:split2]-S2_pha[split1:split2])
TSout['North']['S2']['Amp'] = np.mean(S2_amp_obs[-2:]/S2_amp[-2:])
TSout['North']['S2']['Pha'] = np.mean(S2_pha_obs[-2:]-S2_pha[-2:])
print ("SoG", S2_pha_obs[split1:split2]-S2_pha[split1:split2])
print (TSout['West']['S2']['Amp'])
print (TSout['West']['S2']['Pha'])
print ("North", S2_pha_obs[-2:]-S2_pha[-2:])
print (TSout['North']['S2']['Amp'] )
print (TSout['North']['S2']['Pha'] )
SoG [-- 0.10345625361050992 -- -1.5336799792765277 -- -- -0.19457565898763107 -- -- -- -- -- 2.1199963929421344] 0.997477741519 0.123799252072 North [-- -0.21518163627450804] 1.00745155133 -0.215181636275
In [83]:
#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)
TSout['West']['N2']['Amp'] = np.mean(N2_amp_obs[split1:split2]/N2_amp[split1:split2])
TSout['West']['N2']['Pha'] = np.mean(N2_pha_obs[split1:split2]-N2_pha[split1:split2])
TSout['North']['N2']['Amp'] = np.mean(N2_amp_obs[-2:]/N2_amp[-2:])
TSout['North']['N2']['Pha'] = np.mean(N2_pha_obs[-2:]-N2_pha[-2:])
print ("SoG", N2_pha_obs[split1:split2]-N2_pha[split1:split2])
print (TSout['West']['N2']['Amp'])
print (TSout['West']['N2']['Pha'])
print ("North", N2_pha_obs[-2:]-N2_pha[-2:])
print (TSout['North']['N2']['Amp'] )
print (TSout['North']['N2']['Pha'] )
SoG [-- -3.596496465510313 -- -0.7693696568419748 -- -- -0.11044247201252944 -- -- -- -- -- 3.9300069634948596] 1.0056646519 -0.136575407717 North [-- 0.0740307660430517] 1.00052960163 0.0740307660431
In [84]:
#Plotting - Q1
fig=tidetools.plot_scatter_pha_amp(Q1_amp,Q1_amp_obs,Q1_pha,Q1_pha_obs,'Q1',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)
TSout['West']['Q1']['Amp'] = np.mean(Q1_amp_obs[split1:split2]/Q1_amp[split1:split2])
TSout['West']['Q1']['Pha'] = np.mean(Q1_pha_obs[split1:split2]-Q1_pha[split1:split2])
TSout['North']['Q1']['Amp'] = np.mean(Q1_amp_obs[-2:]/Q1_amp[-2:])
TSout['North']['Q1']['Pha'] = np.mean(Q1_pha_obs[-2:]-Q1_pha[-2:])
print ("SoG", Q1_pha_obs[split1:split2]-Q1_pha[split1:split2])
print (TSout['West']['Q1']['Amp'])
print (TSout['West']['Q1']['Pha'])
print ("North", Q1_pha_obs[-2:]-Q1_pha[-2:])
print (TSout['North']['Q1']['Amp'] )
print (TSout['North']['Q1']['Pha'] )
SoG [-- -2.9419299835951165 -- 1.2049915510871188 -- -- 0.3592464198349603 -- -- -- -- -- 0.4567320990878443] 0.992251988257 -0.230239978396 North [-- 7.218067213692507] 0.955535702159 7.21806721369
The model performs well when the dots are close to the red line.
In [85]:
print (TSout)
{'North': {'S2': {'Amp': 1.007451551325621, 'Pha': -0.21518163627450804}, 'K1': {'Amp': 0.90919240501507015, 'Pha': 6.6534800674565417}, 'N2': {'Amp': 1.0005296016330354, 'Pha': 0.074030766043051699}, 'O1': {'Amp': 0.93156674821510022, 'Pha': 7.0982035048885166}, 'Q1': {'Amp': 0.9555357021589741, 'Pha': 7.2180672136925068}, 'M2': {'Amp': 0.99929495129476509, 'Pha': 0.030673451358239845}}, 'West': {'S2': {'Amp': 0.99747774151918511, 'Pha': 0.12379925207212139}, 'K1': {'Amp': 0.99002509395526239, 'Pha': -0.066133299415329166}, 'N2': {'Amp': 1.0056646519018293, 'Pha': -0.13657540771748944}, 'O1': {'Amp': 0.99381450442468544, 'Pha': -0.10726238789911235}, 'Q1': {'Amp': 0.99225198825699834, 'Pha': -0.23023997839629828}, 'M2': {'Amp': 0.99994657049399938, 'Pha': 0.0049726383435359893}}}
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 [86]:
def mean(diff):
return np.mean(abs(diff))
RMS Error¶
In [87]:
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 [88]:
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 [89]:
#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_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))
D_Q1= complex_diff(np.ma.array(Q1_amp_obs),np.ma.array(Q1_pha_obs), np.ma.array(Q1_amp),np.ma.array(Q1_pha))
print (D_M2[2:6], np.mean(D_M2[2:6]))
[ 0.0379399 0.03605283 0.00978363 0.03761776] 0.0303485306333
Saving the results¶
We will now save these statistics in a spreadsheet
In [90]:
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])])
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-90-4b0dafe57c28> in <module>() 6 'Station Name', 7 'R (M2)', 'Delta phi (M2)', 'D (M2)', ----> 8 'R (K1)', 'Delta phi (K1)', 'D (K1)' 9 ]) 10 for stn in range(numsta): TypeError: a bytes-like object is required, not 'str'
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 [ ]:
plt.figure(figsize=(20,12))
plt.subplot(3,2,1)
plt.plot(np.array(M2_amp), '-bo', label = 'model')
plt.plot(M2_amp_obs, 'r-o', label = 'observation')
plt.ylim((0.3,1.3))
plt.title('M2 Amplitude')
plt.legend( loc='upper left' )
plt.subplot(3,2,2)
plt.plot(np.array(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_uwm = 180./np.pi * np.unwrap((np.array(M2_pha))*np.pi/180.)
plt.plot(pha_uwm, '-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.ylim((245,295))
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.ylim((0,0.5))
plt.title('D error')
plt.legend( loc='upper left' )
In [ ]:
S_M2_amp =[0.81796594981353787, 0.58179663165988404, 0.32721527983734328, 0.33301854629278987, 0.33182009758572201, 0.33335988978569214, 0.37262690513061936, 0.47803724001914821, 0.76169885554259253, 0.70303608264417961, 0.83441811463548532, 0.81038178372584047, 0.88924863927669251, 0.92291668835873131, 0.92528642573620279, 0.92708554337687221, 0.97186745021953891, 0.99008962519318389, 1.0104938653574274, 1.0424872258774132, 1.0339748559161059, 0.81647033002975733, 0.62855943384043089, 0.6281276494314838, 0.59982497040887706, 0.91766195442050613, 0.98470695914773221, 1.1750995194900644, 1.1641418000518871]
S_M2_pha =[231.70087338077767, 246.8346272505064, 280.61345724675436, 294.01113561226316, 298.32945634086707, 302.34910604358527, 349.45735869033183, 9.6187167682703034, 33.026761185393298, 25.25761537536448, 32.142427075739192, 28.724244493671563, 29.854819687552208, 30.589764003661628, 30.716301323815514, 30.845594548277688, 31.352881756715988, 31.583460490483901, 32.581761034889155, 32.88157718881746, 33.349768476521263, 18.907952568145106, 357.37521447974547, 344.966653323424, 351.41260102119486, 310.27266040290533, 299.6131712504702, 274.6100536345358, 273.39389940946035]
S_K1_amp =[0.46620238668937491, 0.54775592624668346, 0.62336982624343718, 0.63973158805461383, 0.65184922440835646, 0.66057847399548664, 0.72550375684688451, 0.75355773256988012, 0.84304869161948126, 0.82304770052036202, 0.8564923283826591, 0.84454355069955411, 0.86725095546444986, 0.87129988589140794, 0.86593626071139251, 0.87399549321899139, 0.88682087783902708, 0.89052058701213388, 0.89968632710580798, 0.90846210405972316, 0.91010616927977983, 0.84667321748426194, 0.79016417661225991, 0.75020026652394067, 0.79220776158699524, 0.66728988733516137, 0.64016020034578169, 0.56781550023732319, 0.5644943600526271]
S_K1_pha =[249.16405591522613, 257.05973209735021, 266.58707764365994, 267.05565730516582, 267.3703139483444, 267.57406182793716, 275.3147133634908, 279.37270299241624, 287.55876409635982, 284.10168038262879, 287.02142153805926, 284.99250055965047, 285.91139251534031, 286.31767919206936, 285.53519011464579, 286.15531294135207, 285.95283693592955, 286.94132412236769, 286.58602905644142, 286.95109254814417, 287.2626660314782, 284.94162283332389, 282.4623568253017, 279.60944497950373, 283.00778612693961, 270.05868447585368, 267.32225157060719, 261.00793896779192, 261.23071856768365]
S_D_M2 =[ 0.16625969, 0.17951281, 0.15906926, 0.14399509, 0.11621831, 0.13389999,
0.09472401, 0.10776225, 0.04434292, 0.02324491, 0.01103787, 0.01309181,
0.02582918, 0.05829261, 0.01065969, 0.02335071, 0.00825068, 0.00585006,
0.03045116, 0.04977764, 0.04220142, 0.01121144, 0.1262386, 0.0603622,
0.53398219, 0.09257765, 0.12175468, 0.03495713, 0.03319541]
S_D_K1 = [ 0.04172159, 0.04244914, 0.02657573, 0.01213842, 0.03215966, 0.03138739,
0.03245773, 0.00826599, 0.03342904, 0.03920169, 0.01186245, 0.01091408,
0.03149171, 0.03627737, 0.0093824, 0.01469177, 0.00721647, 0.01116176,
0.00431927, 0.02450076, 0.00716952, 0.01392518, 0.02894072, 0.02017138,
0.17277305, 0.01172017, 0.03845835, 0.00998477, 0.01369324]
In [ ]:
fig, ax = plt.subplots(3,1,figsize=(10,10))
ax[0].plot(S_M2_amp, '-o', color='grey')
ax[0].plot(S_K1_amp, '-s', color='grey')
ax[0].plot(M2_amp_obs, 'g-o', label = 'M2 observation')
ax[0].plot(np.array(M2_amp), '-ro', label = 'M2 model')
ax[0].plot(K1_amp_obs, 'k-s', label = 'K1 observation')
ax[0].plot(np.array(K1_amp), '-bs', label = 'K1 model')
ax[0].set_title('Barotropic Tides - Sea Surface Height')
ax[0].set_ylabel('Amplitude(m)')
ax[0].legend(loc='upper left')
pha_uw = 180./np.pi * np.unwrap(np.array(S_M2_pha)*np.pi/180.)
ax[1].plot(pha_uw, '-o', color='grey')
pha_uw = 180./np.pi * np.unwrap(np.array(S_K1_pha)*np.pi/180.)
ax[1].plot(pha_uw, '-s', color='grey')
pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha_obs)*np.pi/180.)
ax[1].plot(pha_uw, 'g-o', label = 'M2 observation')
pha_uwm = 180./np.pi * np.unwrap((np.array(M2_pha))*np.pi/180.)
ax[1].plot(pha_uwm, '-ro', label = 'M2 model')
ax[1].plot(K1_pha_obs, 'k-s', label = 'K1 observation')
pha_uw = 180./np.pi * np.unwrap((np.array(K1_pha))*np.pi/180.)
ax[1].plot(pha_uw, '-bs', label = 'K1 model')
ax[1].set_ylabel('Phase (deg UTC)')
ax[1].legend(loc=(0.36, 0.4))
ax[2].plot(S_D_M2, '-o', color = 'grey')
ax[2].plot(S_D_K1, '-s', color = 'grey')
ax[2].plot(D_M2, '-ro', label = 'M2')
ax[2].plot(D_K1, '-bs', label = 'K1')
ax[2].yaxis.grid(True)
ax[2].set_ylim((0,0.6))
ax[2].legend( loc='upper left' )
ax[2].set_ylabel('Complex Difference (m)')
ax[2].set_xlabel('Station Number (see map)')
ax[0].set_xlim((0,28))
ax[2].set_xlim((0,28))
ax[1].set_xlim((0,28))
ax[0].text(15, 0.35, 'Grey is Soontiens et al (2016) Tides')
In [ ]:
plt.figure(figsize=(12,5))
plt.plot(M2_amp_obs, 'r-s', label = 'observation')
plt.plot(M2_amp, '-mo', label='results')
print (np.mean(M2_amp[2:6]), np.mean(M2_amp[2:6])*SoGratio)
print (np.mean(M2_amp_obs[2:6]))
In [ ]:
plt.figure(figsize=(12,5))
pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha_obs)*np.pi/180.)
plt.plot(pha_uw, 'r-s', label = 'observation')
pha_uw_mod = 180./np.pi * np.unwrap(np.array(M2_pha)*np.pi/180.)
plt.plot(pha_uw_mod, 'b*-', label='results')
In [ ]:
plt.plot(pha_uw_mod-pha_uw, 'bo-')
In [ ]:
diffy_obs = pha_uw[1:]-pha_uw[0:-1]
diffy_new = pha_uw_mod[1:]-pha_uw_mod[:-1]
plt.plot(diffy_obs,'r',label='obs')
plt.plot(diffy_new,'g*-', label='new')
plt.xlim((0,15))
plt.ylim((-5,50))
plt.legend()
In [ ]:
je=15
plt.figure(figsize=(14,7))
plt.subplot(1,2,1)
plt.plot(pha_uw[:je], M2_amp_obs[:je],'ro-')
plt.plot(pha_uw_mod[:je],np.array(M2_amp[:je]),'ms-')
plt.arrow(365-20, 0.55+0.2, 20, -0.2)
plt.text(365-100, 0.55+0.2, "too much amplitude drop in Boundary Pass")
plt.subplot(1,2,2)
plt.plot(pha_uw[:je], M2_amp_obs[:je],'ro-')
plt.plot(pha_uw_mod[:je],np.array(M2_amp[:je]),'ms-')
plt.arrow(265+20, 0.55+0.2, -20, -0.2)
plt.text(265+20, 0.55+0.2, "too little phase shift in JdF")
In [ ]:
cmap = plt.get_cmap('PuBu')
cmap.set_bad('burlywood')
fig,axs=plt.subplots(3, 2, figsize=(8,20))
constituent = ('M2', 'K1', 'O1', 'S2', 'N2', 'Q1')
error_D = (D_M2, D_K1, D_O1, D_S2, D_N2, D_Q1)
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])
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 [ ]:
fig, axs = plt.subplots(4,2,figsize=(10,15))
axs[0,0].plot(np.array(O1_amp)/np.array(K1_amp), '-bo', label = 'model')
axs[0,0].plot((0,28),(0.560,0.560), 'r-', label = 'observation')
axs[0,0].set_title('O1/K1 Amplitude')
pha_uw = 180./np.pi * np.unwrap((np.array(O1_pha)-np.array(K1_pha))*np.pi/180.)
axs[0,1].plot(pha_uw, '-bo', label = 'model')
axs[0,1].plot((0,28),(-22.9,-22.9), '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.249,0.249), '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(N2_amp)/np.array(M2_amp), '-bo', label = 'model')
axs[2,0].plot((0,28),(0.200,0.200), 'r-', label = 'observation')
axs[2,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[2,1].plot(pha_uw, '-bo', label = 'model')
axs[2,1].plot((0,28),(-28.3, -28.3), 'r-', label = 'observation')
axs[2,1].set_title('N2-M2 Phase')
axs[3,0].plot(np.array(Q1_amp)/np.array(K1_amp), '-bo', label = 'model')
axs[3,0].plot((0,28),(0.089,0.089), 'r-', label = 'observation')
axs[3,0].set_title('Q1/K1 Amplitude')
pha_uw = 180./np.pi * np.unwrap((np.array(Q1_pha)-np.array(K1_pha))*np.pi/180.)
axs[3,1].plot(pha_uw+360., '-bo', label = 'model')
axs[3,1].plot((0,28),(-27.3+360,-27.3+360), 'r-', label = 'observation')
axs[3,1].set_title('Q1-K1 Phase')
In [ ]:
sample = 17
start = np.zeros(sample)
tend = np.zeros(sample)
for i in range(sample):
start[i] = 196+(480-196)*np.random.rand()
tend[i] = te-(480-196)*np.random.rand()
print start
print tend
timelength = (tend-start)/96.
print np.mean(timelength),2*np.std(timelength)
print time[start[1]:tend[1]]
In [ ]:
#allocate space for our arrays
M2_amp=np.zeros((numsta,sample)); M2_pha=np.zeros((numsta,sample))
K1_amp=np.zeros((numsta,sample)); K1_pha=np.zeros((numsta,sample))
O1_amp=np.zeros((numsta,sample)); O1_pha=np.zeros((numsta,sample))
S2_amp=np.zeros((numsta,sample)); S2_pha=np.zeros((numsta,sample))
N2_amp=np.zeros((numsta,sample)); N2_pha=np.zeros((numsta,sample))
Q1_amp=np.zeros((numsta,sample)); Q1_pha=np.zeros((numsta,sample))
for it,tst,tet in zip(range(sample),start.astype(int),tend.astype(int)):
timeplus = np.zeros([5,tet-tst])
for stn in range(numsta):
fT1 = NC.Dataset(name+stations[stn]+'.nc','r')
time = (fT1.variables["time_counter"][tst:tet]-timeinc)/3600. # want hours not seconds
ssh = fT1.variables["sossheig"][:,0,0]
timeplus[0] = time
timeplus[1] = P1K1_amp[stn]*np.ones_like(time); timeplus[2] = P1K1_pha[stn]*np.ones_like(time)
timeplus[3] = K2S2_amp[stn]*np.ones_like(time); timeplus[4] = K2S2_pha[stn]*np.ones_like(time)
fitted, cov = curve_fit(sixplustwo,timeplus,ssh[tst:tet])
if fitted[0] < 0:
fitted[0] = -fitted[0]
fitted[1] = fitted[1]+180
M2_amp[stn,it] = fitted[0]
pha = fitted[1]
if pha > 360:
pha=pha-360
elif pha < 0:
pha = pha + 360
M2_pha[stn,it] = pha
if fitted[2] < 0:
fitted[2] = -fitted[2]
fitted[3] = fitted[3]+180
K1_amp[stn,it] = fitted[2]
pha= fitted[3]
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]
pha= fitted[5]
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]
pha= fitted[7]
if pha > 360:
pha=pha-360
S2_pha[stn,it]= pha
if fitted[8] < 0:
fitted[8] = -fitted[8]
fitted[9] = fitted[9]+180
N2_amp[stn,it] = fitted[8]
pha= fitted[9]
if pha > 360:
pha=pha-360
N2_pha[stn,it] = pha
if fitted[10] < 0:
fitted[10] = -fitted[10]
fitted[11] = fitted[11]+180
Q1_amp[stn,it] = fitted[10]
pha= fitted[11]
if pha > 360:
pha=pha-360
Q1_pha[stn,it] = pha
In [ ]:
jdef = range(3)
south = range(14,18)
north = range(27,29)
print 'M2'
print ' JdeFuca'
print np.mean(M2_amp[jdef]),2*np.std(np.mean(M2_amp[jdef],axis=0))
print np.mean(M2_amp_obs[jdef]), np.mean(M2_amp_obs[jdef])-np.mean(M2_amp[jdef])
print np.mean(M2_amp_obs[jdef])/np.mean(M2_amp[jdef])
print np.mean(M2_pha[jdef]),2*np.std(np.mean(M2_pha[jdef],axis=0))
print np.mean(M2_pha_obs[jdef]), np.mean(M2_pha_obs[jdef])-np.mean(M2_pha[jdef])
print ' South'
print np.mean(M2_amp[south]),2*np.std(np.mean(M2_amp[south],axis=0))
print np.mean(M2_amp_obs[south]), np.mean(M2_amp_obs[south])-np.mean(M2_amp[south])
print np.mean(M2_amp_obs[south])/np.mean(M2_amp[south])
print np.mean(M2_pha[south]),2*np.std(np.mean(M2_pha[south],axis=0))
print np.mean(M2_pha_obs[south]), np.mean(M2_pha_obs[south])-np.mean(M2_pha[south])
print ' North'
print np.mean(M2_amp[north]),2*np.std(np.mean(M2_amp[north],axis=0))
print np.mean(M2_amp_obs[north]), np.mean(M2_amp_obs[north])-np.mean(M2_amp[north])
print np.mean(M2_amp_obs[north])/np.mean(M2_amp[north])
print np.mean(M2_pha[north]),2*np.std(np.mean(M2_pha[north],axis=0))
print np.mean(M2_pha_obs[north]), np.mean(M2_pha_obs[north])-np.mean(M2_pha[north])
print '==============================================='
print 'K1'
print ' JdeFuca'
print np.mean(K1_amp[jdef]),2*np.std(np.mean(K1_amp[jdef],axis=0))
print np.mean(K1_amp_obs[jdef]), np.mean(K1_amp_obs[jdef])-np.mean(K1_amp[jdef])
print np.mean(K1_amp_obs[jdef])/np.mean(K1_amp[jdef])
print np.mean(K1_pha[jdef]),2*np.std(np.mean(K1_pha[jdef],axis=0))
print np.mean(K1_pha_obs[jdef]), np.mean(K1_pha_obs[jdef])-np.mean(K1_pha[jdef])
print ' South'
print np.mean(K1_amp[south]),2*np.std(np.mean(K1_amp[south],axis=0))
print np.mean(K1_amp_obs[south]), np.mean(K1_amp_obs[south])-np.mean(K1_amp[south])
print np.mean(K1_amp_obs[south])/np.mean(K1_amp[south])
print np.mean(K1_pha[south]),2*np.std(np.mean(K1_pha[south],axis=0))
print np.mean(K1_pha_obs[south]), np.mean(K1_pha_obs[south])-np.mean(K1_pha[south])
print ' North'
print np.mean(K1_amp[north]),2*np.std(np.mean(K1_amp[north],axis=0))
print np.mean(K1_amp_obs[north]), np.mean(K1_amp_obs[north])-np.mean(K1_amp[north])
print np.mean(K1_amp_obs[north])/np.mean(K1_amp[north])
print np.mean(K1_pha[north]),2*np.std(np.mean(K1_pha[north],axis=0))
print np.mean(K1_pha_obs[north]), np.mean(K1_pha_obs[north])-np.mean(K1_pha[north])
print '==============================================='
In [ ]:
print 'O1'
print ' South'
print np.mean(O1_amp[south]/K1_amp[south]),2*np.std(np.mean(O1_amp[south]/K1_amp[south],axis=0))
print np.mean(O1_amp_obs[south]/K1_amp_obs[south]), (np.mean(O1_amp_obs[south]/K1_amp_obs[south])
-np.mean(O1_amp[south]/K1_amp[south]))
print np.mean(O1_amp_obs[south]/K1_amp_obs[south])/np.mean(O1_amp[south]/K1_amp[south])
print np.mean(O1_pha[south]-K1_pha[south]),2*np.std(np.mean(O1_pha[south]-K1_pha[south],axis=0))
print np.mean(O1_pha_obs[south]-K1_pha_obs[south]), (np.mean(O1_pha_obs[south]-K1_pha_obs[south])
-np.mean(O1_pha[south]-K1_pha[south]))
print ' North'
print np.mean(O1_amp[north]/K1_amp[north]),2*np.std(np.mean(O1_amp[north]/K1_amp[north],axis=0))
print np.mean(O1_amp_obs[north]/K1_amp_obs[north]), (np.mean(O1_amp_obs[north]/K1_amp_obs[north])
-np.mean(O1_amp[north]/K1_amp[north]))
print np.mean(O1_amp_obs[north]/K1_amp_obs[north])/np.mean(O1_amp[north]/K1_amp[north])
print np.mean(O1_pha[north]-K1_pha[north]),2*np.std(np.mean(O1_pha[north]-K1_pha[north],axis=0))
print np.mean(O1_pha_obs[north]-K1_pha_obs[north]), (np.mean(O1_pha_obs[north]-K1_pha_obs[north])
-np.mean(O1_pha[north]-K1_pha[north]))
print '==============================================='
print 'S2'
code = ('south','north')
for dir,dire in zip(code,(south,north)):
print dir
print np.mean(S2_amp[dire]/M2_amp[dire]),2*np.std(np.mean(S2_amp[dire]/M2_amp[dire],axis=0))
print np.mean(S2_amp_obs[dire]/M2_amp_obs[dire]), (np.mean(S2_amp_obs[dire]/M2_amp_obs[dire])
-np.mean(S2_amp[dire]/M2_amp[dire]))
print np.mean(S2_amp_obs[dire]/M2_amp_obs[dire])/np.mean(S2_amp[dire]/M2_amp[dire])
unwrap = np.unwrap(np.array(S2_pha)*np.pi/180.)*180./np.pi
M2_un = np.unwrap(np.array(M2_pha)*np.pi/180.)*180./np.pi
plt.plot (unwrap[dire],'r',M2_un[dire],'b')
print np.mean(unwrap[dire]-M2_un[dire])+360.,2*np.std(np.mean(unwrap[dire]-M2_un[dire],axis=0))
print np.mean(S2_pha_obs[dire]-M2_pha_obs[dire]), (np.mean(S2_pha_obs[dire]-M2_pha_obs[dire])
-np.mean(unwrap[dire]-M2_un[dire]))-360.
In [ ]:
const = ('Q1', 'N2')
model_amp = (Q1_amp, N2_amp)
model_pha = ()
for const, model_amp, model_pha, obs_amp, obs_pha in zip(('Q1','N2'),
(Q1_amp, N2_amp),(Q1_pha, N2_pha),
(Q1_amp_obs, N2_amp_obs), (Q1_pha_obs, N2_pha_obs)):
print const
for dir,dire in zip(code,(south,north)):
print dir
print np.mean(model_amp[dire]/K1_amp[dire]),2*np.std(np.mean(model_amp[dire]/K1_amp[dire],axis=0))
print np.mean(obs_amp[dire]/K1_amp_obs[dire]), (np.mean(obs_amp[dire]/K1_amp_obs[dire])
-np.mean(model_amp[dire]/K1_amp[dire]))
print np.mean(obs_amp[dire]/K1_amp_obs[dire])/np.mean(model_amp[dire]/K1_amp[dire])
unwrap = np.unwrap(np.array(model_pha)*np.pi/180.)*180./np.pi
K1_un = np.unwrap(np.array(K1_pha)*np.pi/180.)*180./np.pi
print np.mean(unwrap[dire]-K1_un[dire]),2*np.std(np.mean(unwrap[dire]-K1_un[dire],axis=0))
print np.mean(obs_pha[dire]-K1_pha_obs[dire]), (np.mean(obs_pha[dire]-K1_pha_obs[dire])
-np.mean(unwrap[dire]-K1_un[dire]))
In [ ]:
fig,axs = plt.subplots(6,2,figsize=(15,25))
for i in range(sample):
pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha[:,i])*np.pi/180.)
axs[0,1].plot(pha_uw ,'-ob', 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')
for i in range(sample):
axs[0,0].plot(M2_amp[:,i], '-bo', label = 'model')
axs[0,0].plot(M2_amp_obs, 'r-*', label = 'observation')
axs[0,0].set_title('M2 Amp')
for i in range(sample):
if K1_pha[0,i] < -360:
K1_pha[0,i] = K1_pha[0,i] +720.
elif K1_pha[0,i] < 0:
K1_pha[0,i] = K1_pha[0,i] + 360.
pha_uw = 180./np.pi * np.unwrap(np.array(K1_pha[:,i])*np.pi/180.)
axs[1,1].plot(pha_uw, '-bo', label = 'model')
axs[1,1].plot(K1_pha_obs, 'r-*', label = 'observation')
axs[1,1].set_title('K1 Phase')
for i in range(sample):
axs[1,0].plot(K1_amp[:,i], '-bo', label = 'model')
axs[1,0].plot(K1_amp_obs, 'r-*', label = 'observation')
axs[1,0].set_title('K1 Amp')
for i in range(sample):
if O1_pha[0,i] < 0:
O1_pha[0,i] = O1_pha[0,i] + 360
pha_uw = 180./np.pi * np.unwrap(np.array(O1_pha[:,i])*np.pi/180.)
axs[2,1].plot(pha_uw, '-bo', label = 'model')
axs[2,1].plot(O1_pha_obs, 'r-*', label = 'observation', markersize = 15)
axs[2,1].set_title('O1 Phase')
for i in range(sample):
axs[2,0].plot(O1_amp[:,i], '-bo', label = 'model')
axs[2,0].plot(O1_amp_obs, 'r-*', label = 'observation', markersize = 15)
axs[2,0].set_title('O1 Amp')
for i in range(sample):
if S2_pha[0,i] < 0:
S2_pha[0,i] = S2_pha[0,i] + 360
pha_uw = 180./np.pi * np.unwrap(np.array(S2_pha[:,i])*np.pi/180.)
axs[3,1].plot(pha_uw, '-bo', label = 'model')
pha_uw = 180./np.pi * np.unwrap(np.array(S2_pha_obs)*np.pi/180.)
vsmall = 1e-6
pha_uwm = np.ma.masked_array(pha_uw, mask=(abs(pha_uw-360)<vsmall))
axs[3,1].plot(pha_uwm, 'r-*', label = 'observation', markersize = 15)
axs[3,1].set_title('S2 Phase')
for i in range(sample):
axs[3,0].plot(S2_amp[:,i], '-bo', label = 'model')
axs[3,0].plot(S2_amp_obs, 'r-*', label = 'observation', markersize = 15)
axs[3,0].set_title('S2 Amp')
for i in range(sample):
axs[4,0].plot(N2_amp[:,i], '-bo', label = 'model')
axs[4,0].plot(N2_amp_obs, 'r-*', label = 'observation', markersize = 15)
axs[4,0].set_title('N2 Amp')
for i in range(sample):
pha_uw = 180./np.pi * np.unwrap(np.array(N2_pha[:,i])*np.pi/180.)
axs[4,1].plot(pha_uw, '-bo', label = 'model')
pha_uw = 180./np.pi * np.unwrap(np.array(N2_pha_obs)*np.pi/180.)
pha_uwm = np.ma.masked_array(pha_uw, mask=(abs(pha_uw-360)<vsmall))
axs[4,1].plot(pha_uwm, 'r-*', label = 'observation', markersize = 15)
axs[4,1].set_title('N2 Phase')
for i in range(sample):
axs[5,0].plot(Q1_amp[:,i], '-bo', label = 'model')
axs[5,0].plot(Q1_amp_obs, 'r-*', label = 'observation', markersize = 15)
axs[5,0].set_title('Q1 Amp')
for i in range(sample):
pha_uw = 180./np.pi * np.unwrap(np.array(Q1_pha[:,i])*np.pi/180.)
for j in range(numsta):
if pha_uw[j] < 0:
pha_uw[j] += 360
axs[5,1].plot(pha_uw, '-bo', label = 'model')
axs[5,1].plot(Q1_pha_obs, 'r-*', label = 'observation', markersize = 15)
axs[5,1].set_title('Q1 Phase')
In [ ]: