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.
In [162]:
# 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 [240]:
# 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 = 'corr15'
runname = 'topogbottfric'
#joining the two string together
name = path +runname +'/'
print name
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 [164]:
# 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 [165]:
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
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 [166]:
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
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 [241]:
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 [168]:
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])
Out[168]:
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 [242]:
#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
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
new = 'true'
if new == 'true':
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
In [170]:
# 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 [171]:
# 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 [172]:
# 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 [173]:
# 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.))
In [174]:
# function for fitting 8 frequencies
def octuple(x, M2amp, M2pha, K1amp, K1pha, O1amp, O1pha, S2amp, S2pha,
P1amp, P1pha, N2amp, N2pha, Q1amp, Q1pha, K2amp, K2pha):
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.)+
Q1amp*np.cos(Q1freq*x-Q1pha*np.pi/180.)+
K2amp*np.cos(K2freq*x-K2pha*np.pi/180.))
Now we can apply this fit to our model output.
In [243]:
fig, ax = plt.subplots(1,1,figsize=(12,5))
for stn in (0,4,14,23):
print stations[stn]
fT1 = NC.Dataset(name+stations[stn]+'.nc','r')
time = fT1.variables["time_counter"][:]/3600. # want hours not seconds
ssh = fT1.variables["sossheig"][:,0,0]
ax.plot(time,ssh)
In [244]:
print ssh.shape
print 25*48*2, 30*24*4
In [266]:
#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=[]
Q1_amp=[]; Q1_pha=[]; K2_amp=[]; K2_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)
Q1_amp_obs=np.zeros(numsta); Q1_pha_obs=np.zeros(numsta)
K2_amp_obs=np.zeros(numsta); K2_pha_obs=np.zeros(numsta)
ts = 240
te = ssh.shape[0]
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(octuple,time[ts:te],ssh[ts:te])
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
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]*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)
if fitted[12] < 0:
fitted[12] = -fitted[12]
fitted[13] = fitted[13]+180
Q1_amp.append(fitted[12]*Q1ft)
pha= fitted[13]+Q1uvt
if pha > 360:
pha=pha-360
Q1_pha.append(pha)
if fitted[14] < 0:
fitted[14] = -fitted[14]
fitted[15] = fitted[15]+180
K2_amp.append(fitted[14]*K2ft)
pha= fitted[15]+K2uvt
if pha > 360:
pha = pha-360
if pha < 0:
pha = pha + 360
K2_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]
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]
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]
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]
#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)
Q1_amp_obs =np.ma.masked_values(Q1_amp_obs, 0)
Q1_pha_obs =np.ma.masked_values(Q1_pha_obs, 0)
K2_amp_obs =np.ma.masked_values(K2_amp_obs, 0)
K2_pha_obs =np.ma.masked_values(K2_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 [246]:
#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[246]:
In [247]:
#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[247]:
In [210]:
#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[210]:
In [211]:
#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[211]:
In [212]:
#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[212]:
In [213]:
#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[213]:
In [214]:
#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)
Out[214]:
In [215]:
#Plotting - K2
fig=tidetools.plot_scatter_pha_amp(K2_amp,K2_amp_obs,K2_pha,K2_pha_obs,'K2',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[215]:
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 [186]:
def mean(diff):
return np.mean(abs(diff))
RMS Error¶
In [187]:
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 [188]:
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 [267]:
#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)*1.03,np.array(M2_pha)+2.3)
D_K1= complex_diff(np.array(K1_amp_obs),np.array(K1_pha_obs), np.array(K1_amp)*0.99,np.array(K1_pha)-0.5)
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))
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))
D_K2= complex_diff(np.ma.array(K2_amp_obs),np.ma.array(K2_pha_obs), np.ma.array(K2_amp),np.ma.array(K2_pha))
Saving the results¶
We will now save these statistics in a spreadsheet
In [249]:
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 [269]:
plt.figure(figsize=(20,12))
plt.subplot(3,2,1)
plt.plot(np.array(M2_amp)*1.03, '-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(np.array(K1_amp)*0.99, '-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+2.3, '-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-0.5, '-bo', label = 'model')
plt.plot(K1_pha_obs, 'r-o', label = 'observation')
plt.title('K1 Phase')
plt.subplot(3,2,5)
plt.plot(D_M2, '-bo', label = 'M2')
plt.plot(D_K1, '-go', label = 'K1')
plt.plot((0,30),(0.05,0.05),'k')
plt.plot((0,30),(0.10,0.10),'r')
plt.title('D error')
plt.legend( loc='upper left' )
Out[269]:
In [251]:
M2_amp_topogbottfric = M2_amp
M2_pha_topogbottfric = pha_uwm
In [234]:
M2_amp_rnshlat2 = M2_amp
M2_pha_rnshlat2 = pha_uwm
In [219]:
M2_amp_bot1em2B = M2_amp
M2_pha_bot1em2B = pha_uwm
In [192]:
M2_amp_bot1em2 = M2_amp
M2_pha_bot1em2 = pha_uwm
In [148]:
M2_amp_corr15 = M2_amp
M2_pha_corr15 = pha_uwm
In [104]:
M2_amp_topog = M2_amp
M2_pha_topog = pha_uwm
In [73]:
M2_amp_rnshlat = M2_amp
M2_pha_rnshlat = pha_uwm
In [256]:
plt.figure(figsize=(12,5))
plt.plot(np.array(M2_amp_topog)*0.97, '-bo', label = 'topog')
plt.plot(M2_amp_obs, 'r-s', label = 'observation')
plt.plot(np.array(M2_amp_rnshlat2)*1.08, '-m^')
plt.plot(np.array(M2_amp_bot1em2B)*1.07, '-c^')
plt.plot(M2_amp_corr15, '-go', label='corr15')
plt.plot(np.array(M2_amp_topogbottfric)*1.03, '-yo')
makeit = np.array(M2_amp_corr15) + (np.array(M2_amp_bot1em2B)*1.07-np.array(M2_amp_corr15)) + (np.array(M2_amp_topog)*0.97-np.array(M2_amp_corr15))
plt.plot(makeit, 'k*-')
Out[256]:
In [261]:
plt.figure(figsize=(12,5))
plt.plot(np.array(M2_pha_topog)+2.9, '-bo', label = 'topog')
pha_uw = 180./np.pi * np.unwrap(np.array(M2_pha_obs)*np.pi/180.)
plt.plot(pha_uw, 'r-s', label = 'observation')
plt.plot(np.array(M2_pha_rnshlat2)-0.8, '-m^')
plt.plot(np.array(M2_pha_bot1em2B)-0.5, '-c^')
plt.plot(M2_pha_corr15, '-go', label = 'corr15')
plt.plot(np.array(M2_pha_topogbottfric)+2, '-yo')
makeit = np.array(M2_pha_corr15) + (np.array(M2_pha_bot1em2B)-0.8-np.array(M2_pha_corr15)) + (np.array(M2_pha_topog)+2.9-np.array(M2_pha_corr15))
plt.plot(makeit, 'k*-')
Out[261]:
In [270]:
cmap = plt.get_cmap('PuBu')
cmap.set_bad('burlywood')
fig,axs=plt.subplots(4, 2, figsize=(8,20))
constituent = ('M2', 'K1', 'O1', 'S2', 'P1', 'N2', 'Q1', 'K2')
error_D = (D_M2, D_K1, D_O1, D_S2, D_P1, D_N2, D_Q1, D_K2)
for row in range(4):
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])