In [117]:
# 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 [118]:
# 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
/data/nsoontie/MEOPAR/SalishSea/results/tides/RC5_wO1S2P1N2/
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 [119]:
# 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 [120]:
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.
In [121]:
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
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 Harbour 8.413 -123.399 37.00 247.80 19.70 264.60
4 Port Townsend 8.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 7.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
[32 rows x 15 columns]
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 [122]:
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 [123]:
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 Winchelsea 17 Halfmoon Bay 18 Irvines Landing 19 Powell River 20 Little River 21 Lund 22 Twin Islets 23 Campbell River 24 Maude Island E 25 Nymphe Cove 26 Seymour Narrows 27 Brown Bay 28 Chatham Point 29 Kelsey Bay 30 Yorke Island
Out[123]:
[-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 [124]:
#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
In [125]:
# 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 [126]:
# 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 [127]:
# 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.))
In [128]:
# 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.))
Now we can apply this fit to our model output.
In [129]:
#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)
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 [130]:
#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[130]:
[<matplotlib.lines.Line2D at 0x79b00d0>]
In [131]:
#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[131]:
[<matplotlib.lines.Line2D at 0x13d5bc90>]
In [132]:
#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)
Out[132]:
[<matplotlib.lines.Line2D at 0x14554450>]
In [133]:
#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)
Out[133]:
[<matplotlib.lines.Line2D at 0xd96fc10>]
In [134]:
#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)
Out[134]:
[<matplotlib.lines.Line2D at 0x1be14410>]
In [135]:
#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)
Out[135]:
[<matplotlib.lines.Line2D at 0x18688bd0>]
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 [136]:
def mean(diff):
return np.mean(abs(diff))
RMS Error¶
In [137]:
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 [138]:
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 [176]:
#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))
Saving the results¶
We will now save these statistics in a spreadsheet
In [108]:
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 [109]:
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' )
Out[109]:
<matplotlib.legend.Legend at 0xdc7a850>
In [181]:
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])
Green: D error <= 0.05 m, Yellow: 0.05 m < D error <= 0.1 m, Red: D error > 0.1 m
In [111]:
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')
Out[111]:
<matplotlib.text.Text at 0xe25c090>
In [112]:
#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
In [113]:
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])
M2
1.14308359042 0.0276748723456
1.17 -0.0269164095794
S2
South
0.262078298327 0.0062534109364
22.082293453
North
0.319881225881 0.00354225144697
17.625497454 2.02667467963
N2
South
0.22706545742 0.0134577946162
[[ 29.64920466 28.47956307 28.6954054 ]
[ 30.01288391 28.96340385 29.17468752]
[ 31.83671097 31.27244748 31.02132175]
[ 30.85940277 30.02646649 29.90739474]]
[[-0.02588851 -0.78357794 2.06147684]
[ 0.08708312 -0.41956367 2.42481443]
[ 1.62589875 0.06752611 2.55098799]
[ 1.58874627 0.23379867 2.86045696]]
-28.9689278 2.7870453083
North
0.174558732276 0.00610342281199
-49.1341775611 5.39293575174
In [114]:
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')
Out[114]:
<matplotlib.text.Text at 0xaab7450>
In [115]:
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')
Out[115]:
[<matplotlib.lines.Line2D at 0x1027ae50>]
In [115]:
In [115]:
In [115]:
In [115]:
In [115]:
In [115]: